1. 主文件motionEstAnalysis.m
% This script uses all the Motion Estimation algorithms written for the
% final projectand save their results.
close all
clear all
% imageName =’caltrain.avi’;
VideoName =’shaky_car.avi’;
video =aviread(VideoName);
% movie(video);
mbSize = 16;
p = 7;
for i = 1:6
imgINumber = i;
imgPNumber = i+2;
videoI = video(imgINumber);
videoP = video(imgPNumber);
imgI = double(videoI.cdata);
imgP = double(videoP.cdata);
[row col] = size(imgI);
% Exhaustive Search 基于块的全搜索算法
[BlockCenter, motionVect, computations] =motionEstES(imgP,imgI,mbSize,p);
% P 帧当前重构图像
imgPComp = motionComp(imgI, motionVect, mbSize);
% P 帧当前图像 和 P 帧当前重构图像的PSNR值
ESpsnr(i+1) = imgPSNR(imgP, imgPComp, 255);
EScomputations(i+1) = computations;
% P 帧当前重构误差图像
imagePDiff = imgP – imgPComp;
if i == 4
figure;
imageI = uint8(imgI);
imageP = uint8(imgP);
imagePComp = uint8(imgPComp);
imagePDiff = uint8(imagePDiff);
subplot(221);imshow(imageI);
title(‘I 帧参考图像’);
subplot(222);imshow(imageP);
title(‘P 帧当前图像’);
subplot(223);imshow(imagePComp);
title(‘P 帧当前重构图像’);
subplot(224);imshow(imagePDiff);
title(‘P 帧当前重构误差图像’);
% 画运动矢量图
figure;
quiver( BlockCenter(2,:),BlockCenter(1,:), motionVect(2,:), motionVect(1,:), .2,’r’);
axis([0 320 0 240]);
for i=mbSize:mbSize:col-mbSize
x = [i,i];
y = [0,row];
line(x,y,’LineStyle’,’-‘,’Marker’,’none’);
end
for j=mbSize:mbSize:row-mbSize
x = [0,col];
y = [j,j];
line(x,y,’LineStyle’,’-‘,’Marker’,’none’);
end
xlabel(‘X’);
ylabel(‘Y’);
end
end
2. 文件motionEstES.m
% Computes motion vectors using exhaustive searchmethod（全搜索法计算运动矢量）
%
% Input
% imgP : Theimage for which we want to find motion vectors（当前图像）
% imgI : Thereference image（参考图像）
% mbSize :Size of the macroblock（宏块尺寸）
% p : Searchparameter (read literature to find whatthis means)（搜索参数）
%
% Ouput
% motionVect :the motion vectors for each integral macroblock in imgP（当前图像中每一个积分宏块的运动矢量）
% EScomputations: The average number of points searched for a macroblock（每个宏块搜索的平均点数）
%
% Written by Aroh Barjatya
function [BlockCenter, motionVect, EScomputations] =motionEstES(imgP, imgI, mbSize, p) % 定义函数文件motionEstES.m，imgP、 imgI、 mbSize、 p为传入参数，BlockCenter、motionVect、 EScomputations为返回参数
[row col] = size(imgI); % 将参考图像的行数赋值给row，列数赋值给col
blockcenter = zeros(2,row*col/mbSize^2);
vectors = zeros(2,row*col/mbSize^2); % 定义全0的矢量矩阵的大小
costs = ones(2*p + 1, 2*p +1) * 65537; % 定义最小绝对差矩阵的大小
computations = 0; % 搜索点数赋初值为0
% we start off from the top left of the image（从图像左上角开始）
% we will walk in steps of mbSize（以宏块尺寸为步长）
% for every marcoblock that we look at we will lookfor
% a close match p pixels on the left, right, top andbottom of it （对于每一个宏块，在它的上下左右找到与搜索参数p最匹配的像素）
mbCount = 1; %搜索的宏块数赋初值为1
%1为循环起始值，mbSize为步长值，row-mbSize+1为循环终止值
for i = 1 : mbSize : row-mbSize+1
for j = 1 :mbSize : col-mbSize+1
% theexhaustive search starts here（全搜索开始）
% wewill evaluate cost for (2p + 1) blocksvertically
% and(2p + 1) blocks horizontaly（我们将计算水平方向上(2p + 1)个块的最小绝对差和垂直方向上(2p + 1)个块的最小绝对差）
% m isrow(vertical) index（m为行指数）
% n iscol(horizontal) index（n为列指数）
% thismeans we are scanning in raster order
for m = -p :p %水平方向上位移矢量范围
for n = -p :p %垂直方向上位移矢量范围
% 补充下面程序
% row/Vertco-ordinate for ref block （参考块的行(垂直方向)的范围）
refBlkVer= i+m;
%col/Horizontal co-ordinate（参考块的列(水平方向)的范围）
refBlkHor= j+n;
%如果参考块的行列范围的任意一个在已经搜索过的宏块之外，则继续下一步的搜索
if (refBlkVer < 1 || refBlkVer+mbSize-1 > row …
|| refBlkHor < 1 ||refBlkHor+mbSize-1 > col)
continue;
end
costs(m+p+1,n+p+1) = costFuncMAD(imgP(i:i+mbSize-1,j:j+mbSize-1), …
imgI(refBlkVer:refBlkVer+mbSize-1,refBlkHor:refBlkHor+mbSize-1), mbSize);
% 搜索下一个点
computations= computations + 1;
end
end
% Now wefind the vector where the cost is minimum
% andstore it … this is what will be passed back.（现在找到有最小绝对差的矢量并存储它，这就是将被返回的东西）
% 补充下面程序
blockcenter(1,mbCount) = i+p;
blockcenter(2,mbCount) = j+p;
% findswhich macroblock in imgI gave us min Cost（找到参考图像中最小绝对差的宏块）
[dx, dy,min] = minCost(costs);
% rowco-ordinate for the vector（矢量的行集合）
vectors(1,mbCount) = dx;
% colco-ordinate for the vector（矢量的列集合）
vectors(2,mbCount) = dy;
%搜索下一个宏块
mbCount = mbCount + 1;
costs = ones(2*p + 1, 2*p +1) * 65537;
end
end
BlockCenter = blockcenter;
motionVect =vectors; %返回当前图像中每一个积分宏块的运动矢量
EScomputations = computations/(mbCount – 1); %返回每个宏块搜索的平均点数
3. 文件costFuncMAD.m
% Computes the Mean Absolute Difference (MAD) for thegiven two blocks（对给定的两个块计算最小绝对差）
% Input
% currentBlk : The block for which we are finding the MAD（当前块）
% refBlk :the block w.r.t. which the MAD is being computed（参考块）
% n : theside of the two square blocks
%
% Output
% cost :The MAD for the two blocks（两个块的最小绝对差）
%
% Written by Aroh Barjatya
% 定义函数文件costFuncMAD.m，currentBlk、refBlk、 n为传入参数，cost为返回参数
function cost = costFuncMAD(currentBlk,refBlk, n)
% 补充下面程序
cost=sum(sum(abs(currentBlk-refBlk)))/(n*n);
4. 文件minCost.m
% Finds the indices of the cell that holds the minimumcost（找到拥有最小绝对差的点的指数）
%
% Input
% costs : Thematrix that contains the estimation costs for a
% macroblock（包含宏块的估计代价的矩阵）
%
% Output
% dx : the motionvector component in columns（列方向上运动矢量组成）
% dy : themotion vector component in rows（行方向上运动矢量组成）
%
% Written by Aroh Barjatya
function [dx, dy, min] = minCost(costs)
[row, col] = size(costs);
% we check whether the current value of costs is lessthen the already
% present value in min.
% If its inded smaller then we swap the min value withthe current one and
% note the indices.
% （检测costs的当前值是否比已经出现的最小值小。如果小的话，我们将当前值与最小值对调，并注明指数）
% 补充下面程序
minnum=256;
x=8;
y=8;
for i=1:row
for j=1:col
if(costs(i,j)<minnum)
minnum=costs(i,j);
x=i;
y=j;
end
end
end
dx=x;
dy=y;
min=minnum;
5. 文件motionComp.m
% Computes motion compensated image using the givenmotion vectors（用给定的运动矢量计算运动补偿图像）
%
% Input
% imgI : Thereference image （参考图像）
% motionVect :The motion vectors（运动矢量）
% mbSize :Size of the macroblock（宏块大小）
%
% Ouput
% imgComp :The motion compensated image（运动补偿图像）
%
% Written by Aroh Barjatya
function imgComp = motionComp(imgI, motionVect, mbSize)
[row col] = size(imgI);
% we start off from the top left of the image（从图像左上角开始）
% we will walk in steps of mbSize（以宏块大小为步长）
% for every marcoblock that we look at we will readthe motion vector（对于看到的每一个宏块，读出它的运动矢量）
% and put that macroblock from refernce image in thecompensated image（并将参考图像中的该宏块放到补偿图像中）
mbCount = 1;
for i = 1:mbSize:row-mbSize+1
for j =1:mbSize:col-mbSize+1
% dy isrow(vertical) index（dy为垂直方向上的指数）
% dx iscol(horizontal) index（dx为水平方向上的指数）
% thismeans we are scanning in order
dy =motionVect(1,mbCount);
dx =motionVect(2,mbCount);
refBlkVer = i + dy;
refBlkHor = j + dx;
if (refBlkVer < 1 || refBlkVer+mbSize-1 > row …
|| refBlkHor < 1 || refBlkHor+mbSize-1 > col)
imageComp(i:i+mbSize-1,j:j+mbSize-1)=imgI(i:i+mbSize-1,j:j+mbSize-1);
continue;
end
imageComp(i:i+mbSize-1,j:j+mbSize-1)= imgI(refBlkVer:refBlkVer+mbSize-1, refBlkHor:refBlkHor+mbSize-1);
mbCount= mbCount + 1;
end
end
imgComp = imageComp;
6. 文件imgPSNR.m
% Computes motion compensated image’s PSNR（计算运动补偿图像的峰值信噪比）
%
% Input
% imgP : The original image （原始图像）
% imgComp : The compensated image（补偿图像）
% n : the peak value possible of any pixel inthe images（图像中任何一个像素的可能的峰值）
%
% Ouput
% psnr : The motion compensated image’s PSNR（运动补偿图像的峰值信噪比）
%
% Written by ArohBarjatya
function psnr =imgPSNR(imgP, imgComp, n)
% 补充下面程序
MSE=(1/(n*n))*sum(sum((imgP-imgComp).^2));
PSNR=10*log(max(max(imgComp)).^2/MSE);
psnr=PSNR;