传统图像降噪算法之BM3D原理详解

        图像降噪是一个十分具有实用价值的研究方向，因为噪声总是无处不在的。当处于比较昏暗的环境时，噪声将极大地影响着我们所拍摄的图像。如今，随着深度学习算法以及相关硬件的不断发展，深度卷积网络同样在图像降噪领域占据了主流，并且代表了该领域最优异的成绩。但是，深度神经网络同样有着其缺点，例如模型过于庞大而计算复杂度过高，以及缺乏一些理论上的解释性，当然这些缺点正不断地得到弥补。为了更好地理解图像降噪的基本原理，我们有必要回过头来仔细研读一些传统算法的具体思路，了解其所使用基本理论依据，以及一些巧妙的改进方法。在这些传统降噪算法中，最经典而强大的莫过于 BM3D 了。这篇文章将全面地对其原理进行解读，并且对其论文中一些没有提及的细节进行补充，让各位读者能够更加轻松地理解其算法的内核。在开始这篇文章之前，本人建议大家可以先看一下以下的文章，主要是对本文中一些需要用到但是为了节省篇幅而没有细讲的基本原理进行补充：
1.https://blog.csdn.net/qq_33552519/article/details/108236388
2.https://blog.csdn.net/qq_33552519/article/details/108372176
3.https://blog.csdn.net/qq_33552519/article/details/108589799

1 加性高斯白噪声

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

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

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

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

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

       综上所述，记真实无噪声图像为 u，AWGN 噪声为 n，有噪声图像为 v，那么三者的关系可表示为

        其中 x 代表像素点位置。标准差 σ 反映了噪声值的水平，根据其概率密度函数，99.6% 的噪声值都会在 3 σ 的范围之内。在这里，噪声 n 为随机变量，而真实图像 u 为常量，所以有

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

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

2 基本降噪思路

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

       既然不能方便地通过时域的独立性来消除 AWGN，那我们只能在空域上下手了，因为 AWGN 在空域上也是满足独立同分布假设的。如果我们能够找出所有具有相近真实值 u 的像素，然后求平均，就可以同样得到一个新的像素值，其方差相比于原始像素成倍地缩小。但是，由于我们只能得到有噪声图像 v ，要寻找具有相近真实值u的像素还是十分困难的。从随机变量的角度来讲，对于有噪声图像 v 中任意两个像素，记其误差（差值平方）为

为了简便，有时会省略坐标 x xx。可求得其数学期望和方差（见链接1）分别为

        因此，即便两个像素在真实无噪声图像 u 中的值是一样的，在有噪声图像 v 中也会有很大的差别，特别是当噪声的方差很大，例如低照度拍摄时。如果直接在有噪声图像 v 上选择像素值相近的像素进行叠加，很有可能破坏掉图像的纹理特征，影响降噪的质量。为此，一些解决方法也相继被提出和改进。

2.1 局部滤波

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

2.2 非局部滤波

        既然在有噪声图像 v 上直接比较两个像素的相似性会有较大的误差，那如果是比较两个块的相似性呢？由于噪声 n 在空域上也是满足独立同分布假设的，如果有噪声图像 v 上的两个块没有重叠，也就是每个像素都是独立的，那么它们的均方误差的数学期望和方差分别为

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

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

        其中 Z(x) 为归一化系数，即所有未归一化权值的和，Bv(x)代表以 x 为中心点的在有噪声图像 v 上的邻域块。h 用于控制滤波的强度，h 越大，权值随着块误差增大的衰减速度越慢，也就是接近于等权重叠加；h 越小，权值衰减速度越快，也就是更倾向于只叠加那些误差很小的块。通过这种方法，我们一方面可以获得更低的噪声方差，同时也不会对图像的纹理造成太大的破坏，因为这些叠加像素的相似性是可以得到比较好的保证的。当然，这只是理想的情况，因为均方误差较小的两个块的中心像素其实也有可能差别很大。另外，这种方法的缺点也很明显，因为我们需要为图像中的每个像素寻找相似块，其时间复杂度无疑是非常高的，处理一张低分辨率的图片也可能需要十几秒之久。

