空間域濾波器(實驗報告)

1、數(shù)字圖像處理作業(yè) 空間域濾波器 摘要在圖像處理的過程中,消除圖像的噪聲干擾是一個非常重要的問題。本文利用matlab軟件,采用空域濾波的方式,對圖像進行平滑和銳化處理。平滑空間濾波器用于模糊處理和減小噪聲，經(jīng)常在圖像的預處理中使用；銳化空間濾波器主要用于突出圖像中的細節(jié)或者增強被模糊了的細節(jié)。本文使用的平滑濾波器有中值濾波器和高斯低通濾波器，其中，中值濾波器對去除椒鹽噪聲特別有效，高斯低通濾波器對去除高斯噪聲效果比較好。使用的銳化濾波器有反銳化掩膜濾波、Sobel邊緣檢測、Laplacian邊緣檢測以及Canny算子邊緣檢測濾波器。不同的濾波方式，在特定的圖像處理應用中有著不同的效果和各自的優(yōu)

2、勢。1、分別用高斯濾波器和中值濾波器去平滑測試圖像test1和2，模板大小分別是3x3 ， 5x5 ，7x7；利用固定方差 sigma=1.5產(chǎn)生高斯濾波器. 附件有產(chǎn)生高斯濾波器的方法。實驗原理分析：空域濾波是直接對圖像的數(shù)據(jù)做空間變換達到濾波的目的。它是一種鄰域運算，其機理就是在待處理的圖像中逐點地移動模板，濾波器在該點地響應通過事先定義的濾波器系數(shù)與濾波模板掃過區(qū)域的相應像素值的關(guān)系來計算。如果輸出像素是輸入像素鄰域像素的線性組合則稱為線性濾波（例如最常見的均值濾波和高斯濾波），否則為非線性濾波（中值濾波、邊緣保持濾波等）?？沼驗V波器從處理效果上可以平滑空間濾波器和銳化空間濾波器：平滑空

3、間濾波器用于模糊處理和減小噪聲，經(jīng)常在圖像的預處理中使用；銳化空間濾波器主要用于突出圖像中的細節(jié)或者增強被模糊了的細節(jié)。模板在源圖像中移動的過程中，當模板的一條邊與圖像輪廓重合后，模板中心繼續(xù)向圖像邊緣靠近，那么模板的某一行或列就會處于圖像平面之外，此時最簡單的方法就是將模板中心點的移動范圍限制在距離圖像邊緣不小于（n1）/2個像素處，單處理后的圖像比原始圖像稍小。如果要處理整幅圖像，可以在圖像輪廓邊緣時用全部包含于圖像中的模板部分來濾波所有圖像，或者在圖像邊緣以外再補上一行和一列灰度為零的像素點（或者將邊緣復制補在圖像之外）。中值濾波器的設計：中值濾波器是一種非線性統(tǒng)計濾波器，它的響應基于圖

4、像濾波器包圍的圖像區(qū)域中像素的排序，然后由統(tǒng)計排序的中間值代替中心像素的值。它比小尺寸的線性平滑濾波器的模糊程度明顯要低，對處理脈沖噪聲（椒鹽噪聲）非常有效。中值濾波器的主要功能是使擁有不同灰度的點看起來更接近于它的鄰近值，去除那些相對于其鄰域像素更亮或更暗，并且其區(qū)域小于濾波器區(qū)域一半的孤立像素集。在一維的情況下，中值濾波器是一個含有奇數(shù)個像素的窗口。在處理之后，位于窗口正中的像素的灰度值，用窗口內(nèi)各像素灰度值的中值代替。例如若窗口長度為5，窗口中像素的灰度值為80、90、200、110、120，則中值為110，因為按小到大（或大到?。┡判蚝螅谌坏闹凳?10。于是原理的窗口正中的灰度值2

