基于非线性扩散滤波结构信息的图像去噪方法

基于非线性扩散滤波结构信息的图像去噪方法
张建伟;王译禾;陈允杰
【摘要】To retain more marginal information while denoising,by utilizing the rotational invariance of Laplacian kernel,a Lapacian ker-nel was proposed to fit different angles and the corresponding weighting function was put forward to make the nonlinear differential filter self-
adapted.Experimental results show that the proposed method is able to enhance the edge-preservation capability.%为在去除噪声的同时保持更多的边缘信息，利用拉普拉斯算子的旋转不变性提出一种适用于更多角度的拉普拉斯算子模板，根据模板适应角度的不同提出相应的权重函数，使非线性扩散滤波在原有的去噪基础上能够自适应地选取对应的算子对含噪图像进行扩散处理。

实验结果表明，该方法在去除噪声的同时可以保持更多的边缘信息。

【期刊名称】《计算机工程与设计》
【年(卷),期】2016(037)011
【总页数】8页(P2971-2978)
【关键词】非线性偏微分方程;图像去噪;拉普拉斯算子;梯度信息;边缘信息
【作者】张建伟;王译禾;陈允杰
【作者单位】南京信息工程大学数学与统计学院，江苏南京 210044;南京信息工程大学数学与统计学院，江苏南京 210044;南京信息工程大学数学与统计学院，江苏南京 210044
【正文语种】中文
【中图分类】TP391.41
许多学者在基于非线性偏微分方程的去噪领域取得了重要的研究成果：Perona等提出了经典的基于非线性扩散方程的P-M方程[1]；由于P-M方程的初值问题在
数学意义上是病态的，Catté等提出了改进方法[2]，解决了该方程理论上解的适
定性问题，可是扩散速度仍然较快，在演化过后边缘和细小结构仍会被模糊；Rudin等提出了模型以较好保持图像边缘特征[3]，但是该模型往往会在去噪的同
时使图像在灰度过度区域产生比较明显的“阶梯效应”。

基于非线性偏微分方程的扩散滤波对于边界保持的能力很大一部分取决于扩散方程是否能够完好的检测出边界并且控制扩散过程。

从站在保持边缘和去除噪声的角度，许多研究致力于提出新的扩散系数[4]，或改写扩散系数中的参数[5]，抑或是使用
其它的方法，其基本思想都是在探测到边界时，能够在一定程度上减小扩散项的大小，从而控制在边界处扩散的速度。

然而，这种对边界保持的能力越大就会导致去噪后的图像越容易出现“针孔现象”[6](即有个别噪声点增强)，因为在噪声比较强烈的区域可能会将噪声点误判为边界进行保持甚至于增强，不能达到理想的去噪效果，对视觉判断造成影响。

鉴于上述问题，本文提出一种基于结构信息的边界自适应去噪方法。

该方法利用拉普拉斯算子的旋转不变性提出适应于多种角度的拉普拉斯算子，在扩散过程当中，根据边界的不同角度采用不同的算子进行计算，使用图像的结构信息而不改变传统方法中的扩散系数，在良好去噪的同时，使多种角度的边缘得到与传统方法相比更明显的保护。

1.1 P-M方程[1]简介
经典的基于非线性偏微分方程的滤波器如下所示
式中表示L2范数，div是散度符号，c(·)是扩散系数，代表图像u的梯度。

该种去噪方法的特点是：c(·)为单调递减函数，自变量为0时，函数值取1；自变量趋近于无穷大时，函数值趋近于0。

也就是说，在图像的边界处梯度足够大时，c(·)下降到几乎等于零，扩散几乎被“停止”，边界得以保持；而在图像的同质区域，梯度值很小，c(·)取值接近于1，此时式(1)接近于热扩散方程，扩散相对剧烈，去除噪声。

Perona和Maik在其模型中给出了如下两个扩散系数
其中，K是给定的一个参数，控制扩散的程度，将式(1)展开，可以写成下式
式中表示拉普拉斯算子。

