Non Local Means 图像降噪算法

Non Local Means filter Author: Antoni Buades, Computer Vision and Pattern Recognition, 2, 60-65 (2005).

Antoni Buades 提出指标 method noise 对数字图像降噪方法的性能进行了评价和比较。他首先针对几个被广泛使用的降噪算法计算并分析了降噪性能。同时，基于图像中所有像素的非局部平均，提出了全新的数字图像降噪算法 Non Local means Algorithm，并通过实验比较了新算法与常用的平滑滤波方法的性能。

Introduction

数字图像噪声

数字图像的两个主要限制是模糊及噪声（blur & noise）。模糊是图像采集系统的固有性质，同时数字图像对连续信号离散采样的形式必须遵循 Shannon-Nyquist 采样定律。另一种主要的干扰形式是噪声。数字图像中每个像素点的值 $u (i)$ 都是对局部光强测量的结果，通常通过 CCD 及光学聚焦元件实现。CCD 中每个方形探测单元（captor）都将记录曝光时间内探测区域的入射光子数。在光源强度恒定的情况下，每个探测单元每个曝光周期内采集的光子数的概率分布将遵循中心极限定理，在光强均值周围震荡。另外，如果 CCD 没有经过充分冷却就会接收热光子（heat spurious photons），由此产生的扰动通常称为 obscurity noise。

数字图像降噪方法的目标是从观测到的噪声图像中还原出原始信号，

$v(i)=u(i)+n(i), \tag{1}\label{1}$

其中 $v (i)$ 为实测图像， $u (i)$ 为原始信号， $n (i)$ 则是噪声信号。评估图像中噪声水平通常采用信噪比（signal noise ratio (SNR)）：

$SNR = \frac{\sigma(u)}{\sigma(n)}, \tag{2}\label{2}$

其中， $\sigma (n)$ 为噪声信号标准差， $\sigma (u)$ 表示真实信号的经验标准差，

$\sigma(u) = \sqrt{\frac{1}{|I|}\sum_i(u(i)- \overline{u} )^2}, \tag{3}\label{3}$

$\overline {u}=\frac {1}{|I|}\sum_{i\in I} u (i)$ 为图像的平均灰度值， $|I|$ 指全图像素数。当噪声模型和参数已知时，噪声的标准差也可以用经验测量法或形式化方法得到。

迄今为止，图像处理领域已经提出了诸多抑制噪声、还原真实信号的算法。即便它们通常拥有截然不同的数学形式，但都拥有一个共性：平均。这种平均可以在局部进行：高斯滤波器 (Gabor 1994)，各向异性滤波 (Perona-Malik 1990, Alvarez et al. 1992)，邻域滤波器 (Yaroslavsky 1985, Smith et al. 1997)；也可以通过计算 variations 实现：TV 滤波 (Rudin-Osher-Fatemi 1992)；或是在频域进行：经验维纳滤波 (Yaroslavsky 1985)，小波阈值方法 (Coiffman-Donoho 1995)。

Method noise

不妨令 $u$ 表示实测图像， $D_hu$ 表示降噪方法的输出结果， $h$ 为滤波参数。Antoni Buades 定义 method noise 为降噪前后图像之差：

$u-D_hu. \tag{4}\label{4}$

完美的降噪算法在应用中不应该改变无噪声图像。因此，当图像具有某种规律性时，method noise 理应很小。对理想的降噪算法，Method noise 必须看起来与随机噪声无异，几乎不包含原始信号的结构。因为即便是质量非常高的实测图像，噪声也是不可避免的，计算 method noise 对评估任何降噪算法都是有意义的，而非传统的 “添加噪声，再去除噪声” 的把戏。

局部平均算法

高斯滤波 (Gaussian Filtering)

对数字图像进行各向同性过滤，本质上可以归结为图像与各向同性核的卷积。采用数值呈现出高斯分布的卷积核，既是高斯滤波，是图像处理中最常用的操作之一。通俗的讲，高斯滤波就是对整幅图像进行加权平均的过程，每一个像素点的值，都由其本身和邻域内的其他像素值经过加权平均后得到（所有的局部平滑滤波方法都是如此）。高斯卷积核：

$G_h(x)=\frac{1}{4\pi h^2}e^{-\frac{|x|^2}{4h^2}}. \tag{5}\label{5}$

Theorem 1 (Babor 1960): 当高斯卷积核的特征尺寸 $h$ 极小时，高斯滤波的 method noise 为：

