王 磊，王 敏，张鹏程，任时磊，高晓玲，桂志国

(1. 中北大学 生物医学成像与影像大数据山西省重点实验室，山西 太原 030051；2. 临汾市人民医院，山西 临汾 041000)

0 引 言

数字图像是生活中重要的信息来源，然而，图像在产生、传输、使用的过程中，不可避免地会引入噪声，由于噪声往往与图像高频信息相关，因此，在保留重要特征的同时去除噪声是一项值得研究的工作. 图像去噪在图像增强、图像分割、目标识别等领域扮演着重要的角色[1-3]，如何在去除噪声的同时保留图像边缘和纹理细节，一直是研究的重点. 在图像空间域经典的去噪算法有双边滤波、偏微分方程[4]、非局部均值[5]、导向滤波[6]等. 其中，文献[4]提出的各向异性扩散方法也称为P-M模型，将像素值看作热流，把图像梯度引入扩散系数，在图像边缘处，扩散系数小，在平滑处，扩散系数大，从而有效保护了图像边缘，该模型自提出以来，就受到研究人员的广泛关注.

近年来，国内外众多研究学者基于P-M模型提出了很多改进算法. Ochotorena等[7]提出将各向异性扩散与加权导向滤波[8]结合，利用加权平均实现最大扩散，同时保留图像的强边缘，解决了细节光晕的问题. 方政等[9]提出基于多方向中值滤波的各向异性扩散算法，通过局部方差和图像梯度改进扩散系数，结合中值滤波算法，在滤波过程中注重对边缘细节的保持，提高了图像质量. Chao等[10，11]提出在各向异性扩散模型中针对不同类型图像，在模型中加入灰度方差，进一步保留了图像的纹理区域. 董婵婵等[12]在P-M模型的基础上，将差分曲率引入各向异性扩散方程中，较好区分了图像边缘、平坦区域和噪声，自适应扩散取得了较好的降噪效果. Hossein等[13]分析图像纹理的振荡属性，在残差图像中利用纹理信息构造纹理检测算子，结合梯度模值改进扩散系数，有效保护了图像纹理信息. 这些方法虽然有效改进了P-M模型，但处理结果仍存在受噪声影响大、图像细节易丢失等问题. 近年来，研究学者陆续提出新的去噪方法，非局部均值算法借鉴多幅相同图像叠加取均值的方法可有效去除噪声的原理，利用图像中大量的冗余信息，在图像搜索窗中寻找局部相似块，根据图像块相似程度确定去噪权重，达到较好的去噪效果，但其较长的运行时间和边缘易模糊的缺点会影响其大规模应用. 研究人员从各向异性扩散和非局部均值的优缺点出发，提出了很多改进算法. Yuan[14]提出在非局部均值模型中加入亮度、空间位置和梯度信息减弱P-M模型的阶梯效应，有效保护了图像细节. 陈强等[15]在扩散系数函数中，将像素的片相似代替图像梯度，同时扩展到彩色图像的去噪中，取得了较好效果，但依旧存在去噪不彻底、边缘易模糊等问题. Yang等[16]提出的基于非局部均值理论的各向异性扩散(NL P-M)算法研究P-M模型的不稳定性，在迭代过程中，先将噪声图像使用非局部均值算法处理，然后将处理后图像的梯度值代替噪声图像的梯度值，取得了较好的去噪效果，但其也存在过度平滑和效率低的问题. 为了提高非局部均值算法的效率，研究人员提出了很多解决方案，其中著名的有基于积分图[17-19]的方法，其先计算出像素的积分图，用加减运算代替图像块欧氏距离计算，以存储空间换计算时间的方式，降低算法的复杂度.

本文改进算法受文献[16]的启发，在P-M模型中，将非局部均值引入迭代过程中，利用图像中大量相似冗余信息有效地平滑图像和去除噪声，同时在扩散模型中，引入残差图像局部能量纹理检测算子，更好地保护了图像纹理和细节；针对非局部均值算法效率低的问题，采用积分图的方法降低算法复杂度，减少了算法运行时间.

