图像去噪-传统方法总结

目录

.1图像去噪的重要性... 2

.1图像去噪的分类... 3

噪声的分类... 3

图像去噪方法分类... 3

.2传统方法调研... 4

.2.1空域像素特征去噪算法... 4

算术均值滤波Arithmetic Mean Filter与几何均值滤波Geometric Mean Filter 5

高斯滤波Gaussian filter 6

统计中值滤波Median filtering. 9

双边滤波Bilateral Filter 10

引导滤波Guided Filter 12

NLM非局部均值滤波Non-Local means. 16

.2.2变换域去噪算法... 19

BM3D去噪算法... 19

.1图像去噪的重要性

图像、视频从采集到播放的整个生命周期中会经历各种各样的处理过程，比如采集、剪辑、编码、转码、传输、显示等，每个处理过程都会引入失真。“噪声”就是在信号采集过程中引入的一种普遍失真。许多重要的图像包含一定程度的噪声，即图像强度中无法解释的数据干扰，这些干扰要么无法解释，要么使用者不感兴趣，图像去噪就是去除噪声图像中的噪声，从而恢复真实的图像。降低噪声强度可以使图像主观效果更好。另外，在图像、视频压缩时也不必浪费码率在编码噪声上。同时，会使得视频编码中的运动估计更准确、熵编码速度更快。

图像降噪的英文名称是Image Denoising， 图像处理中的专业术语。是指减少数字图像中噪声的过程，有时候又称为图像去噪。图像的噪声来源相对复杂，搞清楚图像噪声的成因对我们进行图像去噪的工作有帮助。因为对于满足某些数学统计规律的噪声，那么逆向去除的时候就变得好办了。另外，深度学习技术很多也用在图像去噪领域，深度学习依赖数据，明白噪声的分布特点有利于制作数据集。

图像去噪是许多其他图像工作的基础，如边缘提取，图像分割，图像识别等等，都需要先经过图像去噪，去除干扰信息。

.1图像去噪的分类

噪声的分类

我们一般将噪声分为三类：加性噪声，乘性噪声和量化噪声。用

表示给定图像，

表示图像信号，

表示噪声。

1. 加性噪声，此类噪声与输入图像信号无关，含噪图像可表示为

   

   ,

信道噪声及光导摄像管的摄像机扫描图像时产生的噪声就属这类噪声；

2.乘性噪声，此类噪声与图像信号有关，含噪图像可表示为

,

飞点扫描器扫描图像时的噪声，电视图像中的相干噪声，胶片中的颗粒噪声就属于此类噪声；

3.量化噪声，此类噪声与输入图像信号无关，是量化过程存在量化误差，再反映到接收端而产生。

图像去噪方法分类

目前来说图像去噪分为三大类：基于滤波器的方法（Filtering-Based Methods）、基于模型的方法（Model-Based Methods）和基于学习的方法（Learning-Based Methods）。或者分为（更便于理解一些）：基于人工特征的传统去噪方法和基于深度学习的去噪方法。

相对于传统图像去噪方法，基于深度学习的图像去噪方法具有强大的学习能力，不仅可以拟合复杂噪声分布，还节省了计算时间。

.2传统方法调研

.2.1空域像素特征去噪算法

在信号处理教科书中，虽然介绍过很多经典的图像去噪方法，但主要都是针对随机噪声的，对于sensor缺陷导致的一些脉冲噪声（impulse noise）这里我们不考虑。

什么是随机噪声呢？相比于图像的真实信号来说随机噪声就是一种或高或低呈现出不确定变化的一种信号，如下图所示虚线代表真实信号，红蓝线表示的就是随机噪声信号，所有的随机噪声信号求和后结果为0。

由于这个零和特点，目前几乎所有的空域降噪算法都是基于这个理论为出发点来进行降噪处理的。

基于空域像素特征的方法，是通过分析在一定大小的窗口内，中心像素与其他相邻像素之间在灰度空间的直接联系，来获取新的中心像素值的方法，因此往往都会存在一个典型的输入参数，即滤波半径r。此滤波半径可能被用于在该局部窗口内计算像素的相似性，也可能是一些高斯或拉普拉斯算子的计算窗口。最具有代表性的滤波方法有以下几种：

算术均值滤波Arithmetic Mean Filter与几何均值滤波Geometric Mean Filter

这是最简单的均值滤波器，可以去除均匀噪声和高斯噪声，但会对图像造成一定程度的模糊。令

表示中心点在

处，大小为

的滤波器窗口。算术均值滤波器就是简单的计算窗口区域的像素均值，然后将均值赋值给窗口中心点处的像素：

其中，

表示原始图像，

表示均值滤波后得到的图片。

滤波后图像的像素由模板窗口内像素的乘积的

幂给出。和算术均值滤波器相比，几何均值滤波器能够更好的去除高斯噪声，并且能够更多的保留图像的边缘信息。但其对0值是非常敏感的，在滤波器的窗口内只要有一个像素的灰度值为0，就会造成滤波器的输出结果为0。

算数均值滤波在去除椒盐噪声的同时，图像也变得模糊。几何均值滤波去除噪声同时，相对于算数均值滤波锐化。

均值滤波结果随着滤波半径取值的增大而变得越来越模糊，图像对比度越来越小。经过均值处理之后，噪声部分被弱化到周围像素点上，所得到的结果是噪声幅度减小，但是噪声点的颗粒面积同时变大，所以污染面积反而增大。为了解决这个问题，可以通过设定阈值，比较噪声和邻域像素灰度，只有当差值超过一定阈值时，才被认为是噪声。不过阈值的设置需要考虑图像的总体特性和噪声特性，进行统计分析。自适应均值滤波算法通过方向差分来寻找噪声像素，从而赋予噪声像素与非噪声像素不同的权重，并自适应地寻找最优窗口大小，优于一般的均值滤波方法。

高斯滤波Gaussian filter

高斯滤波(Gaussian filter) 包含许多种，包括低通、带通和高通等，我们通常图像上说的高斯滤波，指的是高斯模糊(Gaussian Blur) 是一种高斯低通滤波，其过滤掉图像高频成分（图像细节部分），保留图像低频成分（图像平滑区域），所以对图像进行‘高斯模糊’后，图像会变得模糊。高斯模糊对于抑制高斯噪声 (服从正态分布的噪声) 非常有效。

高斯滤波性质

高斯函数具有五个重要的性质，这些性质使得它在早期图像处理中特别有用．这些性质表明，高斯平滑滤波器无论在空间域还是在频率域都是十分有效的低通滤波器，且在实际图像处理中得到了工程人员的有效使用．高斯函数具有五个十分重要的性质，它们是：

1.二维高斯函数具有旋转对称性，即滤波器在各个方向上的平滑程度是相同的．一般来说，一幅图像的边缘方向是事先不知道的，因此，在滤波前是无法确定一个方向上比另一方向上需要更多的平滑．旋转对称性意味着高斯平滑滤波器在后续边缘检测中不会偏向任一方向．

