图像去噪是非常基础也是非常必要的研究，去噪常常在更高级的图像处理之前进行，是图像处理的基础。可惜的是，目前去噪算法并没有很好的解决方案，实际应用中，更多的是在效果和运算复杂度之间求得一个平衡，再一次验证了我老师的一句话：
好了，废话不多说，来看看效果比较好的去噪算法吧。
图像中噪声的来源有许多种，这些噪声来源于图像采集、传输、压缩等各个方面。噪声的种类也各不相同，比如椒盐噪声，高斯噪声等，针对不同的噪声有不同的处理算法。
对于输入的带有噪声的图像v(x)，其加性噪声可以用一个方程来表示： $$v(x) = u(x) + \eta (x),\quad x \in \Omega ,$$
其中$u(x)$是原来没有噪声的图像。$x$是像素集合，$\eta (x)$是加项噪声项，代表噪声带来的影响。$\Omega$是像素的集合，也就是整幅图像。从这个公式可以看出，噪声是直接叠加在原始图像上的，这个噪声可以是椒盐噪声、高斯噪声。理论上来说，如果能够精确地获得噪声，用输入图像减去噪声就可以恢复出原始图像。但现实往往很骨感，除非明确地知道噪声生成的方式，否则噪声很难单独求出来。
工程上，图像中的噪声常常用高斯噪声${\rm N}(\mu ,{\sigma ^2})$来近似表示，其中$\mu = 0$，${\sigma ^2}$是噪声的方差，${\sigma ^2}$越大，噪声越大。一个有效的去除高斯噪声的方式是图像求平均，对N幅相同的图像求平均的结果将使得高斯噪声的方差降低到原来的N分之一，现在效果比较好的去噪算法都是基于这一思想来进行算法设计。
NL-Means的全称是：Non-Local Means，直译过来是非局部平均，在2005年由Baudes提出，该算法使用自然图像中普遍存在的冗余信息来去噪声。与常用的双线性滤波、中值滤波等利用图像局部信息来滤波不同的是，它利用了整幅图像来进行去噪，以图像块为单位在图像中寻找相似区域，再对这些区域求平均，能够比较好地去掉图像中存在的高斯噪声。NL-Means的滤波过程可以用下面公式来表示： $$\tilde u(x) = \sum\limits_{y \in {\Omega _x}} {w(x,y)v(y)}.$$
在这个公式中，$w(x,y)$是一个权重，表示在原始图像$v$中，像素 $x$ 和像素 $y$ 的相似度。这个权重要大于0，同时，权重的和为1，用公式表示是这样： $$w(x,y) > 0\quad and\quad \sum\limits_{y \in {\Omega _x}} {w(x,y) = 1}, \quad \forall x \in \Omega ,y \in {\Omega _x} .$$
$\Omega _x$是像素 $x$ 的邻域。这个公式可以这样理解：对于图像中的每一个像素 $x$ ，去噪之后的结果等于它邻域中像素 $y$ 的加权和，加权的权重等于 $x$ 和 $y$ 的相似度。这个邻域也称为搜索区域，搜索区域越大，找到相似像素的机会也越大，但同时计算量也是成指数上升。在提出这个算法的文献中，这个区域是整幅图像！导致的结果是处理一幅512x512大小的图像，最少也得几分钟。
衡量像素相似度的方法有很多，最常用的是根据两个像素的亮度值的差的平方来估计。但因为有噪声的存在，单独的一个像素并不可靠。对此解决方法是，考虑它们的邻域，只有邻域相似度高才能说这两个像素的相似度高。衡量两个图像块的相似度最常用的方法是计算他们之间的欧氏距离： $$w(x,y) = {1 \over {n(x)}}\exp ({{\left\| {{\bf{V}}(x) - {\bf{V}}(y)} \right\|_{2,a}^2} \over {{h^2}}}).$$
其中： $n(x)$ 是一个归一化的因子，是所有权重的和，对每个权重除以该因子后，使得权重满足和为1的条件。 $h>0$ 是滤波系数，控制指数函数的衰减从而改变欧氏距离的权重。 ${\bf{V}}(x)$ 和 ${\bf{V}}(y)$ 代表了像素 $x$ 和像素 $y$ 的邻域，这个邻域常称为块(Patch)邻域。块邻域一般要小于搜索区域。${\left\| {{\bf{V}}(x) - {\bf{V}}(y)} \right\|_{2,a}^2}$ 是两个邻域的高斯加权欧式距离。其中 $a>0$ 是高斯核的标准差。在求欧式距离的时候，不同位置的像素的权重是不一样的，距离块的中心越近，权重越大，距离中心越远，权重越小，权重服从高斯分布。实际计算中考虑到计算量的问题，常常采用均匀分布的权重。
讲了这么多，是时候用图来说明问题了：
如上图所示，p为去噪的点，因为q1和q2的邻域与p相似，所以权重$w(p,q1)$和$w(p,q2)$比较大，而邻域相差比较大的点q3的权重值$w(p,q3)$很小。如果用一幅图把所有点的权重表示出来，那就得到下面这些权重图：
这6组图像中，左边是原图，中心的白色色块代表了像素 $x$ 块邻域，右边是计算出来的权重 $w(x,y)$ 图，权重范围从0（黑色）到1（白色）。这个块邻域在整幅图像中移动，计算图像中其他区域跟这个块的相似度，相似度越高，得到的权重越大。最后将这些相似的像素值根据归一化之后的权重加权求和，得到的就是去噪之后的图像了。
这个算法参数的选择也有讲究，一般而言，考虑到算法复杂度，搜索区域大概取21x21，相似度比较的块的可以取7x7。实际中，常常需要根据噪声来选取合适的参数。当高斯噪声的标准差 $\sigma$ 越大时，为了使算法鲁棒性更好，需要增大块区域，块区域增加同样也需要增加搜索区域。同时，滤波系数 $h$ 与 $\sigma$ 正相关：$h=k\sigma $，当块变大时，$k$ 需要适当减小。
NL-Means算法的复杂度跟图像的大小、颜色通道数、相似块的大小和搜索框的大小密切相关，设图像的大小为N$\times$N，颜色通道数为$N_c$，块的大小为k$\times$k，搜索框的大小为n$\times$n，那么算法复杂度为：$O({N^2}{N_c}{k^2}{n^2})$。对512$\times$512的彩色图像而言，设置k=7，n=21，OpenCV在使用了多线程的情况下，处理一幅图像所需要的时间需要几十秒。虽然有人不断基于这个算法进行改进、提速，但离实时处理还是比较远。
最后来看一下这个算法的去噪效果[3]：
BM3D（Block-matching and 3D filtering，3维块匹配滤波）可以说是当前效果最好的算法之一。该算法的思想跟NL-Means有点类似，也是在图像中寻找相似块的方法进行滤波，但是相对于NL-Means要复杂得多，理解了NL-Means有助于理解BM3D算法。BM3D算法总共有两大步骤，分为基础估计（Step1）和最终估计（Step2）：
BM3D算法流程图
在这两大步中，分别又有三小步：相似块分组（Grouping），协同滤波（Collaborative Filtering）和聚合（Aggregation）。上面的算法流程图已经比较好地将这一过程表示出来了，只需要稍加解释。
(1) Grouping：有了NL-Means的基础，寻找相似块的过程很容易理解。首先在噪声图像中选择一些$k \times k$ 大小的参照块（考虑到算法复杂度，不用每个像素点都选参照块，通常隔3个像素为一个不长选取，复杂度降到1/9），在参照块的周围适当大小（$n \times n$）的区域内进行搜索，寻找若干个差异度最小的块，并把这些块整合成一个3维的矩阵，整合的顺序对结果影响不大。同时，参照块自身也要整合进3维矩阵，且差异度为0。寻找相似块这一过程可以用一个公式来表示：$$G(P) = \{ Q:d(P,Q) \le {\tau ^{step1}}\} .$$
d(P,Q)代表两个块之间的欧式距离。最终整合相似块获得的矩阵就是流程图Step1中左下角的蓝色R矩阵。
(2) Collaborative Filtering：形成若干个三维的矩阵之后，首先将每个三维矩阵中的二维的块（即噪声图中的某个块）进行二维变换，可采用小波变换或DCT变换等，通常采用小波BIOR1.5。二维变换结束后，在矩阵的第三个维度进行一维变换，通常为阿达马变换（Hadamard Transform）。变换完成后对三维矩阵进行硬阈值处理，将小于阈值的系数置0，然后通过在第三维的一维反变换和二维反变换得到处理后的图像块。这一过程同样可以用一个公式来表达：$$Q(P) = T_{3Dhard}^{ - 1}(\gamma (T_{3Dhard}(Q(P)))).$$
在这个公式中，二维变换和一维变换用一个$T_{3Dhard}$ 来表示。$\gamma$是一个阈值操作：$$\gamma (x) = \left\{ {\matrix{ {0\quad if{\kern 1pt} \left| x \right| \le {\lambda _{3D}}\sigma } \cr {x\quad \quad \; otherwise} \cr } } \right. .$$
$\sigma$是噪声的标准差，代表噪声的强度。
(3) Aggregation：此时，每个二维块都是对去噪图像的估计。这一步分别将这些块融合到原来的位置，每个像素的灰度值通过每个对应位置的块的值加权平均，权重取决于置0的个数和噪声强度。
(1) Grouping：第二步中的聚合过程与第一步类似，不同的是，这次将会得到两个三维数组：噪声图形成的三维矩阵${Q^{basic}}(P)$和基础估计结果的三维矩阵$Q(P)$。
(2) Collaborative Filtering：两个三维矩阵都进行二维和一维变换，这里的二维变换通常采用DCT变换以得到更好的效果。用维纳滤波（Wiener Filtering）将噪声图形成的三维矩阵进行系数放缩，该系数通过基础估计的三维矩阵的值以及噪声强度得出。这一过程同样可以用一个公式来表达：$$Q(P) = T_{3Dwein}^{ - 1}(w_p \cdot {T_{3Dwein}}(Q(P))).$$
在这个公式中，二维变换和一维变换用一个$T_{3Dwein}$ 来表示。$w_p$是一个维纳滤波的系数：$${w_p}(\xi ) = {{{{\left| {\tau _{3D}^{wien}({Q^{basic}}(P))(\xi )} \right|}^2}} \over {{{\left| {\tau _{3D}^{wien}({Q^{basic}}(P))(\xi )} \right|}^2} + {\sigma ^2}}}$$
$\sigma$是噪声的标准差，代表噪声的强度。
(3) Aggregation：与第一步中一样，这里也是将这些块融合到原来的位置，只是此时加权的权重取决于维纳滤波的系数和噪声强度。
经过最终估计之后，BM3D算法已经将原图的噪声显著地去除。可以来看一组结果：
该算法的主要运算量还是在相似块的搜索与匹配上，在与NL-Means同样大小的相似块和搜索区域的情况下，BM3D的算法复杂度是要高于NL-Means的，应该大概在NL-Means的3倍左右。梦想着实时处理的同学可以死心了。
要比较算法效果，必然离不开评价体系。由于人带有主观因素，每个人的评价可能都不一样，因此有必要用几种客观的评价方法来对结果进行评价。目前，用得比较多的评价方式是MSE(Mean-Squared Error，均方误差)和PSNR(Peak Signal-to-Noise Ratio，峰值信噪比)。
两幅$M \times N$大小的图像$u(x,y)$和$v(x,y)$的MSE计算公式如下：$${e_{MSE}} = {1 \over {MN}}\sum\limits_{n = 1}^N {\sum\limits_{m = 1}^M {{{[u(n,m) - v(n,m)]}^2}} } .$$
在这个公式里没有表现出像素值范围对结果的影响，同样的均方误差8-bit的图像和12-bit的图像显然没有可比性。因此，又引入了峰值信噪比：$$PSNR = - 10{\log _{10}}{{{e_{MSE}}} \over {{s^2}}}.$$
上式中，$s$ 是图像像素最大值，对于8-bit的图像而言 $s$=255，PSNR的单位是分贝(dB)。 通常 PSNR 值越高表示品质越好，一般而言，当 PSNR<30dB 时，代表以人的肉眼看起来是不能容忍的范围。因此大部分PSNR值都要>30dB。但PSNR高，并不代表图像质量一定好，有时候还是必须要靠人的肉眼去辅助判断图像的质量才较为正确。
不同PSNR对应的视觉效果
我对上面两种方法获得的结果针对原图计算了PSNR，结果如下：
| NL-Means | BM3D | |
| PSNR | 32.0913 | 33.6711 |
NL-Means和BM3D可以说是目前效果最好的去噪算法，其中BM3D甚至宣称它可以得到迄今为止最高的PSNR。从最终的结果也可以看出来，BM3D的效果确实要好于NL-Means，噪声更少，能够更好地恢复出图像的细节。在效果这一点上BM3D胜。无愧于State-of-the-art这一称号。当然，这里进行测试的样本比较少，可能还不足以完全说明问题。
这两者可以说是目前最有效的图像去噪算法了，但是都不可避免地要面对一个同问题：尽管计算机性能已经成百上千倍地提高，还是远不能满足很多算法的实时计算的需求，这很大程度上限制了这些算法的使用范围：用户无法处理一张照片需要等待长达几分钟的时间，因此，距离真正意义上的实用还是有一段距离。我们只能期待，要是有一天计算机性能不再是问题，又或者，大牛们能够研究出又快又好的算法吧。
[1] Buades A, Coll B, Morel J M. A non-local algorithm for image denoising[C]//Computer Vision and Pattern Recognition, 2005. CVPR 2005. IEEE Computer Society Conference on. IEEE, 2005, 2: 60-65.
[2] Buades A, Coll B, Morel J M. Nonlocal image and movie denoising[J]. International journal of computer vision, 2008, 76(2): 123-139.
[3] Antoni Buades, Bartomeu Coll, and Jean-Michel Morel, Non-Local Means Denoising, Image Processing On Line, 1 (2011). http://dx.doi.org/10.5201/ipol.2011.bcm_nlm
[4] Jacques Froment, Parameter-Free Fast Pixelwise Non-Local Means Denoising, Image Processing On Line, 4 (2014), pp. 300–326. http://dx.doi.org/10.5201/ipol.2014.120
[5] Dabov K, Foi A, Katkovnik V, et al. Image denoising by sparse 3-D transform-domain collaborative filtering[J]. Image Processing, IEEE Transactions on, 2007, 16(8): 2080-2095.
[6] http://www.cs.tut.fi/~foi/GCF-BM3D/
[7] Marc Lebrun, An Analysis and Implementation of the BM3D Image Denoising Method, Image Processing On Line, 2 (2012), pp. 175–213. http://dx.doi.org/10.5201/ipol.2012.l-bm3d
Last update: 2014-11-21