式(1)常用的离散格式[7]如下
式中：d∈{N,S,W,E}代表上下左右4个方向；τ是步长；n是迭代次数；ui,j代表图像上的各个点。

式(4)进一步展开，还可以写成
式中与uηη分别表示边界沿切向ξ与法向η的分量，对于不同的扩散函数，在切线方向的扩散始终是正向的，而沿法向的扩散有可能是正向的也有可能是逆向的，因为如果当时，c函数选用式(3)，即时，就会发生逆扩散。

因此，选取适当的K 值不仅可以平滑图像，也可能起到边缘增强的效果，但是如果在平滑区域发生了逆扩散，就有可能造成“阶梯效应”。

1.2 使用八方向拉普拉斯算子的去噪方法
在式(4)中，所代表的拉普拉斯变换通常是由Rosenfeld-Kak拉普拉斯算子模板[8]对u做卷积得到的，该拉普拉斯算子的模板表示为
文献[9]中，提及了不同方向的拉普拉斯算子模板，如下所示
为了使用当前点附近更多信息，结合这两种不同的拉普拉斯算子，经过计算和线性组合，得到了适应水平竖直边界，即可以比较完好保存水平竖直边界的拉普拉斯算子[10]，定义如下
起初，该种算子用来估计图像中噪声的大小，而如果将该算子运用到基于偏微分方程的去噪当中，会发现虽然在水平竖直边界处该算子的表现非常出色，可以很好地
保持水平竖直边界，但是在对角边界，会出现很严重的振铃现象，如图1所示。

之所以出现这种现象是因为L3是完全根据水平竖直边界被设计出来的。

L1与L2的线性组合可写为
根据式(8)，不同的α，β可以确定不同的L。

令L对水平边界与对角边界做卷积运算，边界如图2所示，u1，u2代表不同的灰度值，为像素点p编号1-9，假设没有被显示出来的边界与图中所示结构相同。

水平边界计算结果如下。

对角边界计算结果如下
计算的目的是为了在边界处使趋于0，降低扩散幅度。

在对角边界处，的解有且仅有α=β=0，因此不存在对角边界处几乎无扩散的L；在水平边界处，根据，解得α=-β，为计算方便，取α=4，得到L3。

因为L3只在水平竖直边界方向无扩散，所以会出现图1的现象。

针对这个问题Hajiaboli等在文献[11]中提出了适应于水平竖直及对角方向边界的八方向算子(EAPM模型)，其离散格式如下
式中：γ为给定的参数，f是一个用来调整扩散项的函数[11]，定义为
式中是文献[12]中提到的对梯度模值的一种非标准近似，定义为
可以看出，该迭代格式所使用的3×3大小的八方向拉普拉斯算子结合了水平竖直以及对角方向的信息，融合了L1和L2两种模板的优点，有效去除噪声的同时，在一定迭代次数内可以比较完好地保持图像的水平垂直及对角方向的边界。

与L3相比，使用该八方向算子可以降低图2中所出现的振铃现象的影响，将该算子运用到不同的非线性方程滤波当中，会比运用原L1算子的结果有不同程度的改进，处理之后的图片更加清晰，效果更好。

但是该扩散方程只能保持水平竖直及对角方向的边界，对其余边界的保持能力相对较差，并且需要人为赋值参数，无法做到根据图像本身的特点来确定扩散项的大小。

式(11)中提出的八方向的算子只使用了目标点周围的8个像素点信息，只能对水平
竖直及对角方向的边界得到较好的保持，对于其它角度的边界则无法达到如此好的效果，因此我们采用更多的信息，更大的模板来适应更多角度的边界。

与尺寸为
3×3的模板相比，5×5大小的模板可以适应更多的角度，当前点及周围各离散点
如图3所示。

按照定义
仿照L1的构造，使用水平竖直方向的信息，定义模板如下
使用对角方向信息得到的模板定义如下
可以看到，与L1和L2相似，它使用的同样是当前点水平竖直及对角方向的信息，不过与已有模板相比较，新定义的模板考虑的像素点更多。

