图像在采集或传输过程中不可避免地受到噪声污染，造成图像质量的退化。这样的图像不仅影响视觉效果，而且还阻碍了图像视频编码、识别、场景理解等计算机视觉应用。因此，图像去噪一直是图像处理中一个重要的研究课题。高效合理地去除噪声，同时保留图像原本的特性是去噪算法的关键。

针对图像去噪，研究人员提出大量的算法。这些算法主要在空间域或频率域上利用滤波技术重建图像。以Gaussian滤波去噪为主的算法平滑效果好，但是边缘轮廓以及图像结构遭到破坏；另一方面，以各向异性滤波去噪为主的算法具有较好的保边效果，但易引入阶梯人造特征。文献[1]提出了一种异于传统去噪算法的新方法——非局部均值(non-local means，NLM)去噪算法，利用局部图像块(patch)冗余结构信息，将图像中结构相似的像素采用加权平均来获得去噪估计值。由于图像块比单个像素更好地表达图像结构信息，在不损失图像细节的情况下表现出优良的去噪性能，故其性能优于其他的经典去噪算法。NLM虽然性能优越，但效率与参数需要进一步优化。近年，很多学者对NLM提出了改进方法，这些方法主要集中在以下方面：1)降低计算复杂度(K-SVD^[2-3]、DCT^[4]、Walsh-Hadamard投影^[5]、主成分分析^[6-7])。2)优化参数^[8-12](平滑核参数、相似邻域的尺寸和搜索区域的大小)，提高性能。3)以NLM思想为基础，融合其他算法以提升去噪性能^[13-14]。

本文针对NLM算法存在的问题，提出以下改进：1)针对算法计算复杂度高的问题，在搜索窗口内可能有很多不相似的块，利用$ \ell_2$范数逐次消元法，在图像积分图上通过剔除预测结构相似性低的图像块，避免大量的欧氏距离计算，有效地降低算法运行时间。2)针对固定搜索窗口去噪性能的问题，利用图像空间相关性强的特点以及剔除相似性低的图像块比率，提出基于Patch测地线距离寻找同质信息的动态调整搜索范围的方法。从PNSR和SSIM数值和视觉上可以看出，边缘去噪效果与经典的NLM算法相比有较好的提高，同时也节约时间开销。

1 NLM去噪算法描述

一般加性噪声图像的去噪模型可表示为：

$ v(n,m)=u(n,m)+\gamma (n,m) $

式中，n, m代表二维图像空间像素的坐标值；v(n, m)代表坐标为(n, m)的受噪声污染的图像像素值；u(n, m)代表无噪声污染的真实图像值；γ(n, m)代表高斯加性白噪声。对于一张大小N×M的二维图像v(n, m)，在NLM去噪算法中，对u(n, m)的估计值定义为：

$ u[f][n,m]=\sum\limits_{k=0}^{N-1}{\sum\limits_{l=0}^{M-1}{w}}[n,m;k,l]f[k,l] $

(1)

其中

$ \begin{array}{l} w[n,m;k,l] = \frac{1}{{Z[n,m]}}{{\rm{e}}^{ - \frac{{\left\| {N[{{\mathit P}_{n,m}}] - N[{{\mathit P}_{k,l}}]} \right\|_{2,a}^2}}{h}}}\\ \;\;\;\;\;\;Z[n,m] = \sum\limits_{k = 0}^{N - 1} {\sum\limits_{l = 0}^{M - 1} {{{\rm e}^{ - \frac{{\left\| {N[{{\mathit P}_{n,m}}] - N[{{\mathit P}_{k,l}}]} \right\|_{2,a}^2}}{h}}}} } \end{array} $

式中，$w[n,m;k,l]$为$u(k,l)$与$u(n,m)$相似性权重系数；$Z[n,m]$为归一化因子，满足$w[n,m;k,l] \in [0,1]$，且$\sum\limits_{n,m} {w[n,m;k,l]} = 1$；$ P_{n,m} $代表以(n, m)为中心大小为$(2b + 1) \times (2b + 1)$的图像块。h控制指数函数的衰减程度，一般设置成噪声标准方差σ的倍数。各图像块灰度值向量之间的相似性通过高斯加权的欧氏距离衡量。

