數(shù)字圖像處理實(shí)驗(yàn)報(bào)告-基于Matlab

1、-華東師大學(xué)電子工程系2021.6實(shí)驗(yàn)1：圖像灰度級(jí)修正【實(shí)驗(yàn)?zāi)康摹空莆粘Ｓ玫膱D像灰度級(jí)修正方法灰度變換法和直方圖均衡化，加深對(duì)直方圖的理解。觀察圖像的增強(qiáng)效果，對(duì)灰度級(jí)修正前后的圖像加以比較。【實(shí)驗(yàn)容】1編程實(shí)現(xiàn)圖像的灰度變換，改變圖像的輸入、輸出映射參數(shù)圍線性拉伸和反比；2修改參數(shù)gamma值大于、小于、等于1，觀察處理結(jié)果；3對(duì)圖像直方圖作均衡化處理，顯示均衡前后的圖像及其直方圖?！緦?shí)驗(yàn)代碼】original=imread('lena.bmp');linstr=imadjust(original,0.3 0.7,0 1);%線性拉伸opposite=imadjust(or

2、iginal,0 1,1 0); %反比above=imadjust(original,0 1,0 1,2); %gamma>1equal=imadjust(original,0 1,0 1,1); %gamma=1below=imadjust(original,0 1,0 1,0.5); %gamma<1subplot(3,3,1);imshow(original);title('原圖像');subplot(3,3,2);imshow(linstr);title('線性拉伸');subplot(3,3,3);imshow(opposite);tit

3、le('反比');subplot(3,3,4);imshow(above);title('gamma>1');subplot(3,3,5);imshow(equal);title('gamma=1');subplot(3,3,6);imshow(below);title('gamma<1');subplot(3,3,7);imhist(original);title('原圖像直方圖');histequal=histeq(original);%對(duì)圖像均衡化subplot(3,3,8);imshow(his

4、tequal);title('均衡后的圖像');subplot(3,3,9);imhist(histequal);title('均衡圖像的直方圖');a*is(0 256 0 2000);【輸出圖像】【實(shí)驗(yàn)思考】根據(jù)以以下圖片以及實(shí)驗(yàn)結(jié)果可知gamma>1時(shí)圖像整體變暗，灰度級(jí)整體變??；gamma<1時(shí)圖像整體變亮，灰度級(jí)整體變??；而gamma=1時(shí)，圖像維持不變。實(shí)驗(yàn)2：圖像的平滑濾波【實(shí)驗(yàn)?zāi)康摹科交哪康氖菧p少噪聲對(duì)圖像的影響。掌握線性濾波和中值濾波兩種最典型、最常用的圖像平滑方法，對(duì)輸出結(jié)果加以比較、加深理解?！緦?shí)驗(yàn)容】1編寫(xiě)并調(diào)試窗口為3&#

5、215;3、5×5的平滑濾波函數(shù)；如1 1 1;1 1 1 ;1 1 1/9、1 2 1;2 4 2;1 2 1/16等2編寫(xiě)并調(diào)試窗口為3×3、5×5的中值濾波函數(shù)。3比較均值濾波和中值濾波的優(yōu)缺點(diǎn)，分析窗口尺寸對(duì)濾波結(jié)果的影響。附：可供參考的Matlab函數(shù)有imnoise、imfilter、medfilt2【實(shí)驗(yàn)代碼】function fliterI = imread('lena.bmp'); %原始圖像讀取J = imnoise(I,'salt & pepper',0.02); %含噪圖像加椒鹽噪聲subplot(2

6、,3,1);imshow(J);title('含噪圖像'); Newbuf1=AverageFilter(J,256,256,3);%3×3標(biāo)準(zhǔn)平均，調(diào)用均值濾波函數(shù)subplot(2,3,2);imshow(Newbuf1);title('3×3標(biāo)準(zhǔn)平均');Newbuf2=AverageFilter(J,256,256,5);%5×5標(biāo)準(zhǔn)平均，調(diào)用均值濾波函數(shù)subplot(2,3,3);imshow(Newbuf2);title('5×5標(biāo)準(zhǔn)平均');W=1 2 1;2 4 2;1 2 1/16; %

7、設(shè)置加權(quán)平均掩膜Newbuf3=WeighFilter(J,W,256,256,3);%3×3加權(quán)平均，調(diào)用加權(quán)平均函數(shù)subplot(2,3,4);imshow(Newbuf3);title('3×3加權(quán)平均');Newbuf4=MedianFilter(J,256,256,3);%3×3中值濾波，調(diào)用中值濾波函數(shù)subplot(2,3,5);imshow(Newbuf4);title('3×3中值濾波');Newbuf5=MedianFilter(J,256,256,5);%5×5中值濾波，調(diào)用中值濾波函數(shù)s

