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

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

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

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

4、格式：J = imnoise(I,type)J = imnoise(I,type,parameters)參數(shù)Type對(duì)應(yīng)的噪聲類(lèi)型如下：'gaussian'高斯白噪聲'localvar'0均值白噪聲'poisson'泊松噪聲'salt & pepper'鹽椒噪聲'speckle'乘性噪聲(3)處理結(jié)果：1)添加高斯噪聲原始圖像lena添加高斯噪聲(u=0, sA2=)lena添加高斯噪聲(u=0, $人2=)lena添加高斯噪聲(u=0, $人2=)lena添加高斯噪聲（u=, sA2=）lena添加高斯噪

5、聲（u=, $八2=）2）圖像恢復(fù)（選取被均值為0,方差為的高斯噪聲污染的圖像為例）利用算術(shù)均值濾波器恢復(fù)圖像（5x5模板）利用中值濾波器恢復(fù)圖像（5x5模板）lena添加高斯噪聲（u=0, $八2=）中值濾波的結(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í)也模

6、糊了圖像。2、推導(dǎo)維納濾波器并實(shí)現(xiàn)下邊要求；(a) 實(shí)現(xiàn)模糊濾波器如方程Eq. .(b) 模糊 lena 圖像： 45 度方向，T=1；(c) 再模糊的lena 圖像中增加高斯噪聲，均值= 0 ，方差 =10 pixels 以產(chǎn)生模糊圖像；(d)分別利用方程Eq.和，恢復(fù)圖像。( 1 )問(wèn)題分析：1 )維納濾波器的推導(dǎo)：圖像的退化模型為：x(n1,n2) b(n1,n2) s(n1,n2) w(n1,n2)( 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ò)程。圖像

7、的復(fù)原模型為：s?(n1,n2)=h(n1,n2) x(n1,n2)h(l1,l2) x(n1 l1,n2 l2)( 2)l1 l2其中，費(fèi)必，”)為恢復(fù)的圖像，h(n1,n2)為恢復(fù)濾波器。誤差度量為：e2E( s(n1,n2) s?(n1,n2)2( 3)基于正交性原理，若要求誤差最小，則必有下式成立：Ee(n1,n2) x (m1,m2)0( 4)將(3)式帶入(4)式有：Es(n 1, n2) x (m1,m2) Es?(n1,n2) x (m1,m2)( 5)即6)換元得：7)Rsx (n1, n2)h(n1,n2) Rxx(n1,n2)等式兩端同時(shí)取傅里葉變換得:PSx(Wl,W2

8、) H(Wi,W2)PX(Wl,W2)即Psx(Wl,W2)H (w 1, w2)Px(Wi, W2)公式（10）兩端同時(shí)取傅里葉變換得:PSx(Wi,W2)B(w1, W2) PS(w1,W2)(8)(9)(10)(11)公式(8公式,(12)IpXIWillltaiwiiHijkps (Willlhpw(w1,wJ_)_(13).(11)式和8)WS)BJRxTn公式(8)中將符號(hào)化成與書(shū)中一致的表示(15)故表達(dá)式由下式給出(16)2)約束最小二乘方濾波對(duì)于約束最小二乘方濾波，期望是找一個(gè)最小準(zhǔn)則函數(shù)C,定義如下:M 1 N 122C 2 f(x, y)2x 0 y 0(17)其約束為g

9、 Hf?|(18)其中，|wwTw是歐幾里得向量范數(shù)，？是未退化圖像的估計(jì)。這個(gè)最佳問(wèn)題在頻率域中的解決由下面的表達(dá)式給出：F?(u,v) H2-(u2)2G(u, v)H(uv) P(u v)1) imfilter功能：對(duì)任意類(lèi)型數(shù)組或多維圖像進(jìn)行濾波。用法：B = imfilter(A,H)B = imfilter(A,H,option1,option2,)或?qū)懽?g = imfilter(f, w, filtering_mode, boundary_options, size_options)其中，f為輸入圖像，w為濾波掩模，g為濾波后圖像。filtering_mode用于指定在濾波過(guò)程

