運動估計算法MATLAB課程設(shè)計

通信工程應(yīng)用技術(shù)綜合訓(xùn)練與實習(xí)課程設(shè)計課程設(shè)計任務(wù)書學(xué)生姓名： 專業(yè)班級： 指導(dǎo)教師： 工作單位： 題 目: 運動估計算法 初始條件： MATLAB軟件平臺要求完成的主要任務(wù):1. 設(shè)計任務(wù)設(shè)計視頻壓縮系統(tǒng)中的運動估計算法：全搜索法（FS: Full Search）和三步法(TSS: Three Step Search)，比較二種方法的搜索點和每幀的峰值信噪比（PSNR： peak signal to noise ratio）2. 設(shè)計要求 編制算法代碼；對視頻進行運動估計；計算PSNR時間安排： 答辯時間2013年1月24日。指導(dǎo)教師簽名： 年 月 日系主任簽名： 年 月 日目 錄摘 要11 運動估計算法概念21.1 運功估計算法基本思想21.2 運動估計算法實驗原理22 設(shè)計原理和方法32.1 三步法32.2 新三步法32.3 全搜索法42.4 峰值信噪比53 運動估計算法的MATLAB編程63.1全搜索法63.2三步法93.3全搜索法指標(biāo)113.4三步法指標(biāo)113.4仿真結(jié)果分析114 小結(jié)與體會11參考文獻11附錄12摘 要運動估計的基本思想是盡可能準(zhǔn)確地獲得序列圖像幀間的運動位移，即運動矢量。因為運動估計越準(zhǔn)確，預(yù)測補償?shù)膱D像質(zhì)量越高，補償?shù)臍埐罹驮叫。a償編碼所需位數(shù)越少，需要傳輸?shù)谋忍芈示驮叫　＠玫玫降倪\動矢量在幀間進行運動補償。補償殘差經(jīng)過變換、量化、編碼后與運動矢量一起經(jīng)過熵編碼，然后以比特流形式發(fā)送出去。在視頻編碼和處理系統(tǒng)中，運動估計和運動補償技術(shù)對降低視頻序列時間冗余度、提高編碼效率起著非常關(guān)鍵的作用。運動估計的準(zhǔn)確程度將直接決定視頻編碼器的編碼效率。關(guān)鍵詞：運動估計、運動補償技術(shù)、位移（運動）矢量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 accurate 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 transformation, 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 plays 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 運動估計算法概念視頻原始圖像中存在著大量的信息冗余，如時間冗余、空間冗余、信息熵冗余、譜間冗余、幾何結(jié)構(gòu)冗余、視覺冗余和知識冗余等等。運動估計是視頻壓縮編碼中的核心技術(shù)之一，采用運動估計和運動補償技術(shù)可以消除視頻信號的時間冗余以提高編碼效率。如何提高運動估計的效率，使運動估計算法的搜索過程更健壯、更快速、更高效成為目前研究的熱點。運動估計的基本思想是盡可能準(zhǔn)確地獲得序列圖像幀間的運動位移，即運動矢量。因為運動估計越準(zhǔn)確，預(yù)測補償?shù)膱D像質(zhì)量越高，補償?shù)臍埐罹驮叫?，補償編碼所需位數(shù)越少，需要傳輸?shù)谋忍芈示驮叫?。利用得到的運動矢量在幀間進行運動補償。補償殘差經(jīng)過變換、量化、編碼后與運動矢量一起經(jīng)過熵編碼，然后以比特流形式發(fā)送出去。運動估計算法多種多樣，大體上可以把它們分成四類：塊匹配法、遞歸估計法、貝葉斯估計法和光流法。其中塊匹配運動估計算法因其具有算法簡單、便于VLSI實現(xiàn)等優(yōu)點得到廣泛應(yīng)用。所以本文將重點介紹塊匹配運動估計算法，并對各種塊匹配算法在計算速度和估計精度上進行簡單比較。1.1 運功估計算法基本思想運動估計的基本思想是盡可能準(zhǔn)確地獲得序列圖像幀間的運動位移，即運動矢量。運動估計算法的目標(biāo)是效率和準(zhǔn)確性。由于在成象的場景中一般有多個物體作不同的運動，如果直接按照不同類型的運動將圖像分割成復(fù)雜的區(qū)域是比較困難的。最直接和不受約束的方法是在每個像素都指定運動矢量，這就是所謂基于像素表示法。這種表示法是對任何類型圖像都是適用的，但是它需要估計大量的未知量，并且它的解時常在物理上是不正確，除非在估計過程中施加適當(dāng)?shù)奈锢砑s束。這在具體實現(xiàn)時是不可能的，通常采用基于塊的物體運動表示法。因為運動估計越準(zhǔn)確，預(yù)測補償?shù)膱D像質(zhì)量越高，補償?shù)臍埐罹驮叫?，補償編碼所需位數(shù)越少，需要傳輸?shù)谋忍芈示驮叫　?.2 運動估計算法實驗原理在幀間預(yù)測編碼中，由于活動圖像鄰近幀中的景物存在著一定的相關(guān)性。因此，可將活動圖像分成若干塊或宏塊，并設(shè)法搜索出每個塊或宏塊在鄰近幀圖像中的位置，并得出兩者之間的空間位置的相對偏移量，得到的相對偏移量就是通常所指的運動矢量，得到運動矢量的過程被稱為運動估計。 運動矢量和經(jīng)過運動匹配后得到的預(yù)測誤差共同發(fā)送到解碼端，在解碼端按照運動矢量指明的位置，從已經(jīng)解碼的鄰近參考幀圖像中找到相應(yīng)的塊或宏塊，和預(yù)測誤差相加后就得到了塊或宏塊在當(dāng)前幀中的位置。 運動估計的準(zhǔn)確程度往往用補償圖像與原圖像比較的PSNR來衡量表示。2 設(shè)計原理和方法2.1 三步法三步法是應(yīng)用得相當(dāng)廣泛的一種次優(yōu)的運動估計搜索算法它的搜索區(qū)間一般為-7，7 即在候選區(qū)中與編碼塊相同坐標(biāo)位置處為原點，將參考塊在其上下左右距離為7的范圍內(nèi)按照一定規(guī)律移動移到一個位置就做匹配計算它總共進行了三步搜索在下一次搜索時步長減半以前一步搜索得到的最優(yōu)點為中心。下圖為三步法的搜索示意圖。三步法(TSS)搜索示意圖算法的中心思想是，采用一種由粗到細(xì)的搜索模式，從原點開始，按一定步長取周圍8個點構(gòu)成每次搜索的點群，然后進行匹配計算，利用上一步搜索得到的最小塊誤差MBD點作為當(dāng)前搜索的中心位置，每做一步，搜索的步長減1。步搜索算法搜索窗選取(-7，+7)，最多只需要做25個位置的匹配計算，相對于全搜索來比，大大減少了匹配運算的復(fù)雜度，而且數(shù)據(jù)讀取比較規(guī)則。2.2 新三步法TSS 假定運動矢量分布特點是在搜索窗口中均勻分布, 但事實證明運動矢量是偏置中心的，Renxiang Li等人在TSS的基礎(chǔ)上提出了一種增強運動矢量中心偏置搜索和減小補償誤差的新三步法。新三步法(NTSS)搜索示意圖NTSS 是對TSS 的一個改進，對運動量比較小的視頻序列如可視電話序列有比較好的性能。對于絕大多數(shù)的視頻序列，運動矢量的分布都是在中心位置上的概率最大，隨著與中心位置的距離的增大，概率會急劇地下降，這也就是前面所說的運動矢量的中心偏移特性。運動量比較小的視頻序列的這一特性會更加明顯。NTSS算法在最好的情況下只需要做17個點的匹配，在最壞的情況下需要做33個點的匹配，由于運動矢量中心偏置在現(xiàn)實視頻序列中是普遍存在的，在通常情況下，NTSS算法需要做33點匹配的概率比較小，因此，在低速率視頻應(yīng)用中，如視頻電話或視頻會議中，NTSS算法的優(yōu)點可以得到較好的發(fā)揮。2.3 全搜索法全搜索法（Full Search Method, FS）也稱為窮盡搜索法，是對（M 2dx ）（N 2dy ）搜索范圍內(nèi)所有可能的候選位置計算MAD (i,j)值，從中找出最小MAD，其對應(yīng)偏移量即為所求運動矢量。此算法雖計算量大，但最簡單、可靠，找到的必為全局最優(yōu)點。FS算法描述如下：從原點出發(fā)，按順時針螺旋方向由近及遠(yuǎn)，在逐個像素處計算MAD值，直到遍歷搜索范圍內(nèi)聽有的點，然后在計算的所有點的MAD中找到最小值，該點所在位置即對應(yīng)最佳運動矢量。但是正因為它是窮盡搜索因此會產(chǎn)生巨大的計算量如 7, 7的搜索區(qū)間每個宏塊16*16 需計算225個MAD值，這就直接制約了編碼的實時實現(xiàn)?？焖偎惴ū举|(zhì)上是一種窮盡搜索法其計算量仍是相當(dāng)巨大的。全搜索算法是簡單也是效果最好的一種匹配算法，通過的全搜索匹配得到的結(jié)果是全局最優(yōu)的，但由于計算量很大，我們在編解碼中往往不采用這種方法，而只把他作為與其他算法的一種比較。開始獲取圖像分成N*N大小的塊計算MSE進行塊匹配計算MAD比較結(jié)束圖2-1程序流程圖2.4 峰值信噪比峰值信噪比（經(jīng)?？s寫為PSNR）是一個表示信號最大可能功率和影響它的表示精度的破壞性噪聲功率的比值的工程術(shù)語。由于許多信號都有非常寬的動態(tài)范圍，峰值信噪比常用對數(shù)分貝單位來表示。PSNR是“Peak Signal to Noise Ratio”的縮寫。peak的中文意思是頂點。而ratio的意思是比率或比列的。整個意思就是到達噪音比率的頂點信號，psnr是一般是用于最大值信號和背景噪音之間的一個工程項目。通常在經(jīng)過影像壓縮之后，輸出的影像通常都會有某種程度與原始影像不一樣。為了衡量經(jīng)過處理后的影像品質(zhì)，我們通常會參考PSNR 值來認(rèn)定某個處理程序夠不夠令人滿意。PSNR計算公式如下：Peak 就是指8 bits 表示法的最大值255。MSE 指 Mean Square Error（均方誤差，各值相差的n次方和的平均值的n次平方根(這幾個字應(yīng)該沒有)），I(角標(biāo)n)指原始影像第n個pixel 值，P(角標(biāo)n)指經(jīng)處理后的影像第n個pixel 值。PSNR 的單位為dB。所以PSNR值越大，就代表失真越少。PSNR 是最普遍，最廣泛使用的評鑒畫質(zhì)的客觀量測法，不過許多實驗結(jié)果都顯示，PSNR 的分?jǐn)?shù)無法和人眼看到的視覺品質(zhì)完全一致，有可能 PSNR 較高者看起來反而比 PSNR 較低者差。這是因為人眼的視覺對于誤差的敏感度并不是絕對的，其感知結(jié)果會受到許多因素的影響而產(chǎn)生變化（例如：人眼對空間頻率較低的對比差異敏感度較高，人眼對亮度對比差異的敏感度較色度高，人眼對一個區(qū)域的感知結(jié)果會受到其周圍鄰近區(qū)域的影響）。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 恢復(fù)后的第二幀圖像圖3-1-6 第一幀圖像的運動矢量圖圖3-1-7 第二幀圖像的運動矢量圖3.2三步法 圖3-2-1 the first frame 圖3-2-2 the second frame 圖3-2-3 幀間差值 圖3-2-4 DFD圖3-2-5 恢復(fù)后的第二幀圖像圖3-2-6 第一幀圖像的運動矢量圖圖3-2-7 第二幀圖像的運動矢量圖3.3全搜索法指標(biāo)全搜索算法：耗時長為Elapsed time is 1. seconds.3.4三步法指標(biāo)三步法：耗時長為Elapsed time is 0. seconds.3.4仿真結(jié)果分析根據(jù)峰值信噪比可知全搜索算法的效果比三步法的效果要好；但全搜索所消耗的時間要比三步法所消耗的時間長，其運行較慢，因為全搜索算法的搜索點數(shù)要比三步法的搜索點數(shù)多很多。 4 小結(jié)與體會 通過本次對全搜索塊匹配算法課程設(shè)計，讓我更深刻的理解了多媒體信息的獲取、搜索和保存的基本方法?；趬K的運動估計，是先將視頻圖像分成一個個規(guī)則的圖像塊，然后對每個圖像塊估計運動矢量?；趬K的運動估計和運動補償已經(jīng)廣泛應(yīng)用于各種視頻壓縮編碼標(biāo)準(zhǔn)。因此，本次課程設(shè)計是非常有作用的，讓我們清醒的認(rèn)識到自己還有好大的差距，它讓我了解了多媒體處理的基本原理。課程設(shè)計中會遇到很多很多的困難，但在同學(xué)的幫助和自己思考下最終還是沒有解決問題。所以，本次課程設(shè)計我個人的收獲不大，但是讓我知道自己還在圖像編程方面有很多欠缺的地方，不僅對均方誤差(MSE)、絕對值誤差(MAD)更深的理解，并且采用全搜索塊匹配方法得到所有宏塊的運動矢量，把c語言的知識又重新的過了一遍，同時還理解了些多媒體信息處理的基本原理。參考文獻1羅軍輝. MATLAB7.0在圖像處理中的應(yīng)用. 機械出版社.20052王家文. MATLAB7.0編程基礎(chǔ).機械出版社.20053鐘麟. MATLAB仿真技術(shù)與應(yīng)用教程.國防工業(yè)出版社.20044岡薩雷斯. 數(shù)字圖像處理(MATLAB版).電子工業(yè)出版社.2005附錄FS：clear all;I1=imread(claire1.bmp); %read the first frameI2=imread(claire2.bmp); %read the 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: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/blocksize; A=; %為了找到最小的均方誤差，A用于設(shè)定一個很大的初值 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*blocksize+1:(x+1)*blocksize;% tempy=y*blocksize+1:(y+1)*blocksize; for p=-dm:dm for q=-dm:dm %（p,q）表示x,y對應(yīng)子塊在前一幀所的搜索位置 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; yrecord(x+1,y+1)=q; end end end A=; 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; 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(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=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*blocksize+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; end; 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 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); 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=; % Eij=0; xrecord=ones(16,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+blocksize)-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+blocksize,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(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=; 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; 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+blocksize