1 算法原理

1.1 P-M模型及其不稳定性分析

经典的非线性P-M模型根据图像梯度模值自适应扩散，其迭代公式为

t∈(0,T),

(1)

(2)

g(|∇u|)=e-(|∇u|/k)2,

(3)

式中：k为梯度扩散阈值，k选择较大值时，|∇u|/k有较小值，g(|∇u|)趋于1，图像扩散趋于各向同性扩散，图像模糊，相反，k选择较小值时，g(|∇u|)趋于0，抑制图像扩散，但容易造成图像去噪不彻底的问题.

文献[16]研究P-M模型的不稳定性，将图像分解为梯度方向η和切线方向ξ，得到

(4)

(5)

则有

(6)

(7)

式中：uη和uηη分别为图像u在梯度方向的一阶和二阶偏导数；uxx，uyy，uxy分别为图像u(x,y,t)的二阶方向导数.类似得，u(x,y,t)在切线方向的二阶偏导数为

(8)

将式(1)分解为

[g(|∇u|)·ux]x+[g(|∇u|)·uy]y=

g(|∇u|)(uxx+uyy)=

g(|∇u|)(uηη+uζξ)=

[g′(|∇u|)|∇u|+g(|∇u|)]uηη+g(|∇u|)uζξ,

(9)

式中：g(|∇u|)永远为正值，意味着正向扩散和图像平滑.然而，g′(|∇u|)|∇u|+g(|∇u|)可以是正值也可以是负值，导致不稳定的结果. 再加上噪声的影响，各向异性扩散容易出现阶梯效应等问题.

P-M模型扩散的离散形式为

un+1=un+Δtdiv[g(|∇u|)∇u]=

un+Δtg(|∇u|)|∇u|，

(10)

式中：Δt为时间步长，其范围大小为0.05≤Δt≤0.25，∇u梯度可用有限差分表示为

(11)

且有

g(|∇u|)|∇u|=g(|∇Nui,j|)|∇Nui,j|+

g(|∇Sui,j|)|∇Sui,j|+g(|∇Wui,j|)|∇Wui,j|+

g(|∇Eui,j|)|∇Eui,j|.

(12)

1.2 非局部均值算法及积分图加速

非局部均值理论自提出以来，一直受到研究人员的广泛关注，其根据多幅相同图像叠加可有效去除图像中随机噪声的原理，利用图像中像素的邻域相似性，确定图像搜索窗内每个像素的去噪权值，计算公式为

(13)

式中：v为含噪声图像；uNL为去噪后的图像；权值w(i,j)表示像素i和像素j之间的相似性，其由像素i和像素j为中心的邻域像素决定，保证其相似性度量的准确性.

(14)

Z(x)为归一化系数，计算公式为

(15)

式中：h为平滑参数，其值越大图像越平滑，其中

(16)

i为图像邻域内索引.

假设图像中共有N个像素，搜索窗大小为D×D，计算像素相似度的匹配窗口大小为d×d(d=2ds+1)，ds为匹配窗邻域半径，因此，计算像素邻域相似度的时间复杂度为O(d2), 每个像素在搜索窗内共进行D2次计算，所以，非局部均值算法的时间复杂度为O(Nd2D2).由于非局部均值算法的时间复杂度较高，研究人员提出了积分图加速的方法.

在图像邻域求相似度中，定义积分图像为

(17)

其中：

(18)

两个像素领域的相似性可以用下式计算，大大减少运行时间.

St(x1-ds-1,x2-ds-1)-St(x1+ds,x2-ds-1)-

St(x1-ds-1,x2+ds)).

(19)

因此，积分图加速后，非局部均值算法的时间复杂度降为O(ND2).

2 本文算法

(20)

式(10)改进为

(21)