为了对应更多的角度，考虑到拉普拉斯算子模板的旋转不变性，得到如下两个模板
由定义可知，La，Lb,Lc和Ld与式(11)中使用的两种模板相比，可以适应比水平
竖直及对角方向更多的角度。

而且式(11)中，模板的权重无法根据图像本身特点决定，需要人为给出。

因此，这里就需要根据图像本身特征构造恰当的权值函数使扩散项中使用的拉普拉斯算子符合图像边界的角度。

定义图像中每一点的梯度为矩阵A的元素，即
令每一个5×5模板对矩阵A分别做卷积。

本文设计的模板是为适应角度而存在的，图像边界处的梯度较大，当有某角度的边界存在时，那么这个角度所适应的模板对该边界上的点做卷积之后得到的值也应是较大的。

计算各个模板所对应角度的边界在总体边界中所占的比例就可以大致确定每个模板前的权重。

令La，Lb，Lc和
Ld分别与A做卷积运算，即定义
其余元素以此类推。

例如，计算La类型的边界在总体中所占比例，即计算，该值可作为La模板前的权值，其余模板前的系数同样用该方法计算。

将该权重计算式运用到每个拉普拉斯算子模板中，那么扩散式(11)就可以改写为
该扩散式可以从理论上适应16方向上的角度。

使用5×5大小的拉普拉斯模板，
可以尽可能使这些方向上的边界在扩散过程中得以很大程度的保持。

为检测该算子是否可以得到较好的结果，令式(20)与式(5)的P-M方程，式(11)的EAPM方法，分别对含有不同角度边界的图片进行扩散，用峰值信噪比(peak signal to noise ratio，PSNR)[13]来评价边界的完好程度
MSE(mean square error)为均方误差，即各数据误差平方的平均数；Iw为原始图像第w个点的像素值；Pw为经处理后的图像第w个点的像素值；W为图像拥有
的像素点的个数。

PSNR值越大，代表失真越少。

根据模板的旋转不变性，这里只抽样测试0°-45°的边界，3种方法中的扩散系数均采用式(3)的计算方法，取
K2=10，t=0.25，n=40，γ=1.2，测试图片大小均为128×128，PSNR结果保留四位小数，见表1。

由表1可以看出，EAPM方法保边效果在各个角度上均好于最初的P-M方法，而对于适应于更多角度的5×5模板来说，5×5模板的保边效果要好于3×3模板，直至角度减小至10°以下，5×5模板处理之后图像的PSNR数值比EAPM方法的结
果要小，也就是说在接近于水平边界处，5×5模板的处理结果并不十分理想，如
图4所示。

图4中的竖线代表一条竖直方向的边界，假设边界上的梯度大小均为1。

根据前文的估计，5×5大小模板在如图所示位置计算权重时，竖直边界应该更符合La，权
重更大，但实际上La与Lb的权重是相同的均为1/4，与预期结果不同，出现了
偏差。

前文提到，文献[10]中提出的L3可以完好的适应水平竖直边界使其得到最
大程度的保持，虽然使用L3会在处理对角边界时出现严重的振铃现象，但是本文只需要在水平竖直方向使用该算子。

这里需判断所处理的边界是否为水平竖直边界，本文引入结构张量信息[14,15]来判断水平竖直角度。

结构张量定义为
式中：J为一半正定矩阵，用于描述线性结构方向信息，其对应的角度为
这里规定，当边界的角度小于10°或大于80°时即认为该边界属于水平竖直角度一
类，在扩散项中使用L3算子，计算式为
其余角度均可由0°-90°的角度对称，按照前文中提出的5×5大小的模板进行计算。

