0 引言

数字图像在传输过程中容易受到椒盐噪声的污染, 从而导致图像质量的下降, 不利于感兴趣目标的提取及其他相关的应用。因此, 如何高效地滤除椒盐噪声, 同时有效地保护图像的细节信息成为图像处理领域一个重要的基本问题。

为此, 众多降噪方法被提出用于滤除图像中的椒盐噪声。其中经典的中值滤波器(Standard Median Filter, SMF)^[1]具有较好的椒盐噪声滤除效果和较高的执行效率, 但该算法没有区分噪声像素点与非噪声像素点, 对图像中所有的像素点均采用相同的滤波窗进行修正, 其滤波过程存在盲目性, 容易导致结果图像产生模糊效应。针对该问题, 学者们提出了综合噪声检测和图像滤波的先进去噪算法。例如, Zhang等^[2]提出了一种基于卷积的噪声检测方法, 该算法能够有效地去除低密度椒盐噪声, 且随着噪声密度的增大, 其去噪能力也随之急剧下降。Chan等^[3]则提出了使用正则化方法去除椒盐噪声的中值滤波算法, 具有较好的去噪及细节保护的性能, 但算法的时间复杂度较高。而Chang等^[4]采用了自适应中值滤波器(Adaptive Median Filter, AMF)进行去噪, 尽管该方法比SMF方法具备更好的去噪性能, 但不适用于高密度噪声和细节丰富的图像场景。Wang等^[5]通过对一维拉普拉斯算子作卷积运算, 提出了一种新的噪声检测及恢复算法, 适用于包含较多边缘细节的图像, 且能够有效地去除大范围的椒盐噪声, 但需要设置一些外部值。文献 [6]同样提出了包含噪声检测及恢复的自适应模糊开关中值滤波算法, 弥补了中值滤波算法固定滤波窗的不足, 在噪声密度低于50%的场景中具有良好的去噪性能, 但算法的执行时间随着噪声密度的递增而显著增大。Fabijanska等^[7]则在改进噪声检测方法的基础上, 根据噪声密度自适应地选择不同大小的中值滤波器, 能够滤除各种密度的椒盐噪声, 但算法在自适应扩展滤波窗的过程中所引入的非噪声像素点可能远离邻域中心的噪声像素点, 容易产生模糊效应。为解决自适应滤波窗容易产生模糊效应的问题, Lu等^[8]提出了改进型的方向加权中值滤波(Modified Directional-Weighted-Median Filter, MDWMF)算法, 但该算法所判决的最佳方向并不一定是真实的边缘方向, 且简单地将灰度为0或255作为条件来剔除近邻噪声像素点的做法可能会剔除真实边缘像素点, 算法存在一定的局限性, 且不能适应高密度噪声的场景。为去除高密度椒盐噪声, Vijaykumar等^[9]则提出了快速开关中均值滤波(Fast Switching based Median-Mean Filter, FSMMF)算法, 具有较低的时间复杂度和较好的细节保护性能, 但算法仍然存在噪声误检问题以及固定滤波窗的不足之处。尽管上述基于中值滤波的改进型去噪方法能够较好地去除椒盐噪声, 但该类方法在噪声恢复阶段仅利用了局部邻域内像素点的排序信息, 而忽略了图像的局部结构特征, 容易引起图像细节信息的丢失, 保护图像细节的能力仍然有限。为此, Wu等^[10]基于噪声检测器和整体变分修补模型提出了另一种去除椒盐噪声的非中值滤波算法, 能够有效地去除噪声, 且具备一定的细节保护性能, 但算法复杂、执行效率较低。同样, Zhang等^[11]也提出了去除椒盐噪声的图像修补(Image Inpainting, II)算法, 该算法根据图像局部区域的方向特征自适应地使用四个不同方向的滤波窗进行迭代卷积滤波, 提高了去噪结果图像的视觉质量, 但在高密度噪声污染的情况下, 需要迭代执行多次卷积滤波才能实现较好的去噪效果, 算法的执行时间较长。

综合上述算法的优缺点, 针对椒盐噪声的滤除和图像细节保护问题, 本文在引入一种高效脉冲检测器对椒盐噪声进行精确检测的基础上, 应用核回归拟合方法对所检测的噪声像素点进行去噪处理, 提出一种基于核回归拟合的开关去噪算法。实验结果表明, 本文算法不仅能够适应各种噪声密度的应用需求, 具有良好的去噪效果, 而且能够有效地保护图像的细节信息。

