西交大數(shù)字圖像處理第六次作業(yè)

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

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

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

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

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

6、模糊濾波器如方程Eq. (5.6-11).(b) 模糊lena圖像：45度方向，T=1；(c) 再模糊的lena圖像中增加高斯噪聲，均值= 0 ，方差=10 pixels 以產生模糊圖像；(d)分別利用方程 Eq. (5.8-6)和(5.9-4)，恢復圖像。（1）問題分析：1）維納濾波器的推導： 圖像的退化模型為： （1） 其中，s(n1,n2)為原始圖像，b(n1,n2)為退化函數(shù)，w(n1,n2)為噪聲函數(shù)，x(n1,n2)為退化的圖像。并假設s與w不相關，w為0均值的平穩(wěn)隨機過程。圖像的復原模型為： （2） 其中，為恢復的圖像，為恢復濾波器。 誤差度量為： （3） 基于正交性原理，若要求

7、誤差最小，則必有下式成立： （4） 將（3）式帶入（4）式有： （5） 即 （6） 換元得： （7） 等式兩端同時取傅里葉變換得： （8） 即 （9） 公式（8）中（10）公式（10）兩端同時取傅里葉變換得： （11） 公式（8）中（12） 公式（12）兩端同時取傅里葉變換： （13） 將（11）式和（13）式帶入（8）式得 （14） 將符號化成與書中一致的表示 （15） 故表達式由下式給出 （16） 2）約束最小二乘方濾波 對于約束最小二乘方濾波，期望是找一個最小準則函數(shù)C，定義如下： （17） 其約束為 （18） 其中，是歐幾里得向量范數(shù)，是未退化圖像的估計。 這個最佳問題在頻率域中的解決

8、由下面的表達式給出： （19）其中，是一個參數(shù)，必須對它進行調整以滿足式（18）的條件，P(u,v)是函數(shù) （20） 的傅里葉變換。（2）MATLAB函數(shù)：1）imfilter功能：對任意類型數(shù)組或多維圖像進行濾波。用法：B = imfilter(A,H)B = imfilter(A,H,option1,option2,.)或寫做g = imfilter(f, w, filtering_mode, boundary_options, size_options)其中，f為輸入圖像，w為濾波掩模，g為濾波后圖像。filtering_mode用于指定在濾波過程中是使用“相關”還是“卷積”。bound

9、ary_options用于處理邊界充零問題，邊界的大小由濾波器的大小確定。具體參數(shù)選項見下表：選項描述filtering_mode corr通過使用相關來完成，該值為默認。conv通過使用卷積來完成。boundary_options X輸入圖像的邊界通過用值X（無引號）來填充擴展，其默認值為0。replicate圖像大小通過復制外邊界的值來擴展。symmetric圖像大小通過鏡像反射其邊界來擴展。circular圖像大小通過將圖像看成是一個二維周期函數(shù)的一個周期來擴展size_optionsfull輸出圖像的大小與被擴展圖像的大小相同same輸出圖像的大小與輸入圖像的大小相同。這可通過將濾波掩

10、模的中心點的偏移限制到原圖像中包含的點來實現(xiàn)，該值為默認值。2）fspecial功能：fspecial函數(shù)用于建立預定義的濾波算子。用法：h = fspecial(type)h = fspecial(type，para)其中type指定算子的類型，para指定相應的參數(shù)；選項描述typeaverageaveraging filter為均值濾波，參數(shù)為hsize代表模板尺寸，默認值為3 3diskcircular averaging filter為圓形區(qū)域均值濾波，參數(shù)為radius代表區(qū)域半徑，默認值為5.gaussianGaussian lowpass filter為高斯低通濾波，有兩個參數(shù)