2.高斯函数是单值函数．这表明，高斯滤波器用像素邻域的加权均值来代替该点的像素值，而每一邻域像素点权值是随该点与中心点的距离单调增减的．这一性质是很重要的，因为边缘是一种图像局部特征，如果平滑运算对离算子中心很远的像素点仍然有很大作用，则平滑运算会使图像失真．

3.高斯函数的傅立叶变换频谱是单瓣的．正如下面所示，这一性质是高斯函数付立叶变换等于高斯函数本身这一事实的直接推论．图像常被不希望的高频信号所污染(噪声和细纹理)．而所希望的图像特征（如边缘），既含有低频分量，又含有高频分量．高斯函数付立叶变换的单瓣意味着平滑图像不会被不需要的高频信号所污染，同时保留了大部分所需信号．

4.高斯滤波器宽度(决定着平滑程度)是由参数σ表征的，而且σ和平滑程度的关系是非常简单的.σ越大，高斯滤波器的频带就越宽，平滑程度就越好．通过调节平滑程度参数σ，可在图像特征过分模糊(过平滑)与平滑图像中由于噪声和细纹理所引起的过多的不希望突变量(欠平滑)之间取得折衷．

5.由于高斯函数的可分离性，较大尺寸的高斯滤波器可以得以有效地实现．二维高斯函数卷积可以分两步来进行，首先将图像与一维高斯函数进行卷积，然后将卷积结果与方向垂直的相同一维高斯函数卷积．因此，二维高斯滤波的计算量随滤波模板宽度成线性增长而不是成平方增长．

高斯模糊原理

要模糊一张图像，可以这么做：对于每个像素点，以它为中心，取其3x3区域内所有像素灰度值的平均作为中心点的灰度值。可是，如果仅使用简单平均，显然不是很合理，因为图像都是连续的，越靠近的点关系越密切，越远离的点关系越疏远。因此，加权平均更合理，距离越近的点权重越大，距离越远的点权重越小。而正态分布显然是一种可取的权重分配模式。由于图像是二维的，所以需要使用二维的高斯函数。

高斯模糊本质上就是利用高斯函数生成的高斯核(高斯模板)对图像进行卷积操作。

高斯函数

在计算每个像素时，都把当前中心点看作坐标原点，可以使得均值

，简化高斯函数的公式：

一维高斯分布:

二维高斯分布:

高斯核(高斯模板)

理论上，高斯分布在所有定义域上都有非负值，这就需要一个无限大的卷积核。实际上，仅需要取均值周围3倍标准差(即3

)内的值，以外部分直接去掉即可。

高斯滤波的重要两步就是先找到高斯模板然后再进行卷积，模板（有的地方也称作掩膜或者是高斯核）。

假定中心点的坐标是（0,0），那么取距离它最近的8个点坐标，为了计算，需要设定

的值。假定

=1.5，则模糊半径为1的高斯模板就计算如下：

这个时候我们我们还要确保这九个点加起来为1（即归一化，这个是高斯模板的特性）否则的话，使用总值大于1的模板会让图像偏亮，小于1的模板会让图像偏暗。这9个点的权重总和等于0.4787147，因此上面9个值还要分别除以0.4787147，得到最终的高斯模板：

高斯滤波计算

有了高斯模板，就可以利用高斯模板对图像进行卷积了。对图像进行卷积的原理，其实就是将模板作为权值，与对应像素相乘再求和，得到的结果就是中心点卷积后的结果，本质上就是加权平均操作。

假设现有9个像素点，灰度值（0-255）的高斯滤波计算如下：

将这9个值加起来，就是中心点25进行高斯滤波后的值。对所有点重复这个过程，就得到了高斯模糊后的图像。

边界处理

高斯模板在对图像边缘像素进行卷积时，会有一部分权重没有对应像素，因此我们需要在图像的边缘补0 00。这种方法称作Zero Padding。并且权值g（卷积核）要进行归一化操作。若没有进行归一化，滤波后图像边缘像素会比较暗。

优化：高斯函数分离特性

直接进行二维高斯模糊效率较低，实际上高斯模糊也可以在二维图像上对两个独立的一维空间分别进行计算，这叫作线性可分：

上式说明，可以先对图像进行y(或x)轴方向的一维高斯模糊，再在得到的结果上进行x(或y)轴方向的一维高斯模糊，其结果与直接对图像进行二维高斯模糊效果相同。

将二维高斯模糊分离为两个一维高斯模糊后，时间复杂度从

降到

，这里

和

表示图像宽高，

表示图像通道数，一般彩色图是3通道，

表示高斯模板的宽高，一般为奇数，例如上图的高斯核为5x5大小，r=5

统计中值滤波Median filtering

在使用邻域平均法去噪的同时也使得边界变得模糊。而中值滤波是非线性的图像处理方法，在去噪的同时可以兼顾到边界信息的保留。选一个含有奇数点的窗口

，将这个窗口在图像上扫描，把窗口中所含的像素点按灰度级的升或降序排列，取位于中间的灰度值来代替该点的灰度值。计算过程如下图所示：

双边滤波Bilateral Filter

双边滤波（Bilateral filter）是一种非线性的滤波方法，是结合图像的空间邻近度和像素值相似度的一种折衷处理，同时考虑空域信息和灰度相似性，达到保边去噪的目的。具有简单、非迭代、局部的特点。双边滤波器的好处是可以做边缘保存（edge preserving），一般用高斯滤波去降噪，会较明显地模糊边缘，对于高频细节的保护效果并不明显。双边滤波器顾名思义比高斯滤波多了一个高斯方差sigma－d，它是基于空间分布的高斯滤波函数，所以在边缘附近，离的较远的像素不会太多影响到边缘上的像素值，这样就保证了边缘附近像素值的保存。但是由于保存了过多的高频信息，对于彩色图像里的高频噪声，双边滤波器不能够干净的滤掉，只能够对于低频信息进行较好的滤波。

原理：

与高斯滤波的高斯核结构一直保持不变不同，为了使图像的边缘得到保留，就要根据当前被卷积像素的邻域进行观察，“推断”是否是边缘点和接近边缘的点。因此，结构元素就会改变，从而保留边缘点。下面的一组图中，图a是原图像，图c是输出。而中间的图像是什么呢？显然，这是原图中根据某个点的邻域生成的，专属于这个点的结构元素。

可以看到，原图中显然有一个灰度的突变，这就表示是边缘。灰度值高的地方不应该和灰度低的区域进行混合，所以，图像中接近边缘的一个点就会生成图b这样的结构元素。那么这个接近边缘的点在哪里呢？大概就在标出的这个区域。

而生成这样的结构元素的方法，是将我们原本的高斯核，与一个能“推断”出是否在边缘点的结构元素相乘，如下图中间的结构元素。