5、00就由110取代。如果200是一個噪聲的尖峰，則將被濾除。然而，如果它是一個信號，則濾波后就被消除，降低了分辨率。因此中值濾波在某些情況下抑制噪聲，而在另一些情況下卻會抑制信號。 將中值濾波推廣到二維的情況。二維窗口的形式可以是正方形、近似圓形的或十字形等。本次作業(yè)使用正方形模板進行濾波，它的中心一般位于被處理點上。窗口的大小對濾波效果影響較大。根據(jù)上述算法利用MATLAB軟件編程，對源圖像test1和test2進行濾波處理，結(jié)果如下圖： 可見，窗口的大小對濾波效果影響較大。窗口越大，平滑效果越明顯，圖像細節(jié)越模糊。 高斯濾波器的設計：高斯濾波是一種根據(jù)高斯函數(shù)的形狀來選擇模板權(quán)值的線性平滑

6、濾波方法。高斯平滑濾波器對去除服從正態(tài)分布的噪聲是很有效果的。一維零均值高斯函數(shù)為 。其中，高斯分布參數(shù)決定了高斯濾波器的寬度。對圖像來說，常用二維零均值離散高斯函數(shù)作平滑濾波器，函數(shù)表達式如下： （1）高斯函數(shù)具有5個重要性質(zhì)：（1）二維高斯函數(shù)具有旋轉(zhuǎn)對稱性，即濾波器在各個方向上的平滑程度是相同的。一般來說一幅圖像的邊緣方向是不知道的。因此，在濾波之前是無法確定一個方向比另一個方向上要更多的平滑的。旋轉(zhuǎn)對稱性意味著高斯濾波器在后續(xù)的圖像處理中不會偏向任一方向。（2）高斯函數(shù)是單值函數(shù)。高斯濾波器用像素鄰域的加權(quán)均值來代替該點的像素值，而每一鄰域像素點的權(quán)值是隨著該點與中心點距離單調(diào)遞減的。

7、這一性質(zhì)是很重要的，因為邊緣是一種圖像局部特征。如果平滑運算對離算子中心很遠的像素點仍然有很大的作用，則平滑運算會使圖像失真。（3）高斯函數(shù)的傅立葉變換頻譜是單瓣的。這一性質(zhì)是高斯函數(shù)傅立葉變換等于高斯函數(shù)本身這一事實的直接推論。圖像常被不希望的高頻信號所污染，而所希望的圖像特征，既含有低頻分量，又含有高頻分量。高斯函數(shù)傅立葉變換的單瓣意味著平滑圖像不會被不需要的高頻信號所污染，同時保留了大部分所需要的信號。（4）高斯濾波器的寬度(決定著平滑程度)是由參數(shù)表征的，而且和平滑程度的關(guān)系是非常簡單的。越大，高斯濾波器的頻帶就越寬，平滑程度就越好。通過調(diào)節(jié)平滑程度參數(shù)，可在圖像特征分量模糊(過平滑)

8、與平滑圖像中由于噪聲和細紋理所引起的過多的不希望突變量(欠平滑)之間取得折衷。（5）由于高斯函數(shù)的可分離性，大高斯濾波器可以有效實現(xiàn)。通過二維高斯函數(shù)的卷積可以分兩步來進行，首先將圖像與一維高斯函數(shù)進行卷積，然后將卷積的結(jié)果與方向垂直的相同一維高斯函數(shù)進行卷積。因此，二維高斯濾波的計算量隨濾波模板寬度成線性增長而不是成平方增長。這些性質(zhì)使得它在早期的圖像處理中特別有用，表明高斯平滑濾波器無論在空間域還是在頻率域都是十分有效的低通濾波器。根據(jù)上述分析，利用MATLAB軟件設計高斯濾波器，對源圖像test1和test2進行濾波處理，結(jié)果如下圖： 可見，對于高斯濾波器，模板的大小對濾波效果影響不大。

9、高斯濾波雖然能夠在一定程度上去掉噪聲，但也使得圖象變得模糊不清，效果并不能令人滿意。 2、利用高通濾波器濾波測試圖像test3,4：包括unsharp masking, Sobel edge detector, and Laplace edge detection；Canny algorithm.實驗原理分析：銳化濾波能減弱或消除圖像中的低頻率分量，但不影響高頻率分量。因為低頻分量對應圖像中灰度值緩慢變化的區(qū)域，因而與圖像的整體特性，如整體對比度和平均灰度值等有關(guān)。銳化濾波將這些分量濾去可使圖像反差增加，邊緣明顯。在實際應用中，銳化濾波可用于增強被模糊的細節(jié)或者低對比度圖像的目標邊緣。圖像銳化