文献[13]为自适应选择扩散去噪算法，在各向异性扩散中，引入残差局部能量纹理检测算子，可以更准确地检测图像纹理和细节，同时降低阶梯伪影. 将残差图像定义为非局部均值处理结果与再经过P-M模型扩散结果的差值图像，如式(22)所示

(22)

PR=Ps+Pn,

(23)

式中：PR,Ps,Pn分别为uR,us,un的局部能量.因为Pn是近似恒定的，可以将其认为是PR的能量偏移强度，纹理信息Ps在PR中有较大的值.PR的计算公式为

PR(x,y)=

(24)

式中：μ(·)为均值算子；w(x,y)(x,y)为归一化的高斯算子，经过高斯运算，残差图中的噪声部分被平滑，纹理部分更加突出.

如图 1 所示，在残差局部能量图中，纹理区域对应图像中较大值，使纹理区域有了较好的提取，随着迭代次数的增加，局部能量算子可以有效地滤除噪声的影响，从而更准确地检测纹理细节.

图 1 局部能量及第一次和第二次迭代的残差局部能量

由于梯度算子往往较大，为了使局部能量的取值更为合理，使局部能量归一化到梯度范围

(25)

式中：PRScaled为纹理检测算子.

将纹理检测算子与梯度算子结合，即为文献[13]中提到的算法，在图像扩散过程中，可以提升图像保护纹理细节的能力，但容易出现去噪不彻底的问题. 本文结合上述两种算法的优点，在NL P-M算法中，改进扩散函数，添加局部能量纹理检测算子，提出改进的非局部均值各向异性扩散去噪算法，其基本流程为：

(26)

3) 根据式(24)计算残差图像的局部能量，并根据式(25)将局部能量归一化到图像梯度范围，其结果如图 1 所示.

(27)

5) 重复1)～4)，直到获取最终的去噪结果.

3 实验结果与分析

为验证所提算法的有效性，使用大小为512×512的标准灰度图像进行测试. 使用Matlab 2014a进行图像实验，计算机配置为：CPU为Intel(R) Core(TM) i7-4770 @3.50 GHz，内存为8 GB. 因为本文算法是P-M模型结合非局部均值算法的改进，所以采用经典的P-M模型算法、自适应选择扩散去噪算法和NL P-M算法为对比算法.

3.1 实验参数设置

本文算法参数较多，当面对不同类型图像和不同强度噪声时，需要适当手动调整以获取最佳的去噪结果. 对于扩散阈值k，取较大值时，容易使去噪图像平滑模糊，取较小值时，去噪不彻底. 当图像噪声强度较大时，应适当选择较大的扩散阈值. 平滑参数h在非局部均值滤波中发挥作用，其值较大时，容易使图像过度平滑，因此，当图像含有较多纹理细节时，其值应适当减小. 纹理控制阈值β在文献[13]中的范围为β≥1，β值越大，保留的纹理信息越多.高斯核参数σg为文献[13]中的参考值.在所提算法中，非局部均值中平滑参数h和纹理控制阈值β的值对算法有较大影响.

3.2 主观视觉效果分析

本文在标准灰度测试图像中加入均值为0，标准差为20的高斯白噪声，实验效果及局部放大图像如图 2 所示.

图2 添加标准差σ=20时，不同图像去噪结果

仔细观察图 2 及局部放大图像可以看到，经典P-M算法在强边缘处取得较好结果，在图像平滑处，容易受到噪声的影响，出现阶梯效应和少量斑点伪影；文献[13]所提算法中结合纹理检测算子，弱边缘和纹理处的效果都比较好，但在平滑处出现去噪不彻底的问题，视觉效果较差；NL P-M算法中，虽然视觉效果较好，但在纹理处出现过平滑的问题；本文算法有效克服了上述缺点，在平滑图像和保留纹理细节方面都有较好的表现，在去除噪声的同时保留了更多的纹理细节，比如在Lena和Barbara的放大图像中更好保留的纹理，Airplane图像中较为清晰的数字.

