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

   % This script uses all the Motion Estimation algorithms written for the

% final project and 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 search method（全搜索法计算运动矢量）

%

% Input

%   imgP : The image for which we want to find motion vectors（当前图像）

%   imgI : The reference image（参考图像）

%   mbSize : Size of the macroblock（宏块尺寸）

%   p : Search parameter  (read literature to find what this 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 look for

% a close match p pixels on the left, right, top and bottom 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        

        % the exhaustive search starts here（全搜索开始）

        % we will evaluate cost for  (2p + 1) blocks vertically

        % and (2p + 1) blocks horizontaly（我们将计算水平方向上(2p + 1)个块的最小绝对差和垂直方向上(2p + 1)个块的最小绝对差）

        % m is row(vertical) index（m为行指数）

        % n is col(horizontal) index（n为列指数）

        % this means we are scanning in raster order

    for m = -p : p   %水平方向上位移矢量范围

    for n = -p : p %垂直方向上位移矢量范围

   % 补充下面程序

   % row/Vert co-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 we find the vector where the cost is minimum

        % and store it ... this is what will be passed back.（现在找到有最小绝对差的矢量并存储它，这就是将被返回的东西）

        % 补充下面程序

        blockcenter(1,mbCount) = i+p;  

        blockcenter(2,mbCount) = j+p;       

        % finds which macroblock in imgI gave us min Cost（找到参考图像中最小绝对差的宏块）            

        [dx, dy, min] = minCost(costs);         

        % row co-ordinate for the vector（矢量的行集合）

        vectors(1,mbCount) =  dx;  

        % col co-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 the given 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 : the side 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 minimum cost（找到拥有最小绝对差的点的指数）

%

% Input

%   costs : The matrix that contains the estimation costs for a

%   macroblock（包含宏块的估计代价的矩阵）

%

% Output

%   dx : the motion vector component in columns（列方向上运动矢量组成）

%   dy : the motion 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 less then the already

% present value in min.

% If its inded smaller then we swap the min value with the 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=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 given motion vectors（用给定的运动矢量计算运动补偿图像）

%

% Input

%   imgI : The reference 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 read the motion vector（对于看到的每一个宏块，读出它的运动矢量）

% and put that macroblock from refernce image in the compensated image（并将参考图像中的该宏块放到补偿图像中）

mbCount = 1;

for i = 1:mbSize:row-mbSize+1

    for j = 1:mbSize:col-mbSize+1

        % dy is row(vertical) index（dy为垂直方向上的指数）

        % dx is col(horizontal) index（dx为水平方向上的指数）

        % this means 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 in the images（图像中任何一个像素的可能的峰值）

%

% Ouput

%   psnr : The motion compensated image's PSNR（运动补偿图像的峰值信噪比）

%

% Written by Aroh Barjatya

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;