1、數字圖像處理第六次作業姓名：班級：學號：提交日期：2015年4月27日1、在測試圖像上產生高斯噪聲lena圖-需能指定均值和方差；并用濾波器（自選）恢復圖像；問題分析：圖像退化模型：圖像退化過程被建模為一個退化函數和一個加性噪聲項，對一幅輸入圖像f(x,y)進行處理，產生一副退化后的圖像g(x,y)。給定g(x,y)和關于退化函數H的一些知識以及關于加性噪聲項n(x,y)的一些知識后，圖像復原的目的就是獲得原始圖像的一個估計。如果H是一個線性的、位置不變的過程，那么空間域中的退化圖像可由下式給出：g(x,y)=h(x,y)f(x,y)+n(x,y)其中，h(x,y)是退化函數的空間表示；符號”

2、表示空間卷積.等價的頻率域表示：G(u,v)=H(u,v)F(u,v)+N(u,v)高斯噪聲：所謂高斯噪聲是指它的概率密度函數服從高斯分布(即正態分布)的一類噪聲。一個高斯隨機變量z的PDF可表示為：P(z)=其中z代表灰度，u是z的均值，是z的標準差。高斯噪聲的灰度值多集中在均值附近。算術均值濾波器：令Sxy表示中心在點(x,y)處，大小為mXn的矩形子圖像窗口的一組坐標。算術均值濾波器在Sxy定義的區域中計算被污染的圖像g(x,y)的平均值。在點(x,y)處復原圖像的值：1f(x,y)二-工g(s,t)mn(s,t)eSXy這個操作可以使用大小為mXn的一個空間濾波器來實現，其所有系數均為

3、其值的1/mn。均值濾波器平滑一幅圖像中的局部變化，雖然模糊了結果，但降低了噪聲。中值濾波器：統計排序濾波器是空間域濾波器，空間濾波器的響應是基于由該濾波器包圍的圖像區域中的像素值的排序。中值濾波器使用一個像素鄰域中灰度級的中值來替代該像素值，即：f(x,y)二mediang(s,t)(s，t)Sxy編程思路：首先利用imnoise函數給圖像添加高斯噪聲,之后分別利用算術平均值濾波器和中值濾波器進行濾波并觀察效果。MATLAB函數：g=imnoise(f,type,parameters)函數功能：使用函數imnoise來用噪聲污染一幅圖像、調用格式：J=imnoise(I,type)J=imn

4、oise(I,type,parameters)參數Type對應的噪聲類型如下：gaussian高斯白噪聲localvar0均值白噪聲poisson泊松噪聲salt&pepper鹽椒噪聲speckle乘性噪聲處理結果：1)添加高斯噪聲原始圖像lena添加高斯噪聲(u=0,s2=)lena添加高斯噪聲lena添加高斯噪聲(u=0,s2=)lena添加高斯噪聲(u=0,s2=)lena添加高斯噪聲(u=,s2=)lena添加高斯噪聲(u=,s2=)2)圖像恢復(選取被均值為0,方差為的高斯噪聲污染的圖像為例)利用算術均值濾波器恢復圖像(5x5模板)lena添加高斯噪聲(u=0,s2=)算術均值濾波的

5、結果利用中值濾波器恢復圖像(5x5模板)lena添加高斯噪聲(u=0,s2=)中值濾波的結果結果分析及總結：首先通過imnoise函數分別產生了被不同均值和方差的高斯噪聲污染的圖像。當高斯噪聲均值不變為0時,隨著方差增加,圖像噪聲越嚴重；當高斯噪聲方差不變時,均值會影響到整個圖像的灰度值,使整個圖像變亮。與理論上均值和方差對圖像的影響一致。分別使用算術均值濾波器和中值濾波器對加噪圖像進行恢復。兩種方法在一定程度上都可以降低噪聲。算術均值濾波器降低噪聲的同時也模糊了圖像。2、推導維納濾波器并實現下邊要求；實現模糊濾波器如方程Eq.模糊lena圖像：45度方向，T=1；再模糊的lena圖像中增加高