10、的主要目的有兩個：一是增強圖像邊緣，使模糊的圖像變得更加清晰，顏色變得鮮明突出，圖像的質(zhì)量有所改善，產(chǎn)生更適合人眼觀察和識別的圖像；二是希望經(jīng)過銳化處理后，目標物體的邊緣鮮明，以便于提取目標的邊緣、對圖像進行分割、目標區(qū)域識別、區(qū)域形狀提取等，為進一步的圖像理解與分析奠定基礎。由于銳化使噪聲受到比信號還要強的增強，所以要求銳化處理的圖像有較高的信噪比；否則，銳化后圖像的信噪比更低。反銳化掩膜圖像增強（unsharp masking） 圖像的反銳化掩蔽算法可以表示為： (1)其中fs(x,y)表示經(jīng)過反銳化掩蔽得到的銳化圖像，是f(x,y)的模糊形式。反銳化掩蔽進一步的普遍形式稱為高提升濾波。在

11、圖像中任何一點(x,y)處，高提升濾波后的圖像可定義如下： (2)其中A1，與前式一樣，是的模糊形式，此式也可以寫成：(3)結(jié)合式(1)，可以得到： (4)這一表達式可計算高提升濾波圖像。如果選擇拉普拉斯變換，式(4)變成： (5)高提升濾波處理可以通過任何一個圖1所示的掩模得以實現(xiàn)。當A=1時，高提升濾波處理就是標準的拉普拉斯變換。隨著A超過l不斷增大，銳化處理的效果越來越不明顯。最終，當A足夠大時，高提升圖像將近似等于經(jīng)常數(shù)調(diào)制的圖像。 圖1 高頻提升濾波技術(shù)可以用其中一種掩膜來實現(xiàn)（） 本文采用的反銳化掩膜濾波器中A=1。根據(jù)以上分析，利用MATLAB軟件設計反銳化掩膜濾波器對test3

12、和test4進行濾波，結(jié)果如下圖： 可見，經(jīng)過反銳化掩膜濾波后，圖像的邊緣得到了增強，細節(jié)更加明顯，但同時也帶來了高頻噪聲的影響。索貝爾邊緣檢測（Sobel edge detector） 索貝爾算子（Sobel operater）主要用作邊緣檢測，在技術(shù)上，它是一離散性差分算子，用來運算圖像亮度函數(shù)的灰度之近似值。在圖像的任何一點使用此算子，將會產(chǎn)生對應的灰度矢量或其法矢量。該算子包含兩組3x3的矩陣，分別為橫向及縱向，將之與圖像作平面卷積，即可分別得出橫向及縱向的亮度差分近似值。如果以A代表原始圖像，Gx及Gy分別代表經(jīng)橫向及縱向邊緣檢測的圖像灰度值，其公式如下： 圖像的每一個像素的橫向及縱

13、向梯度近似值可用以下的公式結(jié)合，來計算梯度的大小。然后可用以下公式計算梯度方向。如果角度等于零，即代表圖像在該處擁有縱向邊緣，左方較右方暗。 根據(jù)以上算法分析，利用MATLAB軟件設計Sobel邊緣檢測濾波器，對源圖像test3和test4進行濾波，結(jié)果如下圖示： Laplacian邊緣檢測（Laplacian edge detector）拉普拉斯算子是最簡單的各向同性微分算子，具有旋轉(zhuǎn)不變性。一個二維圖像函數(shù)的拉普拉斯變換是各向同性的二階導數(shù)，定義為： （6）為了更適合于數(shù)字圖像處理，將該方程表示為離散形式： （7）另外，拉普拉斯算子還可以表示成模板的形式，如圖2所示。圖2（a）表示離散拉普

