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

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

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

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

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

5、結(jié)果利用中值濾波器恢復(fù)圖像(5x5模板)lena添加高斯噪聲(u=0,s2=)中值濾波的結(jié)果結(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)下邊要求；實(shí)現(xiàn)模糊濾波器如方程Eq.模糊lena圖像：45度方向，T=1；再模糊的lena圖像中增加高

6、斯噪聲，均值=0，方差=10pixels以產(chǎn)生模糊圖像；分別利用方程Eq.和，恢復(fù)圖像。問(wèn)題分析：1)維納濾波器的推導(dǎo)：圖像的退化模型為：x(n,n)=b(n,n)*s(n,n)+w(n,n)(1)12121212其中，s(n1,n2)為原始圖像，b(n1,n2)為退化函數(shù)，w(nl,n2)為噪聲函數(shù)，x(n1,n2)為退化的圖像。并假設(shè)s與w不相關(guān)，w為0均值的平穩(wěn)隨機(jī)過(guò)程。圖像的復(fù)原模型為：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)為恢復(fù)的圖像，h(n,n)為恢復(fù)濾波器。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等式兩端同時(shí)取傅里葉變換得：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)兩端同時(shí)取傅里葉變換得：P(w,w)=B*(w,w)xP(w,w)TOC o 1-5 h zsx1212s12公式(8)中公式（12）兩端同時(shí)取傅里葉變換：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將符號(hào)化成與書(shū)中一致的表示15)故表達(dá)式由下式給出(16)2)約束最小二乘方濾波對(duì)于約束最小二乘方濾波，期望是找一個(gè)最小準(zhǔn)則函數(shù)C,定義如下:17)C=峻丈NTWf(x,y)217)x=0y=0其約束為18

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

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

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

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

13、motionfilter為運(yùn)動(dòng)模糊算子，有兩個(gè)參數(shù)，表示攝像物體逆時(shí)針?lè)较蛞詔heta角度運(yùn)動(dòng)了len個(gè)像素，len的默認(rèn)值為9，theta的默認(rèn)值為0;prewittPrewitthorizontaledge-emphasizingfilter用于邊緣增強(qiáng)，大小為33,無(wú)參數(shù)sobelSobelhorizontaledge-emphasizingfilter用于邊緣提取，無(wú)參數(shù)unsharpunsharpcontrastenhancementfilter為對(duì)比度增強(qiáng)濾波器。參數(shù)alpha用于控制濾波器的形狀，范圍為0,1,默認(rèn)值為.3）處理結(jié)果：（a）實(shí)現(xiàn)模糊濾波器如方程Eq.模糊濾波器的頻

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

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

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

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

18、北京：電子工業(yè)出版社，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算術(shù)均值濾波恢復(fù)圖像)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(算術(shù)均值濾波的結(jié)果(5x5);imwrite(I2,算術(shù)均值濾波的結(jié)果(5x5).bmp);(題目1中值濾波恢復(fù)圖像)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(中值濾波的結(jié)果(5x5);imwrite(I2,中值濾波的結(jié)果(5x5).bmp);(題目2運(yùn)動(dòng)模糊濾波器)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(運(yùn)動(dòng)模糊化);imwrite(g,lena運(yùn)

22、動(dòng)模糊的結(jié)果.bmp);(題目2維納濾波器自編版)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=;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(維納濾波的將結(jié)果);imwrite(f,維納濾波的結(jié)果(K=.bmp);(題目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=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(維納濾波的結(jié)果)；e=f-uint8(g);SNR(i)=sum(sum(g.2)/sum(sum(e.2);i=i+1;endid