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

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

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

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

4、參數(shù)Type對(duì)應(yīng)的噪聲類型如下： gaussian高斯白噪聲 localvar0均值白噪聲 poisson泊松噪聲 salt & pepper鹽椒噪聲 speckle乘性噪聲（3） 處理結(jié)果：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）圖像恢復(fù)（選取被均值為0，方差為0.01的高斯噪聲污染的圖像為例）利用算術(shù)均值濾波器恢復(fù)圖像(5x5模板) len

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

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

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

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

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

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

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

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

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

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

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

16、達(dá)到更好的效果可能需要尋找更加合適的K值。最后采用MATLAAB提供的deconvreg函數(shù)進(jìn)行約束最小二乘方濾波。從濾波后的結(jié)果看，約束最小二乘方濾波得到了比維納濾波更好的結(jié)果，尤其是對(duì)噪聲的濾除。由于實(shí)在是沒看懂書上的公式應(yīng)該怎么實(shí)現(xiàn)，所以只好調(diào)用函數(shù)了。尤其實(shí)P(u,v)的維數(shù)和其他幾個(gè)矩陣的維數(shù)不同，怎么遍歷啊？附錄【參考文獻(xiàn)】1 岡薩雷斯.數(shù)字圖像處理（第三版）北京：電子工業(yè)出版社，20112 周品.MATLAB數(shù)字圖像處理 北京：清華大學(xué)出版社，20123 楊杰.數(shù)字圖像處理及MATLAB實(shí)現(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算術(shù)均值濾波恢復(fù)圖像)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(算術(shù)均值濾波的結(jié)果（5x5）);imwrite(I2,算術(shù)均值濾波的結(jié)果（5x5）.bmp);median_filter.m(題目1中值濾波恢復(fù)圖像)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(中值濾波的結(jié)果（5x5）);imwrite(I2,中值濾波的結(jié)果（5x5）.bmp);mohu_filter.m(題目2運(yùn)動(dòng)模糊濾波器)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(運(yùn)動(dòng)模糊化lena.bmp);imwrite(g,lena運(yùn)動(dòng)模糊的結(jié)果.bmp); Wiener_filter.m(題目2維納濾波器自編版)I=imread(lena運(yùn)動(dòng)模糊+高斯噪聲.bmp);figure(1);imshow(I);title(lena運(yùn)動(dòng)模糊+高斯噪聲);imw

22、rite(I,lena運(yùn)動(dòng)模糊+高斯噪聲.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(維納濾波的將結(jié)果);imwrite(f,維納濾波的結(jié)果(K=0.0259).bmp);wiener_filter_k.m(題目2維納濾波器尋找信噪比最大的K值)I=imread(lena運(yùn)動(dòng)模糊+高斯噪聲.bmp);figure(1);imshow(I);title(lena運(yùn)動(dòng)模糊+高斯噪聲);imwrite(I,lena運(yùn)動(dòng)模糊+高斯噪聲.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(維納濾波的結(jié)果); 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