14、拉斯算子的模板，圖2（b）表示其擴展模板，圖2（c）和（d）則分別表示其他兩種拉普拉斯的實現(xiàn)模板。從模板形式容易看出，如果在圖像中一個較暗的區(qū)域中出現(xiàn)了一個亮點，那么用拉普拉斯運算就會使這個亮點變得更亮。因為圖像中的邊緣就是那些灰度發(fā)生跳變的區(qū)域，所以拉普拉斯銳化模板在邊緣檢測中很有用。一般增強技術(shù)對于陡峭的邊緣和緩慢變化的邊緣很難確定其邊緣線的位置。但該算子卻可用二次微分正峰和負峰之間的過零點來確定，對孤立點或端點更為敏感，因此特別適用于以突出圖像中的孤立點、孤立線或線端點為目的的場合。同梯度算子一樣，拉普拉斯算子也會增強圖像中的噪聲，有時用拉普拉斯算子進行邊緣檢測時，可將圖像先進行平滑處理

15、。0101111-411-81010111（a）拉普拉斯運算模板 （b）拉普拉斯運算擴展模板0-10-11-1-14-118-10-10-11-1（c）、（d）拉普拉斯其他兩種模板圖2 拉普拉斯的4種模板拉普拉斯銳化的基本方法可以由下式表示：這種簡單的銳化方法既可以產(chǎn)生拉普拉斯銳化處理的效果，同時又能保留背景信息，將原始圖像疊加到拉普拉斯變換的處理結(jié)果中去，可以使圖像中的各灰度值得到保留，使灰度突變處的對比度得到增強，最終結(jié)果是在保留圖像背景的前提下，突現(xiàn)出圖像中小的細節(jié)信息。 根據(jù)以上算法分析，利用MATLAB軟件設計Laplacian邊緣檢測濾波器，對源圖像test3和test4進行濾波，

16、結(jié)果如下圖示：由上圖可知，運用Laplacian算子對test4 copy.bmp進行邊緣檢測后，圖像的邊緣信息沒有被很好的檢測出來。這是由于拉普拉斯梯算子會增強圖像中的噪聲，因此本文再用拉普拉斯算子對test4 copy.bmp進行邊緣檢測時，先將圖像進行高斯平滑濾波處理。處理結(jié)果如下圖： Canny算子邊緣檢測（Canny algorithm）在圖像邊緣檢測中，抑制噪聲和邊緣精確定位是無法同時滿足的。邊緣檢測算法通過平滑濾波去除圖像噪聲的同時，也增加了邊緣定位的不確定性；反之，提高邊緣檢測算子對邊緣敏感性的同時，也提高了對噪聲的敏感性。Canny算子力圖在抗噪聲干擾和精確定位邊緣之間尋求最

17、佳折中方案。用Canny算子檢測圖像邊緣的步驟如下： step1:用高斯濾波器平滑圖象； step2:用一階偏導的有限差分來計算濾波后圖像梯度的幅值和方向； step3:對梯度幅值進行非極大值抑制，其過程為找出圖像梯度中的局部極大值點，把其他非局部極大值點置零以得到細化的邊緣。 step4:用雙閾值算法檢測和連接邊緣。使用兩個閾值T1和T2（T1T2），T1用來找到每條線段，T2用來在這些線段的兩個方向上延伸尋找邊緣的斷裂處，并連接這些邊緣。利用上述原理和MATLAB軟件設計濾波器，對test3和test4進行邊緣檢測，結(jié)果如下圖所示： 可見，經(jīng)過Canny算子較好的檢測出了圖像的所有邊緣信息

18、，具有較強的去噪能力。附錄一、 參考文獻1 岡薩雷斯著.數(shù)字圖像處理(第三版).北京：電子工業(yè)出版社，20102 楊杰 李慶著.數(shù)字圖像處理及MATLAB實現(xiàn)學習與實驗指導.北京：電子工業(yè)出版社，20103 蘇金明 王永利著.MATLAB圖形圖像. 北京：電子工業(yè)出版社，20054 朱習軍 隋思漣等著.MATLAB在信號與圖像處理中的應用. 北京：電子工業(yè)出版社，20095 張秀蘭著.基于MATLAB的數(shù)字圖像的邊緣檢測.吉林化工學院學報，20106 楊先花 黎粵華著.基于MATLAB圖像邊緣檢測算法效果對比.機電產(chǎn)品開發(fā)與創(chuàng)新，2010二、源代碼：第一題(1) 平滑濾波器1. 中值濾波器（以