10、 中是使用"相關(guān)"還是"卷積"。boundary_options用于處理邊界充零問(wèn)題，邊界的大小由濾波器的大小確定。具體參數(shù)選項(xiàng)見(jiàn)下表：選項(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

11、 '圖像大小通過(guò)將圖像看成是一個(gè)二維周期函數(shù) 的一個(gè)周期來(lái)擴(kuò)展size_options'full '輸出圖像的大小與被擴(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)描述type'average'averaging filter 為均值

12、濾波， 參數(shù)為hsize代表模板尺寸， 默認(rèn)值為3 3'disk'circular averaging filter 為圓形 區(qū)域均值濾波，參數(shù)為radius 代表區(qū)域半徑，默認(rèn)值為 5.'gaussian'Gaussian lowpass filter 為高斯 低通濾波，后兩個(gè)參數(shù)，hsize 表示模板尺寸，默認(rèn)值為33, sigma為濾波器的標(biāo)準(zhǔn)值， 單位為像素，默認(rèn)值為.'laplacian'laplacian filter 為拉普拉斯算 子，參數(shù)alpha用丁控制算子 形狀，取值范圍為0, 1,默 認(rèn)值為。'log'Lap

13、lacian of Gaussian filter 為拉普拉斯高斯算子，有兩個(gè)參數(shù)，hsize表示模板尺寸，默認(rèn)值為3 3, sigma為濾波 器的標(biāo)準(zhǔn)差，單位為像素，默 認(rèn)值為.'motion'motion filter為運(yùn)動(dòng)模糊算子，有兩個(gè)參數(shù)，表小攝像物 體逆時(shí)針?lè)较蛞詔heta角度運(yùn) 動(dòng)了 len個(gè)像素，len的默認(rèn) 值為9, theta的默認(rèn)值為0;'prewitt'Prewit thorizontaledge-emphasizing filter用于邊緣增強(qiáng)，大小為3 3, 無(wú)參數(shù)'sobel'Sobel horizontal edg

14、e-emphasizing filter用于邊緣提取，無(wú)參數(shù)'unsharp'unsharpcontrastenhancement filter 為對(duì)比度增強(qiáng)濾波器。參數(shù)alpha用丁 控制濾波器的形狀，范圍為 0, 1,默認(rèn)值為.(3)處理結(jié)果：(a)實(shí)現(xiàn)模糊濾波器如方程Eq.模糊濾波器的頻域表達(dá)式為：H(u,v)sin (ua vb)e j ()(ua vb)故實(shí)現(xiàn)該濾波器，只需先將輸入圖像進(jìn)行傅里葉變換并移至圖像中心，之后將圖像的傅里葉變換和模糊濾波器的傅里葉變換進(jìn)行陣列相乘，將得到的結(jié)果經(jīng)過(guò)傅里葉反變換返回到空間域即可實(shí)現(xiàn)該濾波器。具體程序見(jiàn)附件(b)模糊 lena 圖

15、像：45 度方向，T=1; ( a=, b=; T=1)原始圖像lena運(yùn)動(dòng)模糊的結(jié)果(a=, b=; T=1)原始圖像lena運(yùn)動(dòng)模糊的結(jié)果(調(diào)用MATLAB中的函數(shù))(c)再模糊的lena圖像中增加高斯噪聲，均值=0 ,方差=10 pixels以產(chǎn)生模糊圖像;lena運(yùn)動(dòng)模糊的結(jié)果(a=, b=; T=1) 添加高斯噪聲的結(jié)果(均值為0,方差為)lena運(yùn)動(dòng)模糊的結(jié)果（MATLAB版）添加高斯噪聲的結(jié)果（MATLAB版）維納濾波結(jié)果（K=）（d）分別利用方程 Eq.和，恢復(fù)圖像。維納濾波的結(jié)果：lena運(yùn)動(dòng)模糊+高斯噪聲.bmp維納濾波結(jié)果（K=）維納濾波結(jié)果（K=）lena運(yùn)動(dòng)模糊+高斯