对图像每个像素点进行$u[f][n,m]$估计，最终获得去噪结果。从NLM算法公式中可以看出，图像块(n, m)与图像块(k, l)越相似，用于去噪估计$u[f][n,m]$的权重越高。计算图像任意像素点的$u(n,m)$值需要图像中所有点构成的图像块计算权重。因此，图像块大小为$B \times B$，总时间复杂度为$O({N^2}{M^2}{B^2})$，由于图像块的结构相似度具有空间集中性，为了减少计算复杂度，将搜索范围限定在以(n, m)为中心的S区域内，时间复杂度降为$O(NM{S^2}{B^2})$。在去噪过程中，每个像素所在的图像块需要和搜索窗口的其他所有图像块，通过高斯加权的欧氏距离计算相似性获取权重，导致算法的计算复杂度高，难以应用到实际工程中。

2 选择性计算及自适应搜索区域的NLM算法

由于经典的NLM算法的去噪效果好，大量文献都希望通过降低计算复杂度来改进算法，以此达到推广应用的目的，其中降低度量像素间的相似性计算量是主要的改进方向。文献[2]提出利用K-SVD(singular value decomposition)分解技术降低维数从而减少计算量。文献[4]在离散余弦变换的低频系数子空间内度量像素间的相似性。文献[5]利用Walsh-Hadamard变换的快速运算和能量集中特性，通过投影到低维子空间进行度量像素点之间的相似性，减少高维计算。文献[6-7]通过主成分分析方法将图像块投影到低维子空间，并在该低维子空间中度量像素点之间的相似性。文献[15]利用图像块的均值和方差聚类图像块避免权值低的图像块的欧氏距离计算。这些方法主要是通过变换在子空间里度量相似性，减少计算时间，但同时会影响去噪效果。

本文的思路与文献[15]具有一定的相似，利用选择性计算欧氏距离降低运行时间。通过逐次消元法判定剔除相似性低的图像块。文献[15]需要计算图像块的均值和方差并聚类，本文的方法仅需简单的加法运算，从而大幅度地降低了计算量。

2.1 基于逐次消元法的加速策略

逐次消元法^[16]本质上是一种利用范数不等式淘汰冗余点的快速全局搜索算法，减少平均绝对差(mean absolute difference，MAD)和计算量，有效地降低算法复杂性，同时可以提供与全搜索算法相同的结果。逐次消元法在视频目标运动估计等应用中取得良好的效果。

X表示以(n，m)为中心的图像块像素灰度值的向量，Y(k, l)表示搜素窗口内的以(k, l)为中心的图像块像素灰度值的向量。首先考虑L1范数相似度度量方式应用不等式$\left| {{{\left\| \mathit{\boldsymbol{A}} \right\|}_1} - {{\left\| \mathit{\boldsymbol{B}} \right\|}_1}} \right| \leqslant {\left\| {\mathit{\boldsymbol{A}} - \mathit{\boldsymbol{B}}} \right\|_1}$，其中A=X，B=Y(k, l)，$\left| {{{\left\| \mathit{\boldsymbol{X}} \right\|}_1} - {{\left\| {\mathit{\boldsymbol{Y}}(k,l)} \right\|}_1}} \right| \leqslant {\left\| {\mathit{\boldsymbol{X}} - \mathit{\boldsymbol{Y}}(k,l)} \right\|_1}$，${\left\| \mathit{\boldsymbol{X}} \right\|_1} = $$\sum\limits_{n,m \in P} {{x_{n,m}}} $。

根据L1范数与L2范数的边界定义可知，${\left\| {\mathit{\boldsymbol{X}} - \mathit{\boldsymbol{Y}}(k,l)} \right\|_1} \leqslant \sqrt {\left| \mathit{\boldsymbol{X}} \right|} \times {\left\| {\mathit{\boldsymbol{X}} - \mathit{\boldsymbol{Y}}(k,l)} \right\|_2}$，其中$\left| \mathit{\boldsymbol{X}} \right|$为X向量的维数，则$\frac{{\left| {{{\left\| \mathit{\boldsymbol{X}} \right\|}_1} - {{\left\| {\mathit{\boldsymbol{Y}}(k,l)} \right\|}_1}} \right|}}{{2b + 1}} \leqslant {\left\| {\mathit{\boldsymbol{X}} - \mathit{\boldsymbol{Y}}(k,l)} \right\|_2}$，在不等式两边同时乘上高斯加权矩阵，则有：

