运动估计算法MATLAB课程设计

1、通信工程应用技术综合训练与实习课程设计课程设计任务书学生姓名： 专业班级： 指导教师： 工作单位： 题 目: 运动估计算法 初始条件： MATLAB软件平台要求完成的主要任务:1. 设计任务设计视频压缩系统中的运动估计算法：全搜索法（FS: Full Search）和三步法(TSS: Three Step Search)，比较二种方法的搜索点和每帧的峰值信噪比（PSNR： peak signal to noise ratio）2. 设计要求 编制算法代码；对视频进行运动估计；计算PSNR时间安排： 答辩时间2013年1月24日。指导教师签名： 年 月 日系主任签名： 年 月 日目 录摘 要11

2、 运动估计算法概念21.1 运功估计算法基本思想21.2 运动估计算法实验原理22 设计原理和方法32.1 三步法32.2 新三步法32.3 全搜索法42.4 峰值信噪比53 运动估计算法的MATLAB编程63.1全搜索法63.2三步法93.3全搜索法指标113.4三步法指标113.4仿真结果分析114 小结与体会11参考文献11附录12摘 要运动估计的基本思想是尽可能准确地获得序列图像帧间的运动位移，即运动矢量。因为运动估计越准确，预测补偿的图像质量越高，补偿的残差就越小，补偿编码所需位数越少，需要传输的比特率就越小。利用得到的运动矢量在帧间进行运动补偿。补偿残差经过变换、量化、编码后与运动

3、矢量一起经过熵编码，然后以比特流形式发送出去。在视频编码和处理系统中，运动估计和运动补偿技术对降低视频序列时间冗余度、提高编码效率起着非常关键的作用。运动估计的准确程度将直接决定视频编码器的编码效率。关键词：运动估计、运动补偿技术、位移（运动）矢量AbstractThe basic idea is that the motion estimation as accurate as possible the image sequence interframe motion displacement, i.e. the motion vector. Motion estimation more a

4、ccurate prediction compensation, the higher the image quality is compensated residuals is smaller, less compensation coding bits required, the smaller the transmission bit rate. Performing motion compensation using the motion vector obtained in the interframe. Compensation residuals through transfor

5、mation, quantization, entropy-coded together with the motion vector is encoded, and then sent out in the form of a bit stream.In video coding and processing system, the motion estimation and motion compensation to reduce the temporal redundancy of video sequence to improve the coding efficiency play

6、s a crucial role. The degree of accuracy of the motion estimation will directly determine the encoding efficiency of the video encoder.Keywords: Motion estimation Motion compensation techniques The vector of displacement (movement)多媒体信息处理1 运动估计算法概念视频原始图像中存在着大量的信息冗余，如时间冗余、空间冗余、信息熵冗余、谱间冗余、几何结构冗余、视觉冗余和

7、知识冗余等等。运动估计是视频压缩编码中的核心技术之一，采用运动估计和运动补偿技术可以消除视频信号的时间冗余以提高编码效率。如何提高运动估计的效率，使运动估计算法的搜索过程更健壮、更快速、更高效成为目前研究的热点。运动估计的基本思想是尽可能准确地获得序列图像帧间的运动位移，即运动矢量。因为运动估计越准确，预测补偿的图像质量越高，补偿的残差就越小，补偿编码所需位数越少，需要传输的比特率就越小。利用得到的运动矢量在帧间进行运动补偿。补偿残差经过变换、量化、编码后与运动矢量一起经过熵编码，然后以比特流形式发送出去。运动估计算法多种多样，大体上可以把它们分成四类：块匹配法、递归估计法、贝叶斯估计法和光流

8、法。其中块匹配运动估计算法因其具有算法简单、便于VLSI实现等优点得到广泛应用。所以本文将重点介绍块匹配运动估计算法，并对各种块匹配算法在计算速度和估计精度上进行简单比较。1.1 运功估计算法基本思想运动估计的基本思想是尽可能准确地获得序列图像帧间的运动位移，即运动矢量。运动估计算法的目标是效率和准确性。由于在成象的场景中一般有多个物体作不同的运动，如果直接按照不同类型的运动将图像分割成复杂的区域是比较困难的。最直接和不受约束的方法是在每个像素都指定运动矢量，这就是所谓基于像素表示法。这种表示法是对任何类型图像都是适用的，但是它需要估计大量的未知量，并且它的解时常在物理上是不正确，除非在估计过

