数字图像处理 实验八：运动估计

DIP实验8：运动估计

实验目的

熟悉运动估计的块匹配（BMA）算法原理，编程实现全搜索算法（三步搜索或钻石搜索算法），了解运动估计在混合编码器中的作用。

实验内容

1）编写全搜索算法函数，将运动矢量叠加到当前帧上并显示输出；
2）显示输出预测帧、残差帧和重建图像，计算预测帧的PSNR。

附：可供参考的Matlab函数有hold、quiver

输出图像排列格式如下：

参考帧	当前帧	运动矢量图
预测帧	残差帧	重建帧

需要编写的函数：
1）SAD计算函数：
对给定的两个块计算SAD值

function cost = costSAD(currentBlk,refBlk,mbSize)
% Input
%       currentBlk : 当前块
%       refBlk : 参考块
%       mbSize : MB尺寸
% Output
%       cost : 两个块之间的误差代价（SAD）

2）求具有最小SAD值的函数：
找出具有最小代价的块的下标

function [dx, dy, minc] = minCost(costs)
% Input
%       costs : 包含当前宏块所有运动估计误差代价的SAD矩阵
% Output
%       dx : MV的垂直分量（行位移）
%       dy : MV的水平分量（列位移）

3）求预测帧的函数：
由给定的运动矢量进行运动补偿计算预测帧

function imgComp = motionComp(imgI, motionVect, mbSize)
% Input
%       imgI : 参考帧
%       motionVect : MV（dx为垂直分量，dy为水平分量）
%   	mbSize : MB尺寸
% Ouput
%   	imgComp : 运动补偿后的图像

4）PSNR函数：
计算尺寸相同的两幅图像之间的峰值信噪比

function psnr = imgPSNR(img, imgComp)
% Input
%   	img : 原始帧
%   	imgComp : 预测帧
% Ouput
%   	psnr : 运动补偿后的图像的PSNR

5）FS算法：
全搜索（Full Search/Exhaustive Search）

function [motionVect, EScomputations] = ME_ES(img, imgI, mbSize, dm)
% Input
%   	img : 当前帧
%   	imgI : 参考帧
%   	mbSize : MB尺寸
%   	dm : 搜索窗口大小（2dm+1）×（2dm+1）
% Ouput
%   	motionVect : 整像素精度MV
%   	EScomputations: 搜索每个宏块所需的平均点数

6）默认参数：

mbSize=16;      % 块尺寸为16*16
dm=7;		    % 匹配次数为(2dm+1)*(2dm+1)

---

参考代码

mbSize = 16;      % 块尺寸为16*16
dm = 7;			  % 匹配次数为(2dm+1)*(2dm+1)

imgI = imread('girl0.bmp');     % 定义参考帧
img = imread('girl1.bmp');      % 定义当前帧
subplot(231); imshow(imgI); title('参考帧');
subplot(232); imshow(img); title('当前帧');
imgI = double(imgI);
img = double(img);

[motionVect, blk_center,costs] = ME_ES(img, imgI, mbSize, dm); %基于块的全搜索算法  

subplot(233); imshow(uint8(img)); title('运动矢量图'); hold on
% 参考帧指向当前帧
y = blk_center(1, : );
x = blk_center(2, : );
v = -motionVect(1, : );
u = -motionVect(2, : ); 
quiver(x, y, u, v, 'g');
hold on

%运动补偿后的图像：根据运动矢量计算预测帧，并传输残差帧
imgComp = motionComp(imgI, motionVect, mbSize); 
%计算当前帧和预测帧的峰值信噪比
psnr = imgPSNR(img, imgComp) 
subplot(234); imshow(uint8(imgComp));
title('预测帧');

 
imgErr = img - imgComp; %残差帧
cal = Calibration(imgErr); %标定
subplot(235); imshow(cal); title('残差帧');
 
%根据运动矢量指明的位置及残差帧重建图像
rebuild = imgComp + imgErr; 
subplot(236); imshow(uint8(rebuild)); title('重建帧');