$u-G_h* u=-h^2\Delta u+o(h^2). \tag{6}\label{6}$

高斯滤波的 method noise 在图像谐波部分几乎为零，而在边缘、纹理区域非常大。因此，高斯滤波在图像的平坦区域相对表现优秀，但在边缘、纹理区域较为模糊（blurred）。

各向异性滤波（Anisotropic Filtering, AF）

各向异性滤波提出之初旨在解决高斯滤波在边缘及纹理区域的模糊问题。该算法通过只在 $Du (\vec x)$ 正交方向上计算图像 $u$ 在 $\vec x$ 处的卷积。这种想法可以追溯到 Perona & Malik。各向异性滤波算法的定义为：

$AF_hu(\vec x)=\int G_h(t)u(\vec x+t\frac{Du(\vec x)^\perp}{|Du(\vec x)|})dt, \tag{7}\label{7}$

在 $\vec x$ 处，当 $Du (\vec x)\neq0$ 时成立。 $(x,y)^\perp=(-y,x)$ ，且 $G_h$ 代指方差 $h^2$ 的一维高斯函数。假设原始图像 $u$ 在 $\vec x$ 处二次连续可微（twice continuously differentiable ( $C^2$ )）, 将式 \eqref {5} 二次泰勒展开（second order Taylor expansion）可以推导出：

Theorem 2: 当 $Du (\vec x)\neq0$ 时，各向异性滤波 $AF_h$ 的 method noise 为：

$u(\vec x)-AF_hu(\vec x)=-\frac{1}{2}h^2|Du|curv(u)(\vec x)+o(h^2), \tag{8}\label{8}$

$curv (u)(\vec x)$ 指局部曲率（此处存疑），即经过 $\vec x$ 点的水平直线上曲率半径的逆（signed inverse）。在图像 $u$ 中局部几乎为一条直线的区域，method noise 几乎为零，而在弯曲的边缘或纹理区域较大。因而，各向异性滤波对直边区域得以保持，而平坦、纹理区域图像精度有所退化。

TV 滤波（Total Variation Minimization）

全变分图像去噪算法最早由 Rudin, Osher and Fatemi 提出。给定一幅实测图像 $v (\vec x)$ ，该算法将恢复原始信号 $u (\vec x)$ 的问题转化为式 \eqref {7} 的最小化问题：

$TVF_\lambda(v)={\rm arg}\,\underset{u}{\rm \mathop {min}}\,TV(u)+\lambda\int|v(\vec x)-u(\vec x)|^2d\vec x, \tag{9}\label{9}$

其中 $TV (u)$ 为图像 $u$ 的全变分， $\lambda$ 为给定的拉格朗日乘子（Lagrange multiplier）。上述最小化问题的最小值存在且唯一。参数 $\lambda$ 与噪声的统计信息相关，并控制了滤波程度。

Theorem 3: TV 滤波的 method noise 为：

$u(\vec x)-TVF_\lambda(u)(\vec x)=-\frac{1}{2\lambda}curv(TVF_\lambda(u))(\vec x). \tag{10}\label{10}$

在各向异性的情况下，直边由于曲率小而得以保持。但 $\lambda$ 过小时细节及纹理会被过度平滑。

邻域滤波（Neighborhood Filtering）

我们称邻域滤波为将临近区域内具有相似灰度值的像素取平均以期恢复原始信号的滤波器。Yaroslavsky(1985) 首次提出的方法通过计算空间邻域 $B_\rho (\vec x)$ 内具有相似灰度值像素的平均来恢复信号：

$YNF_{h,\rho}u(\vec x)=\frac{1}{C(\vec x)}\int_{B_\rho (\vec x)} u(\vec y)e^{-\frac{|u(\vec y)-u(\vec x)|^2}{h^2}}d\vec y, \tag{11}\label{11}$

其中， $\vec x\in\Omega$ ， $C (\vec x)=\int_{B_\rho (\vec x)} e^{-\frac {|u (\vec y)-u (\vec x)|^2}{h^2}} d\vec y$ 为归一化常量， $h$ 为滤波参数。

较晚提出的邻域滤波算法相较 Yaroslavsky 滤波要更广为人知一些，既 SUSAN 滤波（1995）及双边滤波（1998）。这些滤波算法都引入到参考像素 $\vec x$ 的距离作为权重因子，而非单纯的考虑一个固定范围的邻域，