9、程中施加适当的物理约束。这在具体实现时是不可能的，通常采用基于块的物体运动表示法。因为运动估计越准确，预测补偿的图像质量越高，补偿的残差就越小，补偿编码所需位数越少，需要传输的比特率就越小。1.2 运动估计算法实验原理在帧间预测编码中，由于活动图像邻近帧中的景物存在着一定的相关性。因此，可将活动图像分成若干块或宏块，并设法搜索出每个块或宏块在邻近帧图像中的位置，并得出两者之间的空间位置的相对偏移量，得到的相对偏移量就是通常所指的运动矢量，得到运动矢量的过程被称为运动估计。 运动矢量和经过运动匹配后得到的预测误差共同发送到解码端，在解码端按照运动矢量指明的位置，从已经解码的邻近参考帧图像中找到相

10、应的块或宏块，和预测误差相加后就得到了块或宏块在当前帧中的位置。 运动估计的准确程度往往用补偿图像与原图像比较的PSNR来衡量表示。2 设计原理和方法2.1 三步法三步法是应用得相当广泛的一种次优的运动估计搜索算法它的搜索区间一般为-7，7 即在候选区中与编码块相同坐标位置处为原点，将参考块在其上下左右距离为7的范围内按照一定规律移动移到一个位置就做匹配计算它总共进行了三步搜索在下一次搜索时步长减半以前一步搜索得到的最优点为中心。下图为三步法的搜索示意图。三步法(TSS)搜索示意图算法的中心思想是，采用一种由粗到细的搜索模式，从原点开始，按一定步长取周围8个点构成每次搜索的点群，然后进行匹配计

11、算，利用上一步搜索得到的最小块误差MBD点作为当前搜索的中心位置，每做一步，搜索的步长减1。步搜索算法搜索窗选取(-7，+7)，最多只需要做25个位置的匹配计算，相对于全搜索来比，大大减少了匹配运算的复杂度，而且数据读取比较规则。2.2 新三步法TSS 假定运动矢量分布特点是在搜索窗口中均匀分布, 但事实证明运动矢量是偏置中心的，Renxiang Li等人在TSS的基础上提出了一种增强运动矢量中心偏置搜索和减小补偿误差的新三步法。新三步法(NTSS)搜索示意图NTSS 是对TSS 的一个改进，对运动量比较小的视频序列如可视电话序列有比较好的性能。对于绝大多数的视频序列，运动矢量的分布都是在中心

12、位置上的概率最大，随着与中心位置的距离的增大，概率会急剧地下降，这也就是前面所说的运动矢量的中心偏移特性。运动量比较小的视频序列的这一特性会更加明显。NTSS算法在最好的情况下只需要做17个点的匹配，在最坏的情况下需要做33个点的匹配，由于运动矢量中心偏置在现实视频序列中是普遍存在的，在通常情况下，NTSS算法需要做33点匹配的概率比较小，因此，在低速率视频应用中，如视频电话或视频会议中，NTSS算法的优点可以得到较好的发挥。2.3 全搜索法全搜索法（Full Search Method, FS）也称为穷尽搜索法，是对（M 2dx ）（N 2dy ）搜索范围内所有可能的候选位置计算MAD (i

13、,j)值，从中找出最小MAD，其对应偏移量即为所求运动矢量。此算法虽计算量大，但最简单、可靠，找到的必为全局最优点。FS算法描述如下：从原点出发，按顺时针螺旋方向由近及远，在逐个像素处计算MAD值，直到遍历搜索范围内听有的点，然后在计算的所有点的MAD中找到最小值，该点所在位置即对应最佳运动矢量。但是正因为它是穷尽搜索因此会产生巨大的计算量如 7, 7的搜索区间每个宏块16*16 需计算225个MAD值，这就直接制约了编码的实时实现。快速算法本质上是一种穷尽搜索法其计算量仍是相当巨大的。全搜索算法是简单也是效果最好的一种匹配算法，通过的全搜索匹配得到的结果是全局最优的，但由于计算量很大，我们在