示意图如下：

双边滤波器的输出像素依赖于当前被卷积像素的邻域。

和

是当前被卷积像素的坐标点，

和

是邻域像素的坐标点:

可以看到，它 取决于被卷积像素的灰度值和邻域像素的灰度值的差。我们知道，边缘会有较大的灰度变化，而这个公式就会使边缘和边缘另一边的区域生成比较小的权值，与被卷积像素的灰度值类似的区域会生成比较大的权值，就像之前图中的一个“断崖”。

相乘就得到加权系数

，即新的核:

引导滤波Guided Filter

高斯滤波等线性滤波算法所用的核函数相对于待处理的图像是独立无关的，这里的独立无关也就意味着，对任意图像都是采用相同的操作。

引导滤波就是在滤波过程中加入引导图像中的信息，这里的引导图可以是单独的图像也可以是输入图像，当引导图为输入图像时，引导滤波就成为了一个可以保持边缘的去噪滤波操作。

除了速度优势以外，导向滤波的一个很好的性能就是可以保持梯度，这是bilateral做不到的，因为会有梯度翻转现象。

基本原理：

第一步：假设该引导滤波函数的输出与输入在一个二维窗口内满足线性关系，如下：

其中，

为输入图像，

为导向图，

为输出图像。在这里我们认为输出图像可以看成导向图

的一个局部线性变换，其中

是局部化的窗口的中点，因此属于窗口

的pixel，都可以用导向图对应的pixel通过

的系数进行变换计算出来。同时，我们认为输入图像

是由

加上我们不希望的噪声或纹理得到的，因此有

和

是像素索引，

和

是当窗口中心位于

时该线性函数的系数。（当引导图为输入图像时，引导滤波就成为一个保持边缘的滤波操作，即

对上式两边取梯度可得

即当输入图

有梯度时，输出

也有类似的梯度，这也就可以解释为什么引导滤波有边缘保持特性了。

第二步：

接下来就是解出这样的系数，使得

和

的差别尽量小，而且还可以保持局部线性模型。这里利用了带有正则项的 linear ridge regression（岭回归）

求解以上方程得到

和

在局部的值，对于一个要求的pixel可能含在多个window中，因此平均后得到：

在这里，

和

表示

在局部窗口

中的均值和方差。

是窗口内的所有像素数，

表示

在窗口

中的均值，

就是规整化参数，当

时，上面公式即可简化为：

• 如果

  

  =0，显然

  

  =1,

  

  =0是

  

  为最小值的解，从上式可以看出，这时的滤波器没有任何作用，将输入原封不动地输出。
• 如果

  

  >0，在像素强度变化小的区域（方差不大），即图像

  

  在窗口

  

  中基本保持固定，此时有

  

  <<

  

  ，于是有

  

  ≈0和

  

  ≈

  

  即做了一个加权均值滤波，而在高方差区域，即表示图像

  

  在窗口

  

  中变化比较大，此时我们有

  

  >>

  

  ，于是有

  

  ≈1和

  

  ≈0，对图像的滤波效果很弱，有助于保持边缘。
• 在窗口大小不变的情况下，随着

  

  的增大，滤波效果越明显。

第三步：

在计算每个窗口的线性系数时，我们可以发现一个像素会被多个窗口包含，也就是说，每个像素都由多个线性函数所描述。因此，如之前所说，要具体求某一点的输出值

时，只需将所有包含该点的线性函数值平均即可，如下：

其中，输出值

又与两个均值有关，分别为

和

在窗口

中的均值，我们将上一步得到两个图像

和

都进行盒式滤波，得到两个新图：

和

.然后用

乘以引导图像

，再加上

，即得最终滤波之后的输出图像

.

NLM非局部均值滤波Non-Local means

前面基于邻域像素的滤波方法，基本上只考虑了有限窗口范围内的像素灰度值信息，没有考虑该窗口范围内像素的统计信息如方差，也没有考虑整个图像的像素分布特性，和噪声的先验知识。针对其局限性，NLM算法被提出，该算法使用自然图像中普遍存在的冗余信息来去噪声。与常用的双线性滤波、中值滤波等利用图像局部信息来滤波不同的是，它利用了整幅图像来进行去噪，以图像块为单位在图像中寻找相似区域，再对这些区域求平均，能够比较好地去掉图像中存在的高斯噪声。这里直接用图片来更好地说明问题：

如上图所示，其中p为去噪的点，从图中看出q1和q2的邻域与p相似，所以权重和比较大，而q3因为与q邻域相差比较大所以赋予的权重值就很小。NLM就是将一幅图像中所有点的权重都表示出来，那就得到下面这些权重图：

上面权值图像中，左边是原图，中心的白色色块代表了像素块邻域，右边是计算出来的权重图，权重范围从0（黑色）到1（白色）。

这个块邻域在整幅图像中移动，计算图像中其他区域跟这个块的相似度，相似度越高，得到的权重越大。最后将这些相似的像素值根据归一化之后的权重加权求和，得到的就是去噪之后的图像了。