1 开关核回归拟合去噪算法

噪声检测作为采用开关策略进行去噪处理的关键技术, 其误检率和漏检率直接影响去噪的效果。基于椒盐噪声点通常取值为某个灰度动态范围的最大值或最小值的思想, 在噪声检测阶段引入一种高效的脉冲检测器^[12], 对图像中的椒盐噪声点进行精确检测, 有效地降低噪声检测过程中的误检率和漏检率。而在噪声恢复阶段, 首先将所检测到的噪声像素点当作缺失数据, 应用核回归方法对以噪声像素点为中心的邻域内的非噪声像素点进行拟合, 得到符合图像的局部结构特征的核回归拟合曲面;最后, 以噪声像素点的空间坐标对核回归拟合曲面进行重采样, 获得噪声像素点恢复后的灰度值, 实现对椒盐噪声的有效滤除。

1.1 椒盐噪声的精确检测

对于图像中坐标为(x, y)、灰度值为f (x, y)的任意像素点p, 采用高效脉冲检测器对其进行噪声检测的过程具体如下:

首先, 令W_{x, y}表示以像素点p为中心、大小为3×3的检测窗, 如式(1) 所示:

${{W}_{x,y}}=\left\{ {{p}_{i,j}}\left| x-1\le i\le x+1,y-1\le j\le y+1 \right. \right\}$

(1)

其次，假设MaxinW_{x, y}、MininW_{x, y}分别为当前检测窗W_{x, y}中像素点的最大、最小灰度值, 则从第1个检测窗W_{0, 0}至当前检测窗W_{x, y}的最大灰度值Max_{x, y}、最小灰度值Min_{x, y}可分别采用式(2) 、(3) 进行定义:

$\begin{eqnarray*}Ma{{x}_{x,y}}=\left\{ \begin{array}{*{35}{l}} Ma{{x}_{x,y-1}}\text{, }&Ma{{x}_{x,y-1}}\ge Maxin{{W}_{x,y}} \\ Maxin{{W}_{x,y}},&\text{ 其他} \\ \end{array} \right.\end{eqnarray*}$

(2)

$\begin{eqnarray*}Mi{{n}_{x,y}}=\left\{ \begin{array}{*{35}{l}} Mi{{n}_{x,y-1}}\text{, }&Mi{{n}_{x,y-1}}\le Minin{{W}_{x,y}} \\ Minin{{W}_{x,y}},&\text{ 其他} \\ \end{array} \right.\end{eqnarray*}$

(3)

再者, 由于椒盐噪声像素点通常取值为某个灰度动态范围的最大值或最小值, 因此可采用如式(4) ~(5) 所示的N_max、N_min分别作为当前检测窗W_{x, y}内的盐噪声灰度值、椒噪声灰度值的估计:

$\begin{eqnarray*}{{N}_{\max }}=\left\{ \begin{array}{*{35}{l}} Ma{{x}_{x,y}}\text{, }&Ma{{x}_{x,y}}=Ma{{x}_{x,y-1}} \\ 255,&\text{ 其他} \\ \end{array} \right.\end{eqnarray*}$

(4)

$\begin{eqnarray*}{{N}_{\min }}=\left\{ \begin{array}{*{35}{l}} Mi{{n}_{x,y}}\text{, }&Mi{{n}_{x,y}}=Mi{{n}_{x,y-1}} \\ 255,&\text{ 其他} \\ \end{array} \right.\end{eqnarray*}$

(5)

最后, 依据式(4) ~(5) 所确定椒盐噪声的灰度估计值对当前像素点p进行判别。即, 如果当前像素点p的灰度值f (x, y)等于N_max或N_min, 则p即为椒盐噪声点, 从而在后续的噪声恢复阶段需要对其执行去噪处理;否则, 意味着当前像素点p未受到噪声的污染, 无需对其执行去噪处理。

1.2 基于核回归拟合的噪声恢复

核回归作为一种非线性方法, 已广泛应用于图像处理的多个领域^[13-18]。应用该方法计算椒盐噪声像素点的恢复值的过程具体如下:

假设由上述检测窗W_{x, y}内的非噪声像素点的灰度值拟合的函数r (x)属于再生核Hilbert空间。定义r_m为检测窗W_{x, y}内各个非噪声像素点的灰度值所构成的列向量, m为非噪声像素点的数量, 且第i个非噪声像素点的灰度值对应的核函数如下:

${\psi _i}(\mathit{\boldsymbol{x}}) = k(\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}_i}),{\rm{ }}i = 1,2, \cdots ,m$

(6)

其中：k (·, ·)为再生核Hilbert空间的核函数；x为非噪声像素点在检测窗中对应的位置；x_i表示第i个非噪声像素点在检测窗中对应的位置。则拟合函数r (x)由多个核函数线性叠加而成, 即

$r(\mathit{\boldsymbol{x}}) = \sum\limits_{i = 1}^m {{a_i}k(\mathit{\boldsymbol{x}},{\mathit{\boldsymbol{x}}_i})} $

(7)

其中：a_i为未知系数, 由其所构成的系数向量a_m可以使用最小二乘准则估计得到。a_m的计算公式如下:

${\mathit{\boldsymbol{a}}_m} = {\left[ {{a_1},{a_{\rm{2}}}, \cdots ,{a_m}} \right]^{\rm{T}}} = \mathit{\boldsymbol{K}}_m^ + {\mathit{\boldsymbol{r}}_m}$

(8)

其中：符号“+”表示矩阵的Moore-Penrose广义逆；K_m为m行m列的核矩阵。第s行t列的元素如式(9) 所示:

${\left( {{\mathit{\boldsymbol{K}}_m}} \right)_{st}} = k\left( {{\mathit{\boldsymbol{x}}_s},{\mathit{\boldsymbol{x}}_t}} \right)$

(9)

上述模型也称为基于最小二乘的核拟合(LS-based kernel fitting)模型^[16]。

在式(8) 中, 当核矩阵K_m的维数过大时, 其求逆的过程对计算机内存要求较高。为了避免这一现象, 采用递归模型的增量KNR (Kernel Nonlinear Representor)算法^[13]进行处理。首先, 定义参量τ_m+1、α_m+1, 如式(10) ~(11) 所示:

${\mathit{\boldsymbol{\tau }}_{m + 1}} = \mathit{\boldsymbol{K}}_m^ + {\mathit{\boldsymbol{\zeta }}_{m + 1}}$

(10)

${\alpha _{m + 1}} = {\left\| {{\psi _{m + 1}}} \right\|^2} \ge 0$

(11)

其中：ζ_m+1表示核矩阵K_m中第m+1列对角线以上的m个元素所构成的列向量。因此, 可采用如下递归式(12) 定义K_m⁺:

$\begin{array}{l} \mathit{\boldsymbol{K}}_{m + 1}^ + = \left\{ {\begin{array}{*{20}{c}} {\left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{K}}_m^ + + \frac{{{\mathit{\boldsymbol{\tau }}_{m + 1}}\mathit{\boldsymbol{\tau }}_{_{m + 1}}^ * }}{{{\alpha _{m + 1}}}}}&{ - \frac{{{\mathit{\boldsymbol{\tau }}_{m + 1}}}}{{{\alpha _{m + 1}}}}}\\ { - \frac{{\mathit{\boldsymbol{\tau }}_{_{m + 1}}^ * }}{{{\alpha _{m + 1}}}}}&{\frac{1}{{{\alpha _{m + 1}}}}} \end{array}} \right],}&{{\rm{ }}{\alpha _{m + 1}} > 0}\\ {\left[ {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{K}}_m^ + }&{{{\bf{0}}_{m \times 1}}}\\ {{{\bf{0}}_{1 \times m}}}&0 \end{array}} \right],}&{{\rm{ }}{\alpha _{m + 1}} = 0} \end{array}} \right.\\ \end{array}$

(12)

其次, 由式(8) 可得到包含m+1项系数的向量a_m+1的表达式如下:

${\mathit{\boldsymbol{a}}_{m + 1}} = \mathit{\boldsymbol{K}}_{m + 1}^ + {\mathit{\boldsymbol{r}}_{m + 1}}$

(13)