$ {G_{\rm{ \mathsf{ σ} }} } \times \frac{{\left| {\sum\limits_{i = - b}^{i < b} {\sum\limits_{j = - b}^{j < b} {I(n + i,m + j)} } - \sum\limits_{i = - b}^{i < b} {\sum\limits_{j = - b}^{j < b} {I(k + i,l + j)} } } \right|}}{{2b + 1}} \leqslant\\ \ \ \ {G_{\rm{ \mathsf{ σ} }} } \times \sqrt {\sum\limits_{i = - b}^{i < b} {\sum\limits_{j = - b}^{j < b} {{{(I(n + i,m + j) - I(k + i,l + j))}^2}} } } $

(2)

式(2)中的高斯加权计算可以忽略。在搜索窗口S中可能有很多不相似的图像块，在去噪估计时，这些图像块对应的像素分配的权重越小越好。由于NLM算法使用指数形式的权重函数，当不相似的图像块较多时，这些小的权重值总和将较大，这对去噪过程将产生不利影响。因此，本文的NLM算法不再用搜索窗口S中所有的像素点，而是在搜索窗口S中选取q个权值$w[n,m;k,l]$最高的图像块对应像素点进行加权去噪，一般设置为$q = \alpha {(2s + 1)^2}$，$\alpha $为阈值。假设坐标为(k, l)的像素点为当前第q个候选参与坐标为(n, m)的像素点去噪估计点，有：

$ \frac{{\left| {{{\left\| \mathit{\boldsymbol{X}} \right\|}_1} - {{\left\| {\mathit{\boldsymbol{Y}}(i,j)} \right\|}_1}} \right|}}{{2b + 1}} \leqslant {\left\| {\mathit{\boldsymbol{X}} - \mathit{\boldsymbol{Y}}(k,l)} \right\|_2} $

(3)

如果式(3)满足，需要计算${\left\| {\mathit{\boldsymbol{X}} - \mathit{\boldsymbol{Y}}(i,j)} \right\|_2}$；如果不满足，则有${\left\| {\mathit{\boldsymbol{X}} - \mathit{\boldsymbol{Y}}(i,j)} \right\|_2} > {\left\| {\mathit{\boldsymbol{X}} - \mathit{\boldsymbol{Y}}(k,l)} \right\|_2}$，说明像素点(i, j)权值较低，无需计算欧氏距离。

为了比较当前像素点是否满足式(3)，需要计算$\left| {{{\left\| \mathit{\boldsymbol{X}} \right\|}_1} - {{\left\| {\mathit{\boldsymbol{Y}}(i,j)} \right\|}_1}} \right|$，即计算图像块P_{n, m}与P_{i, j}之差的绝对值。为了得到${\left\| {{\mathit{\boldsymbol{P}}_{n,m}}} \right\|_1}$和${\left\| {{\mathit{\boldsymbol{P}}_{i,j}}} \right\|_1}$(图像块P_{n, m}与P_{i, j}的像素灰度值积分)，本文利用图像积分图来加速计算。积分图^[17]概念用于计算窗口图像的Harr特征。在计算机视觉领域里，图像积分图主要用于规则区域特征值的计算，有着广泛的应用。图像积分图算法的基本原理为：

$ II(x,y) = \sum\limits_{x' < x,y' < y} {I(x',y')} $

(4)

式中，$II(x,y)$表示矩阵元素的积分图；$I(x',y')$表示图像矩形区域元素的像素值。

临时变量$X(x,y)$表示点$(x,y)$所在列纵坐标不大于y的所有像素灰度值的积分，其定义为$X(x,y) = \sum\limits_{l = 1}^y {I(x,l)} $，根据定义$X(x,0) = 0$，$X(0,y) = $ $0$，有：

$ X(x,y) = X(x,y - 1) + I(x,y) $

(5)

$ II(x,y) = II(x - 1,y) + X(x,y) $

(6)