由于原始NLM方法需要用图像中所有的像素来估计每一个像素的值，因此计算量非常大，研究者不断对该方法进行了几点改进。(a)采用一定的搜索窗口代替所有的像素，使用相似度阈值，对于相似度低于某一阈值的像素，不加入到权重的计算(即不考虑其相对影响，这些都可以降低计算复杂度。(b)使用块之间的显著特征，如纹理特征等代替灰度值的欧氏距离来计算相似度，在计算上更加有优势，应用上也更加灵活。

基本原理：

考虑到了图像的自相似性质，它充分利用了图像中的冗余信息，在去噪的同时能够最大程度的保持图像的细节特征。

该算法需要计算图像中所有像素与当前像素之间的相似性，考虑到这个计算量与效率的问题，一般会设定两个固定大小的窗口，一个大的搜索窗口(D×D)和一个小的邻域窗口(d×d)，邻域窗口在搜索窗口中进行滑动，根据邻域间的相似性来确定对应中心像素对当前像素的影响度，也就是权值。

下图是NLM算法执行过程，大窗口是以目标像素为中心的搜索窗口，两个灰色小窗口分别是以x，y为中心的邻域窗口。其中以y为中心的邻域窗口在搜索窗口中滑动，通过计算两个邻域窗口间的相似程度为y赋以权值w(x,y)。

假设含噪声图像为

，去噪图像为

，

中

处的像素点灰度值计算如下：

其中权值

表示像素点

和

间的相似度，它的值由以

、

为中心的矩形领域

、

间的距离

决定：

其中，

为归一化系数，

为平滑参数，

为小搜索窗区域，控制高斯函数的衰减程度。

越大高斯函数变化越平缓，去噪水平越高，但同时也会导致图像越模糊；

越小，边缘细节成分保持得越多，但会残留过多的噪声点；

的具体取值应当以图像中的噪声水平为依据。

.2.2变换域去噪算法

BM3D去噪算法

原理详解

1 加性高斯白噪声

       BM3D主要用于去除图像中的加性高斯白噪声（Additive White Gaussian Noise, AWGN）。这里主要涉及三个概念：

     1) 加性，即噪声其对原始信号的影响可表示为线性叠加，对于AWGN，我们一般假设其与原始信号无关，在同一幅图像中，无论各像素点的明暗如何，其噪声的分布函数都是一致的，即独立同分布（i.i.d.）。另一种常见的加性噪声是散粒噪声，主要由光子的量子特性所造成，像素单元实际接收到的光子数量符合泊松分布，是一种与信号相关的噪声，像素本身的亮度越高，噪声越大，但整体的信噪比却是在提高的，这里不详细讨论。

     2) 高斯，即噪声瞬时值服从高斯分布。之所以为高斯分布，是因为图像成像过程中会受到很多因素的影响，例如电子的热扰动，信号的放大与量化误差，光子量子特性所带来的散粒噪声，黑体辐射，太阳或其他宇宙天体的活动等等，根据中心极限定理，大量独立随机事件的影响的总和将会趋近于高斯分布，这里通常假设其均值为 0。高斯分布的方差与所处的光照环境、相机传感器质量、快门速度、感光度等等相关，在相同条件下一般认为是一样的。

     3) 白噪声，即噪声功率谱密度服从均匀分布。类似于白光包含了可见光中的所有频率，白噪声也包含了所有频率的波，并且在各频率上的功率谱密度都是一样的，这表明 AWGN 的强度并不会与像素的颜色相关（但实际上图像像素的颜色是由多种频率的光混合而成的，并不是严格的单频率光，所以也不是严格的无关）。另外，因为频域中的均匀分布对应着时域中的冲激函数，所以当前时刻的噪声不会与其他时刻相关，例如前一时刻的噪声值很大，并不意味当前时刻的噪声值也会很大，也就是说在相同拍摄条件下不同时刻拍摄的图像也是符合独立同分布假设的。注意，这里所说的时刻是指拍摄单张图像的整个过程，而不是单纯的某一瞬间。例如相机内部传感器以及其他一些信号处理器件的温度会随着拍摄时间的延长而升高，从而带来更严重的热噪声，所以在这个拍摄过程中不同瞬间的噪声的强度是有一定相关性的。

       显然，AWGN 是一种理想化的模型，真实图像的噪声远要比以上的假设复杂。但大量的实验与经验也表明，在低照度以及高温的情况下，高斯噪声占据了图像噪声的主要部分，这时使用 AWGN 一方面能够较好地反映图像的真实噪声，另一方面也方便了理论上的分析以及各种算法实验的展开。例如，想要获取噪声图像对应的原始无噪声图像并不是一件易事，而在有监督深度学习中，我们需要大量的样本来训练我们的网络，这时使用人工生成的 AWGN 就能大大地减轻数据收集的工作压力，同时也能在一定程度上保证该降噪算法对于真实噪声图像的适用性。

       综上所述，记真实无噪声图像为

，AWGN 噪声为

，有噪声图像为

，那么三者的关系可表示为

  （1.1）

其中

代表像素点位置。标准差

反映了噪声值的水平，根据其概率密度函数，99.6% 的噪声值都会在3

的范围之内。在这里，噪声

为随机变量，而真实图像

为常量，所以有

  （1.2）

那么，对于有限个相互独立的像素点的线性叠加，有

 （1.3）

这也是大多数 AWGN 降噪算法的基本理论依据，即多个独立同分布的高斯分布的加权平均会使得新像素的方差相对于原始噪声方差成比例地缩小，相应地我们就可以得到更加纯净的图像。

2 基本降噪思路

根据式 (1.3)，我们可以找到一种消除 AWGN 的有效方法，就是在同一拍摄条件下重复拍摄多幅图像并进行平均，利用噪声在时间上的独立性将其方差缩小到可以接受的范围。这也是一些训练数据集中真实无噪声图像的常用采集方法。但这种方法的局限性也十分明显，即重复拍摄需要保证被拍摄物的始终静止，否则图像在叠加时就会产生拖影。这导致我们无法对那些包含运动物体的图像进行降噪，自然地也不可能处理普遍包含运动物体的视频。

     既然不能方便地通过时域的独立性来消除 AWGN，那我们只能在空域上下手了，因为 AWGN 在空域上也是满足独立同分布假设的。如果我们能够找出所有具有相近真实值

的像素，然后求平均，就可以同样得到一个新的像素值，其方差相比于原始像素成倍地缩小。但是，由于我们只能得到有噪声图像

，要寻找具有相近真实值

的像素还是十分困难的。从随机变量的角度来讲，对于有噪声图像

中任意两个像素，记其误差（差值平方）为

  （2.1）

为了简便，有时会省略坐标

 。可求得其数学期望和方差分别为

  （2.2）

因此，即便两个像素在真实无噪声图像

中的值是一样的，在有噪声图像

中也会有很大的差别，特别是当噪声的方差很大，例如低照度拍摄时。如果直接在有噪声图像

上选择像素值相近的像素进行叠加，很有可能破坏掉图像的纹理特征，影响降噪的质量。为此，一些解决方法也相继被提出和改进。

2.1 局部滤波

因为图像的内容大多数都是低频的，相邻像素之间差值通常不会很大，我们可以认为当前像素的一定邻域内的所有像素都是相似的，那么把这些像素直接加起来求平均就可以了。然而问题在于，相似不等于相同，把所有像素都认为同等重要会导致图像变得非常模糊，丢失掉很多细节。所以一个直观的解决办法是使用不等权重的窗函数，常用的就是高斯窗，从而赋予中心像素更高的权重，离得越远的像素其重要性衰减得越多，这样就可以更好地保留一些图像的细节。除此以外，前面所说的窗函数的权重分配通常只是考虑空间上的像素距离，这样在跨越一些边缘的时候会使得边缘变得模糊。为了解决这种问题，一些滤波算法，比如双边滤波，在计算权重的时候不仅考虑像素之间的空间距离，还会考虑值域上的距离，即如果两个像素的值相差很多，即便两者空间距离比较小，也只会获得较小的权重，从而更好地保护边缘。由于局部滤波需要考虑图像纹理的变化，一般滤波窗口不能设置得太大，从而限制了参与滤波的像素个数，在噪声方差较大时对噪声的衰减并不是很理想。

2.2 非局部滤波

   既然在有噪声图像

上直接比较两个像素的相似性会有较大的误差，那如果是比较两个块的相似性呢？由于噪声

在空域上也是满足独立同分布假设的，如果有噪声图像

上的两个块没有重叠，也就是每个像素都是独立的，那么它们的均方误差的数学期望和方差分别为

  （2.3）

  （2.4）