8、ubplot(2,3,6);imshow(Newbuf5);title('5×5中值濾波');%標(biāo)準(zhǔn)平均濾波函數(shù)function Newbuf=AverageFilter(Oldbuf,M,N,m)%Newbuf 濾波后圖像矩陣%Oldbuf 含噪圖像矩陣%M、N 含噪圖像像素矩陣行、列%m 均值濾波窗口大小f=zeros(M+m-1,N+m-1); %將原圖像像素復(fù)制到f矩陣上，空出(m-1)/2大小的邊界f(m-1)/2+1:M+(m-1)/2,(m-1)/2+1:N+(m-1)/2)=Oldbuf(1:M,1:N);%將與邊界相鄰的(m-1)/2行或列的像素值復(fù)

9、制到邊界，以填充邊界f(m-1)/2+1:M+(m-1)/2,1:(m-1)/2)=Oldbuf( : ,1:(m-1)/2);f(m-1)/2+1:M+(m-1)/2,N+(m-1)/2:N+m-1)=Oldbuf( : ,N-(m-1)/2:N);f(1:(m-1)/2,(m-1)/2+1:N+(m-1)/2)=Oldbuf(1:(m-1)/2, : );f(M+(m-1)/2:M+m-1,(m-1)/2+1:N+(m-1)/2)=Oldbuf(M-(m-1)/2:M, : );g=zeros(M+m-1,N+m-1);Im=zeros(M,N);%根據(jù)公式計(jì)算出處理后g(*,y)的像素值

10、for *=(m-1)/2+1:M+(m-1)/2for y=(m-1)/2+1:N+(m-1)/2for s=-(m-1)/2:(m-1)/2for t=-(m-1)/2:(m-1)/2 g(*,y)=g(*,y)+f(*+s,y+t)*1/(m*m);endendendendIm(1:M,1:N)=g(m-1)/2+1:M+(m-1)/2,(m-1)/2+1:N+(m-1)/2);%將double型轉(zhuǎn)換為uint8型才可以用imshow正常顯示Newbuf=uint8(Im);%加權(quán)平均濾波函數(shù)function Newbuf=WeighFilter(Oldbuf,W,M,N,m)%Newb

11、uf 濾波后圖像矩陣%Oldbuf 含噪圖像矩陣%W 掩模%M、N 含噪圖像像素矩陣行、列%m 掩模模板窗口大小f=zeros(M+m-1,N+m-1);%將原圖像像素復(fù)制到f矩陣上，空出(m-1)/2大小的邊界f(m-1)/2+1:M+(m-1)/2,(m-1)/2+1:N+(m-1)/2)=Oldbuf(1:M,1:N);%將與邊界相鄰的(m-1)/2行或列的像素值復(fù)制到邊界，以填充邊界f(m-1)/2+1:M+(m-1)/2,1:(m-1)/2)=Oldbuf( : ,1:(m-1)/2);f(m-1)/2+1:M+(m-1)/2,N+(m-1)/2:N+m-1)=Oldbuf( : ,

12、N-(m-1)/2:N);f(1:(m-1)/2,(m-1)/2+1:N+(m-1)/2)=Oldbuf(1:(m-1)/2, : );f(M+(m-1)/2:M+m-1,(m-1)/2+1:N+(m-1)/2)=Oldbuf(M-(m-1)/2:M, : );g=zeros(M+m-1,N+m-1);Im=zeros(M,N);%根據(jù)公式計(jì)算出處理后g(*,y)的像素值for *=(m-1)/2+1:M+(m-1)/2for y=(m-1)/2+1:N+(m-1)/2for s=-(m-1)/2:(m-1)/2for t=-(m-1)/2:(m-1)/2 g(*,y)=g(*,y)+W(s+

13、(m+1)/2,t+(m+1)/2)*f(*+s,y+t);endendendendIm(1:M,1:N)=g(m-1)/2+1:M+(m-1)/2,(m-1)/2+1:N+(m-1)/2);Newbuf=uint8(Im);%中值濾波函數(shù)function Newbuf=MedianFilter(Oldbuf,M,N,m)%Newbuf 濾波后圖像矩陣%Oldbuf 含噪圖像矩陣%M、N 含噪圖像矩陣像素行、列%m 中值濾波窗口大小f=zeros(M+m-1,N+m-1);%將原圖像像素復(fù)制到f矩陣上，空出(m-1)/2大小的邊界f(m-1)/2+1:M+(m-1)/2,(m-1)/2+1:N