2.3 变换域稀疏表示

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

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

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

        正交变换在图像变换中十分普遍，例如 DCT, DST 以及 Haar 变换等等。对于正交变换，我们可以使用一个正交矩阵 A 来表示其变换核，并且有

        其中 v , u , n  分别代表有噪声信号 v ，无噪声信号 u 和噪声信号 n 的一个向量（例如某一列/行），Tv, Tu, Tn为相应的正交变换系数。注意，虽然这里只表示了一维的变换，但对于可分离多维正交变换，其相当于先对某一维（例如每一列）进行正交变换，然后对所得的结果进行其他维（例如每一行）的正交变换，比如对于二维的变换，有

        所以，对于可分离的多维正交变换，有噪声信号 v 的变换系数总可以表示为无噪声信号 u 的变换系数与噪声信号 n 的变换系数的和。而前面已经提到，AWGN 经过多次正交之后仍然是 AWGN，所以有噪声信号 v 的变换系数其实就是无噪声信号 u 的变换系数加上与原始噪声同样方差的 AWGN 的结果。

       因此，对于方差较小的噪声，因为无噪声图像本身的变换系数主要集中在低频区域，并且具有较大的幅度，而噪声在变换域中由于方差较小的缘故，通常只有较小的值（例如 95.4% 的值都在两个标准差的范围内），我们可以通过设定阈值的方法（一般为 2~3 个噪声标准差），只保留有噪声图像v的变换系数中那些具有较大幅度的值，其他的则置零。这样，绝大部分的噪声的变换系数就得以消除，只有少数还残留在低频的系数中，从而经过反变换后能够极大地减少噪声，同时还能较好地保留图像的纹理。当然，那些物体边缘的区域可能会由于高频信息的丢失而稍微显得模糊，或者有少量的振铃效应等等。

        然而，当噪声方差本身很大时，噪声的能量和无噪声图像本身的能量级别相当，那么在变换域中无噪声图像部分次低频的变换系数反而有可能被噪声给掩盖掉。例如，图 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 展示了在不同数量的相似块下有噪声图像的变换系数幅度的变化情况。注意，为了方便作图，这里只使用了图 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 块的协同变换系数进行收缩，然后对收缩后的系数进行反变换以及整合，直接得到最终所要的降噪结果。所以在这个过程中，我们实际处理的还是输入的原始有噪声图像，而基础估计只是作为一种辅助，用于更准确地确定原始图像中相似块的位置，以及相应的系数收缩比例。

        为了便于理解，我这里重新画了一个流程图，如图 6 所示：

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

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

        其实，经过 Step1 后，我们得到的基础估计图像已经能够极大地去除噪声的能量。虽然前面也说了，简单的硬阈值操作并不能去除残留在低频系数上噪声，并且很多的细节信息也有可能被一同抹杀掉，不过经过反变换之后这部分噪声的能量会被分散稀释到各个像素中，另外通过整合不同的协同滤波的结果也可以在一定程度上弥补这两个内在的缺陷。这时我们可能会想，一些细节信息会被抹杀掉，其实一部分原因是 3D 组合内的块不够相似，从而该细节部分的能量不够集中而无法与噪声区分开来，而块与块不够相似的一部分原因是输入图像中的噪声过于强大，无法对块匹配误差进行准确的估计。因为通过 Step1 我们已经能够得到一张比较纯净的图像，那作为一种改善手段是不是可以通过在基础估计图像上寻找相似块，然后重复 Step1 的操作呢？这似乎可以形成一种迭代的过程。但实际上，这种方法的性能提升并不是很明显，因为我们始终无法避免硬阈值滤波的缺点。为此，不同于硬阈值滤波直接将小于阈值的系数置零而完全保留那些较大的系数，BM3D 引入了一种称为经验维纳收缩（Emprical Wiener Shrinkage）的方法，其实际上是根据基础估计图像上的变换系数的功率谱来对原始有噪声图像的变换系数进行收缩，从而能够更好地去除包括低频系数中的噪声能量，同时尽可能地保留原本属于无噪声图像的那部分高频系数。

4 算法的一些扩展

        BM3D 原本主要用于 2D 灰度图像的降噪，但我们接触到绝大多数都是 RGB 彩色图像，所以有必要对其进行扩展。在加性高斯白噪声的定义中，白噪声表明其功率谱密度是服从均匀分布的，也就是不同颜色像素上噪声分布应该是一样的，那么直观的解决办法是分别对 R, G, B 分量进行降噪。但实际上，由于光线的能量是有限的，当我们把一束光线的能量分散到 R, G, B 分量上时，其信噪比相应地就会有所降低，这不利于保证相似块匹配的准确度。并且，R, G, B 分量在空间上具有很强的相关性，它们的纹理细节基本是一致的，但可能对比度有所不同。因此，我们可以通过将其转换到其他一些色彩空间来综合这些信息，例如常用的 YUV 等等。YUV 包括一个亮度分量 Y 和两个色度分量 UV，一般认为 Y 分量综合了 RGB 之间的空间相关性，只保留了图像的亮度信息，包括图像边缘纹理等等都能够很好地表现，同时使得光线能量更加集中，从而具有更高的信噪比；UV 分量分别表示 B, R 分量与 Y 分量之间的差异，从而去掉了大部分的纹理细节，主要反映了图像的色彩信息，其变化一般是比较平滑的，所以上面包含的主要为低频信息，人眼对这些色度的变化相对不会那么敏感，但由于光线能量主要集中在 Y 分量上，所以也具有较低的信噪比。所以，为了达到更好的降噪效果，论文提出在 YUV 空间进行处理，并且相似块的匹配只在 Y 分量上进行，而 UV 分量则直接复用 Y 分量的位置信息。这是因为 YUV 三个分量本身具有相似的轮廓信息，但是 UV 分量具有较低的信噪比，如果单独在 UV 分量上寻找相似块容易造成错误的匹配，而 Y 分量的信噪比相对较高，同时包含较多的细节信息，可以降低匹配的误差。另外，只在 Y 分量上进行相似块匹配也可以节省掉 2/3 的时间。

       除了彩色图像以外，BM3D 还可以进一步扩展到视频的降噪，目前主要有 VBM3D，VBM4D 等等。VBM3D 与 VBM4D 的区别在于，VBM3D 还是使用 2D 参考块，但是相似块可以跨越不同的帧，最终这些相似块还是组成一个 3D 组合；而 VBM4D 则直接使用 3D 参考块，即把多帧图像堆叠作为一幅 3D 图像，所有相似块也是这幅 3D 图像上某一部分，最终我们得到的是一个 4D 组合。在 VBM4D 中，为了解决不同帧之间的运动，一些论文也提出使用光流信息对 3D 块内不同的帧进行补偿。除此以外，还有一些比如自适应的块大小以及变换类型等等的改进，总之，这些扩展的基本原理都是相似的，即寻找相似块，然后进行协同滤波，最后把同一个像素的不同滤波结果整合。这里不再详述。