14、编解码中往往不采用这种方法，而只把他作为与其他算法的一种比较。开始获取图像分成N*N大小的块计算MSE进行块匹配计算MAD比较结束图2-1程序流程图2.4 峰值信噪比峰值信噪比（经常缩写为PSNR）是一个表示信号最大可能功率和影响它的表示精度的破坏性噪声功率的比值的工程术语。由于许多信号都有非常宽的动态范围，峰值信噪比常用对数分贝单位来表示。PSNR是“Peak Signal to Noise Ratio”的缩写。peak的中文意思是顶点。而ratio的意思是比率或比列的。整个意思就是到达噪音比率的顶点信号，psnr是一般是用于最大值信号和背景噪音之间的一个工程项目。通常在经过影像压缩之后，输

15、出的影像通常都会有某种程度与原始影像不一样。为了衡量经过处理后的影像品质，我们通常会参考PSNR 值来认定某个处理程序够不够令人满意。PSNR计算公式如下：Peak 就是指8 bits 表示法的最大值255。MSE 指 Mean Square Error（均方误差，各值相差的n次方和的平均值的n次平方根(这几个字应该没有)），I(角标n)指原始影像第n个pixel 值，P(角标n)指经处理后的影像第n个pixel 值。PSNR 的单位为dB。所以PSNR值越大，就代表失真越少。PSNR 是最普遍，最广泛使用的评鉴画质的客观量测法，不过许多实验结果都显示，PSNR 的分数无法和人眼看到的视觉品质

16、完全一致，有可能 PSNR 较高者看起来反而比 PSNR 较低者差。这是因为人眼的视觉对于误差的敏感度并不是绝对的，其感知结果会受到许多因素的影响而产生变化（例如：人眼对空间频率较低的对比差异敏感度较高，人眼对亮度对比差异的敏感度较色度高，人眼对一个区域的感知结果会受到其周围邻近区域的影响）。3 运动估计算法的MATLAB编程3.1全搜索法 图3-1-1 the first frame 图3-1-2 the second frame 图3-1-3 帧间差值 图3-1-4 DFD图3-1-5 恢复后的第二帧图像图3-1-6 第一帧图像的运动矢量图图3-1-7 第二帧图像的运动矢量图3.2三步法

17、图3-2-1 the first frame 图3-2-2 the second frame 图3-2-3 帧间差值 图3-2-4 DFD图3-2-5 恢复后的第二帧图像图3-2-6 第一帧图像的运动矢量图图3-2-7 第二帧图像的运动矢量图3.3全搜索法指标全搜索算法：耗时长为Elapsed time is 1.053320 seconds.3.4三步法指标三步法：耗时长为Elapsed time is 0.881349 seconds.3.4仿真结果分析根据峰值信噪比可知全搜索算法的效果比三步法的效果要好；但全搜索所消耗的时间要比三步法所消耗的时间长，其运行较慢，因为全搜索算法的搜索点数要

18、比三步法的搜索点数多很多。 4 小结与体会 通过本次对全搜索块匹配算法课程设计，让我更深刻的理解了多媒体信息的获取、搜索和保存的基本方法。基于块的运动估计，是先将视频图像分成一个个规则的图像块，然后对每个图像块估计运动矢量。基于块的运动估计和运动补偿已经广泛应用于各种视频压缩编码标准。因此，本次课程设计是非常有作用的，让我们清醒的认识到自己还有好大的差距，它让我了解了多媒体处理的基本原理。课程设计中会遇到很多很多的困难，但在同学的帮助和自己思考下最终还是没有解决问题。所以，本次课程设计我个人的收获不大，但是让我知道自己还在图像编程方面有很多欠缺的地方，不仅对均方误差(MSE)、绝对值误差(MA

19、D)更深的理解，并且采用全搜索块匹配方法得到所有宏块的运动矢量，把c语言的知识又重新的过了一遍，同时还理解了些多媒体信息处理的基本原理。参考文献1罗军辉. MATLAB7.0在图像处理中的应用. 机械出版社.20052王家文. MATLAB7.0编程基础.机械出版社.20053钟麟. MATLAB仿真技术与应用教程.国防工业出版社.20044冈萨雷斯. 数字图像处理(MATLAB版).电子工业出版社.2005附录FS：clear all;I1=imread(claire1.bmp); %read the first frameI2=imread(claire2.bmp); %read the