当图像中的边界接近于水平竖直时，根据式(26)进行计算，使该种边界得到极大程度的保持；而当边界属于其它角度时，则使用适用于多角度的5×5算子模板式(20)进行计算，该方法可以自适应的调整拉普拉斯算子前的权值大小，根据边界的角度决定使用的算子，使拉普拉斯算子自适应实际边界，以期达到较好的保持效果。

图5为人工合成图像去噪实验结果。

图5(a)为513×469大小的人工合成图像加入均值为0，方差为200的高斯噪声后，使用P-M方程的去噪结果；图5(b)为该高斯噪声图使用EAPM的去噪结果；图5(c)为该高斯噪声图使用本文方法的去噪结果；扩散系数均选用式(3)，参数均取K2=10，t=0.25，n=80，γ=1.4；图
5(d)～图5(f)分别为3种方法处理结果图中的三角形右下角黑色框内角点放大图。

使用PSNR评价去噪效果，保留四位小数，见表2。

可以看到，在一定的迭代次数后，使用P-M方法处理并不能很好的保持图像的边
缘信息，导致了图像边缘的模糊；EAPM方法可以在一定迭代次数内较完好的保
持图像的边缘信息，可是在迭代次数逐渐增加时，有可能产生比较严重的马赛克现象，造成图像失真，如图5(e)所示；本文方法不仅没有快速的模糊边缘，反而在
边界处有所增强，如图5(f)所示。

在一定的迭代次数内，使用本文方法处理后结果的可以保持比EAPM方法处理更多的边界信息，PSNR数值也相对较高，见表2。

图6为Barbara图像去噪实验结果。

图6(a)是强纹理结构，大小为450×450的Barbara图添加均值为0，方差为60的高斯噪声后，使用P-M方程的去噪结果；图6(b)为该高斯噪声图使用EAPM方法的去噪结果；图6(c)为该高斯噪声图使用
本文方法的去噪结果图，扩散系数选用式(3)，参数均取K2=10，t=0.25，n=40，γ=1.3；图6(d)为EAPM的处理结果与原图像的差图；图6(e)为本文方法与原图
像的差图；图6(f)为本文方法与EAPM方法处理结果的差图。

使用PSNR来评价
去噪的效果，保留四位小数，见表3。

可以看出，虽然对于强纹理图片来说3种方法对纹理结构的保持都不尽如人意，但是与P-M和EAPM方法相比，本文方法所能保持的边缘信息还是要远远多于这两种方法的，图6(f)可以更加直观的看出EAPM与本文方法处理结果的差别，本文方法表现还是比较出色的，PSNR数值也较高，见表3。

图7～图10为辣椒图像去噪实验结果。

图7(a)是大小512×512，含有均值为0，方差为20高斯噪声的常用辣椒测试；图7(b)～图7(d)为P-M方程，EAPM及本文方法对图7(a)的处理结果；图8(a)含有均值为0，方差为50的高斯噪声；图
8(b)～图8(d)为P-M方程，EAPM及本文方法对图8(a)的处理结果；图9(a)为含有均值为0，方差为100的高斯噪声；图9(b)～图9(d)为P-M方程，EAPM及本文方法对9(a)的处理结果。

扩散系数选用式(3)，参数均取K2=10，t=0.25，
n=40，γ=1.3；图10(a)与图10(b)分别为图9(c)与图9(d)
与无噪原始图像的差；图10(c)与图10(d)分别为图10(a)与图10(b)的右下角白色方框内部分放大。

使用PSNR来评价去噪的效果，保留四位小数。

对比图10(c)与图10(d)，可以清晰的看到图10(c)中有很明显的细线结构，而图10(d)中的细线结构则较少且不明显，这就可以说明，使用EAPM方法所模糊掉的边界，使用本文方法处理时可以在一定的迭代次数内起到保护作用，而且本文方法的PSNR数值也相对较高，见表4。

本文根据非线性偏微分方程与结构信息，提出了一种基于非线性扩散滤波的去噪方法，该方法可以根据边界的结构信息确定用于扩散的拉普拉斯算子，以确保可以保持尽可能多的边界信息。

