一種基于局部方向尺度的各向異性擴散圖像濾波方法*

章東平，陳非予，童 超，何海波

(中國計量學院信息工程學院，杭州310018)

在圖像處理與分析中，如何在除去噪聲的同時更好的保護圖像的邊緣及紋理信息一直是極為重要的技術問題。由于Perona[1]提出的各向異性擴散濾波算法的性能優于那些傳統的圖像去噪算法[2]，導致了大量的學者研究圖像濾波的PDE方法。如Parker和Schnable[3]建議用一個基于熵的統計值與文獻[1]中提出的擴散傳導函數相乘。Edoardo[4]提出了使用邊緣驅動的自適應各向異性擴散濾波，指出在邊緣區域各向異性擴散應考慮邊緣方向的影響。Hye[5]為了濾除噪聲的同時更好的保持邊緣，提出了基于邊緣方向的各向異性擴散濾波法。他利用Prewitt算子得到邊緣方向，并且自適應的根據邊緣方向選擇擴散的權值。Chao[6]利用圖像中細小區域比噪聲背景擁有更高灰度級方差的特性，提出了一種基于修正的各向異性擴散邊緣保持平滑的方法。

此外研究者們在各向異性擴散基礎上還結合其他濾波去噪技術以獲得更好的濾波效果。由于小波變換的多分辨表示能力在圖像去噪的同時也能有效保留圖像細節[7]，Zhong[8]利用絕對小波系數的總和把小波系數分為規則和不規則的區域，然后在小波域內對不同層級使用不同的參數進行各向異性擴散。與之類似，Wang[9]提出了基于小波域內自適應歸一化濾波大小和步長的各向異性擴散。此外在Black[10]指出各向異性擴散可以視為魯棒估算的最小化能量過程之后，Couprise[11]等首次將用于圖像分割的能量分水嶺原理用于各向異性的能量最優化過程中來，代替之前的魯棒估算濾波模型。

最初由Perona提出的各向異性擴散濾波方法沒有使用任何形態或結構信息來控制不同區域擴散的范圍，從而使細小的圖像結構常會被除去，模糊的邊緣會更模糊。為了克服這些問題，Saha[12]引入了“局部尺度”的概念，在擴散濾波中，用局部尺度控制圖像中不同區域平滑的程度。根據Saha對局部尺度的定義，具有細小結構的區域和邊界附近的局部尺度較小，因此，在小區域內(對應于具有細小結構的區域或邊界附近)可采用小的參數進行濾波，而在大尺度區域內(對應于大的均勻區域的內部)使用大的參數濾波。

Saha算法由于利用了每個象素的局部結構信息而在性能上優于Perona算法，但Saha并沒有充分利用每個象素的局部結構信息。Saha在對每個象素尺度的定義中，每個象素的尺度信息只由區域一致性所決定的，實際上，有的象素在各個方向上形態和結構是不同的。根據圖像這一特征，本文重新定義每個象素的局部尺度，我們稱之為“局部方向尺度”。本文定義的局部方向尺度與Saha定義的局部尺度的主要不同之處是:本文定義的象素的局部方向尺度是有方向性的，即:每個象素在每個方向都具有一個局部尺度。對于二維圖像我們使用東、南、西、北、東南、西南、東北和西北八個方向，也就是說每個象素具有東、南、西、北、西南、東南、東北和西北八個方向的“局部方向尺度”。在擴散濾波中，本文采用局部方向尺度來控制擴散傳導函數在各個方向擴散的參數。

1 象素局部方向尺度計算

Saha算法的主要思想是在一個區域一致性標準下決定每個象素局部結構的尺寸，然后使用這個參數控制濾波的廣度。Saha的尺度概念最初在文獻[13]中引入。

(1)Saha局部尺度的計算:在二維圖像C中，中心點位于c(i，j)，半徑為k的圓Bk(c)定義為:

根據圓的定義，Saha定義一個函數FOk(c)，表示圓的邊界與c是否具有一致屬性(例如強度相似):