14、+(m-1)/2)=Oldbuf(1:M,1:N);%將與邊界相鄰的(m-1)/2行或列的像素值復(fù)制到邊界，以填充邊界f(m-1)/2+1:M+(m-1)/2,1:(m-1)/2)=Oldbuf( : ,1:(m-1)/2);f(m-1)/2+1:M+(m-1)/2,N+(m-1)/2:N+m-1)=Oldbuf( : ,N-(m-1)/2:N);f(1:(m-1)/2,(m-1)/2+1:N+(m-1)/2)=Oldbuf(1:(m-1)/2, : );f(M+(m-1)/2:M+m-1,(m-1)/2+1:N+(m-1)/2)=Oldbuf(M-(m-1)/2:M, : );g=zeros

15、(M+m-1,N+m-1);Im=zeros(M,N);for *=(m-1)/2+1:M+(m-1)/2for y=(m-1)/2+1:N+(m-1)/2 j=1;for s=-(m-1)/2:(m-1)/2for t=-(m-1)/2:(m-1)/2 a(j)=f(*+s,y+t);%將窗口里的二維元素變成一維元素 j=j+1;endend g(*,y)=SeekMid(a,m);endendIm(1:M,1:N)=g(m-1)/2+1:M+(m-1)/2,(m-1)/2+1:N+(m-1)/2);Newbuf=uint8(Im);%找出中值函數(shù)function mid=SeekMid(w

16、inbuf,m)%mid 排序后的中值%winbuf 待排序窗口%m 窗口大小%采用冒泡排序方法將窗口像素值從小到大排列，返回中間像素值for i=1:m*m-1for k=1:m*m-iif winbuf(k)>winbuf(k+1) temp=winbuf(k); winbuf(k)=winbuf(k+1); winbuf(k+1)=temp;endendendmid=winbuf(m*m+1)/2);【輸出圖像】【實(shí)驗(yàn)思考】1. 比較均值濾波和中值濾波的優(yōu)缺點(diǎn)均值濾波可以減小圖像灰度級(jí)的“sharp變化，可以降低噪聲，但是降噪的同時(shí)也使邊緣局部變得模糊，還可以平滑偽輪廓，去除圖像中

17、的不相關(guān)的小于掩模尺寸的細(xì)節(jié)。中值濾波器的主要功能是使具有不同灰度的點(diǎn)看起來(lái)更接近它的相鄰點(diǎn)，去除那些相對(duì)于其鄰域像素更亮或更暗、且區(qū)域小于n2/2的孤立像素集。中值濾波對(duì)降低*些類(lèi)型的隨機(jī)噪聲性能優(yōu)異，模糊程度低。在處理椒鹽噪聲時(shí)，均值濾波使圖像變得模糊，并且噪聲去除性能很差，而中值濾波的效果卻很好。顯然，中值濾波比均值濾波更適合去除椒鹽噪聲。2. 分析窗口尺寸對(duì)濾波結(jié)果的影響窗口尺寸越大，圖像越模糊，圖像邊緣和與掩膜大小接近的細(xì)節(jié)受到的影響也越大實(shí)驗(yàn)3：圖像的銳化處理【實(shí)驗(yàn)?zāi)康摹夸J化的目的是加強(qiáng)圖像的邊界和細(xì)節(jié)，熟悉Robert、Sobel和Laplace算子進(jìn)展檢測(cè)，使圖像特征如邊緣、輪

18、廓等進(jìn)一步增強(qiáng)并突出?！緦?shí)驗(yàn)容】1編寫(xiě)Robert算子濾波函數(shù)；2編寫(xiě)Sobel算子濾波函數(shù)；3編寫(xiě)Laplace算子濾波函數(shù)；4編寫(xiě)限幅和標(biāo)定函數(shù)，給出增強(qiáng)后的圖像?！緦?shí)驗(yàn)代碼】function E*3I=imread('rice.bmp');subplot(2,4,1);imshow(I);title('原始圖像');rob=RobertFilter(I);subplot(2,4,2);imshow(rob);title('Robert算子濾波結(jié)果');R1=I+rob;la1=LimitAmplitude(R1);subplot(2,4,6

19、);imshow(la1);title('Robert算子增強(qiáng)結(jié)果');a2=-1 -2 -1;0 0 0;1 2 1;b2=-1 0 1;-2 0 2;-1 0 1; sob=SobelFilter(I,a2,b2);subplot(2,4,3);imshow(sob);title('Sobel算子濾波結(jié)果');R2=I+sob;la2=LimitAmplitude(R2);subplot(2,4,7);imshow(la2);title('Sobel算子增強(qiáng)結(jié)果');% s=0 1 0;1 -4 1;0 1 0;s=1 1 1;1 -8 1;