$SNF_{h,\rho}u(\vec x)=\frac{1}{C(\vec x)}\int_{\Omega} u(\vec y)e^{-\frac{|y-x|^2}{\rho^2}}e^{-\frac{|u(\vec y)-u(\vec x)|^2}{h^2}}d\vec y, \tag{12}\label{12}$

这里 $C (\vec x)=\int_{\Omega} e^{-\frac {|y-x|^2}{\rho^2}} e^{-\frac {|u (\vec y)-u (\vec x)|^2}{h^2}} d\vec y$ 为归一化常量， $\rho$ 为空间滤波参数（spatial filtering parameter）。事实上， $YNF_{h,\rho}$ 与 $SNF_{h,\rho}$ 之间并没有本质区别（此处存疑）。如果两个区域的灰度值差异大于 $h$ ，这些算法都将计算来自于同一区域的像素灰度平均值来恢复参考点的原始信号。因而该算法不会模糊边界区域，这是正是该类算法最核心的用途。

然而这类算法的问题是，只将单个像素作为参考点，而如若该参考像素恰好被噪声干扰严重，滤波效果将不够鲁棒（robust）。同时，邻域滤波器也会制造人为干扰（artificial shocks），这将会在它的 method noise 中展示出来。

非局部平均算法（Non Local Means Algorithm）

Antoni Buades 于 2005 年提出非局部平均数字图像降噪算法（Non Local Algorithm）。给定一幅实测图像 $v=\lbrace v (i)|i\in I\rbrace$ ，像素 $i$ 处的估计值 $NL [v](i)$ 是该图像上所有像素点的加权平均值，

$NL[v](i)=\sum_{j\in I}\omega(i,j)v(j), \tag{13}\label{13}$

这里的权重系数 $\lbrace\omega (i,j)\rbrace_j$ 取决于像素 $i$ 与 $j$ 之间的相似程度，且始终满足如下标准： $0\leq\omega (i,j)\leq 1$ 且 $\sum_{j} \omega (i,j)=1$ （等价于上文介绍的滤波算法中归一化常量）。

俩个像素 $(i,j)$ 之间的相似程度取决于邻域灰度矩阵 $v (N_i)$ 及 $v (N_j)$ （intensity gray level vectors），这里 $N_k$ 指以像素 $k$ 为中心，给定尺寸的方形邻域。这种相似性被定义为加权欧式距离的递减函数， $\parallel v (N_i)-v (N_j)\parallel_{2,a}^2$ ，其中 $a>0$ 是高斯卷积核的标准差。欧式距离对噪声邻域的应用引入了如下等式：

$E\parallel v(N_i)-v(N_j)\parallel_{2,a}^2 = \parallel u(N_i)-u(N_j) \parallel_{2,a}^2+2\sigma^2. \tag{14}\label{14}$

这个等式证明了该算法的鲁棒性（robustness）。因为含噪声的实测图像 $v$ 的欧式距离期望恰恰遵循真正的原始信号之间的相似性。

与 $v (N_i)$ 具有相似灰度邻域的像素在计算平均时的权重因子更大，

$\omega(i,j)=\frac{1}{Z(i)}e^{-\cfrac{\parallel v(N_i)-v(N_j) \parallel_{2,a}^2}{h^2}}, \tag{15}\label{15}$

这里， $Z (i)$ 为归一化常量，

$Z(i) = \sum_j e^{-\cfrac{\parallel v(N_i)-v(N_j) \parallel_{2,a}^2}{h^2}}, \tag{16}\label{16}$

参数 $h$ 控制滤波程度，它直接影响了指数函数的衰减趋势，进而控制欧式距离对权重因子衰减速度的影响。

NL-means 算法不仅仅考虑单个像素的灰度值，而是考虑该像素整个邻域的几何构型，这正是 NL-means 算法比邻域滤波更鲁棒的原因。图 $(1)$ 也说明了这个问题，像素 $q3$ 与 $p$ 具有完全一致的灰度值，而邻域的几何构型完全不同，导致 NL-means 中的权重因子 $\omega (p,q3)$ 几乎为零。

$图1. NL-means 算法方案。相似的像素邻域将提供更大的权重，如 $\omega(p,q1), \omega(p,q2)$，而不相似的邻域提供的权重几乎为零，如 $\omega(p,q3)$。$