同时, 定义如式(14) ~(15) 所示的参量β_m+1、γ_m+1:

${\beta _{m + {\rm{1}}}} = {\rm{1}} + < {\mathit{\boldsymbol{\tau }}_{m + {\rm{1}}}},{\mathit{\boldsymbol{\tau }}_{m + {\rm{1}}}} > $

(14)

${{\gamma }_{m+1}}={\left( {{\mathit{\boldsymbol{r}}}_{m+1}}-<{{\mathit{\boldsymbol{r}}}_{m+1}},{{\mathit{\boldsymbol{ }}\!\!\tau\!\!\rm{ }}_{m+1}}> \right)}/{{{\beta }_{m+1}}}\;$

(15)

其中〈·, ·〉表示向量的内积。则由式(12) 和(13) 可得到求解系数向量a_m+1的递归公式如下:

${\mathit{\boldsymbol{a}}_{m + 1}} = \left\{ {\begin{array}{*{20}{c}} {\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{a}}_m} - \frac{{{\gamma _{m + 1}}{\beta _{m + 1}}}}{{{\alpha _{m + 1}}}}{\tau _{m + 1}}}\\ {\frac{{{\gamma _{m + 1}}{\beta _{m + 1}}}}{{{\alpha _{m + 1}}}}} \end{array}} \right],}&{{\rm{ }}{\alpha _{m + 1}} > 0}\\ {\left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{a}}_m}}\\ 0 \end{array}} \right],}&{{\rm{ }}{\alpha _{m + 1}} = 0} \end{array}} \right.$

(16)

由上述模型可见, 通过递归式(16) 估计出系数向量a_m+1之后, 按照式(7) 可生成检测窗内非噪声像素点的灰度值拟合曲面r (x), 且该拟合曲面符合检测窗所对应的图像局部区域的结构特征;最后, 以噪声像素点的空间坐标对该拟合函数进行重采样, 即可获得噪声像素点恢复后的灰度值。

1.3 算法步骤

本文算法的实现步骤归纳如下:

1) 应用1.1节所述的脉冲检测器对图像中的椒盐噪声像素点进行精确的检测。

2) 将步骤1) 所检测到的噪声像素点p当作缺失数据, 并采用1.2节所述的增量KNR模型对以噪声像素点p为中心的检测窗W_{x, y}内的非噪声像素点进行拟合, 得到其核回归拟合函数r (x)。

3) 以噪声像素点p的空间坐标对步骤2) 所得到的核回归拟合曲面r (x)进行重采样, 获得噪声像素点p恢复后的灰度值。

2 实验与分析

本章对如图 1所示的大小为256×256的6幅测试图像分别叠加密度为10%、50%、90%的椒盐噪声, 并采用本文算法以及近年具有的代表性的算法SMF、AMF、MDWMF、FSMMF、II进行去噪处理，以验证本文算法在低密度、中等密度、高密度椒盐噪声污染场景下的去噪性能。

不同密度噪声污染的Lena测试图像如图 2所示, 其相应的去噪结果图像及局部细节图像分别如图 3~6所示。实验结果的客观评价指标包括峰值信噪比(Peak Signal-to-Noise Ratio, PSNR)、平均绝对误差(Mean Absolute Error, MAE), 其计算公式分别如式(17) ~(18) 所示:

$PSNR = 10 \ \lg \frac{{{{255}^2}}}{{\frac{1}{{MN}}\sum\limits_{i = {\rm{1}}}^M {\sum\limits_{j = {\rm{1}}}^M {{{\left[ {f\left( {i,j} \right) - g\left( {i,j} \right)} \right]}^2}} } }}$

(17)

$MAE = \frac{1}{{MN}}\sum\limits_{i = 1}^M {\sum\limits_{i = 1}^N {\left| {f\left( {i,j} \right) - g\left( {i,j} \right)} \right|} } $

(18)

其中：f表示未受噪声污染、大小为M×N的原始图像；g则表示相应噪声污染图像的去噪结果图像; f (i, j)、g (i, j)分别表示原始图像、去噪结果图像中空间坐标为(i, j)的像素点的灰度值。各种算法的去噪结果的客观评价指标对比如表 1所示。

对比图 3、5、6所示的去噪结果图像可知, 上述各种算法具备不同程度的去噪性能, 且随着噪声密度的增加, 其去噪结果的主观视觉质量均有所降低。