20、1 1 1;lap=LapFilter(I,s);cal=Calibration(lap);subplot(2,4,4);imshow(cal);title('Laplace算子濾波結(jié)果');lap=uint8(lap);lapr=I-lap;lapr3=LimitAmplitude(lapr);subplot(2,4,8);imshow(lapr3);title('Laplace算子增強(qiáng)結(jié)果');%Robert算子濾波function rob=RobertFilter(F)a1=-1 0;0 1;b1=0 -1;1 0;%Robert算子模板M,N=size(

21、F);f=zeros(M+1,N+1);f(1:M,1:N)=F(1:M,1:N);f(1:M,N+1:N+1)=F( : ,N:N);f(M+1:M+1,1:N)=F(M:M, : );%邊界填充g=zeros(M+1,N+1);for *=1:Mfor y=1:N mod=f(*,y) f(*,y+1);f(*+1,y) f(*+1,y+1); gs*=a1.*mod; gsy=b1.*mod; g(*,y)=abs(sum(gs*(:)+abs(sum(gsy(:);endendIm=zeros(M,N);Im(1:M,1:N)=g(1:M,1:N);rob=uint8(Im);%Sob

22、el算子濾波function sob=SobelFilter(F,s*,sy)%s*,sy為Sobel算子模板M,N=size(F);m,n=size(s*);f=zeros(M+m-1,N+n-1);f(m-1)/2+1:M+(m-1)/2,(n-1)/2+1:N+(n-1)/2)=F(1:M,1:N);f(m-1)/2+1:M+(m-1)/2,1:(n-1)/2)=F( : ,1:(n-1)/2);f(m-1)/2+1:M+(m-1)/2,N+(n-1)/2:N+m-1)=F( : ,N-(n-1)/2:N);f(1:(m-1)/2,(n-1)/2+1:N+(n-1)/2)=F(1:(m-

23、1)/2, : );f(M+(m-1)/2:M+m-1,(n-1)/2+1:N+(n-1)/2)=F(M-(m-1)/2:M, : );%邊界填充g=zeros(M+m-1,N+n-1);for *=(m-1)/2+1:M+(m-1)/2for y=(n-1)/2+1:N+(n-1)/2 mod=f(*-1,y-1) f(*-1,y) f(*-1,y+1); f(*,y-1) f(*,y) f(*,y+1);f(*+1,y-1) f(*+1,y) f(*+1,y+1); gs*=s*.*mod; gsy=sy.*mod; g(*,y)=abs(sum(gs*(:)+abs(sum(gsy(:)

24、;endendIm=zeros(M,N);Im(1:M,1:N)=g(m-1)/2+1:M+(m-1)/2,(n-1)/2+1:N+(n-1)/2);sob=uint8(Im);%Laplace算子濾波function lap=LapFilter(F,S)M,N=size(F);m,n=size(S);f=zeros(M+m-1,N+n-1);f(m-1)/2+1:M+(m-1)/2,(n-1)/2+1:N+(n-1)/2)=F(1:M,1:N);f(m-1)/2+1:M+(m-1)/2,1:(n-1)/2)=F( : ,1:(n-1)/2);f(m-1)/2+1:M+(m-1)/2,N+(n

25、-1)/2:N+m-1)=F( : ,N-(n-1)/2:N);f(1:(m-1)/2,(n-1)/2+1:N+(n-1)/2)=F(1:(m-1)/2, : );f(M+(m-1)/2:M+m-1,(n-1)/2+1:N+(n-1)/2)=F(M-(m-1)/2:M, : );g=zeros(M+m-1,N+n-1);for *=(m-1)/2+1:M+(m-1)/2for y=(n-1)/2+1:N+(n-1)/2 mod=f(*-1,y-1) f(*-1,y) f(*-1,y+1); f(*,y-1) f(*,y) f(*,y+1);f(*+1,y-1) f(*+1,y) f(*+1,y

26、+1); gs=S.*mod; g(*,y)=sum(gs(:);endendIm=zeros(M,N);Im(1:M,1:N)=g(m-1)/2+1:M+(m-1)/2,(n-1)/2+1:N+(n-1)/2);lap=Im;%限幅函數(shù)function la=LimitAmplitude(F)f=uint8(F);M,N=size(f);for *=1:Mfor y=1:Nif f(*,y)>=255; f(*,y)=255;elseif f(*,y)<=0 f(*,y)=0;else f(*,y)=f(*,y);%將灰度值限定在0到255之間endendendendla=f;%

27、標(biāo)定函數(shù)function cal=Calibration(F)F=double(F);M,N=size(F);m1=min(min(F);for *=1:Mfor y=1:N fm(*,y)=F(*,y)-m1;endendm2=ma*(ma*(fm);Fm=double(fm);for *=1:Mfor y=1:N fs(*,y)=255*(Fm(*,y)/m2);endendcal=uint8(fs);【輸出圖像】實(shí)驗(yàn)4：圖像的統(tǒng)計(jì)特性【實(shí)驗(yàn)?zāi)康摹坑^察序列圖像幀、幀間差值信號(hào)的分布曲線，理解圖像在空間域和頻率域上的統(tǒng)計(jì)特性及其在壓縮中的重要性?！緦?shí)驗(yàn)容】1編寫(xiě)幀統(tǒng)計(jì)函數(shù)，計(jì)算差值圖像同一行

28、差值、同一列差值，觀察統(tǒng)計(jì)分布曲線；2編寫(xiě)幀間統(tǒng)計(jì)函數(shù)，計(jì)算差值圖像相鄰幀的差值，觀察統(tǒng)計(jì)分布曲線cla0/1或girl0/1。附：可供參考的Matlab函數(shù)有sum、cat、plot【實(shí)驗(yàn)代碼】function E*4oldbuf=imread('rice.bmp');I1=imread('CLA1.bmp');I2=imread('CLA2.bmp');newbuf1=Intrah(oldbuf,1);%幀水平差值統(tǒng)計(jì)特性newbuf2=Intrah(oldbuf,0);%幀垂直差值統(tǒng)計(jì)特性newbuf3=Inter(I1,I2);幀間統(tǒng)計(jì)特

29、性subplot(2,3,1);imshow(oldbuf);title('原始圖像');subplot(2,3,2);draw(newbuf1);title('水平差值統(tǒng)計(jì)特性');subplot(2,3,3);draw(newbuf2);title('垂直差值統(tǒng)計(jì)特性');subplot(2,3,4);imshow(I1);title('CLA1');subplot(2,3,5);imshow(I2);title('CLA2');subplot(2,3,6);draw(newbuf3);title('幀

30、間統(tǒng)計(jì)特性');function newbuf=Intrah(oldbuf,pop2) %幀統(tǒng)計(jì)函數(shù)oldbuf=double(oldbuf);M,N=size(oldbuf); %防止溢出將數(shù)據(jù)類(lèi)型從uint8型轉(zhuǎn)換為double型newbuf=zeros(1,511);ifpop2=1for i=1:Mfor j=1:N-1 dH=oldbuf(i,j)-oldbuf(i,j+1);%幀水平灰度差值 newbuf(dH+256)=newbuf(dH+256)+1;endendelsefor i=1:M-1for j=1:N dV=oldbuf(i,j)-oldbuf(i+1,j);

31、 newbuf(dV+256)=newbuf(dV+256)+1;endendendfunction newbuf=Inter(oldbuf,oldbuf1) %幀間統(tǒng)計(jì)函數(shù)oldbuf=double(oldbuf);oldbuf1=double(oldbuf1);M,N=size(oldbuf); newbuf=zeros(1,511);for i=1:Mfor j=1:N dt=oldbuf(i,j)-oldbuf1(i,j);%計(jì)算幀間差值 newbuf(dt+256)=newbuf(dt+256)+1;endendfunction draw(D)D=D/sum(D);*=-255:25

32、5;plot(*,D);a*is(-100 100 0 0.5);%為了顯示效果好縮小坐標(biāo)軸圍【輸出圖像】實(shí)驗(yàn)6：方塊編碼【實(shí)驗(yàn)?zāi)康摹空莆辗綁K編碼的根本方法及壓縮性能?！緦?shí)驗(yàn)容】1編程實(shí)現(xiàn)子塊為n×n的方塊編碼算法；2分別取n4和8的方塊尺寸進(jìn)展實(shí)驗(yàn)，計(jì)算重建圖像的PSNR和壓縮比?！緦?shí)驗(yàn)代碼】1.主程序I=imread('lena.bmp');M,N=size(I);subplot(1,3,1);imshow(I);title('原圖像');%顯示原圖像I=double(I);newbuf1=BtcCode(I,4);PSNR1,Cr1=Analyz