其中，|Bk(c)－Bk－1(c)|為 Bk(c)－Bk－1(c)中象素的數目，Wψ是一個一致性函數，根據圓的邊界與c具有一致屬性函數公式，可以確定二維圖像C中的任何象素c(c∈C)的局部尺度。象素c(c∈C)的局部尺度為當圓的邊界與c具有一致屬性時最大圓的半徑r(c)。

(2)局部方向尺度的計算:在Saha提出的算法中，使用圓內的所有象素來計算圓的邊界像素與像素c是否具有一致屬性。實際上，在圖像的邊緣附近，象素c在不同方向鄰域內的象素具有不同的屬性，這樣用圓內的所有象素來計算圓的邊界與c是否具有一致屬性就不很恰當。為了克服這一缺陷，本文提出用象素c鄰域內不同方向扇面內的象素來分別來計算不同方向扇面的邊界與c是否具有一致屬性。然后，根據各個方向一致屬性來計算象素在各個方向的局部尺度。

在二維圖像C中，中心點位于c(i，j)，半徑為k方向角為u的“扇面”Bk，u(c)定義為半徑為k，兩條半徑的方向分別為u+π/4和u－π/4的扇形內的區域，如圖1所示，圖中陰影部分是半徑為k方向角u為∠xOA的扇面，其中

圖1 扇面

根據“扇面”的定義，分別定義函數FOk，u(c)和函數FO'k，u(c)，來表示“扇面”的弧形邊界與c具有一致屬性的程度:

其中，|Bk，u(c)－Bk－1，u(c)|為 Bk，u(c)－Bk－1，u(c)中象素的數目，|Bk，u(c)－Bk－2，u(c)|為 Bk，u(c)－Bk－2，u(c)中象素的數目，Wψ是一個一致性函數，與Saha一樣我們使用高斯函數。

計算象素c的局部方向尺度的算法如下:

算法從k=1開始，逐步增加扇面的半徑，直到FOk，u(c)和 FO'k，u(c)小于閾值 ts(k)。FOk，u(c)和FO'k，u(c)小于閾值ts(k)時，認為u方向的扇面內含有了一個不同于c的對象，其中ts(k)是k的函數，它隨著k的增大而增大，它的取值范圍為[0.75，1]。ts(k)的下限取為0.75的理由是:在 c的3×3鄰域內，當象素c的四個鄰域象素中在某個方向的有一個屬性不同的點，但c與它的鄰域還是有可能屬于同一個對象。ts(k)的上限取為1的理由是:在k較大時，只有“扇面”的弧形邊界上的所在象素點的屬性都與象素點c相同時，才把半徑為k的“扇面”與c歸于同一個對象。本文自適應的取。

2 基于局部方向尺度的擴散濾波

在Perona提出的二個擴散傳導函數中，σ是決定濾波范圍的控制參數。Perona僅利用了圖像的全局噪聲估計來確定σ的值，沒有利用圖像中每個對象的局部結構，因而算法的自適應性差。針對Perona算法的這一缺陷，Saha引入局部尺度的概念，并用它來控制擴散傳導函數的參數σ，從而使算法能夠根據圖像的局部結構信息自適應地調整擴散的范圍。但是，Saha并沒有充分利用像素的局部結構信息，僅僅用一個參數來描述像素的局部形狀和結構，這是很不夠的。根據Saha算法的這一不足之處，本文提出了局部方向尺度的概念，用八個方向尺度來描述圖像的局部形狀和結構，并用這八個方向尺度來控制圖像強度在八個方向的擴散范圍。

我們的擴散傳導函數為:

式中，σs，u(c)，是用于控制象素 c 流向東、東北、北、西北、西、西南、南和東南八個方向流量的自適應參數，u∈{0，π/4，π/2，3π/4，π，5π/4，3π/2，7π/4}，σs，u(c)可通過下列公式求取:

式(7)中，rmax是圖像中局部尺度的上限，取為5，σψ是區域一致性函數Wψ的參數，reffu是用于控制象素c流向u方向流量的有效局部方向尺度，它由下面的算式求取:

式(8)中，e是在u方向上與象素c最相鄰的象素，象素d是在－u方向上與象素c最相鄰的象素。

根據上面的討論，可推導出，在第t次迭代中，象素c的流向u方向的流量Vu(c)為:

式(9)中，Ft，u(c)表示沿著u方向的梯度向量。

根據上面這些公式可推導出基于局部對象尺度擴散濾波的迭代處理可由下面的算式表示:

式(10)中KD是一個常量，它具有上限:

式(11)中，μα，u(c)是用于加權的常量。

3 實驗結果及分析

我們使用Canny直方圖估計方法選取參數σψ。Canny直方圖估計方法是先計算整幅圖像的梯度直方圖，σψ選取為從0開始累積梯度直方圖面積直到累積的面積與總面積的比值為閾值Tr時結束，取此時的梯度值為參數σψ的值，Tr=0.9是一個經驗值。

圖像中，同一象素沿不同方向的到邊緣的距離是不同的，采用基于局部方向尺度的濾波擴散可以在對象內更大的范圍內進行。因此，在邊緣附近本文的算法具有更強的除噪聲能力。

為了定性的評價本文提出的算法優越性，我們在標準測試圖像中加入了不同的噪聲，將本文提出的濾波方法與Perona各向異性擴散濾波方法和Saha基于尺度的擴散濾波方法進行比較。由圖2可見，基于尺度(Saha)和基于方向尺度(本文提出的算法)的濾波方法比Perona擴散濾波方法更能保持圖像中的細小結構。基于方向尺度(本文提出的算法)的擴散濾波方法與基于尺度(Saha)的擴散濾波方法相比較，前一種方法對噪聲的平滑能力和保持圖像中細小結構方面的能力都比后面一種方法更強。

圖2 三種濾波方法試驗結果的比較

圖2(a)為加入未歸一化的方差為σ2=100的高斯噪聲圖，圖2(b)為Perona算法5次迭代濾波的結果，圖2(c)為Saha算法5次迭代濾波的結果，圖2(d)為本文提出的算法5次迭代濾波的結果。圖2(e)為加入未歸一化的方差為σ2=1 000的高斯噪聲圖，圖2(f)為Perona算法5次迭代濾波的結果，圖2(g)為Saha算法5次迭代濾波的結果，圖2(h)為本文提出的算法5次迭代濾波的結果。

為定量地評價本文算法，我們對圖3加入具有不同方差的噪聲，其中，方差 σ2=100，110，120，…，200。

圖3 未加噪聲的原始測試圖像

實驗中，我們采用圖像的峰值信噪比作為評價標準，定義如下:其中，I(i，j)為未加入噪聲的原始圖像，ρ(i，j)為去噪后的圖像，圖像大小為M×N，圖像灰度級為0～255，PSNR越大，濾波去噪能力越強。實驗發現，隨著迭代次數的增加，各向異性平滑的輸出圖像的信噪比會逐步達到一個最大值，然后逐步下降(如圖4所示)。因此在評價中，當峰值信噪比達到最大值的時候，停止迭代.評價步驟如下:

(1)分別采用Perona算法，Saha算法和本文提出的算法處理噪聲污染圖像，得到輸出圖像ρPeronaM(i，j)，ρSahaM(i，j)和 ρZhangM(i，j);

(2)根據公式 12 計算 ρPeronaM(i，j)，ρSahaM(i，j)和 ρZhangM(i，j)的峰值信噪比;

(3)如果 PSNR(t)≥PSNR(t－1)，重復步驟(1)和(2);

(4)輸出PSNR=PSNR(t)。

圖4 隨著迭代次數的增加，各向異性平滑的輸出圖像的信噪比會到達一個最大值后逐漸下降