$ SII(x,y) = II(x + i,y + j) - II(x + i,y) - II(x,y + j) + II(x,y) $

(7)

通过式(5)、式(6)对原图像每个像素一次访问便可计算出图像积分图，用式(7)可以快速计算出图 1中窗口$S(x,y)$的像素值和。

积分图创建以后，只需要进行7次加法和1次除法运算，可以得出逐次消元法淘汰冗余点的条件。由于需要寻找前q个权值高的像素点，在极端情况下，如果每次检查点都满足式(3)，那么采用逐次消元法并不能降低其计算量。为了减少该情况发生，在计算图像块间距离前，先对搜索窗口内所有图像块与以(n, m)点构成的图像块以均值为度量进行类聚，选出q个与(n, m)点均值最近的像素点，作为起始选择点。图像块的均值利用积分图可以在初始时一次性计算得到，因此，时间复杂度并不会增加太多。

2.2 动态搜索区域

传统的NLM算法为了减少计算量，搜索范围由整个图像区域压缩到中心像素周围一定大小的区域$(2s + 1) \times (2s + 1)$，称为搜索窗口。固定搜索窗口对像素之间的结构信息有较好的保留，然而没有考虑到图像的同质信息，不能充分利用相似性高的像素点。由于同质信息能较好地反映图像的相似信息，对去噪更有利。根据图像的同质信息，动态调整搜索窗口，获得每个像素点的自适应搜索区域，去噪效果将更加理想。针对该问题，本文提出一种基于Patch测地线(patch geodesic path，PatchGP)距离寻找同质信息动态调整搜索窗口的方法。

文献[18]提出以近似估计PatchGP为图像块之间相似性的NLM去噪，对边缘信息去噪效果好以及图像细节保留较完整，但对低频噪声去除不够理想，平滑区域易出现块状斑迹。由于PatchGP既可以反映图像相似结构，又能利用测地线距离找到同质信息。因此，本文根据PatchGP为线索搜寻同质信息来动态调整搜索窗口的方向与大小。

在图像I中的图像块c和d的PatchGP定义为：

$ {D^{{\rm PatchGP}}}(\mathit{\boldsymbol{c}},\mathit{\boldsymbol{d}}) = \mathop {\min }\limits_\mathit{\Gamma } \sum\limits_{i = 1}^{{N_r} - 1} {\left\| {{N_I}(\mathit{\Gamma } ({p_{i + 1}})) - {N_I}(\mathit{\Gamma } ({p_i}))} \right\|} $

式中，$\mathit{\Gamma } $表示c和d构成的路径集合；${N_r}$表示路径长度；$\mathit{\Gamma } ({p_i})$表示路径上的像素点位置；${D^{{\rm PatchGP}}}(\mathit{\boldsymbol{c}},\mathit{\boldsymbol{d}})$表示路径上所有相邻像素点构成的图像块差异累积和，其实质是以图像块的特征点为度量方式来寻找c和d间的最短路径。基于像素的测地线距离是累积梯度信息，易受噪声的干扰；而PatchGP利用图像块较好地保留图像结构冗余信息，且不易受噪声影响，更能反映图像的本质，是本文选用PatchGP的理由。文献[18]中的PatchGP是以路径上相邻图像块的相似度为度量方式，得到是路径梯度的最小值。然而，本文利用PatchGP是为了找到与中心像素更接近的同质区域，应该寻找与中心像素相似度更高的路径，因此，采用式(8)修改PatchGP度量方式，c为中心像素的图像块，有：