6、斯噪聲，均值=0，方差=10pixels以產生模糊圖像；分別利用方程Eq.和，恢復圖像。問題分析：1)維納濾波器的推導：圖像的退化模型為：x(n,n)=b(n,n)*s(n,n)+w(n,n)(1)12121212其中，s(n1,n2)為原始圖像，b(n1,n2)為退化函數，w(nl,n2)為噪聲函數，x(n1,n2)為退化的圖像。并假設s與w不相關，w為0均值的平穩隨機過程。圖像的復原模型為：s(n,n)=h(n,n)*x(n,n)=Hh(l,l)xx(n-1,n-1)(2)121212121122（3）（3）（14）（14）4）5）4）5）（7）（8）（9）（11）（12）（13）H(w,

7、w)=12其中，s(n,n)為恢復的圖像，h(n,n)為恢復濾波器。1212誤差度量為：e2=E(s(n,n)-S(n,n)21212基于正交性原理，若要求誤差最小，則必有下式成立：Ee(n,n)xx*(m,m)二01212將(3)式帶入(4)式有：Es(n,n)xx*(m,m)=Es(n,n)xx*(m,m)12121212即(6)換元得：R(n,n)=h(n,n)*R(n,n)TOC o 1-5 h zsx1212xx12等式兩端同時取傅里葉變換得：P(w,w)=H(w,w)xP(w,w)sx1212x12即(w,w)P(w,w)H(w,w)=sx12-12P(w,w)x12公式(8)中(

8、10)公式(10)兩端同時取傅里葉變換得：P(w,w)=B*(w,w)xP(w,w)TOC o 1-5 h zsx1212s12公式(8)中公式（12）兩端同時取傅里葉變換：P(w,w)=B(w,w)|2xP(w,w)+P(w,w)x1212丨s12w12將(11)式和(13)式帶入(8)式得B*(w,w)xP(w,w).12s12B(w,w)|2xP(w,w)+P(w,w)12s12w12將符號化成與書中一致的表示15)故表達式由下式給出(16)2)約束最小二乘方濾波對于約束最小二乘方濾波，期望是找一個最小準則函數C,定義如下:17)C=峻丈NTWf(x,y)217)x=0y=0其約束為18

9、)|g-Hf|2=|h|218)其中，|W|2全WTW是歐幾里得向量范數，f是未退化圖像的估計。這個最佳問題在頻率域中的解決由下面的表達式給出：F(u,v)=F(u,v)=H*(u,v)|H(u,v)|2+丫|P(u,v)|2G(u,v)19)其中，丫是一個參數，必須對它進行調整以滿足式(18)的條件，P(u,v)是函數0-10p(x,y)=-14-1(20)0-10的傅里葉變換。MATLAB函數：1)imfilter功能：對任意類型數組或多維圖像進行濾波。用法：B=imfilter(A,H)B=imfilter(A,H,option1,option2,.)或寫做g=imfilter(f,w,

10、filtering_mode,boundary_options,size_options)其中，f為輸入圖像，w為濾波掩模，g為濾波后圖像。filtering_mode用于指定在濾波過程中是使用“相關”還是“卷積”。boundary_options用于處理邊界充零問題，邊界的大小由濾波器的大小確定。具體參數選項見下表：選項描述filteringmodecorr通過使用相關來完成，該值為默認。conv通過使用卷積來完成。boundary_optionsX輸入圖像的邊界通過用值X(無引號)來填充擴展，其默認值為0。replicate圖像大小通過復制外邊界的值來擴展。symmetrie圖像大小通過鏡

11、像反射其邊界來擴展。circular圖像大小通過將圖像看成是一個二維周期函數的一個周期來擴展sizeoptionsfull輸出圖像的大小與被擴展圖像的大小相同same輸出圖像的大小與輸入圖像的大小相同。這可通過將濾波掩模的中心點的偏移限制到原圖像中包含的點來實現，該值為默認值。2)fspecial功能：fspecial函數用于建立預定義的濾波算子。用法：h=fspecial(type)h=fspecial(type，para)其中type指定算子的類型，para指定相應的參數;選項描述typeaverageaveragingfilter為均值濾波,參數為hsize代表模板尺寸，默認值為33di