所以，两个块的误差的方差也是随着块的增大而比例地缩小的，如果我们使用较大的块去比较，它们的相似性就可以得到比较好的保证。因为图像中总是有很多相似的区域，比如天空、道路、建筑砖块等等，如果不考虑计算的复杂度，我们通常可以找到足够多的相似块。将这些相似块叠加起来，噪声的方差就能够衰减得足够小。当然，每个块的匹配误差是不一样的，我们可以赋予那些误差较小的块较大的权重，而误差较大的块较小的权重，至于那些误差超过一定阈值的块则直接舍弃。

    不过，由于不同的块还是具有一定的差别的，直接叠加可能会使得图像变得十分杂乱，所以在非局部均值（Non-Local Means）算法中，只有每个块的中心像素才会被用于滤波。具体地，对于图像中的每个像素，以其一定大小的邻域作为参考块，在图像中的其他位置找到足够多的相似块，那么，该像素的滤波结果就是这些相似块（包括参考块）中心像素的加权平均，即

  (2.5)

  (2.6)

其中

为归一化系数，即所有未归一化权值的和，

代表以

为中心点的在有噪声图像

上的邻域块。

用于控制滤波的强度，

越大，权值随着块误差增大的衰减速度越慢，也就是接近于等权重叠加；

越小，权值衰减速度越快，也就是更倾向于只叠加那些误差很小的块。通过这种方法，我们一方面可以获得更低的噪声方差，同时也不会对图像的纹理造成太大的破坏，因为这些叠加像素的相似性是可以得到比较好的保证的。当然，这只是理想的情况，因为均方误差较小的两个块的中心像素其实也有可能差别很大。另外，这种方法的缺点也很明显，因为我们需要为图像中的每个像素寻找相似块，其时间复杂度无疑是非常高的，处理一张低分辨率的图片也可能需要十几秒之久。

2.3 变换域稀疏表示

图像的变换域一般是相对于空间域来说的，即像素值在水平或垂直方向的波动快慢。很明显，图像的内容通常以低频为主，例如随处可见的大片纯色或者渐变等等，只有少数的细节部分才会包含一些较高频的信息，而特别高频的如密集的马赛克等等反而会引起观察者生理和心理上的不适。所以，图像变换到频域之后通常只有少数几个低频点会有较大的系数或者能量，而绝大部分高频点上的系数或能量则几乎可以忽略，这就是图像在变换域中的稀疏性。

     图 1 展示了一幅自然无噪声图像（来自 SIDD 图像降噪数据库，下同）的 DCT-II 变换结果，右侧为相应的 DCT 系数的绝对值，为了方便作图已经把最低频的 4 × 4的系数置零，因为这部分的系数本身很大，例如 (0,0) 这个点实际上就是图像所有像素值的和除以边长，如果把它们放在图中其他值几乎可以忽略。可以看到，整个变换结果中只有低频区域很小一部分有较大的值，也就是说低频分量占据了图像的绝大部分能量。

图 1无噪声图像与其DCT变换系数

接下来，我们再看看添加了 AWGN 的图像的 DCT 变换系数是什么样子的。图 2 左侧为图 1 无噪声图像添加了标准差为 σ = 20 的高斯白噪声的结果，右侧为相应的 DCT 变换系数。可以看到，虽然低频分量还是占据了绝大部分的能量，但由于噪声的缘故，其他的频率点的系数已经远远没有图 1 那样平滑。甚至我们可以发现，即便是最高频的地方的变换系数的量级竟然和那些接近低频的地方是差不多的，或者说，除了最低频的地方，其他的变换系数就像是一些独立同分布的噪声。实际上，不仅是 DCT，可以证明 AWGN 经过任意次数的正交变换之后仍然是 AWGN ，并且其分布和原始噪声的分布一致，即具有相同的方差。

图 2添加了AWGN的图像与其DCT变换系数

正交变换在图像变换中十分普遍，例如 DCT, DST 以及 Haar 变换等等。对于正交变换，我们可以使用一个正交矩阵

来表示其变换核，并且有

  (2.7)

其中

分别代表有噪声信号

，无噪声信号

和噪声信号

的一个向量（例如某一列/行），

为相应的正交变换系数。注意，虽然这里只表示了一维的变换，但对于可分离多维正交变换，其相当于先对某一维（例如每一列）进行正交变换，然后对所得的结果进行其他维（例如每一行）的正交变换，比如对于二维的变换，有

  （2.8）

所以，对于可分离的多维正交变换，有噪声信号

的变换系数总可以表示为无噪声信号

的变换系数与噪声信号

的变换系数的和。而前面已经提到，AWGN 经过多次正交之后仍然是 AWGN，所以有噪声信号

的变换系数其实就是无噪声信号

的变换系数加上与原始噪声同样方差的 AWGN 的结果。

因此，对于方差较小的噪声，因为无噪声图像本身的变换系数主要集中在低频区域，并且具有较大的幅度，而噪声在变换域中由于方差较小的缘故，通常只有较小的值（例如 95.4% 的值都在两个标准差的范围内），我们可以通过设定阈值的方法（一般为 2~3 个噪声标准差），只保留有噪声图像

的变换系数中那些具有较大幅度的值，其他的则置零。这样，绝大部分的噪声的变换系数就得以消除，只有少数还残留在低频的系数中，从而经过反变换后能够极大地减少噪声，同时还能较好地保留图像的纹理。当然，那些物体边缘的区域可能会由于高频信息的丢失而稍微显得模糊，或者有少量的振铃效应等等。

图 3真实强噪声图像与其DCT变换系数

然而，当噪声方差本身很大时，噪声的能量和无噪声图像本身的能量级别相当，那么在变换域中无噪声图像部分次低频的变换系数反而有可能被噪声给掩盖掉。例如，图 3 是一幅在低照度环境下真实拍摄的图像，其实际上是用于叠加生成图 1 中的无噪声图像的上百张副本之一。可以看到，在空间域中，图像本身的信息已经极大地被噪声破坏，而在变换域中，我们也可以看到与图 2 人工噪声非常相似的杂乱的噪声，这也在一定程度上表明图像在低照度环境下的噪声确实以 AWGN 为主。另外，相比于图 2，图 3 变换域中的噪声幅度明显要大得多，普遍可以达到 80-100，所以其噪声的标准差应该大概在 40-50 左右。回顾图 1 无噪声图像的变换系数，其实有很大一部分的低频系数的幅度都在 100 以内，如果这时我们直接通过设置阈值的方法去除那些幅度低于 2~3 标准差的变换系数，就会把很多本身有用的信息给丢弃掉，最终我们得到的可能只是一幅过度平滑的图像，没有任何的细节，这无疑是不可接受的。况且，因为 AWGN 在变换域中仍然是 AWGN，所以不仅是高频的系数，我们所保留的低频系数本身也已经是被噪声污染了的。例如在图 3 中，最大的系数幅度（注意这里不包含最低频的 4 × 4系数）只在 300-400 之间，而在图 1 中其应该在 400-500 的范围，所以即便我们只保留低频的信息还是会有比较大的失真的。

 为此，类似于非局部均值（NLM）算法，一种被称之为协同滤波（Collaborative Filtering）的方法被提出。其原理也比较简单，即把那些相似的块叠加起来形成第 3 维，在前两维的变换基础上，再进行第 3 维的正交变换，之后就可以通过设置阈值的方法把那些低于一定幅度的变换系数去掉。那么，这种方法相比于前面的二维变换又有怎样的改进呢？其实，其基本的依据还是在于 AWGN 经过正交变换之后仍然是与原始噪声同分布的 AWGN，也就是说，无论我们进行多少次正交变换，噪声的强度都是不会发生改变的。然而，因为我们要处理的是相似的块，所以这些块本身无噪声的部分在第 3 维上是几乎一样的，也就是低频分量占据了主要部分，从而经过第 3 维的变换之后我们可以使得无噪声信号的能量更加集中，具体表现在属于无噪声信号的那部分系数幅度将变得更加的大。而因为噪声信号经过第 3 维的变换之后幅度没有变化，所以两者的差距将进一步地拉大，这时进行阈值置零操作就能够减少误杀的情况，而在那些没有被置零的系数中噪声信号的能量比例也进一步地缩小，从而更好地保留图像本身的高频信息。