16、噪聲（MATLAB版）維納濾波結(jié)果（MATLAB版）約束最小二乘濾波結(jié)果（MATLAB版）約束最小二乘方濾波的結(jié)果：lena運(yùn)動(dòng)模糊+高斯噪聲（MATLAB版）（3）結(jié)果分析及總結(jié)：首先分別通過(guò)自己編寫(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ù)

17、進(jìn)行維納濾波。調(diào)用MATLAB中函數(shù)濾波后的圖像得到了一定的改善，運(yùn)動(dòng)模糊的影響基本被消除，但噪聲的影 響仍然較大，導(dǎo)致圖像質(zhì)量下降；對(duì)于自己編寫(xiě)的維納濾波函數(shù)，難點(diǎn)在于尋找令信噪比最大的K值，報(bào)告中顯示了部分 K值對(duì)應(yīng)的濾波Z果，其中 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ù)

18、和其他幾個(gè)矩陣的維數(shù)不同，怎么遍歷啊附錄【參考文獻(xiàn)】1 岡薩雷斯.數(shù)字圖像處理(第三版)北京：電子工業(yè)出版社，20112周品.MATLAB數(shù)字圖像處理 北京：清華大學(xué)出版社，20123楊杰.數(shù)字圖像處理及 MATLAB實(shí)現(xiàn) 北京：電子工業(yè)出版社，2010【源代碼】(題目1 添加高斯噪聲)I=imread('');figure(1);imshow(I);title( '源圖像 ');imwrite(I, 'lena 原始圖像.bmp' );I2=imnoise(I, 'gaussian',;figure(2);imshow(I2);

19、title( ' 加入 gaussian 噪聲后的');imwrite(I2, 'lena 力口入 gaussian 噪聲后(u=, sA2=.bmp');(題目1 算術(shù)均值濾波恢復(fù)圖像)I=imread('');figure(1);imshow(I);title( 'lena 加入 gaussian 噪聲后(u=0， sA2=.bmp' );imwrite(I, 'lena 加入 gaussian 噪聲后(u=0， sA2=.bmp' );n=5;h=1/nA2.*ones(n,n);I2=conv2(I,h,

20、'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 , sA2=.bmp');imwrite(I, 'lena 力口入 gaussian 噪聲后(u=0 , sA2=.bmp');figu

21、re(2);n=5;a=ones(n,n);p=size(I);x1=double(I);x2=x1;for i=1:p(1)-n+1for j=1:p(2)-n+1c=x1(i:i+(n-1),j:j+(n-1);e=c(1,:);for u=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

22、=imread('');figure(1);imshow(I);title( '原始圖像');imwrite(I, 'lena 原始圖像.bmp' );f=double(I);F=fft2(f);F=fftshift(F);M,N=size(F);a=;b=;T=1;for u=1:Mfor v=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=2

23、56.*g./max(max(g);g=uint8(real(g);figure(2);imshow(g);title( '運(yùn)動(dòng)模糊化');imwrite(g, 'lena 運(yùn)動(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=fftshif

24、t(G);M,N=size(G);a=;b=;T=1;K=;for u=1:Mfor v=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)=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=uint8(real(f);figure(2);imshow(f);title( '維納濾波的將結(jié)果');imwrite(f, '

25、;維納濾波的結(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=fftshift(G);M,N=size(G);a=;b=;T=1;i=1;format longfor k=:for u=1:Mfor v=1:NH(u,v)=(T/(pi*(u*a+v

26、*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=uint8(real(f);figure;imshow(f);title( '維納濾波的結(jié)果');e=f-uint8(g);SNR(i)=sum(sum(g.A2)/sum(sum(e-2);i=i+1;endidx=find(max(SNR)(題目2 維納濾波器MATLAB 版 )I=imread('');H=fspecial('motion' ,50,45);I1=imfilter(I,H, 'circular' ,'conv');figure(1);imshow(I1);title(&