12、skcircularaveragingfilter為圓形區域均值濾波，參數為radius代表區域半徑，默認值為5.gaussianGaussianlowpassfilter為咼斯低通濾波，有兩個參數，hsize表示模板尺寸，默認值為33，sigma為濾波器的標準值，單位為像素，默認值為.laplacianlaplacianfilter為拉普拉斯算子，參數alpha用于控制算子形狀，取值范圍為0,1,默認值為。logLaplacianofGaussianfilter為拉普拉斯咼斯算子，有兩個參數，hsize表示模板尺寸，默認值為33，sigma為濾波器的標準差，單位為像素，默認值為.motion

13、motionfilter為運動模糊算子，有兩個參數，表示攝像物體逆時針方向以theta角度運動了len個像素，len的默認值為9，theta的默認值為0;prewittPrewitthorizontaledge-emphasizingfilter用于邊緣增強，大小為33,無參數sobelSobelhorizontaledge-emphasizingfilter用于邊緣提取，無參數unsharpunsharpcontrastenhancementfilter為對比度增強濾波器。參數alpha用于控制濾波器的形狀，范圍為0,1,默認值為.3）處理結果：（a）實現模糊濾波器如方程Eq.模糊濾波器的頻

14、域表達式為：TH（u,v）二sin兀（ua+vb）e-丿兀（ua+vb）兀（ua+vb）故實現該濾波器，只需先將輸入圖像進行傅里葉變換并移至圖像中心，之后將圖像的傅里葉變換和模糊濾波器的傅里葉變換進行陣列相乘，將得到的結果經過傅里葉反變換返回到空間域即可實現該濾波器。具體程序見附件（b）模糊lena圖像：45度方向，T=l；（a二，b二；T=1）原始圖像lena運動模糊的結果（a二，b二；T=l）原始圖像lena運動模糊的結果（調用MATLAB中的函數）原始圖像（c）再模糊的lena圖像中增加高斯噪聲，均值=0，方差=10pixels以產生模糊圖像;lena運動模糊的結果lena運動模糊的結果

15、（a二,b二;T=1）添加高斯噪聲的結果（均值為0，方差為）lena運動模糊的結果（lena運動模糊的結果（MATLAB版）添加高斯噪聲的結果（MATLAB版）（d）分別利用方程（d）分別利用方程Eq.和，恢復圖像。維納濾波的結果：lena運動模糊+高斯噪聲.bmp維納濾波結果（K=）維納濾波結果（K=）維納濾波結果（K=）維納濾波結果（K=）lena運動模糊lena運動模糊+高斯噪聲（MATLAB版）維納濾波結果（MATLAB版）約束最小二乘方濾波的結果：約束最小二乘濾波結果（MATLAB版）lena運動模糊+高斯噪聲（約束最小二乘濾波結果（MATLAB版）（3）結果分析及總結：首先分別通過

16、自己編寫的模糊函數和MATLAB中提供的imfilter和fspecial函數配合使用對圖像lena進行了模糊濾波。發現套用書上的公式圖像是斜向下45度運動模糊，而MATLAB中函數模糊的結果是斜向上45度運動模糊，不過這不是重點可以接受，模糊的基本效果還是一致的。之后調用imnoise函數對兩幅圖像加入高斯噪聲，得到第二問的結果。之后分別使用自己編寫的函數和MATLAB中提供的deconvwnr函數進行維納濾波。調用MATLAB中函數濾波后的圖像得到了一定的改善，運動模糊的影響基本被消除，但噪聲的影響仍然較大，導致圖像質量下降；對于自己編寫的維納濾波函數，難點在于尋找令信噪比最大的K值，報告

17、中顯示了部分K值對應的濾波結果，其中K=，為信噪比最大時的濾波結果，從結果看，視覺上的效果并不是很理想，要想達到更好的效果可能需要尋找更加合適的K值。最后采用MATLAAB提供的deconvreg函數進行約束最小二乘方濾波。從濾波后的結果看，約束最小二乘方濾波得到了比維納濾波更好的結果，尤其是對噪聲的濾除。由于實在是沒看懂書上的公式應該怎么實現，所以只好調用函數了。尤其實P(u,v)的維數和其他幾個矩陣的維數不同，怎么遍歷啊附錄【參考文獻】岡薩雷斯.數字圖像處理(第三版)北京：電子工業出版社，2011周品.MATLAB數字圖像處理北京：清華大學出版社，2012楊杰數字圖像處理及MATLAB實現