19、3x3的模板大小為例） I=imread( E:大三下圖像處理英文課件作業(yè)第四次作業(yè)test1.pgm,pgm);n=3; 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; endendI2=uint8(x2);imshow(I2)title(中值濾波后的test1.pgm（3x3）)I=

20、imread(E:大三下圖像處理英文課件作業(yè)第四次作業(yè)test2.tif,tif);n=3;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; endendI2=uint8(x2);imshow(I2)title(中值濾波后的test2.tif（3x3）)2、高斯濾波器（以3x3的模板大小為

21、例） n1=3;sigma1=1.5;n2=3;sigma2=1.5;theta=0; I,map=imread( E:大三下圖像處理英文課件作業(yè)第四次作業(yè)test1.pgm,pgm); r=cos(theta) -sin(theta); sin(theta) cos(theta); for i = 1 : n2 for j = 1 : n1 u = r*j-(n1+1)/2 i-(n2+1)/2; h(i,j)=exp(-u(1)2/(2*sigma12)/(sigma1*sqrt(2*pi)*exp(-u(2)2/(2*sigma22)/(sigma2*sqrt(2*pi); end en

22、d h = h / sqrt(sum(sum(h.*h); f1=conv2(I,h,same); figure(1); imagesc(I);title(test1.pgm); colormap(gray); figure(2);imagesc(f1);title( 高斯濾波平滑后的test1.pgm(3x3); colormap(gray); n1=3;sigma1=1.5;n2=3;sigma2=1.5;theta=0; I,map=imread( E:大三下圖像處理英文課件作業(yè)第四次作業(yè)test2.tif,tif); r=cos(theta) -sin(theta); sin(thet

23、a) cos(theta); for i = 1 : n2 for j = 1 : n1 u = r*j-(n1+1)/2 i-(n2+1)/2; h(i,j)=exp(-u(1)2/(2*sigma12)/(sigma1*sqrt(2*pi)*exp(-u(2)2/(2*sigma22)/(sigma2*sqrt(2*pi); end end h = h / sqrt(sum(sum(h.*h); f1=conv2(I,h,same); figure(1); imagesc(I);title(test2.tif); colormap(gray); figure(2);imagesc(f1);

24、title( 高斯濾波平滑后的test2.tif(3x3); colormap(gray); （2）銳化濾波器：反銳化掩膜（以處理test3_corrupt.pgm為例）pic = imread(E:大三下圖像處理英文課件作業(yè)第四次作業(yè)test4 copy.bmp,bmp);I=pic;picSize = size(pic);pic = cast(pic,int32);t = zeros(picSize(1) + 2 , picSize(2) + 2);t = cast(t,int32); t(2:picSize(1) + 1 , 2:picSize(2) + 1) = pic;t(: , 1

25、) = t(: , 2);t(: , picSize(2) + 2) = t(: , picSize(2) + 1);t(1 , :) = t(2 , :);t(picSize(1) + 2 , :) = t(picSize(1) + 1 , :);A = 1.0;for i=2:1:picSize(1)+1 for j=2:1:picSize(2)+1 pic(i-1,j-1) = t(i,j)*A - (t(i,j)*(-8) + t(i,j-1) + t(i-1,j-1) + t(i-1,j) + t(i,j+1) + t(i+1,j) + t(i+1,j+1) + t(i+1,j-1)

26、 + t(i-1,j+1); endendpic = cast(pic,uint8);imwrite(pic,UnsharpMasking_A_, num2Str(10*A),.png,png);figure;subplot(1,2,1)imshow(I);title(原圖像test4 copy.bmp);subplot(1,2,2)imshow(pic);title(反銳化掩膜后的test4 copy.bmp);Sobel 邊緣檢測（以處理test4_corrupt.pgm為例）I= imread(E:大三下圖像處理英文課件作業(yè)第四次作業(yè)test4 copy.bmp,bmp); N,M=si

27、ze(I)I=double(I);h1=-1,0,1;-2,0,2;-1,0,1;h2=-1,-2,-1;0,0,0;1,2,1;Gx=conv2(I,h1,same);Gy=conv2(I,h2,same);F=abs(Gx)+abs(Gy);for i=1:N for j=1:M I(i,j)=F(i,j); endendI=uint8(I);figure;imshow(I);title(Sobel 邊緣檢測test4 copy.bmp);Laplace邊緣檢測a）Laplace邊緣檢測test3_corrupt.pgmI = imread(E:大三下圖像處理英文課件作業(yè)第四次作業(yè)test

28、3_corrupt.pgm,pgm);N,M=size(I);figure;subplot(1,2,1);imshow(I);title(test3_corrupt.pgm);t=ones(N,M);t=I;for i=2:N-1 for j=2:M-1 I(i,j)=t(i+1,j)+t(i-1,j)+t(i,j+1)+t(i,j-1)-4*t(i,j); endendsubplot(1,2,2);imshow(I);title(Laplacian邊緣檢測-test3_corrupt.pgm);b）Laplace邊緣檢測test4 copy.bmp未進行高斯濾波：I= imread(E:大三

29、下圖像處理英文課件作業(yè)第四次作業(yè)test4 copy.bmp,bmp); N,M=size(I);figure;subplot(1,2,1);imshow(I);title(源圖像test4 copy.bmp);t=ones(N,M);t=I;for i=2:N-1 for j=2:M-1 I(i,j)=t(i+1,j)+t(i-1,j)+t(i,j+1)+t(i,j-1)-4*t(i,j); endendsubplot(1,2,2);imshow(I);title(Laplacian邊緣檢測(未濾波)-test4 copy.bmp);加入高斯濾波：I= imread(E:大三下圖像處理英文課

30、件作業(yè)第四次作業(yè)test4 copy.bmp,bmp); N,M=size(I);figure;subplot(1,2,1);imshow(I);title(源圖像test4 copy.bmp );n1=3;sigma1=0.5;n2=3;sigma2=0.5;theta=0; r=cos(theta) -sin(theta); sin(theta) cos(theta); for i = 1 : n2 for j = 1 : n1 u = r*j-(n1+1)/2 i-(n2+1)/2; h(i,j)=exp(-u(1)2/(2*sigma12)/(sigma1*sqrt(2*pi)*exp

31、(-u(2)2/(2*sigma22)/(sigma2*sqrt(2*pi); end end h = h / sqrt(sum(sum(h.*h); f1=conv2(I,h,same); t=ones(N,M);t=f1;for i=2:N-1 for j=2:M-1 f1(i,j)=t(i+1,j)+t(i-1,j)+t(i,j+1)+t(i,j-1)-4*t(i,j); endendsubplot(1,2,2);imshow(f1);title( Laplacian邊緣檢測(高斯濾波后)-test4 copy.bmp ); Canny算子邊緣檢測（以處理test4 copy.bmp為例

32、）I= imread(E:大三下圖像處理英文課件作業(yè)第四次作業(yè)test4 copy.bmp,bmp); a = im2single(I);m,n = size(a);e = false(m,n);GaussianDieOff = 0.0001; PercentOfPixelsNotEdges = 0.7;ThresholdRatio = 0.4; sigma = 1;thresh = ;pw = 1:30; ssq = sigma2;width = find(exp(-(pw.*pw)/(2*ssq)GaussianDieOff,1,last);if isempty(width) width

33、= 1; end t = (-width:width); gau = exp(-(t.*t)/(2*ssq)/(2*pi*ssq); x,y=meshgrid(-width:width,-width:width); dgau2D=-x.*exp(-(x.*x+y.*y)/(2*ssq)/(pi*ssq); aSmooth=imfilter(a,gau,conv,replicate); aSmooth=imfilter(aSmooth,gau,conv,replicate); ax = imfilter(aSmooth, dgau2D, conv,replicate); ay = imfilte

34、r(aSmooth, dgau2D, conv,replicate); mag = sqrt(ax.*ax) + (ay.*ay); magmax = max(mag(:); if magmax0 mag = mag / magmax; end if isempty(thresh) counts=imhist(mag, 64); highThresh = find(cumsum(counts) PercentOfPixelsNotEdges*m*n,. 1,first) / 64; lowThresh = ThresholdRatio*highThresh; thresh = lowThresh highThresh; elseif length(thresh)=1 highThresh = thresh; if thresh=1 eid = sprintf(Im