33、e(I,newbuf1,M,N,4);subplot(1,3,2);imshow(uint8(newbuf1);title('4*4BTC重建圖像,PSNR=',num2str(PSNR1),'壓縮比=',num2str(Cr1);newbuf2=BtcCode(I,8);PSNR2,Cr2=Analyze(I,newbuf2,M,N,8);subplot(1,3,3);imshow(uint8(newbuf2);title('8*8BTC重建圖像,PSNR=',num2str(PSNR2),'壓縮比=',num2str(Cr2)

34、;function outbuf=BtcBlock(inbuf,n)%btc 方塊編碼算法函數(shù)%inbuf 子塊數(shù)組%n 方塊尺寸%對(duì)每個(gè)子塊的圖像數(shù)據(jù)分別計(jì)算*t、a0、a1值，再用分辨率分量%(a0,a1)替代方塊原來(lái)的數(shù)據(jù)，最后放入方塊圖像數(shù)組中并返回該數(shù)組inbuf=double(inbuf);temp=0;%總的像素值temp0=0;%小于閥值的總像素temp1=0;%大于閥值的總像素q=0;%大于閥值的像素的個(gè)數(shù)m=n*n;for i=1:nfor j=1:n temp=temp+inbuf(i,j);endend*t=temp/m;%平均像素值即閥值for i=1:nfor j=

35、1:n if inbuf(i,j)<*t temp0=temp0+inbuf(i,j);%得出小于閥值的總像素else temp1=temp1+inbuf(i,j);%得出大于閥值的總像素 q=q+1;%大于閥值的像素個(gè)數(shù)endendendif q=m a0=uint8(temp0/(m-q);%得出小于閥值的像素值endif q=0 a1=uint8(temp1/q); %得出大于閥值的像素值endfor i=1:nfor j=1:nif inbuf(i,j)<*t outbuf(i,j)=a0;else outbuf(i,j)=a1;endendendfunction newb

36、uf=BtcCode(oldbuf,n)%調(diào)用方塊編碼算法函數(shù)，輸出編碼后的圖像M,N=size(oldbuf);row_num=M/n;%子塊行數(shù)col_num=N/n;%子塊列數(shù)row_start=(0:row_num)*n+1;%子塊起始行row_end=(1:row_num)*n;%子塊終止行col_start=(0:col_num-1)*n+1;%子塊起始列col_end=(1:row_num)*n;%子塊終止列for i=1:row_numfor j=1:col_num f=oldbuf(row_start(i):row_end(i),col_start(j):col_end(j)

37、;%此式太長(zhǎng)為方便書(shū)寫(xiě)定義f oldbuf(row_start(i):row_end(i),col_start(j):col_end(j)=BtcBlock(f,n);%將原圖像分成一個(gè)個(gè)子塊，在原圖像里一個(gè)個(gè)對(duì)這些子塊進(jìn)展編碼，編碼后的結(jié)果保存原圖像里endendnewbuf=oldbuf;%編碼后的圖像4. Analyze.mfunction PSNR,Cr=Analyze(I1,I2,M,N,n)%計(jì)算重建圖像的PSNR和壓縮比m=n*n;mse=sum(sum(I1-I2).2)/(M*N);PSNR=10*log10(2552)/mse);Cr=8/(1+2*8/m);end【輸出圖

38、像】實(shí)驗(yàn)7：JPEG壓縮編碼【實(shí)驗(yàn)?zāi)康摹空莆課×n塊的DCT圖像變換及頻譜特點(diǎn)。熟悉JPEG根本系統(tǒng)的圖像編解碼方法?！緦?shí)驗(yàn)容】 1編程實(shí)現(xiàn)n×n塊DCT變換的圖像頻譜顯示，塊DCT系數(shù)按照Z(yǔ)ig-Zag掃描并取局部進(jìn)展圖像重建，計(jì)算圖像的均方根誤差RMSE，顯示誤差圖像和誤差直方圖。 2對(duì)8×8塊的DCT系數(shù)，采用JPEG默認(rèn)的量化矩陣進(jìn)展量化和反量化，計(jì)算原圖像與重建圖像之間的均方根誤差RMSE、并顯示誤差圖像?！緦?shí)驗(yàn)代碼】1.主程序F=imread('lena.bmp');subplot(231);imshow(F);title('&

39、#212;­Í¼Ïñ');%顯示原圖像F=double(F);F=F-128;%將原圖像減小一半便于處理%計(jì)算原圖像的8×8塊的DCT系數(shù)，并轉(zhuǎn)換為可視頻譜圖以便觀察dctfre=DctBlock(F,8);subplot(232);imshow(log(abs(dctfre)*5+1),);title('8*8DCT頻譜顯示');%表示將原圖像的最大最小值之間的圍整體映射到0255之間，即做限幅DCTch=10;n=8;I,e,rmse1=ZigIDCT(F,dctfre,DCTch,n);subplot(2

40、33);imshow(uint8(I);title('取',num2str(DCTch),'個(gè)DCT系數(shù)時(shí)的重建圖像');subplot(234);imhist(uint8(abs(e);title('差值直方圖，RMSE=',num2str(rmse1);scale=4;newbuf,err,rmse2=QuanIQuan(F,dctfre,n,scale);subplot(235);imshow(uint8(newbuf);title('scale為',num2str(scale),'時(shí)的重建圖像');subp

41、lot(236);imshow(uint8(abs(err),);title('量化誤差圖像，RMSE=',num2str(rmse2);2. ZigIDCT.mfunction I,e,rmse1=ZigIDCT(oldbuf,dctfre,DCTch,n)%oldbuf:原始圖像數(shù)據(jù)%dctfre:DCT系數(shù)矩陣%DCTch:每個(gè)分塊中需要保存的DCT系數(shù)個(gè)數(shù)%n:分塊的大小%e:原圖像與保存局部DCT系數(shù)后的重建圖像之間的誤差矩陣% 按Zig-Zag掃描順序，根據(jù)DCTch參數(shù)，只保存64個(gè)% DCT系數(shù)中的前DCTch個(gè)系數(shù)，對(duì)修改后的DCT系數(shù)用逆DCT變換重建圖像，

42、得到DCT變% 換的壓縮圖像。計(jì)算重建圖像的均方根誤差RMSE ；顯示誤差圖像和誤差直方圖。zigzag = 1 2 6 7 15 16 28 29 3 5 8 14 17 27 30 43 4 9 13 18 26 31 42 44 10 12 19 25 32 41 45 54 11 20 24 33 40 46 53 55 21 23 34 39 47 52 56 61 22 35 38 48 51 57 60 62 36 37 49 50 58 59 63 64;%設(shè)置z掃描順序mask=zigzag<=DCTch; %根據(jù)當(dāng)前DCTch值得到“Z字掃描的邏輯值,mask為log

43、ic類(lèi)型%對(duì)修改后的DCT系數(shù)用逆DCT變換重建圖像，得到DCT變換的壓縮圖像D=dctmt*(n);I=blkproc(dctfre,n n,'P1*(*.*P2)*P3',D',maskbuf,D);%I為重建的壓縮圖像矩陣e=oldbuf-I;%e:原圖像與保存局部DCT系數(shù)后的重建圖像之間的誤差矩陣I=I+128;rmse1=RMSE(e);endfunction dctfre = DctBlock(oldbuf,n)%分塊DCT函數(shù)：根據(jù)給定的n值，計(jì)算原圖像的n×n塊的DCT系數(shù)，并轉(zhuǎn)換為可視頻譜圖以便觀察 % oldbuf 原始圖像數(shù)據(jù) % n

44、分塊的大小 % dctfre DCT系數(shù)矩陣D=dctmt*(n);%D是返回N×N的DCT變換矩陣，矩陣A的DCT變換可用D×A×D來(lái)計(jì)算dctfre=blkproc(oldbuf,n,n,'P1*P2',D,D'); %D'為D的轉(zhuǎn)置end4. QuanIQuan.mfunction newbuf,e,rmse2=QuanIQuan(oldbuf,dctfre,n,scale)%量化和反量化函數(shù)：根據(jù)給定的默認(rèn)JPEG量化表，%對(duì)每個(gè)n×n塊的DCT系數(shù)進(jìn)展量化和反量化，顯示量化誤差圖像及其直方圖。%oldbuf:原始

45、圖像數(shù)據(jù)%dctfre:DCT系數(shù)矩陣%n:分塊的大小%scale;量化系數(shù)z= 16 11 10 16 24 40 51 61 12 12 14 19 26 58 60 55 14 13 16 24 40 57 69 56 14 17 22 29 51 87 80 62 18 22 37 56 68 109 103 77 24 35 55 64 81 104 113 92 49 64 78 87 103 121 120 101 72 92 95 98 112 100 103 99;%默認(rèn)JPEG量化表 Qvalue=blkproc(dctfre,n n,'round(*./P1)&#

46、39;,scale*z);%量化 IQvalue=blkproc(Qvalue,n n,'*.*P1',scale*z);%反量化%對(duì)經(jīng)過(guò)量化和反量化后的矩陣進(jìn)展逆DCT變換得到重建圖像矩陣 D=dctmt*(n); newbuf=blkproc(IQvalue,n n,'P1*P2',D',D); e=newbuf-oldbuf;%e為量化誤差矩陣 rmse2=RMSE(e);%求均方根誤差 newbuf=newbuf+128;end5. RMSE.mfunction rmse=RMSE(oldbuf)%求均方根誤差M,N=size(oldbuf);

47、e=oldbuf.2; rmse=sqrt(sum(e(:)/(M*N);end【輸出圖像】實(shí)驗(yàn)8：運(yùn)動(dòng)估計(jì)【實(shí)驗(yàn)?zāi)康摹渴煜み\(yùn)動(dòng)估計(jì)的塊匹配BMA算法原理，編程實(shí)現(xiàn)全搜索算法三步搜索或鉆石搜索算法，了解運(yùn)動(dòng)估計(jì)在混合編碼器中的作用?！緦?shí)驗(yàn)容】1編寫(xiě)全搜索算法函數(shù)，將運(yùn)動(dòng)矢量疊加到當(dāng)前幀上并顯示輸出；2顯示輸出預(yù)測(cè)幀、殘差幀和重建圖像，計(jì)算預(yù)測(cè)幀的PSNR。附：可供參考的Matlab函數(shù)有hold、quiver【實(shí)驗(yàn)代碼】1.主程序imgI=imread('CLA1.bmp');%定義參考幀imgP=imread('CLA2.bmp');%定義當(dāng)前幀subplot

48、(231);imshow(imgI);title('參考幀');subplot(232);imshow(imgP);title('當(dāng)前幀');imgI=double(imgI);imgP=double(imgP);mbSize=16;% 塊尺寸為16*16p=7;%搜索窗口為(2p+1)*(2p+1)motionVect,ESputations,blk_center,costs=ME_ES(imgP,imgI,mbSize,p);%基于塊的全搜索算法imgMV(motionVect,imgP,blk_center);%畫(huà)運(yùn)動(dòng)矢量圖imgp=motionp(img

49、I,motionVect,mbSize);%根據(jù)運(yùn)動(dòng)矢量計(jì)算預(yù)測(cè)幀，并傳輸殘差幀psnr=imgPSNR(imgP,imgp);%計(jì)算峰值信噪比subplot(234);imshow(uint8(imgp);title('預(yù)測(cè)幀,PSNR=',num2str(psnr);imgErr=imgP-imgp;%殘差幀cal=Calibration(imgErr);%標(biāo)定，顯示更好效果subplot(235);imshow(cal);title('殘差幀');ChongJian=imgp+imgErr;%根據(jù)運(yùn)動(dòng)矢量指明的位置及殘差幀重建圖像subplot(236);

50、imshow(uint8(ChongJian);title('重建幀');2.ME_ES.mfunction motionVect,ESputations,blk_center,costs = ME_ES(imgP, imgI, mbSize, p)%function:FS算法：全搜索Full Search/E*haustive Search%img:當(dāng)前幀%imgI:參考幀%mbSize:MB尺寸%p:搜索窗口大小2p+1×2p+1%motionVect:整像素精度MV%ESputations:搜索每個(gè)宏塊所需的平均點(diǎn)數(shù)row,col=size(imgP);blk_

51、center=zeros(2,row*col/(mbSize2);%定義每個(gè)宏塊中心點(diǎn)位置motionVect=zeros(2,row*col/(mbSize2);%定義每個(gè)宏塊運(yùn)動(dòng)矢量costs=ones(2*p+1,2*p+1)*65537;putations=0;%搜索的點(diǎn)數(shù)之和mbCount=1;for i = 1:mbSize:row-mbSize+1 %當(dāng)前幀起始行搜索圍，步長(zhǎng)是塊數(shù)for j = 1:mbSize:col-mbSize+1 %當(dāng)前幀起始列搜索圍，步長(zhǎng)是塊數(shù)for m=-p:pfor n=-p:p ref_blk_row=i+m;%參考幀搜索框起始行 ref_blk

52、_col=j+n;%參考幀搜索框起始列%如果參考?jí)K的行列圍的任意一點(diǎn)在已經(jīng)搜索過(guò)的宏塊之外，則跳過(guò)，搜索下一點(diǎn)if (ref_blk_row<1|ref_blk_row+mbSize-1>row|ref_blk_col<1|ref_blk_col+mbSize-1>col) continue; end%否則計(jì)算該點(diǎn)SAD值 costs(m+p+1,n+p+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); putations=putations+1;endend blk_center(1,mbCount) = i+ mbSize/2-1;%記錄中心點(diǎn)行坐標(biāo) blk_center(2,mbCount) = j+ mbSize/2-1;%記錄中心點(diǎn)列坐標(biāo) minc,d*,dy=minCost(c