%% FS算法：全搜索（Full Search/Exhaustive Search）
function [motionVect,blk_center,costs] = ME_ES(imgP, imgI, mbSize, dm)
% Input
%   	img : 当前帧
%   	imgI : 参考帧
%   	mbSize : MB尺寸
%   	dm : 搜索窗口大小（2dm+1）×（2dm+1）
% Ouput
%   	motionVect : 整像素精度MV
%   	EScomputations: 搜索每个宏块所需的平均点数
    [row, col] = size(imgP);
    blk_center = zeros(2, row*col/(mbSize^2)); 
    %定义每个宏块中心点位置
    motionVect = zeros(2,row*col/(mbSize^2)); %定义每个宏块运动矢量
    costs = ones(2*dm+1,2*dm+1)*100000;
    computations = 0; %搜索的点数之和
    mb_cnt= 1;
    for i = 1:mbSize:row-mbSize+1     %当前帧起始行搜索范围，步长是块数 
        for j = 1:mbSize:col-mbSize+1 %当前帧起始列搜索范围，步长是块数
            for m= -dm: dm
                for n= -dm: dm
                    ref_blk_row = i+m; %参考帧搜索框起始行
                    ref_blk_col = j+n; %参考帧搜索框起始列
                    %超出搜索范围
                    if (ref_blk_row<1||ref_blk_row+mbSize-1>row||ref_blk_col<1||ref_blk_col+mbSize-1>col)                     
                        continue;                                                            
                    end
                    %计算SAD
                    costs(m+dm+1,n+dm+1) =...
                        costSAD(imgP(i:i+mbSize-1,j:j+mbSize-1),imgI(ref_blk_row:ref_blk_row+mbSize-1,ref_blk_col:ref_blk_col+mbSize-1));
                    computations = computations+1;
                end
            end
            blk_center(1,mb_cnt) = i+ mbSize/2-1; %记录中心点行坐标                     
            blk_center(2,mb_cnt) = j+ mbSize/2-1; %记录中心点列坐标
            [dx,dy,~]=minCost(costs); %找出有最小代价的块的下标
            motionVect(1,mb_cnt) = dx-dm-1; %垂直运动矢量
            motionVect(2,mb_cnt) = dy-dm-1; %水平运动矢量
            mb_cnt = mb_cnt+1;
            costs = ones(2*dm+1,2*dm+1)*100000; %重新赋值
         end
    end 
end

%% 求预测帧的函数：由给定的运动矢量进行运动补偿计算预测帧
function imgComp = motionComp(imgI, motionVect, mbSize)
% Input
%       imgI : 参考帧
%       motionVect : MV（dx为垂直分量，dy为水平分量）
%   	mbSize : MB尺寸
% Ouput
%   	imgComp : 运动补偿后的图像
    [row,col]=size(imgI);
    mb_cnt=1;
    for i = 1:mbSize:row-mbSize+1                
        for j = 1:mbSize:col-mbSize+1 
             ref_blk_row=i+motionVect(1,mb_cnt); %参考帧搜索块起始行
             ref_blk_col=j+motionVect(2,mb_cnt); %参考帧搜索块起始列
             imgComp(i:i+mbSize-1,j:j+mbSize-1)=imgI(ref_blk_row:ref_blk_row+mbSize-1,ref_blk_col:ref_blk_col+mbSize-1);   
             mb_cnt=mb_cnt+1;
        end
    end
end

%% 求具有最小SAD值的函数：找出具有最小代价的块的下标
function [dx,dy,minc] = minCost(costs)
% Input
%       costs : 包含当前宏块所有运动估计误差代价的SAD矩阵
% Output
%       dx : MV的垂直分量（行位移）
%       dy : MV的水平分量（列位移）
    minc = min(min(costs));
    [dx, dy] = find(costs == minc);
    dx = dx(1);
    dy = dy(1);    
end

%% SAD计算函数：对给定的两个块计算SAD值
function cost = costSAD(currentBlk,refBlk)
% Input
%       currentBlk : 当前块
%       refBlk : 参考块
%       mbSize : MB尺寸
% Output
%       cost : 两个块之间的误差代价（SAD）
    cost=sum(sum(abs(currentBlk-refBlk))); 
end

%% PSNR函数：计算尺寸相同的两幅图像之间的峰值信噪比
function psnr = imgPSNR(img, imgComp)
% Input
%   	img : 原始帧
%   	imgComp : 预测帧
% Ouput
%   	psnr : 运动补偿后的图像的PSNR
    [M,N]=size(img);
    MSE=sum(sum((img-imgComp).^2))/(M*N);
    psnr=10*log10(255^2/MSE);
end

%% 标定
function cal = Calibration(img)
    img = double(img);
    [M, N] = size(img);
    fmin = min(min(img));
    fm = img - fmin * ones(M, N);
    fmmax = max(max(fm));
    fs = 255 * fm ./ fmmax;
    cal = uint8(fs);
end

实验结果