20、second framedm=7;%给图像扩边，每个边都扩大dm大小% |-|-|-|% | B | C | D |% |-|-|-|% | | | |% | | | |% | E | A | F | - 图像% | | | |% | | | |% | | | |% |-|-|-|% | G | H | I |% |-|-|-|% I1=double(I1);I2=double(I2);rownum colnum = size(I1); II=zeros(rownum+2*dm,colnum+2*dm);II(dm+1:dm+rownum,dm+1:dm+colnum)=I1;for i=1:

21、dm II(i,dm+1:dm+colnum)=II(dm+1,dm+1:dm+colnum); II(rownum+dm+i,dm+1:dm+colnum)=II(dm+rownum,dm+1:dm+colnum);endfor j=1:dm II(1:rownum+2*dm,j)=II(1:rownum+2*dm,dm+1); II(1:rownum+2*dm,colnum+dm+j)=II(1:rownum+2*dm,dm+colnum);end %下面进行全搜索算法 blocksize=16; rowblocks =rownum/blocksize; colblocks =colnum

22、/blocksize; A=99999999999999999999; %为了找到最小的均方误差，A用于设定一个很大的初值 Eij=0; xrecord=ones(16,16); %xrecord，yrecord用于存放匹配快的块号，即运动矢量 yrecord=ones(16,16); diff=zeros(256,256); %这幅图的大小为256*256 tic for x=0:(rowblocks-1) %x表示行中第几个子块 row=x*blocksize; for y=0:(colblocks-1) %y表示列中第几个子块 col=y*blocksize;% tempx=x*bloc