18、北京：電子工業出版社，2010【源代碼】(題目1添加高斯噪聲)I=imread();figure(1);imshow(I);title(源圖像);imwrite(I,lena原始圖像.bmp);I2=imnoise(I,gaussian,;figure(2);imshow(I2);title(加入gaussian噪聲后的);imwrite(12,lena力口入gaussian噪聲后(u=，s2=.bmp);(題目1算術均值濾波恢復圖像)I=imread();figure(1);imshow(I);title(lena力口入gaussian噪聲后(u=0,s2=.bmp);imwrite(I,l

19、ena力口入gaussian噪聲后(u=0,s2=.bmp);n=5;h=l/n2.*ones(n,n);I2=conv2(I,h,same);I2=uint8(I2);figure(2);imshow(I2);title(算術均值濾波的結果(5x5);imwrite(I2,算術均值濾波的結果(5x5).bmp);(題目1中值濾波恢復圖像)I=imread();figure(1);imshow(I);title(lena力口入gaussian噪聲后(u=0,s2=.bmp);imwrite(I,lena力口入gaussian噪聲后(u=0,s2=.bmp);figure(2);n=5;a=on

20、es(n,n);p=size(I);x1=double(I);x2=x1;fori=1:p(1)-n+1forj=1:p(2)-n+1c=x1(i:i+(n-1),j:j+(n-1);e=c(1,:);foru=2:ne=e,c(u,:);endmm=median(e);x2(i+(n-1)/2,j+(n-1)/2)=mm;endendI2=uint8(x2);imshow(I2);title(中值濾波的結果(5x5);imwrite(I2,中值濾波的結果(5x5).bmp);(題目2運動模糊濾波器)I=imread();figure(1);imshow(I);title(原始圖像)；imwr

21、ite(I,lena原始圖像.bmp);f=double(I)；F=fft2(f);F=fftshift(F);M,N=size(F);a=;b=;T=1;foru=1:Mforv=1:NH(u,v)=(T/(pi*(u*a+v*b)*sin(pi*(u*a+v*b)*exp(-sqrt(-1)*pi*(u*a+v*b);G(u,v)=H(u,v)*F(u,v);endendG=ifftshift(G);g=ifft2(G);g=256.*g./max(max(g);g=uint8(real(g);figure(2);imshow(g);title(運動模糊化);imwrite(g,lena運

22、動模糊的結果.bmp);(題目2維納濾波器自編版)I=imread(lena運動模糊+高斯噪聲.bmp);figure(1);imshow(I);title(lena運動模糊+高斯噪聲);imwrite(I,lena運動模糊+高斯噪聲.bmp);g=double(I);G=fft2(g);G=fftshift(G);M,N=size(G);a=;b=;T=1;K=;foru=1:Mforv=1:NH(u,v)=(T/(pi*(u*a+v*b)*sin(pi*(u*a+v*b)*exp(-sqrt(-1)*pi*(u*a+v*b);F(u,v)=l/H(u,v)*(abs(H(u,v)2/(ab

23、s(H(u,v)2+K)*G(u,v);endendF=ifftshift(F);f=ifft2(F);f=256.*f./max(max(f);f=uint8(real(f);figure(2);imshow(f);title(維納濾波的將結果);imwrite(f,維納濾波的結果(K=.bmp);(題目2維納濾波器尋找信噪比最大的K值)I=imread(lena運動模糊+高斯噪聲.bmp);figure(1);imshow(I);title(lena運動模糊+高斯噪聲);imwrite(I,lena運動模糊+高斯噪聲.bmp);g=double(I);G=fft2(g);G=fftshif

24、t(G);M,N=size(G);a=;b=;T=1;i=1;formatlongfork=:foru=1:Mforv=1:NH(u,v)=(T/(pi*(u*a+v*b)*sin(pi*(u*a+v*b)*exp(-sqrt(-1)*pi*(u*a+v*b);F(u,v)=(l/H(u,v)*(abs(H(u,v)2/(abs(H(u,v)2+k)*G(u,v);endendF=ifftshift(F);f=ifft2(F);f=256.*f./max(max(f);f=uint8(real(f);figure;imshow(f);title(維納濾波的結果)；e=f-uint8(g);SNR(i)=sum(sum(g.2)/sum(sum(e.2);i=i+1;endid