图 4不同数量相似块对变换系数幅度的影响

为了便于理解，图 4 展示了在不同数量的相似块下有噪声图像的变换系数幅度的变化情况。注意，为了方便作图，这里只使用了图 1 中无噪声图像的某一列，其变换系数如左上角所示，并且直流系数已经被置零。为了模拟不同的相似块，这里只是简单重复地为该列加上标准差为 σ = 50的高斯白噪声，然后把它们按列组合成一张二维图像，这时进行行方向的变换也就相当于我们前面所说的协同变换。因为经过行方向的协同变换后，变换系数的能量主要集中到第 1 列，所以图 4 只给出第 1 列的变换系数的绝对值，并且同样把直流分量置零。

可以看到，当只有单纯一个添加了噪声的块时，如图 4 右上角所示，由于噪声的方差很大，甚至达到了无噪声图像本身的量级，经过变换之后其系数和原始的相比除了直流分量已经完全无法分辨了，自然也不可能保留太多的细节。而随着相似块数量的不断增加，如图 4 下半部分，噪声的变换系数的量级并没有发生改变，始终保持在 100 以下，而其他原本属于无噪声图像的变换系数则随着能量的集中而不断地增大，从而与噪声的变换系数分离开来，这样我们就可以方便地抑制属于噪声信号的那部分能量，而不会过多地影响原本属于无噪声信号的那一部分，更好地保留图像的细节。

这时，你可能会想，如果我们直接把这些相似块的变换系数加起来求平均不是更好吗？其实，对于线性变换来说，变换域的相加和原始空间域的相加是等价的。我们已经说过，由于不同的相似块并不是完全一样的，所以直接把它们叠加起来只会让图像更加糟糕。图 4 只是展示了一种理想的情况，即所有相似块对应的无噪声图像都是一样的，所以左下角和右下角的变换系数才会和左上角那么相似。另外，这里也存在一个矛盾，即越多的相似块能够让低频分量的能量更加集中，但也会导致块与块之间的差异点越来越多，从而引入一部分的较高频的信息，但这部分系数的能量通常又不足以与噪声系数的能量区分开来，所以这部分信息也很有可能和噪声一起被抹去，经过反变换之后丢失更多的细节。

相对于 NLM 算法来说，协同滤波同样需要寻找足够的相似块。但是，NLM 是以每个像素点为操作单位的，也就是需要为每个像素找到足够多的相似块，然后把这些相似块的中心像素进行加权叠加，从而把该像素的噪声方差降低到足够小。所以，虽然其时间复杂度很高，但能够比较好地保留图像的细节。相反，协同滤波可以一次性处理所有块的所有像素，这样速度无疑要快很多，但是在细节保留方面则要逊色不少，因为该算法本身并没有把噪声的方差缩小，而无论我们怎样把低频分量的能量集中，总会有一些细节部分的的能量会被淹没在噪声当中，无法恢复。

另外，为了寻找相似块，我们总不能把整幅图像作为一个参考块，所以协同滤波本身也只是处理了图像的一小部分区域。为了保证图像中所有像素都能被覆盖到，最保险的方法是将图像拆分为规整的矩形块，然后使用栅格扫描的方式，对每一个区域都进行一次相似块匹配以及协同滤波。但这时又有一个问题，因为每个区域的相似块都是在整幅图像上寻找的，所以不同区域的相似块总可能包含同样的像素，也就是同一个像素通常会被进行多次的协同滤波。所以，一个很直观的改进方法就是把这些不同的协同滤波的结果进行整合，从而对当前像素的无噪声值做出更恰当的估计，相应地在使得图像更加纯净的同时能够更好地保留自身的细节。这就是下面我们所要说的 BM3D 降噪算法的基本思路。实际上，为了达到更好的降噪效果，通常来说这些拆分的区域本身也应该是有重叠的，这样可以进一步提高某个像素被多次滤波的可能，这种方法称之为过完备（Overcompleteness），并且被包括 BM3D 等多种算法所采用。

3 BM3D算法原理

在第 2 节中，我们其实已经对 BM3D 算法的基本理论依据做了详细的描述，其主要可分为三个步骤，即首先对每个参考块进行相似块匹配（Block-Matching）并分别得到一个三维的组合，然后对其进行协同变换和滤波（3D-Transform），最后对各个参考块对应组合的滤波结果进行整合（Aggregation），从而得到最终的降噪结果。在此之上，为了进一步改善图像的质量，BM3D 实际进行了两次降噪，即将以上三个步骤再重复了一遍，但是具体的块匹配标准、滤波方式以及整合权重等等会有一些区别，其基本流程如图 5 所示。