11、，hsize表示模板尺寸，默認值為3 3，sigma為濾波器的標準值，單位為像素，默認值為0.5.laplacianlaplacian filter為拉普拉斯算子，參數(shù)alpha用于控制算子形狀，取值范圍為0，1，默認值為0.2。logLaplacian of Gaussian filter為拉普拉斯高斯算子，有兩個參數(shù)，hsize表示模板尺寸，默認值為3 3，sigma為濾波器的標準差，單位為像素，默認值為0.5.motionmotion filter為運動模糊算子，有兩個參數(shù)，表示攝像物體逆時針方向以theta角度運動了len個像素，len的默認值為9，theta的默認值為0；prewit

12、tPrewit thorizontal edge-emphasizing filter用于邊緣增強，大小為3 3，無參數(shù)sobelSobel horizontal edge-emphasizing filter用于邊緣提取，無參數(shù)unsharpunsharp contrast enhancement filter為對比度增強濾波器。參數(shù)alpha用于控制濾波器的形狀，范圍為0，1，默認值為0.2.（3）處理結果：（a）實現(xiàn)模糊濾波器如方程Eq. (5.6-11). 模糊濾波器的頻域表達式為： 故實現(xiàn)該濾波器，只需先將輸入圖像進行傅里葉變換并移至圖像中心，之后將圖像的傅里葉變換和模糊濾波器的傅里

13、葉變換進行陣列相乘，將得到的結果經過傅里葉反變換返回到空間域即可實現(xiàn)該濾波器。具體程序見附件mohu_filter.m（b）模糊lena圖像：45度方向，T=1；（a=0.1，b=0.1；T=1） lena.bmp原始圖像 lena運動模糊的結果（a=0.1，b=0.1；T=1） lena.bmp原始圖像 lena運動模糊的結果（調用MATLAB中的函數(shù)） （c）再模糊的lena圖像中增加高斯噪聲，均值= 0 ，方差=10 pixels 以產生模糊圖像；lena運動模糊的結果（a=0.1，b=0.1；T=1） 添加高斯噪聲的結果（均值為0，方差為0.01） lena運動模糊的結果（MATLAB

14、版） 添加高斯噪聲的結果（MATLAB版） （d）分別利用方程 Eq. (5.8-6)和(5.9-4)，恢復圖像。維納濾波的結果： lena運動模糊+高斯噪聲.bmp 維納濾波結果（K=0.06） 維納濾波結果（K=0.01） 維納濾波結果（K=0.5） lena運動模糊+高斯噪聲（MATLAB版） 維納濾波結果（MATLAB版） 約束最小二乘方濾波的結果： lena運動模糊+高斯噪聲（MATLAB版） 約束最小二乘濾波結果（MATLAB版） （3） 結果分析及總結：首先分別通過自己編寫的模糊函數(shù)和MATLAB中提供的imfilter和fspecial函數(shù)配合使用對圖像lena進行了模糊濾波。

15、發(fā)現(xiàn)套用書上的公式圖像是斜向下45度運動模糊，而MATLAB中函數(shù)模糊的結果是斜向上45度運動模糊，不過這不是重點可以接受，模糊的基本效果還是一致的。之后調用imnoise函數(shù)對兩幅圖像加入高斯噪聲，得到第二問的結果。之后分別使用自己編寫的函數(shù)和MATLAB中提供的deconvwnr函數(shù)進行維納濾波。調用MATLAB中函數(shù)濾波后的圖像得到了一定的改善，運動模糊的影響基本被消除，但噪聲的影響仍然較大，導致圖像質量下降；對于自己編寫的維納濾波函數(shù)，難點在于尋找令信噪比最大的K值，報告中顯示了部分K值對應的濾波結果，其中K=0.06，為信噪比最大時的濾波結果，從結果看，視覺上的效果并不是很理想，要想