23、ksize+1:(x+1)*blocksize;% tempy=y*blocksize+1:(y+1)*blocksize; for p=-dm:dm for q=-dm:dm %（p,q）表示x,y对应子块在前一帧所的搜索位置 Eij=0; Eij=sum(sum(I2(row+1:row+blocksize,col+1:col+blocksize)-II(row+dm+p+1:row+dm+p+blocksize,col+dm+q+1:col+dm+q+blocksize).2)/(blocksize2); if EijA A=Eij; xrecord(x+1,y+1)=p; yrecor

24、d(x+1,y+1)=q; end end end A=999999999999999999; for mx=1:blocksize for ny=1:blocksize diff(row+mx,col+ny)=I2(row+mx,col+ny)-II(row+mx+dm+xrecord(x+1,y+1),col+ny+dm+yrecord(x+1,y+1); end end end endtoc figure,imshow(I1,); title(the first frame); figure,imshow(I2,); title(the second frame); IIII=I2-I1

25、; figure,imshow(IIII,); title(帧间差值); figure,imshow(diff,); title(DFD); %title(利用全搜索算法匹配后的帧间差); for x=0:(rowblocks-1) % row=x*blocksize; for y=0:(colblocks-1) % col=y*blocksize; III(row+1:row+blocksize,col+1:col+blocksize)=II(row+dm+xrecord(x+1,y+1)+1:row+dm+xrecord(x+1,y+1)+blocksize,col+dm+yrecord(

26、x+1,y+1)+1:col+dm+yrecord(x+1,y+1)+blocksize)+diff(row+1:row+blocksize,col+1:col+blocksize); end end %III=I1+abs(diff); figure,imshow(III,); title(); ERR=diff;% figure,imshow(ERR,);% title(DFD); numberarray=0:1:255; for m=1:255 numberarray(m+1)=0; end; zeronumber=0; for n=1:rownum for m=1:colnum dif

27、=abs(ERR(m,n); if(dif=0) temp=zeronumber; zeronumber=temp+1; else numberarray(dif)=numberarray(dif)+1; end; end; end; figure;plot(0,zeronumber,k*);hold on;plot(numberarray,r*),title(DFD distribution);hold off; ERR1=zeros(16,16); for i=0:15 for j=0:15 ERR1(i+1,j+1)=round(sum(sum(ERR(i*blocksize+1:i*b

28、locksize+blocksize,j*blocksize+1:j*blocksize+blocksize)/(blocksize*blocksize); end end numberarray=0:1:255; for m=1:255 numberarray(m+1)=0; end; zeronumber=0; for n=1:16 for m=1:16 dif=abs(ERR1(m,n); if(dif=0) temp=zeronumber; zeronumber=temp+1; else numberarray(dif)=numberarray(dif)+1; end; end; en

29、d; figure;plot(0,zeronumber,k*);hold on;plot(numberarray,r*),title(DFD(block average) distribution);hold off; %figure;mesh(diff);figure;contour(diff,15); figure; for i=1:16 for j=1:16 quiver(i,j,xrecord(i,j)/16,yrecord(i,j)/16); hold on; end end grid on;figure;quiver(1:16,1:16,yrecord,xrecord); grid

30、 on; S_3SS：clear all; I1=imread(claire1.bmp); %read the first frameI2=imread(claire2.bmp); %read the second frame dm=7;I1=double(I1);I2=double(I2);rownum colnum = size(I1);II=zeros(rownum+2*dm,colnum+2*dm);II(dm+1:dm+rownum,dm+1:dm+colnum)=I1;for i=1:dm II(i,dm+1:dm+colnum)=II(dm+1,dm+1:dm+colnum);

31、II(rownum+dm+i,dm+1:dm+colnum)=II(dm+rownum,dm+1:dm+colnum);endfor j=1:dm II(1:rownum+2*dm,j)=II(1:rownum+2*dm,dm+1); II(1:rownum+2*dm,colnum+dm+j)=II(1:rownum+2*dm,dm+colnum);end blocksize=16; rowblocks =rownum/blocksize; colblocks =colnum/blocksize; A=99999999999999999999; % Eij=0; xrecord=ones(16

32、,16); % yrecord=ones(16,16); diff=zeros(256,256); % tic for x=0:(rowblocks-1) % row=x*blocksize; for y=0:(colblocks-1) % col=y*blocksize;% tempx=x*blocksize+1:(x+1)*blocksize;% tempy=y*blocksize+1:(y+1)*blocksize; for p1=-4:4:4 % for q1=-4:4:4 % Eij=0; Eij=sum(sum(I2(row+1:row+blocksize,col+1:col+bl

33、ocksize)-II(row+dm+p1+1:row+dm+p1+blocksize,col+dm+q1+1:col+dm+q1+blocksize).2)/(blocksize2); if EijA A=Eij; xrecord(x+1,y+1)=p1; yrecord(x+1,y+1)=q1; end end end p1=xrecord(x+1,y+1); q1=yrecord(x+1,y+1); for p2=p1-2:2:p1+2 % for q2=q1-2:2:q1+2 if p2=p1 | q2=q1 Eij=0; Eij=sum(sum(I2(row+1:row+blocks

34、ize,col+1:col+blocksize)-II(row+dm+p2+1:row+dm+p2+blocksize,col+dm+q2+1:col+dm+q2+blocksize).2)/(blocksize2); if EijA A=Eij; xrecord(x+1,y+1)=p2; yrecord(x+1,y+1)=q2; end end end end p2=xrecord(x+1,y+1); q2=yrecord(x+1,y+1); for p3=p2-1:1:p2+1 % for q3=q2-1:1:q2+1 if p3=p2 | q3=q2 Eij=0; Eij=sum(sum

35、(I2(row+1:row+blocksize,col+1:col+blocksize)-II(row+dm+p3+1:row+dm+p3+blocksize,col+dm+q3+1:col+dm+q3+blocksize).2)/(blocksize2); if EijA A=Eij; xrecord(x+1,y+1)=p3; yrecord(x+1,y+1)=q3; end end end end A=999999999999999999; for mx=1:blocksize for ny=1:blocksize diff(row+mx,col+ny)=I2(row+mx,col+ny)

36、-II(row+mx+dm+xrecord(x+1,y+1),col+ny+dm+yrecord(x+1,y+1); end end end endtoc figure,imshow(I1,); title(the first frame); figure,imshow(I2,); title(the second frame); IIII=I2-I1; figure,imshow(IIII,); title(); figure,imshow(diff,); title(DFD); for x=0:(rowblocks-1) row=x*blocksize; for y=0:(colblocks-1) col=y*blocksize; III(row+1:row+blocksize,col+1:col+blocksize)=II(row+dm+xrecord(x+1,y+1)+1:row+dm+xrecord(x+1,y+1)+blocksize,col+dm+yrecord(x+1,y+1)+1:col+dm+yrecord(x+1,y+1)+blocksize)+diff(row+1:row+b