基于块匹配算法的运动估计

基于块匹配算法的运动估计摘要：本文首先介绍了运动估计和块匹配的概念和思想，然后详细介绍了块匹配的搜索算法。最后根据块匹配算法的运动估计给出了一个设计。关键字：运动估计，块匹配，算法1研究背景随着现代信息社会对通信业务要求的不断增长，图像通信与通信网容量的矛盾日益突出，尤其是具有庞大数据量的视频图像通信，更是很难传输和存储，极大的制约了图像通信的开展。例如，按CCIR601建议，普通质量的电视信号数字视频的码率约为216Mbit/s，而高清晰度电视HDTV那么在1.2Gb/s以上，如果没有高效率的压缩技术，那么难以传输和存储。而运动检测与估计不仅是军事领域中目标探测与跟踪的有效技术之一，同时也是视频图象编码与压缩的常用方法与核心技术。因此，该项研究对于复杂背景下的目标探测以及视频处理均具有重要的应用价值。运动估计和运动补偿是紧密联系的，它是视频图像压缩编码中使用的一项核心技术，很好的解决了视频图像中时间冗余的问题，经验说明，实用化的压缩方法可以将运动图像数据压缩30倍而不失真。运动估计技术主要分为两大类:象素递归法和块匹配法。考虑到计算复杂度和实时性要求，块匹配法已成为目前最常用的方法。显然，要想获得好的运动补偿，关键是要有准确的运动估计，因此运动估计算法的研究成为视频压缩算法的研究重点。运动估计研究的主要内容就是如何快速、有效的获得有足够精度的运动矢量。即把前一帧的运动局部根据运动矢量补过来，同时用其它方法得到其剩余的不同局部的过程称为运动补偿(MotionCompensation,MC)。就这样，采用运动估计和补偿可以有效地去除视频信号在时间方向的重复信息，到达压缩的目的。其中，在这方面的一种有效方法就是块匹配运动估计BMME(Block-MatchingMotionEstimation),它目前已被许多视频编码标准所采纳。为了提高搜索速度和效率，目前研究最多的是基于块匹配(BlockMatching,BM)的快速搜索算法，例如三步法，四步法，二维对数法，菱形法等。随着计算机网络的普及和开展，很多信息可以通过网络实现共享。形色兼具的视频信息在网络传输的信息中所占据的比例也越来越高。但由于视频信息本身十分庞大，限制了其在网络中的传播速度。于是视频信息的压缩成为视频传输的一个重要环节。数字视频信息是由数字图像的时间序列构成的，每一幅数字图像称为一帧。视频编码的一个主(frame)要目的就是在保证一定重构质量的前提下，以尽量少的比特数来表征视频信息。视频编码是视频压缩的关键技术。视频编码与的原理是：由于表示图像和视频信息所需的大量的数据往往是高度相关的，这些相关性会引起信息的冗余，因此可以通过去除这些冗余信息来实现对视频数据的压缩。静止图像的压缩是在保持重建图像质量可以接收的同时，尽量去除图像本身存在的空间冗余，而视频信号的压缩，除了去除空间冗余之外，还可以通过去除时间冗余到达较高的压缩比。运动估计技术是视频图像压缩编码中使用的一项核心技术，很好的解决了视频图像中时间冗余的问题，运动估计技术主要分为两大类:象素递归法和块匹配算法。考虑到计算复杂度和实时性要求，块匹配法已成为目前最常用的方法。2实验原理如图1所示,运动估计的根本思想是将图像序列的每一帧图像分成许多互不重叠的宏块,并假设块内各像素只作相等的平移,然后对于当前帧中的每一块到前一帧或后一帧某一给定搜索范围内根据一定的匹配准那么找出与当前块最相似的块,即匹配块,由匹配块与当前块的相对位置计算出运动位移,所得运动位移即为当前块的运动矢量。宏块大小为M×N,一般取16×16。搜索范围一般由最大偏移矢量来决定,设可能的最大偏移矢量为(dxmax,dymax),那么搜索范围为(M+2dxmax)×(N+2dymax)。图1运动估计根本原理2.1块匹配介绍按照一般的想法，运动估计应当首先将图像中静止背景和运动物体区分开来，然后对运动物体的实际位移进行估计。但块匹配方法却不是这样，它的思想是将图像划分为许多互不重叠的子块(例如16×16)，并认为子块内所有像素的位移量都相同。这意味着每个子块被视为运动物体。假设在图像序列中，t时刻对应于第k帧图像，t-τ时刻对应于第k-1帧图像。对于k帧中的一个子块，在k-1帧中寻找与其最相似的子块，这个过程称为寻找匹配块，并认为该匹配块在k-1帧中所处的位置就是k帧子块位移前的位置，这种位置的变化用运动矢量D来表示。将图像分割成M×N的小块，并假设块内象素作相同的运动，且只作平移运动。虽然实际上块内各点运动不一定相同，也不定只有平移运动，但当M×N较小时，上述假设可近似成立。这样做的目的只是为了简化运算。块匹配法对当前帧图像的每一块，在上一帧的一定范围内搜索最优匹配，并认为本块就是从上一帧最优匹配块位置处平移过来的。设可能的最大偏移矢量为(r,r)，那么搜索范围为(p+2r)×(q+2r)。图2示出了待匹配块与搜索区的几何位置关系。随着图像压缩编码技术的开展和对压缩效率愈来愈高的要求，在很多场合均要求运动矢量精确到亚象素元级。图2待匹配块与搜索区的几何位置关系块的大小受到两个矛盾的约束:块大时，块内各像素作平移运动的假设易被破坏，影响估计的精度；块小时，那么易受噪声影响，估计不够可靠，而且运算量增加，所需传输的附加信息也增加了。因此必须恰到好处地选择块的大小，以做到两者兼顾。目前的视频压缩标准，如H.26x和MPEG等，一般均以16×16大小的块作为块匹配单元，这是一个已为实践证明的较好的折衷结果。2.2离散余弦变换(DCT)离散余弦变换〔DCT〕是利用傅立叶变换的对称性,将图象描述为不同幅值和频率的正弦值之和的形式；是图象压缩JPEG压缩算法的根底和核心。是一个无信号损失的双向数学过程。通过DCT变换能去除视频信号的空间冗余。一个N×N矩阵的二维DCT定义如下：(1)正变换〔2.1〕(2)反变换〔2.2〕其中2.3块匹配的准那么运动估计算法中常用的匹配准那么有三种，即最小绝对差(拟D)、最小均方误差(MSE)和归一化互相关函数(NCCF)。分别定义如下：1、互相关函数(Cross-CorrelationFunction,简称CCF)〔2.3〕其中，。(2.3)式是计算当前帧中(X×Y)矩阵域与前帧相对应的(X+2P,Y+2P)矩阵区域互相关函数。2、均方误差函数(MSE)〔2.4〕其中，。(2.4)式是计算相邻帧相对应(X×Y)矩形区域的最小均方误差。这是一种非线性测量，能较好地跟踪图像的协方差模型。3、绝对平均误差函数(MAD)〔2.5〕(2.5)式是最简单的匹配函数，用它计算相邻帧的绝对平均误差。在(2.3)，(2.4)和(2.5)式中都在寻找(x,y)，得到一个最小失真矢量。(2.3)式是计算一个最大的自相关函数CCF(x,y)，而(2.4)和(2.5)式是计算最小的均方误差和绝对平均误差函数MSE(x,y)和MAD(x,y)。由于块匹配算法计算简单，能够实时处理，近几年获得广泛的应用。4、最大像素匹配统计(MaximumPixelsCounting，MPC)还有一种匹配准那么叫做最大匹配像素数(MPC：Matching-PixelCount)准那么。首先根据下式将当前块中的像素分成匹配像素和不匹配像素：〔2.6〕〔2.7〕，那么位置的像素为匹配像素，否那么为不匹配像素。2.4块匹配的搜索算法前两节的分析可以发现，在块匹配方法中最重要的两个问题是如何确定:(1)判别两个子块匹配的准那么:(2)计算量最小的搜索方法。对这两个问题的不同解决方案构成了不同的搜索算法。下面将对一些典型的块匹配快速搜索算法进行逐次介绍。2.4.1完全搜索法〔1〕算法思想全搜索法(FullSearchMethod,FS)也称为穷尽搜索法，是对(p+2r)×(q+2r)搜索范围内所有可能的候选位置计算MAD(i,j)值，从中找出最小MAD，其对应偏移量即为所求运动矢量。此算法虽计算量大，但最简单、可靠，找到的必为全局最优点。〔2〕FS算法描述FS算法描述如下:Step1:从原点出发，按顺时针方向由近及远，在逐个像素处计算MAD值，直到遍历搜索范围内所有的点.Step2:在所有点MAD中找到最小值，该点所在位置即对应最正确运动矢量。〔3〕FS算法的分析FS算法是最简单、最原始的块匹配算法，由于可靠，且能够得到全局最优的结果，通常是其它算法性能比拟的标准，但它的计算量确实很大，这就限制了在需要实时压缩场合的应用，所以有必要进一步研究其它快速算法。2.4.2二维对数法二维对数(Two-DimensionalLogarithmicTDL)搜索法由J.R.Jain和A.K.Jain提出，它开创了快速算法的先例，分多个阶段搜索，逐次减小搜索范围直到不能再小而结束。〔1〕TDL算法描述TDL算法的根本思想是从原点开始，以“十〞字形分布的五个点构成每次搜索的点群，通过快速搜索跟踪最小块误差MBD(MininumBlockDistortion)点(MAD值最小的点)，算法具体描述如下:Step1:从原点开始，选取一定的步长，在以十字形分布的五个点处进行块匹配计算并比拟。Step2:假设MBD点在边缘四个点处，那么以该点作为中心点，保持步长不变，重新搜索十字形分布的五个点；假设MBD点位于中心点，那么保持中心点位置不变，将步长减半，构成十字形点群，在五个点处计算。Step3:在中心及周围8个点处找出MBD点，假设步长为1，该点所在位置即对应最正确运动矢量，算法结束；否那么重复Step2。具体的一个搜索例子请参考图3。图中每个点上的数字说明了每个阶段搜索时计算的候选块的位置。图3TDL搜索过程〔2〕TDL算法的分析TDL算法搜索时，最大搜索点数为2+71og2r，假设发现新的十字形点群的中心点位于搜索区的边缘，那么步长也减半，后来有人提出应该在搜索的每个阶段都将步长减半,所有这些改动都是为了使算法搜索范围很快变小，提高收敛速度。TDL算法的前提是假设搜索区内只有一个谷点，如果搜索区内存在多个谷点时，该方法找到的可能是局部最小点。2.4.3三步搜索法三步搜索(ThreeStepSearch,TSS)法与二维对数法类似，是T.KOGA等人提出的，由于简单、健壮、性能良好的特点，为人们所重视。假设最大搜索长度为7，搜索精度取1个像素，那么步长为4,2,1，共需三步即可满足要求，因此而得名三步法。TSS算法描述TSS算法的根本思想是采用一种由粗到细的搜索模式，从原点开始，按一定步长取周围8个点构成每次搜索的点群，然后进行匹配计算，跟踪最小块误差MBD点算法具体描述如下:Step1:从原点开始，选取最大搜索长度的一半为步长，在周围距离步长的8个点处进行块匹配计算并比拟。Step2:将步长减半，中心点移到上一步的MBD点，重新在周围距离步长的8个点处进行块匹配计算并比拟。Step3:在中心及周围8个点处找出MBD点，假设步长为1，该点所在位置即对应最正确运动矢量，算法结束；否那么重复Step2。一个可能的搜索过程如图4所示，图中每个点上的数字说明了每个阶段搜索时计算的候选块的位置。图4三步搜索步骤〔2〕TSS算法的分析TSS算法搜索时，整个过程采用了统一的搜索模板(SearchPattern)，使得第一步的步长过大，容易引起误导，从而对小运动效率较低。最大搜索点数为1＋1og2r，当搜索范围大于7时，仅用3步是不够的，搜索步数的一般表达式为log2(p+1)。总体说来，三步法是一种较典型的快速搜索算法，所以被研究的较多，后来又相继有许多改良的新三步法出现，改良了它对小运动的估计性能。直角搜索法直角搜索法(OrthogonalSearchAlgorithm，OSA)是二维对数法和三步法的一种混合，由Purl提出[is]，为了找到最正确匹配点搜索过程由垂直和水平两个阶段交替进行，逐步减小步长。〔1〕OSA算法描述TSS算法实际上仍是TDL和TSS的模式，但更加简化，使每次只沿水平或垂直方向在三个点处进行块匹配计算，算法具体描述如下:Step1:初始化选取最大搜索长度的一半为步长，先在水平方向距离原点步长的位置处进行块匹配计算并比拟，然后将中心点移到MBD点。Step2:在垂直方向距离当前中心点步长处选取两点进行块匹配，并找到MBD点。Step3:假设步长为1，该点即为最优点，算法结束;否那么，将步长减半，改变中心点位置，重复Step1。图4是一个具体的搜索例子，图中每个点上的数字说明了每个阶段搜索时计算的候选块的位置。图5OSA搜索过程〔2〕OSA算法分析OSA算法每一步可以看成是由水平和垂直两个阶段组成，直到步长变为1时视为算法结束。由于每次计算的匹配点数太少，不能顾及各个方向，虽然速度较快，但得到局部最优的可能性进一步扩大。交叉法1990年，Ghanbari提出了交叉搜索算法(CrossSearchAlgorithm,CSA),它也是在TDL和TSS根底上为进一步减少计算量而开展起来的快速搜索法。〔1〕CSA描述CSA的根本思想是从原点开始，以“X〞字形分布的五个点构成每次搜索的点群，以TDL的搜索方法检测MBD点，仅在最后一步采用“十〞字形点群。算法具体描述如下:Step1:从原点开始，选取最大搜索长度的一半为步长，在以“X〞字形分布的五个点处进行块匹配计算并比拟，然后移动中心点。Step2:以上一步的MBD点为中心，步长减半，继续做“X〞字形的S点搜索。假设步长大于1，那么重复Step2;假设步长为1，那么进行Step3。Step3:最后一步根据MBD点的位置，分别做“+〞字形和“X〞字形搜索:假设上一步MBD点处于中心点、左下角或右上角，做“十〞字形搜索；假设上一步MBD点处于左上角或右下角，做“X〞字形搜索。由当前MBD点得到最正确运动矢量，算法结束。图6是CSA搜索的一个具体实例，图中每个点上的数字说明了每个阶段搜索时计算的候选块的位置，第三步箭头说明了两种不同的搜索模式。图6CSA搜索过程〔2〕CSA分析CSA的最大搜索点数为5+41og2r，搜索速度很快，但是运动补偿的效果不算太好。如图6右图所示，有些点根本就不可能成为候选点，或者说，在搜索区域的边界上有四分之一的点CSA没有考虑，因此它不适用于较复杂的运动。除此以外，还有菱形搜索法和四步搜索法，就不一一介绍了。3设计实现f2=imread('G:\壁纸\car1.BMP')f1=imread('G:\壁纸\car2.BMP')fp=0figure,imshow(f2),title('target')figure,imshow(f1),title('anchor')N=16;R=16height=256width=256fori=1:N:height-N+1forj=1:N:width-N+1MAD_min=256*N*Ndy=0;dx=0;fork=-R:1:Rforl=-R:1:Rifi+k<1MAD=256*N*Nelseifi+k>height-NMAD=256*N*Nelseifj+l<1MAD=256*N*Nelseifj+l>width-NMAD=256*N*NelseMAD=sum(sum(abs(double(f1(i:i+N-1,j:j+N-1))-double(f2(i+k:i+k+N-1,j+l:j+l+N-1)))))endifMAD<MAD_minMAD_min=MADdy=k;dx=l;end;end;end;fp(i:i+N-1,j:j+N-1)=f2(i+dy:i+dy+N-1,j+dx:j+dx+N-1)iblk=floor((i-1)/N+1);jblk=floor((j-1)/N+1);mvx(iblk,jblk)=dx;mvy(iblk,jblk)=dy;end;end;figure,imshow(uint8(fp)),title('predict');[X,Y]=meshgrid(N/2:N:256-N/2);Y=256-Y;figure,quiver(X,Y,mvx,mvy),title('motionvector');diff=abs(double(f1)-fp);figure,imshow(uint8(d