图 1. NL-means 算法方案。相似的像素邻域将提供更大的权重，如 $\omega (p,q1), \omega (p,q2)$ ，而不相似的邻域提供的权重几乎为零，如 $\omega (p,q3)$ 。

NL-means 算法最终的数学形式：

$NL[v](x)=\frac{1}{C(x)} \int_{\Omega} e^{-\cfrac{(G_a * |v(x+.)-v(y+.)|^2)(0)}{h^2}}v(y)dy, \tag{17}\label{17}$

$x\in\Omega$ ， $C (x)=\int_{\Omega} {\rm exp}\lbrack -\frac {(G_a* |u (x+.)-u (z+.)|^2)(0)}{h^2}\rbrack dz$ 为归一化常量， $G_a$ 为高斯核， $h$ 控制过滤程度。

NL-means 算法的中心思想是：像素 $x$ 处的信息恢复，是由整幅图像内所有邻域与像素 $x$ 邻域相似的点取平均得到的。与局部滤波算法或频域滤波算法相比，NL-means 算法的主要区别在于可以系统地运用整幅图像中所有可能自预测局部结构的信息。

实验及讨论

本节将针对以下三点比较 Non Local Means 算法与局部平滑滤波的性能： method noise，视觉效果，均方差（mean square error）。这里均方差指修复图像与真实图像的欧几里得差异（Euclidean difference）。

在实现 NL-means 算法时，我们将相似邻域的搜索范围约束在一个较大的 $S\times S$ pixels 区域内。在所有的实验中均固定该搜索范围为 $21\times 21$ pixels，并指定邻域 $N_i$ 范围为 $7\times 7$ pixels。不妨设输入图像的像素数为 $N^2$ ，那么该算法的时间复杂度约为 $21^2\times 7^2 \times N^2$ 。

图 2. NL-means 权重分布。每组中左图中心像素为参考像素点 $x$ ，每组中右图为权重因子分布，由黑至白对应从 0 到 1。

$7\times 7$ pixels 的邻域范围已经足够大，能够确保算法对噪声是鲁棒的；也足够小，让细节及纹理得以保持。滤波参数 $h$ 被设置为 $10\sigma$ ，这里 $\sigma$ 为人工添加高斯白噪声的标准差。由于指数权重的快速衰减，距离中心较远像素的权重几乎为零，这发挥了自动剔除远处像素的作用。Non Local means 算法的权重分布见图 $(2)$ 。

图 3. 对自然图像的降噪实验。从左至右：含噪声的实测图像（噪声标准差 35），高斯滤波，各向异性滤波，TV 滤波，邻域滤波及 NL-means 算法。

在前面的段落中，我们明确的计算了各种局部平滑滤波方法的 method noise 理论值。图 $(4)$ 的可视化实验证明了前文中公式的正确性。图 $(4)$ 对比了数种降噪方法对含噪声 lena 图像运算的 method noise，噪声为高斯白噪声，标准差 2.5，滤波参数 $h$ 均相同。Method noise 能更好的协助判断降噪算法的性能及局限，因为被移除或改变的纹理、细节将被 method noise 醒目的展示出来。对比图 $(4)$ 中数种降噪算法，NL-means 算法的 method noise 几乎无法察觉任何几何纹理。图 $(2)$ 可以解释这一现象，因为 NL-means 算法选取的权重因子完全适应图像局部及非局部的几何结构。

图 4. 降噪算法的 method noise。从左到右，由上至下：噪声图像（标准差 20），高斯滤波，各向异性滤波，TV 滤波，邻域滤波及 NL-means 算法。视觉实验验证了第二节的计算公式。

由于算法本身的性质，纹理及周期性结构是 NL-means 算法最适用的情况。因为对任意像素 $i$ ，纹理图像或周期性图像中将可以找到大量与该像素具有相似邻域的点，如图 $(2.e)$ 。图 $(3)$ 展示了局部平滑滤波器及 NL-means 算法对自然纹理的平滑效果。

自然图像同样含有足够的信息冗余会被 NL-means 恢复。平坦区域内部将呈现出大量相似的几何结构，见图 $(2.a)$ 。直边界、或弯曲边界将被筛选出一条结构相似的像素线，见图 $(2. b)，(2.c)$ 。并且，NL-means 将会在很远的位置寻找到与参考点相似的结构，如图 $(2.f)$ 。图 $(5)$ 展示了一次对自然图像的可视化实验，这组结果与图 $(4)$ 是相对应的。