在 Step1 中，协同滤波主要通过设置硬阈值（Hard-thresholding）的方法来去除噪声的能量，这也就是我们在第 2 节中所介绍的方法。这种方法相对来说比较简单，但是一方面可能误杀掉部分细节的信息，另一方面也无法去除噪声在低频系数中的能量。另外，在最初的有噪声图像中，块匹配的误差其实是一个随机变量，根据式 (2.4) 可知其方差大概正比于噪声方差的平方，当噪声方差较大时块匹配的质量并不是很理想，也就可能引入更多块与块之间不同的细节信息，并且被硬阈值滤波给误杀掉。所以，Step1 通常只是作为一种基础估计（Basic Estimate），但由于过完备以及块整合的缘故，其结果也是也是可以接受的。

 作为改善手段，Step2 使用了 Step1 中基础估计的结果来进行相似块的匹配，因为噪声已经极大地被消除了，所以我们认为这时的匹配结果是更加准确的。注意，这里并不是直接复用 Step1 中已经找到的相似块的位置，而是利用基础估计图像重新搜寻一遍，所以 Step2 和 Step1 所需的时间其实是差不多的。但是，不同于 Step1 使用硬阈值的方法来进行协同滤波，对基础估计图像再进行一次硬阈值滤波其实是没有多大用处的，因为这样只会让得到的图像更加模糊。当然，这时我们可能会想使用在基础估计图像上找到的相似块位置对原始有噪声图像进行硬阈值滤波，但这种方法始终不能避免前面 Step1 所说的缺点，即误杀细节，以及无法去除低频系数的噪声。实际上，Step2 采用了经验维纳收缩（Emprical Wiener Shrinkage）的方法，其根据基础估计图像上的协同变换系数的功率谱以及噪声的强度，对原始有噪声图像同样位置 3D 块的协同变换系数进行收缩，然后对收缩后的系数进行反变换以及整合，直接得到最终所要的降噪结果。所以在这个过程中，我们实际处理的还是输入的原始有噪声图像，而基础估计只是作为一种辅助，用于更准确地确定原始图像中相似块的位置，以及相应的系数收缩比例。具体的原理将会在后面进行介绍。

       为了方便阅读原论文，这里的所用到符号会与其保持一致。对于一幅添加了 AWGN 的图像，其可表示为

  (3.1)

注意，这里所说的图像都是 2D 单通道图像，对于彩色图像的降噪将在后面进行介绍。我们用

​来代表有噪声图像

上的一个左上角坐标为

的大小为

的正方形块，而

为一个 3D 块，所包含的每个 2D 块的左上角坐标由集合

来定义。为了表示每个 3D 块中的参考块，我们会在相应的坐标下添加

 (Reference) 的标识，即

。同样，对于无噪声图像

上的一个块，我们也可以使用

来表示。因为 Step1 和 Step2 的步骤基本一致，所以它们所用到的一些参数的意义也是一样的，为了区分两者，我们会分别加上 “

” 和 “

” 的标识。对于 Step1 的基础估计图像，我们用

来表示，Step2 所得的最终估计图像则为

。

3.1 Step1

接下来我们将对 BM3D 算法的 Step1 进行详细的介绍。首先有必要再次说明参考块的具体意义，不同于 NLM 算法会把每个像素的邻域都作为一个参考块，为了降低计算的复杂度，一般在算法实现中，我们会从图像的第一个

块开始，采用栅格扫描的方式，每次平移几列（行），例如用

来表示，一般为 3~6，直至铺满整幅图像。因为我们必须保证所有的像素都至少参与一次协同滤波，所以所有的块都会被选择一次作为参考块，而参考块的重叠也是过完备（Overcompleteness）的一种体现。在第 2 节中我们已经提到，对于无噪声图像

上某个大小为

的参考块

和其他位置同样大小的块

，记其匹配误差或者块距离为

  (3.2)

如果这两个块是没有重叠的，也就是每个像素都是独立的，那么在有噪声图像

上对应位置的两个块的的匹配误差的数学期望以及方差分别为

  (3.3)

  (3.4)

所以，当噪声方差很大或者块比较小的时候，直接在有噪声图像

上寻找相似块是很不准确的，有可能无噪声图像

上两个差距很大的块在有噪声图像

上由于方差很大的原因差距显得比较小，从而被判定为相似，造成错误的匹配。

    为了解决这个问题，论文提出先对有噪声图像

上的两个块进行可分离的二维正交变换，然后将那些幅度小于一定阈值的系数置零，这样两个块的匹配误差就可表示为这些系数的均方误差，即

  (3.5)

其中，

为硬阈值操作，相应的阈值为

，一般为两个标准差。

为相应的可分离二维正交变换。我们前面已经有提到，AWGN 经过正交变换之后还是 AWGN，而无噪声图像由于其在变换域中的稀疏特性，一般只有在少数低频的位置上有较大的值。所以经过二维正交变换之后，通过设置硬阈值的方法，我们可以去除绝大部分噪声的能量。虽然也有一部分图像本身的能量被去除，但其影响基本可以忽略，因此这时我们可以使得块匹配更加准确。注意，我们不需要对系数进行反变换再在空间域上计算两个块的匹配误差，因为对于正交变换以及其反变换来说，变换域和空间域的距离其实是一样的。

     实际上，并不是所有的情况都需要先对变换系数进行硬阈值操作再计算块距离，论文中也指出一般在噪声的标准差

>40 时才需要设置一个非零的阈值，对于阈值为零的情况，式 (3.5) 所得的块距离与直接在空间域计算是等价的。另外，我们可以先计算所有可能的候选块的 2D 变换系数并缓存起来，因为相邻的参考块的需要搜索的候选块大多数是一样的，所以这部分结果也可以进行复用。当然，我们不需要把整幅图像的候选块的变换都存储起来，因为全搜索不可能应用于整幅图像，而是只有某个较小的搜索窗口。当搜索窗口移出某个候选块时，我们就可以删除其变换系数以节省存储空间了。另一方面，因为我们已经计算得到所有相似块的 2D 变换系数，如果协同变换采用的是可分离的 3D 正交变换，那么这些 2D 变换的结果就可以进行复用，我们只需再进行第 3 维的 1D 变换即可，从而进一步降低其复杂度。

       回到 Step1 的块匹配，根据式 (3.5) 计算得到当前参考块与其他所有块（实际上只有某个邻域内的块会被考虑）的匹配误差后，我们只保留那些误差小于一定阈值的块，并得到相应的坐标集合

  (3.6)

其中

为判定两个块相似的最大误差阈值，一般根据经验来设置。很明显当前参考块肯定会被判定为相似块，因为其匹配误差为 0，所以坐标集合中至少有 1 个元素。这时我们把所有的相似块堆叠起来，就可以得到一个形状为

的 3D 数组

，其中

代表该集合元素的个数。注意，通常来说该 3D 数组中块的顺序并不是很重要，但有时候因为第 3 维的正交变换可能需要输入的长度为 2 的幂，例如 Haar 变换，这时我们只保留集合中那些匹配误差最小的块以满足长度的要求。

    得到参考块对应的 3D 数组后，我们就可以进行 3D 的协同变换和滤波，形式化地可表示为

  （3.7）

其中

为相应的 3D 正交变换，一般来说为了降低计算的复杂度我们会选用可分离的变换，因为前面在计算块匹配误差时二维的变换结果已经存在了，我们完全可以对其进行复用。所以，

可以表示为前两维的

和第 3 维的

的组合。具体的变换类型没有太多的限制，一般前两维的变换会选用傅里叶变换类型的如 DCT 等，可以更好地捕捉图像块中周期性的信息，而第 3 维的变换则会选用小波变换类型的如 Haar 等，可以通过尺度缩放等更好地捕捉块与块之间的局部相似性，但实验数据也表明其性能差别并不是很大。