首先, 在低密度噪声污染场景下, 由图 3所示的去噪结果图像可知, 本文算法及所比较算法均能够有效去除椒盐噪声, 但算法的细节保护性能各有差异。如图 4所示的去噪结果图像的局部细节对比, SMF算法由于对图像中所有的像素点均采用相同的中值滤波窗进行修正, 导致去噪结果图像整体产生明显的模糊效应;而AMF、MDWMF、FSMMF、II算法的去噪结果图像的主观视觉质量相当, 且优于SMF算法, 但图像细节区域仍然存在较为明显的模糊效应和锯齿现象;相对于上述算法, 本文算法的去噪结果图像的灰度分布更为平滑, 且边缘和纹理等细节更为锐化和完整, 整体视觉质量最接近原始图像, 其原因在于本文算法应用核回归方法对以噪声点为中心的邻域内的非噪声点进行拟合, 通过对拟合曲面进行重采样所获得的噪声像素点的恢复值更加符合图像的局部结构特征。

其次, 在中等密度噪声污染的场景下, 由图 5所示的去噪结果图像可知, 除算法SMF之外, 所比较算法均能实现较好的去噪效果。其中:AMF、MDWMF、FSMMF、II等算法的去噪结果图像的视觉质量相当, 均能够较好去除噪声, 但模糊效应和细节区域的锯齿现象仍然较明显;而本文算法由于噪声检测阶段的漏检和误检因素, 导致去噪结果图像在不同区域均存在少量的细小斑点, 但其主观视觉整体显得更为平滑, 且细节区域的完整性优于上述算法。

再者, 在高密度噪声污染的场景下, 如图 6的去噪结果所示, 除本文算法之外, 所比较算法的去噪结果图像的视觉质量均急剧下降;AMF算法和MDWMF算法的去噪结果图像模糊效应明显, 不能完整地体现原始图像的内容;而算法FSMMF、II的去噪结果图像的视觉质量尽管较接近于本文算法, 但FSMMF算法的去噪结果图像在细节区域产生严重的毛刺现象, II算法的去噪结果图像的对比度则明显下降;相反, 本文算法的去噪结果不仅完整地体现了原始图像的内容, 且有效地保持了边缘和纹理等细节的完整性。

最后, 通过对比分析表 1去噪结果的客观评价指标数据可知, 随着噪声密度的增加, 上述各种算法的去噪结果的PSNR均有所降低, 但在各种密度噪声污染的场景下, 本文算法对6幅测试图像的去噪结果的PSNR均高于上述所比较算法, 客观反映了本文算法对噪声像素点的恢复能力优于所比较的算法;其中, 在低密度和中等密度噪声污染的场景下, 本文算法的PSNR分别平均提高了6.02 dB、6.33 dB, 而在高密度噪声污染的场景下则平均提高了5.58 dB。同时, 由表 1所示的各项MAE可知, 本文算法在各种密度噪声污染的场景下对6幅测试图像的去噪结果的MAE均为最低, 说明本文算法的去噪结果图像整体上最接近相应的未受噪声污染的原始图像;其中, 在低密度和中等密度噪声污染的场景下, 本文算法的MAE分别平均降低了0.90、5.84, 而在高密度噪声污染的场景下则平均降低了25.29。

3 结语

针对任意密度椒盐噪声的去噪和图像细节保护问题, 本文提出了一种基于核回归拟合的开关去噪算法。算法在引入高效脉冲检测器对椒盐噪声像素点进行精确检测的基础上, 将所检测的噪声像素点当作缺失数据, 应用核回归方法构造以噪声像素点为中心的邻域内的非噪声像素点拟合曲面, 并通过重采样获得符合图像局部结构特征的噪声像素点的恢复值。实验结果表明, 本文算法不仅能够有效去除各种密度的椒盐噪声, 去噪性能明显优于新近的代表性算法, 且具备良好的图像细节保护能力。另一方面, 本文算法在中等密度和高密度噪声污染的场景下, 仍然存在较为明显的噪声漏检和误检问题, 容易导致去噪结果图像产生斑点, 影响了算法的去噪性能;因此, 进一步提高算法检测噪声的准确性是后续研究工作的重点。