一种基于局部方向尺度的各向异性扩散图像滤波方法*

章东平，陈非予，童 超，何海波

(中国计量学院信息工程学院，杭州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.