$ {D^{{\rm PGP'}}}(\mathit{\boldsymbol{c}},\mathit{\boldsymbol{d}}) = \mathop {\min }\limits_\mathit{\Gamma } \sum\limits_{i = 1}^{{N_r} - 1} {\left\| {{N_I}(\mathit{\Gamma } ({p_i})) - {N_I}(\mathit{\boldsymbol{c}})} \right\|} $

(8)

$ {C_{{\rm ratio}}} = \frac{{\sum\limits_{i = - r}^r {\sum\limits_{j = - r}^r {\mathit{\boldsymbol{I}}(i,j) \in {\rm NN}(q)} } }}{{(2r + 1) \times (2r + 1)}} $

(9)

式中，${\rm NN}(q)$为当前最高权重的$q$个像素点；${C_{{\rm ratio}}}$为在大小为$(2r + 1) \times (2r + 1)$的区域R内像素点在${\rm NN}(q)$中的比值。

图 2b所示为图 2a标准图像Barbara的手臂部分，图像块P与P₁、P₂、P₃相比，P₁、P₃与P结构更相似。在搜索窗口S里，P₃被漏掉，没有充分利用图像结构冗余的特点寻找到同质区域。采用以中心像素的S区域计算相似度，然后在S区域边缘像素和中心像素使用式(8)的${D^{{\rm PGP'}}}$(c, d)计算相似度，采用式(10)得到最小PatchGP的边缘点，以该边缘点为边缘中心向外扩展区域R，再次计算R区域与中心像素的图像块的相似度，并将R区域合并到S区域中，用式(9)的权值替换比率控制迭代执行次数。为了控制图像块引起的噪声，当替换比率大于$\delta $，才能再次增加搜索范围，有：

$ \partial s_i^{\min } = \arg \min \left\{ {{D^{{\rm PGP'}}}(\mathit{\boldsymbol{c}},\partial {s_i}),{\rm{ }}{D^{{\rm PGP'}}}(\mathit{\boldsymbol{c}},\partial {s_i}) \ne \infty } \right\} $

(10)

式中，$\partial {s_i}$表示以区域边缘点为中心的图像块。直接计算PatchGP的计算量是巨大的，文中利用中心像素与搜索区域S内其他像素图像块的欧氏距离为度量，被逐次消元算法剔除的像素点设置为$\infty $。由于计算邻域权值时，其相似度距离已经获得，${D^{{\rm PGP'}}}(\mathit{\boldsymbol{c}},\partial {s_i})$只是由中心像素出发采用八邻域矩阵寻找一条最短路径，如果邻域被剔除无法到达某边缘点，设置该边缘点的${D^{{\rm PGP'}}}(\mathit{\boldsymbol{c}},\partial {s_i}) = \infty $。由于设定S区域初始值较小，其计算量可忽略不计。

像素点I(n, m)的邻域搜索范围确定流程如下：

1) 初始化。以I(n, m)中心点的(2s+1)×(2s+1)搜索区域S；利用逐次消元法，计算区域S与中心像素点之间的相似度矩阵。

2) 用式(10)找到PatchGP最小的边缘点，以该边缘点为矩形边中心扩展大小为(2r+1)×(2r+1)的区域R，$\mathit{\boldsymbol{S}} = \mathit{\boldsymbol{S}} \cup \mathit{\boldsymbol{R}}$，利用逐次消元法计算区域R与中心像素点之间的相似度矩阵，并计算替换比率。

3) 如果${C_{{\rm ratio}}} \geqslant \delta $，跳转步骤2)。

4) 最终区域S为I(n, m)的邻域搜索范围，并得到与中心点的相似度。

3 实验结果与分析

本文采用4个标准的图像Lena、Pepper、Barbara和Baboon，并加入多个取值的标准差σ的零均值高斯白噪声。在实验中，本文的方法分别与K-SVD- NLM^[2]、DCT-NLM^[4]、FM-PatchGP^[18]、cNLM^[19]和BM-3D^[20]共5种基于块相似性的去噪算法从速度和性能两方面做定性比较与分析。其中，cNLM代表经典优化参数的NLM算法，DCT-NLM代表变换子空间去噪法、FM-PatchGP、K-SVD-NLM、BM-3D代表去噪性能最佳的方法。为了比较的公平，所有比较算法的搜索窗口和图像块大小参考文献[15]分别设置为21×21和7×7。本文的算法参数设置$\delta = 0.08$，$\alpha = 0.6$，初始搜索窗口S为11×11，扩展窗口R为9×9，图像块大小为7×7，h=0.7×σ。

为了验证本文的方法加速和自适应搜索范围的有效性，本文的算法和DCT-NLM均使用C++语言编程实现，cNLM和K-SVD-NLM源代码来自文献[19]，FM-PatchGP和BM-3D源代码来自作者网站，测试环境为visual C++ 2008、Dell Inspiron、Core i3、2 GB内存。采用PSNR(peak signal-to-noise ratio)和SSIM(structural similarity index)两个指标对去噪结果进行定量比较。SSIM是一种以度量两幅图像间的结构性失真来评估图像质量评价度量的新指标，取值范围[0~1]，其值越大越接近，当图像完全一致时，SSIM为1。

图 3和图 4所示分别为6种算法在$\sigma $={10, 15, 20, 25, 30，35, 40}取值下的PSNR、SSIM值比较曲线图。从图 3中可以看出，DCT-NLM整体指标较好，随着$\sigma $增大，本文的算法数值基本接近DCT-NLM。从图 4可以看出，K-SVD-NLM随着噪声水平的增加，SSIM值下降较快。

从PSNR、SSIM两个指标数值上可得，本文的算法都比cNLM算法有较大的提高。本文的算法虽然当$\sigma $较小时，PSNR值略低于DCT-NLM和FM-PatchGP，接近BM-3D，且高于K-SVD-NLM。但随着$\sigma $增大时，图像Lena、Pepper和Barbara的PSNR值和SSIM值均不同程度地超过了其他5种算法。同时随着噪声水平的增加，该算法性能下降较缓慢。该算法具有良好的性能指标值，调整搜索区域的改进NLM算法通过寻找同质信息，有效地阻止了噪声信息的干扰，增强了去噪能力。

由于图 3d的测试图像Baboon存在细腻且不规则的纹理，本文的算法降低了自适应搜索区域能力，退化为较小固定搜索窗口的NLM且噪声干扰无法辨别，导致PSNR数值下降较快，性能略低于cNLM，但图 4d说明图像整体结构保留较好，同样由于不规则细腻的纹理导致PatchGP不够准确，FM-PatchGP的PSNR和SSIM值都下降较快。

图 5所示为Barbara图像部分细节原图以及6种算法去噪后的对比图(σ=25)。图 5a包含了左上角的平滑区域和头巾部分的纹理边缘区域。对比可见，除了FM-PatchGP，其他的5种算法都可以较好地保留图像的结构与细节。cNLM、K-SVD-NLM和BM-3D算法去噪结果的平滑区域部分不同程度出现明显的伪影。本文的算法将结构相似性和同质信息相结合动态调整搜索窗口，有效地去除了噪声干扰，较好地利用了图像冗余相似性，其边缘区更整齐，图像结构信息保留较完整。由于选取相似性最接近的q个像素点，对于平滑区域，该算法的去噪能力也明显强于cNLM，避免了伪影现象。与DCT-NLM算法相比，本文的算法对图像的边界细节纹理信息保留得稍好，图像细节更清晰。

逐次消元法剔除相似性低的图像点和优化的搜索窗口都有效地降低了计算复杂度。表 1记录了本文的算法与K-SVD-NLM、DCT-NLM、FM-PatchGP、cNLM和BM-3D对本文测试的4幅图像的平均耗时比较。从表中可以看出，对于大小为128×128和256×256的图像，FM-PatchGP耗时最少，但其去噪性能较低以及图像细节丢失严重。对于大小为512×512和1 024×768的图像，由于本文的算法使用逐次消元法剔除像素点，图像越大，相似度低的像素点被丢弃得越多。该算法计算量增幅较小，与其他几种算法相比，具有良好的加速效果。

综上分析，从总体性能指标来看，本文的算法具有一定的优势。同时，该方法也存在如下两点不足：1)不太适应含有较多的细腻且不规则的纹理图像。2)从去噪指标和运行时间比较看，对于噪声σ＜20和低于256×256大小的图像，该算法的整体性能略低于DCT-NLM算法。

4 结论

NLM算法具有较好的去噪性能，但存在运行效率低的问题，本文利用L2范数的逐次消元法和积分图，有效地降低了像素间相似性计算强度，并根据图像空间相关性强的特点，提出一种基于Patch测地线距离寻找同质信息的自适应调整搜索窗口的方法。在实验中，本文的算法与当前性能最佳的5个基于图像块去噪算法相比，从PSNR、SSIM两个指标和运行时间比较上可知，该方法不但获得了较好的加速，而且去噪效果也有一定的提升。