基于块的全搜索运动估计算法实现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. 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;