在该拉普拉斯算子被应用到非线性扩散滤波中时，某些特定角度给予特殊考虑，改善了该方法的一些不足之处。

实验结果表明本文提出的方法可以在有效去除噪声的同时保留尽可能多的边界。

同时，计算时间开销也随之增大，对于一些实时要求不高的应用场合，不失为一种选择。

本文实验中，使用该方
法并不能完整地保留所有的角点，并且在灰度过渡区域有可能会出现阶梯效应，这些都是今后需要研究和改进的重点。

【相关文献】
[1]Kamalaveni V,Anitha Rajalakshmi R,Narayanankutty KA.Image denoising using variations of Perona-Malik model with different edge stopping functions[C]//Procedia Computer Science,2015:673-682.
[2]Wang Y,Ren W,Wang H.Anisotropic second and fourth order diffusion models based on convolutional virtual electric field for image denoising[J].Computers & Mathematics with Applications,2013,66(10):1729-1742.
[3]Jiang L,Huang J,Lyu XG.Alternating direction method for the high-order total variation-based Poisson noise removal problem[J].Numerical Algorithms,2014,69(3):495-516. [4]Guo ZC,Sun JB,Zhang DZ,et al.Adaptive Perona-Malik model based on the variable exponent for image denoising[J].IEEE Transactions Image Processing,2012,21(3):958-967.
[5]Tsiotsios C,Petrou M.On the choice of the parameters for anisotropic diffusion in image processing[J].Pattern Recognition,2013,46(5):1369-1381.
[6]Srivastava A,Bhateja V,Tiwari H.Modified anisotropic diffusion filtering algorithm for MRI[C]//2nd International Conference on Computing for Sustainable Global Development,2015:1885-1890.
[7]Maiseli B,Elisha O,Mei J,et al.Edge preservation image enlargement and enhancement method based on the adaptive Perona-Malik nonlinear diffusion model[J].Iet Image Proces-sing,2014,8(12):753-760.
[8]Cheng FC,Huang SC,Ruan SJ.Scene analysis for object detection in advanced surveillance systems using Laplacian distribution model[J].IEEE Transactions on Systems Man & Cybernetics Part C Applications & Reviews,2011,41(5):589-598.
[9]Brox T,Rousson M,Deriche R,et al.Colour,texture,and motion in level set based segmentation and tracking[J].Image and Vision Computing,2010,28(3):376-390.
[10]James AP,Kavitha AP.Mean-variance blind noise estimation for CT
images[M].Switzerland:Advances in Signal Processing and Intelligent Recognition Systems.Springer International Publishing,2014:235-243.
[11]Hajiaboli MR,Omair Ahmad M,Wang CY.An edge-adapting Laplacian kernel for nonlinear diffusion filters[J].IEEE Transactions on Image Processing,2012,21(4):1561-1572.
[12]Wang Y,Niu R,Zhang L,et al.A scale-based forward-backward diffusion process for adaptive image enhancement and denoising[J].Eurasip Journal on Advances in Signal
Proces-sing,2011(2):1-19.
[13]Hore A,Ziou D.Image quality matrics:PSNR vs. SSIM[C]//IEEE International Conference on Pattern Recognition,2010:2366-2369.
[14]ZHONG Yong,WANG Lan.Method of 3D seismic image enhancement based on structure tensor[J].Computer Engineering and Design,2011,32(3):1010-1013(in Chinese).[钟勇,王岚.基于结构张量的三维地震图像增强方法[J].计算机工程与设
计,2011,32(3):1010-1013.]
[15]CHEN Yunjie,ZHANG Jianwei,WANG Shunfeng,et al.Brain MR image segmentation based on anisotropic Wells model[J].Journal of Computer Research and Development,2010,47(11):1878-1885(in Chinese).[陈允杰,张建伟,王顺凤,等.一种各向异性Wells算法脑核磁共振图像分割模型[J].计算机研究与发展,2010,47(11):1878-1885.]。