为另一个硬阈值操作，相应的阈值为

，通常为 2～3 个标准差。

包含了

个块

，表示经过式 (3.7) 的协同滤波所得的结果，对应的参考块的坐标为

。注意，因为不同的参考块可能包含同样的相似块，所以一个块会对应多次的协同滤波结果，我们有必要通过其对应的参考块坐标进行区分，这些不同的协同滤波结果将会在后面的操作中进行整合，从而为每个像素获得 Step1 的基础估计。

       得到每个参考块以及其对应的3D组的协同滤波结果后，因为参考块本身就是有重叠的，而且不同的参考块也可能包含重叠的相似块，为了得到每个像素的基础估计，我们有必要对这些来自不同组合的协同滤波结果进行整合，例如赋予这些不同的组合以不同的权值。为了获得一个比较好的结果，我们当然是希望每个组合的协同滤波都能够完美地去除噪声，但这明显是不可能的，因为经过多次的正交变换以后，AWGN 还是会以 AWGN 的形式残留在那些没有被去除的低频系数中，没有被置零的系数越多，残留的噪声能量就越多。另外，如果一个组合内的所有块都足够相似，那么其变换系数也本应该是非常稀疏的。所以，我们可以根据协同滤波后剩余的非零系数的个数来进行权值的分配，假设该 3D 组合内所有的像素都是独立的，记硬阈值操作后剩下的系数个数为

，则有

  （3.8）

实际上，一个 3D 组合内的相似块通常都会有重叠，也就是说所有像素之间并不会完全独立，但为了简单起见我们还是会以式 (3.8) 来进行权值的计算。注意，这里的权值对于该 3D 组合内所有的像素都是一样的，但为了减少边界效应，通常都会为每个块加上一个窗函数，从而赋予中心像素更高的权值。论文中使用了 Kaiser 窗，其一维定义为

  (3.9)

其中

为零阶修正的贝塞尔函数，一般在 MATLAB 或者 Python.numpy 中可以直接调用。

为相对于中心点的距离，所以对于偶数长度的离散窗其应该为 1/2 的奇数倍。二维的 Kaiser 窗只需把列向量的一维 Kaiser 窗乘以其转置即可，在 BM3D 中其大小应该为

。注意这个窗是对单个块而言的，同一个 3D 组的各个相似块使用同样的窗函数。因为同一个像素可能包含在不同的相似块中，通过添加 Kaiser 窗也可以使得属于不同块不同位置的同一个像素获得不同的权值。最后，通过对某个像素所属的所有参考块以及对应的相似块的加权叠加，我们就可以得到其基础估计，即

其中我们假设每个块都已经通过补零的方法填充到和原图一样的大小，并且通过

来判断某个像素

是否在块

上。从式中也可以看出，一个像素既可以属于不同的参考块，也可以属于同一个参考块的不同相似块，而同一个块也可以是不同的参考块的相似块。通过这种过完备的协同滤波方法，我们一方面可以尽可能地去除噪声的能量，同时也能较好地保持图像中的细节，并且可以保持较低的计算复杂度，因为这里主要是以块为操作单位的，相比于 NLM 算法可以减少很大一部分寻找相似块的时间。

        以上是 BM3D 论文中对 Step1 流程比较公式化的描述，对于部分刚接触 BM3D 算法的读者来说可能不那么直观，特别是聚集，即式 (3.10)，所表达的意思一时不能理解。为此，我这里也重新画了一个流程图，如图 6 所示，这个流程图对于 Step2 也是适用的，只是具体的实现细节会有所差别。图 6 中最大的矩形框内部展示了图像中某个参考块（即橙色块）的处理流程，在 BM3D 算法开始之前，需要首先把图像划分为一系列的参考块，具体的参考块划分方式并不固定，即参考块可以具有不同的大小和形状，只需保证所有参考块所在图像区域的并集能够完整覆盖图像的所有像素即可，在实际算法实现中为了方便通常会固定参考块的大小，并通过固定步长的光栅扫描顺序来实现图像所有像素的遍历覆盖。

对于当前所遍历到的参考块（即橙色块），我们会在其周围一定区域内（即搜索窗口，图中带斜线阴影的虚线框）进行相似块搜索，假设我们所搜索到的最相似的块如图中灰色块所示，具体相似块的个数由算法配置决定，这里加上参考块本身（匹配误差为 0）总共可得到 4 个块。注意，尽管图中所示的几个相似块是不重叠的，但在实际的运行环境中，相似块之间存在重叠的情况非常常见。我们把这 4 个相似块堆叠起来，即可得到一个三维的相似块组合，之后对其进行三维的空域-频域转换得到相应的三维的变换系数，并对变换系数进行硬阈值滤波，对滤波后的变换系数进行三维的频域-空域转换，即可得到滤波后的三维相似块组合，即 4 个滤波后的二维相似图像块。

    理论上，我们把这 4 个滤波后的图像块放回图像原来的位置，即用滤波后的图像块的像素值替换掉输入图像对应位置图像块的像素值，就相当于对图像的这 4 个区域进行了滤波降噪操作。因为我们前面强调了图像在划分参考块时需要保证所有参考块的并集需要覆盖图像的所有像素，所以当我们遍历所有的参考块后，必然能够保证图像的每个像素都能进行至少一次滤波。但问题也出在这里，即图像同一个像素位置可能会进行多次的滤波，从而产生多个不同的滤波后的像素值。这是因为参考块之间可能会存在重叠，而参考块所搜索到的相似块之间也会存在重叠，使得图像中同一个像素会被多个图像块所包含。更糟糕的是，因为相似块匹配本身几乎是完全随机的，在完成图像所有参考块的遍历之前，我们不可能预知某个像素会被多少个相似块所包含，从而产生多少个不同的滤波后的像素值。为此，最简单的做法是，分别创建一个分子（Numerator）缓冲区和一个分母（Denominator）缓冲区，如图 6 所示，两者的大小和输入图像一样。分子缓冲区负责记录这些不同的滤波后的像素值的总和，而分母缓冲区则负责记录总共有多少个不同的滤波后的像素值。当图像中所有的参考块都完成遍历后，把分子缓冲区的不同滤波后的像素值的累加值除以分母缓冲区所记录的滤波次数，即可得到某个像素的不同滤波结果的平均值，该平均值即可作为最终的输出图像的像素值。但正如式 (3.10) 所示，在滤波的过程中，我们还能给某个像素每次滤波所得的像素值赋予不同的权重，这时我们只需让分子缓冲区负责记录这些不同的滤波后的像素值的加权和，而分母缓冲区则负责记录不同的滤波结果所对应的权重和，在图像所有参考块遍历结束后将分子缓冲区的像素值加权和除以分母缓冲区的权重和即可得到相应的输出像素值。

3.2 Step2

3.3 高效的实现方法

3.4 算法的一些扩展