按照以上步驟，分別對圖3中加入不同方差的噪聲圖像進行處理，我們可以得到模型輸出圖像的峰值信噪比隨噪聲方差的變化曲線，如圖5所示。

圖5 輸出圖像的峰值信噪比隨噪聲方差變化曲線

為了較清晰顯示本文算法在紋理圖像的去噪效果，我們分別對圖6(a)原始紋理圖像，進行三種濾波方法的信噪比測試，結果如圖6(b)所示，本文方法在各測試方差噪聲下信噪比都是最好的。

圖6 三種濾波方法對紋理圖像除噪的結果比較

4 結論

本文針對在去噪聲的同時保護圖像邊緣及紋理信息的問題，提出了基于局部方向尺度的概念的擴散濾波方法。該方法首先本文根據象素鄰域內不同方向的象素來分別來計算邊界與所求像素是否具有一致屬性，然后，根據各個方向一致屬性來計算象素在各個方向的局部尺度，并用各個方向尺度來控制圖像強度在方向上的擴散濾波。在Lena測試圖像及紋理圖像上的實驗結果表明，所提出的算法的去噪效果不僅在圖像視覺效果上比較好，而且在輸出圖像信噪比參數上也取得了較好的結果，從而證明了本文所提出算法的有效性。

[1]Perona P，Malik J.Scale-Space and Edge Detection Using Anisotropic Diffusion[J].IEEE Transaction on Pattern Analysis and Machine Intelligence，1990，12(7):629－639.

[2]盧京晶，方中華，孫勝利.一種自適應的加權均值濾波器[J].傳感技術學報，2005，18(4):880－882.

[3]Parker G J M，Schnabel J A.Enhancement of Anisotropic Diffusive Filtering of MR Images Using Approximate Entropy[C]//Processing of International Society for Magnetic Resonance in Medicine，1999，175.

[4]Ardizzone E，Pirrone R，Gallea R，et al.Noise Filtering Using Edge-Driven Adaptive Anisotropic Diffusion[J].Computer-Based Medical Systems，2008:29－34.

[5]Hye Suk Kim，Yeong Joo Cho.Anisotropic Diffusion for Preserving Boundary-Edge[J].Proceedings ofthe12th International Conference on Advanced Communication Technology[C]//USA，2010，1693－1698.

[6]Shin-Min Chao，Du-Ming Tsai，Wei-Yao Chiu，et al.Anisotropic Diffusion-Based Detail-Preserving Smoothing for Image Restoration[C]//Proceedings of IEEE International Conference on Image Processing，2010.4145－4148.

[7]徐英，洪治.結合小波去噪的THz圖像多尺度增強算法研究[J].傳感技術學報，2011，24(3):398－401.

[8]Zhong J，Sun H.Wavelet-Based Multiscale Anisotropic Diffusion with Adaptive Statistical Analysis for Image Restoration[J].IEEE Transactions on Circuits and Systems—Ⅰ:Regular Papers，2008，55(9):2716－2725.

[9]Wang Zhiming，Bao Hong，Zhang Li.Image Denoising by Anisotropic Diffusion in Wavelet domain[C]//International Conference on Measuring Technology and Mechatronics Automation，2011，359－362.

[10]Black M J，Sapiro G，Marimont D H，et al.Robust Anisotropic Diffusion[J].IEEE Transaction on Image Processing，1998，7(3):421－432.

[11]Couprie C，Grady L，Najman L，et al.Anisotropic Diffusion Using Power Watersheds[C]//IEEE 17th International Conference on Image Processing.September 26－29，2010，Hong Kong，4153－4156.

[12]Saha P K，Udupa J K.Scale-Based Diffusive Image Filtering Preserving Boundary Sharpness and Fine Structure[J].IEEE Transaction on Medical Image，2001，20(11):1140－1155.

[13]Saha P K，Udupa J K，Odhner D.Scale-Based Fuzzy Connected Image Segmentation:Theory，Algorithms，and Validation[J].Computer Vision Image Understanding，2000，77(2):145－174.