基于块的全搜索运动估计算法实现matlab代码

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. whi
97d4
ch 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;