16、達到更好的效果可能需要尋找更加合適的K值。最后采用MATLAAB提供的deconvreg函數(shù)進行約束最小二乘方濾波。從濾波后的結果看，約束最小二乘方濾波得到了比維納濾波更好的結果，尤其是對噪聲的濾除。由于實在是沒看懂書上的公式應該怎么實現(xiàn)，所以只好調用函數(shù)了。尤其實P(u,v)的維數(shù)和其他幾個矩陣的維數(shù)不同，怎么遍歷??？附錄【參考文獻】1 岡薩雷斯.數(shù)字圖像處理（第三版）北京：電子工業(yè)出版社，20112 周品.MATLAB數(shù)字圖像處理 北京：清華大學出版社，20123 楊杰.數(shù)字圖像處理及MATLAB實現(xiàn) 北京：電子工業(yè)出版社，2010【源代碼】add_gaussian_noise.m(題目1

17、添加高斯噪聲)I=imread(lena.bmp);figure(1);imshow(I);title(源圖像lena.bmp);imwrite(I,lena原始圖像.bmp);I2=imnoise(I,gaussian,0.5,0.01);figure(2);imshow(I2);title(加入gaussian噪聲后的lena.bmp);imwrite(I2,lena加入gaussian噪聲后(u=0.5，s2=0.01).bmp);suanshujunzhi_filter.m(題目1算術均值濾波恢復圖像)I=imread(lena_gaussian_noise.bmp);figure(1

18、);imshow(I);title(lena加入gaussian噪聲后(u=0，s2=0.01).bmp);imwrite(I,lena加入gaussian噪聲后(u=0，s2=0.01).bmp);n=5; h=1/n2.*ones(n,n);I2=conv2(I,h,same); I2=uint8(I2);figure(2);imshow(I2); title(算術均值濾波的結果（5x5）);imwrite(I2,算術均值濾波的結果（5x5）.bmp);median_filter.m(題目1中值濾波恢復圖像)I=imread(lena_gaussian_noise.bmp);figure(

19、1);imshow(I);title(lena加入gaussian噪聲后(u=0，s2=0.01).bmp);imwrite(I,lena加入gaussian噪聲后(u=0，s2=0.01).bmp);figure(2);n=5; a=ones(n,n);p=size(I); x1=double(I);x2=x1;for i=1:p(1)-n+1 for j=1:p(2)-n+1 c=x1(i:i+(n-1),j:j+(n-1); e=c(1,:); for u=2:n e=e,c(u,:); end mm=median(e); x2(i+(n-1)/2,j+(n-1)/2)=mm; ende

20、ndI2=uint8(x2);imshow(I2); title(中值濾波的結果（5x5）);imwrite(I2,中值濾波的結果（5x5）.bmp);mohu_filter.m(題目2運動模糊濾波器)I=imread(lena.bmp);figure(1);imshow(I);title(lena.bmp原始圖像);imwrite(I,lena原始圖像.bmp);f=double(I); F=fft2(f); F=fftshift(F); M,N=size(F);a=0.1;b=0.1;T=1;for u=1:M for v=1:N H(u,v)=(T/(pi*(u*a+v*b)*sin(p

21、i*(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(運動模糊化lena.bmp);imwrite(g,lena運動模糊的結果.bmp); Wiener_filter.m(題目2維納濾波器自編版)I=imread(lena運動模糊+高斯噪聲.bmp);figure(1);imshow(I);title(lena運動模糊+高斯噪聲);imw

22、rite(I,lena運動模糊+高斯噪聲.bmp);g=double(I); G=fft2(g); G=fftshift(G); M,N=size(G);a=0.1;b=0.1;T=1;K=0.0259;for u=1:M for v=1:N H(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)=1/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=u

23、int8(real(f);figure(2);imshow(f);title(維納濾波的將結果);imwrite(f,維納濾波的結果(K=0.0259).bmp);wiener_filter_k.m(題目2維納濾波器尋找信噪比最大的K值)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=0.1;b=0.1;T=1;i=1;format longf

24、or k=0.01:0.01:0.11 for u=1:M for v=1:N H(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)=(1/H(u,v)*(abs(H(u,v)2/(abs(H(u,v)2+k)*G(u,v); end end F=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;endidx=find(max(SNR)wiener_filter_matlab.m(題目2維納濾波器MATLAB版)I=imread(lena.bmp);H=fspecial(motion,50,45);I1=imfi