3.3 客观评价指标分析

在图像去噪领域，常用的客观评价标准有峰值信噪比(PSNR)和平均结构相似度(MSSIM)[21].

PSNR的定义如下

(28)

式中：m和n分别为图像行数和列数；u为恢复图像；u0为原始噪声图像.PSNR值越大，代表恢复结果越好.

结合图像亮度、对比度和结构信息，具有符合人眼视觉系统特性的结构相似度(SSIM)定义为

(29)

在实际评估中，常采用MSSIM来估计整幅图像的结构相似度，即为

(30)

式中：μx和μy分别为x和y的均值；σx和σy为标准差；σxy为协方差；J为局部窗口的总个数.MSSIM值越大，表示去噪结果越接近原始图像.

因主观视觉评价受观测者个体差异性较大，本节从PSNR和MSSIM两个客观指标来分析所提算法的去噪效果，结果如表1，表2 所示.

表1 添加高斯白噪声，经不同算法处理的PSNR

从表1 和表2 可以看出，在标准测试图像中添加标准差为20和30的高斯白噪声，本文所提算法相对P-M算法、文献[13]算法和NL P-M算法，都有最好的峰值信噪比和平均结构相似度，因此，从视觉效果和客观评价指标都验证了所提算法的有效性. 通过对不同噪声强度的去噪处理和选择最为常用的Lnea，Barbara，Airplane标准测试图像，验证了所提算法的有效性.

表2 添加高斯白噪声后经不同算法处理的MSSIM

3.4 算法时间复杂度分析

基于各向异性扩散的方法具有运行速度快的优势，但由于结合了非局部均值算法，使本文算法时间复杂度大大增加，因此，采用基于积分图方法的算法加速. 设图像共有N个像素，迭代次数为iter，下面分析对比算法与所提算法的时间复杂度.

P-M算法中用有限差分的方式求每个像素的梯度，并计算散度算子，因此，其时间复杂度为

T(n)=O(4iterN).

(31)

文献[13]在扩散函数中加入纹理检测算子，增加了纹理检测算子的计算时间，其时间复杂度为

T(n)=O(4iterN+N).

(32)

NL P-M算法在每次迭代中都有非局部均值算法，根据前面的分析，非局部均值的时间复杂度为O(Nd2D2)，因此，NL P-M的时间复杂度为

T(n)=O(iterND2d2+4iterN).

(33)

本文算法的时间复杂度和NL P-M算法的时间复杂度基本相同，但是由于采用了积分图加速方法，时间复杂度降低，此外，本算法和NL P-M算法迭代次数较少，算法耗时主要为非局部均值运行时长，时间复杂度为

T(n)=O(iterND2+4iterN).

(34)

所提算法结合矩阵操作，优化代码后可进一步减少运行时间. 为了更直观地观察各种算法的运算效率，以大小为512×512的Lena图像为例，将各种算法的运行时间记录如表3 所示，可见本文所提算法相比NL P-M算法，运行时间大幅度减少.

表3 Lena图像在标准差为20时，各种算法的运行时长

4 结 论

本文算法在P-M模型算法的基础上，将非局部均值算法引入迭代过程中，即在各向异性扩散前进行预处理，减弱了P-M模型的不稳定性；同时在扩散系数函数中融合残差图像局部能量算子，在扩散过程中，更准确地保护了纹理细节，提出一种改进的非局部均值各向异性扩散去噪算法，并且在非局部算法中，运用积分图的方法和使用矩阵运算代替传统像素遍历循环，极大减少了算法耗时，解决了NL P-M算法过度平滑和运行时间较长的问题. 从视觉效果和客观评价指标可以看出，所提算法优于P-M模型算法、文献[13]算法和NL P-M算法，在有效去除噪声的同时，更多地保留了纹理细节信息. 但是，本文主要针对在图像中最容易引入的高斯白噪声，对于其他类型的噪声去除，将是后继研究的重点.