大家好，又见面了，我是你们的朋友全栈君。
图像降噪是一个十分具有实用价值的研究方向，因为噪声总是无处不在的。当处于比较昏暗的环境时，噪声将极大地影响着我们所拍摄的图像。如今，随着深度学习算法以及相关硬件的不断发展，深度卷积网络同样在图像降噪领域占据了主流，并且代表了该领域最优异的成绩。但是，深度神经网络同样有着其缺点，例如模型过于庞大而计算复杂度过高，以及缺乏一些理论上的解释性，当然这些缺点正不断地得到弥补。为了更好地理解图像降噪的基本原理，我们有必要回过头来仔细研读一些传统算法的具体思路，了解其所使用基本理论依据，以及一些巧妙的改进方法。在这些传统降噪算法中，最经典而强大的莫过于 BM3D 了。这篇文章将全面地对其原理进行解读，并且对其论文中一些没有提及的细节进行补充，让各位读者能够更加轻松地理解其算法的内核。在开始这篇文章之前，本人建议大家可以先看一下以下的文章，主要是对本文中一些需要用到但是为了节省篇幅而没有细讲的基本原理进行补充：
- https://blog.csdn.net/qq_33552519/article/details/108236388
- https://blog.csdn.net/qq_33552519/article/details/108372176
- 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 u u，AWGN 噪声为 n n n，有噪声图像为 v v v，那么三者的关系可表示为
v ( x ) = u ( x ) n ( x ) , n ( x ) ∼ N ( 0 , σ 2 ) . (1.1) v(x) = u(x) n(x),;n(x) sim mathcal{N}(0,;{sigma ^2}).tag{1.1} v(x)=u(x) n(x),n(x)∼N(0,σ2).(1.1)
其中 x x x 代表像素点位置。标准差 σ sigma σ 反映了噪声值的水平，根据其概率密度函数，99.6% 的噪声值都会在 3 σ 3sigma 3σ 的范围之内。在这里，噪声 n n n 为随机变量，而真实图像 u u u 为常量，所以有
v ( x ) ∼ N ( u ( x ) , σ 2 ) . (1.2) v(x) sim mathcal{N}left( {u(x),;{sigma ^2}} right).tag{1.2} v(x)∼N(u(x),σ2).(1.2)
那么，对于有限个相互独立的像素点的线性叠加，有
∑ i a i v ( x i ) ∼ N ( ∑ i a i u ( x i ) , ∑ i a i 2 σ 2 ) . (1.3) sumlimits_i { {a_i}v({x_i})} sim mathcal{N}left( {sumlimits_i { {a_i}u({x_i})} ,;sumlimits_i {a_i^2} {sigma ^2}} right).tag{1.3} i∑aiv(xi)∼N(i∑aiu(xi),i∑ai2σ2).(1.3)
这也是大多数 AWGN 降噪算法的基本理论依据，即多个独立同分布的高斯分布的加权平均会使得新像素的方差相对于原始噪声方差成比例地缩小，相应地我们就可以得到更加纯净的图像。
2 基本降噪思路
根据式 (1.3)，我们可以找到一种消除 AWGN 的有效方法，就是在同一拍摄条件下重复拍摄多幅图像并进行平均，利用噪声在时间上的独立性将其方差缩小到可以接受的范围。这也是一些训练数据集中真实无噪声图像的常用采集方法。但这种方法的局限性也十分明显，即重复拍摄需要保证被拍摄物的始终静止，否则图像在叠加时就会产生拖影。这导致我们无法对那些包含运动物体的图像进行降噪，自然地也不可能处理普遍包含运动物体的视频。
既然不能方便地通过时域的独立性来消除 AWGN，那我们只能在空域上下手了，因为 AWGN 在空域上也是满足独立同分布假设的。如果我们能够找出所有具有相近真实值 u u u 的像素，然后求平均，就可以同样得到一个新的像素值，其方差相比于原始像素成倍地缩小。但是，由于我们只能得到有噪声图像 v v v，要寻找具有相近真实值u的像素还是十分困难的。从随机变量的角度来讲，对于有噪声图像 v v v 中任意两个像素，记其误差（差值平方）为
d v ( x 1 , x 2 ) = ( v ( x 1 ) − v ( x 2 ) ) 2 = [ ( u ( x 1 ) − u ( x 2 ) ) ( n ( x 1 ) − n ( x 2 ) ) ] 2 . (2.1) {d_v}left( { {x_1},;{x_2}} right) = {left( {vleft( { {x_1}} right) – vleft( { {x_2}} right)} right)^2} = {left[ {left( {uleft( { {x_1}} right) – uleft( { {x_2}} right)} right) left( {nleft( { {x_1}} right) – nleft( { {x_2}} right)} right)} right]^2}.tag{2.1} dv(x1,x2)=(v(x1)−v(x2))2=[(u(x1)−u(x2)) (n(x1)−n(x2))]2.(2.1)
为了简便，有时会省略坐标 x x x。可求得其数学期望和方差（见链接1）分别为
E ( d v ) = d u 2 σ 2 , D ( d v ) = 8 σ 2 d u 8 σ 4 . (2.2) Eleft( { {d_v}} right) = {d_u} 2{sigma ^2},quad Dleft( { {d_v}} right) = 8{sigma ^2}{d_u} 8{sigma ^4}.tag{2.2} E(dv)=du 2σ2,D(dv)=8σ2du 8σ4.(2.2)
因此，即便两个像素在真实无噪声图像 u u u 中的值是一样的，在有噪声图像 v v v 中也会有很大的差别，特别是当噪声的方差很大，例如低照度拍摄时。如果直接在有噪声图像 v v v 上选择像素值相近的像素进行叠加，很有可能破坏掉图像的纹理特征，影响降噪的质量。为此，一些解决方法也相继被提出和改进。
2.1 局部滤波
因为图像的内容大多数都是低频的，相邻像素之间差值通常不会很大，我们可以认为当前像素的一定邻域内的所有像素都是相似的，那么把这些像素直接加起来求平均就可以了。然而问题在于，相似不等于相同，把所有像素都认为同等重要会导致图像变得非常模糊，丢失掉很多细节。所以一个直观的解决办法是使用不等权重的窗函数，常用的就是高斯窗，从而赋予中心像素更高的权重，离得越远的像素其重要性衰减得越多，这样就可以更好地保留一些图像的细节。除此以外，前面所说的窗函数的权重分配通常只是考虑空间上的像素距离，这样在跨越一些边缘的时候会使得边缘变得模糊。为了解决这种问题，一些滤波算法，比如双边滤波，在计算权重的时候不仅考虑像素之间的空间距离，还会考虑值域上的距离，即如果两个像素的值相差很多，即便两者空间距离比较小，也只会获得较小的权重，从而更好地保护边缘。由于局部滤波需要考虑图像纹理的变化，一般滤波窗口不能设置得太大，从而限制了参与滤波的像素个数，在噪声方差较大时对噪声的衰减并不是很理想。
2.2 非局部滤波
既然在有噪声图像 v v v 上直接比较两个像素的相似性会有较大的误差，那如果是比较两个块的相似性呢？由于噪声 n n n 在空域上也是满足独立同分布假设的，如果有噪声图像 v v v 上的两个块没有重叠，也就是每个像素都是独立的，那么它们的均方误差的数学期望和方差分别为
E ( 1 N 2 ∑ i d v , i ) = 1 N 2 ∑ i d u , i 2 σ 2 . (2.3) Eleft( {frac{1}{ { {N^2}}}sumlimits_i { {d_{v,i}}} } right) = frac{1}{ { {N^2}}}sumlimits_i { {d_{u,i}}} 2{sigma ^2}.tag{2.3} E(N21i∑dv,i)=N21i∑du,i 2σ2.(2.3)
D ( 1 N 2 ∑ i d v , i ) = 8 σ 2 N 4 ∑ i d u , i 8 σ 4 N 2 . (2.4) Dleft( {frac{1}{ { {N^2}}}sumlimits_i { {d_{v,i}}} } right) = frac{ {8{sigma ^2}}}{ { {N^4}}}sumlimits_i { {d_{u,i}}} frac{ {8{sigma ^4}}}{ { {N^2}}}.tag{2.4} D(N21i∑dv,i)=N48σ2i∑du,i N28σ4.(2.4)
所以，两个块的误差的方差也是随着块的增大而比例地缩小的，如果我们使用较大的块去比较，它们的相似性就可以得到比较好的保证。因为图像中总是有很多相似的区域，比如天空、道路、建筑砖块等等，如果不考虑计算的复杂度，我们通常可以找到足够多的相似块。将这些相似块叠加起来，噪声的方差就能够衰减得足够小。当然，每个块的匹配误差是不一样的，我们可以赋予那些误差较小的块较大的权重，而误差较大的块较小的权重，至于那些误差超过一定阈值的块则直接舍弃。
不过，由于不同的块还是具有一定的差别的，直接叠加可能会使得图像变得十分杂乱，所以在非局部均值（Non-Local Means）算法中，只有每个块的中心像素才会被用于滤波。具体地，对于图像中的每个像素，以其一定大小的邻域作为参考块，在图像中的其他位置找到足够多的相似块，那么，该像素的滤波结果就是这些相似块（包括参考块）中心像素的加权平均，即
u ~ ( x ) = ∑ i w ( x , x i ) v ( x i ) , (2.5) tilde u(x) = sumlimits_i {w(x,;{x_i})v({x_i}),} tag{2.5} u~(x)=i∑w(x,xi)v(xi),(2.5)
w ( x , x i ) = 1 Z ( x ) exp ( − ∥ B v ( x ) − B v ( x i ) ∥ 2 2 h 2 ) . (2.6) w(x,;{x_i}) = frac{1}{ {Z(x)}}exp left( { – frac{ {left| { {B_v}(x) – {B_v}({x_i})} right|_2^2}}{ { {h^2}}}} right).tag{2.6} w(x,xi)=Z(x)1exp(−h2∥Bv(x)−Bv(xi)∥22).(2.6)
其中 Z ( x ) Z(x) Z(x) 为归一化系数，即所有未归一化权值的和， B v ( x ) B_v(x) Bv(x) 代表以 x x x 为中心点的在有噪声图像 v v v 上的邻域块。 h h h 用于控制滤波的强度， h h h 越大，权值随着块误差增大的衰减速度越慢，也就是接近于等权重叠加； h h h 越小，权值衰减速度越快，也就是更倾向于只叠加那些误差很小的块。通过这种方法，我们一方面可以获得更低的噪声方差，同时也不会对图像的纹理造成太大的破坏，因为这些叠加像素的相似性是可以得到比较好的保证的。当然，这只是理想的情况，因为均方误差较小的两个块的中心像素其实也有可能差别很大。另外，这种方法的缺点也很明显，因为我们需要为图像中的每个像素寻找相似块，其时间复杂度无疑是非常高的，处理一张低分辨率的图片也可能需要十几秒之久。
2.3 变换域稀疏表示
图像的变换域一般是相对于空间域来说的，即像素值在水平或垂直方向的波动快慢。很明显，图像的内容通常以低频为主，例如随处可见的大片纯色或者渐变等等，只有少数的细节部分才会包含一些较高频的信息，而特别高频的如密集的马赛克等等反而会引起观察者生理和心理上的不适。所以，图像变换到频域之后通常只有少数几个低频点会有较大的系数或者能量，而绝大部分高频点上的系数或能量则几乎可以忽略，这就是图像在变换域中的稀疏性。
图 1 展示了一幅自然无噪声图像（来自 SIDD 图像降噪数据库，下同）的 DCT-II 变换结果，右侧为相应的 DCT 系数的绝对值，为了方便作图已经把最低频的 4 × 4 4times4 4×4 的系数置零，因为这部分的系数本身很大，例如 ( 0 , 0 ) (0, 0) (0,0) 这个点实际上就是图像所有像素值的和除以边长，如果把它们放在图中其他值几乎可以忽略。可以看到，整个变换结果中只有低频区域很小一部分有较大的值，也就是说低频分量占据了图像的绝大部分能量。
图1 无噪声图像与其DCT变换系数
接下来，我们再看看添加了 AWGN 的图像的 DCT 变换系数是什么样子的。图 2 左侧为图 1 无噪声图像添加了标准差为 σ = 20 sigma = 20 σ=20 的高斯白噪声的结果，右侧为相应的 DCT 变换系数。可以看到，虽然低频分量还是占据了绝大部分的能量，但由于噪声的缘故，其他的频率点的系数已经远远没有图 1 那样平滑。甚至我们可以发现，即便是最高频的地方的变换系数的量级竟然和那些接近低频的地方是差不多的，或者说，除了最低频的地方，其他的变换系数就像是一些独立同分布的噪声。实际上，不仅是 DCT，可以证明 AWGN 经过任意次数的正交变换之后仍然是 AWGN ，并且其分布和原始噪声的分布一致，即具有相同的方差（见链接2）。
图2 添加了AWGN的图像与其DCT变换系数
正交变换在图像变换中十分普遍，例如 DCT, DST 以及 Haar 变换等等。对于正交变换，我们可以使用一个正交矩阵 A bf{A} A 来表示其变换核，并且有
T v = A v = A u A n = T u T n . (2.7) { {mathbf{T}}_v} = {mathbf{Av}} = {mathbf{Au}} {mathbf{An}} = { {mathbf{T}}_u} { {mathbf{T}}_n}.tag{2.7} Tv=Av=Au An=Tu Tn.(2.7)
其中 v bf{v} v, u bf{u} u, n bf{n} n 分别代表有噪声信号 v v v，无噪声信号 u u u 和噪声信号 n n n 的一个向量（例如某一列/行）， T v {bf{T}}_v Tv, T u {bf{T}}_u Tu, T n {bf{T}}_n Tn 为相应的正交变换系数。注意，虽然这里只表示了一维的变换，但对于可分离多维正交变换，其相当于先对某一维（例如每一列）进行正交变换，然后对所得的结果进行其他维（例如每一行）的正交变换，比如对于二维的变换，有
T v = ( A ( A V ) T ) T = ( A ( A U ) T ) T ( A ( A N ) T ) T = T u T n . (2.8) { {mathbf{T}}_v} = {left( { {mathbf{A}}{ {({mathbf{AV}})}^T}} right)^T} = {left( { {mathbf{A}}{ {({mathbf{AU}})}^T}} right)^T} {left( { {mathbf{A}}{ {({mathbf{AN}})}^T}} right)^T} = { {mathbf{T}}_u} { {mathbf{T}}_n}.tag{2.8} Tv=(A(AV)T)T=(A(AU)T)T (A(AN)T)T=Tu Tn.(2.8)
所以，对于可分离的多维正交变换，有噪声信号 v v v 的变换系数总可以表示为无噪声信号 u u u 的变换系数与噪声信号 n n n 的变换系数的和。而前面已经提到，AWGN 经过多次正交之后仍然是 AWGN，所以有噪声信号 v v v 的变换系数其实就是无噪声信号 u u u 的变换系数加上与原始噪声同样方差的 AWGN 的结果。
因此，对于方差较小的噪声，因为无噪声图像本身的变换系数主要集中在低频区域，并且具有较大的幅度，而噪声在变换域中由于方差较小的缘故，通常只有较小的值（例如 95.4% 的值都在两个标准差的范围内），我们可以通过设定阈值的方法（一般为 2~3 个噪声标准差），只保留有噪声图像v的变换系数中那些具有较大幅度的值，其他的则置零。这样，绝大部分的噪声的变换系数就得以消除，只有少数还残留在低频的系数中，从而经过反变换后能够极大地减少噪声，同时还能较好地保留图像的纹理。当然，那些物体边缘的区域可能会由于高频信息的丢失而稍微显得模糊，或者有少量的振铃效应等等。
图3 真实强噪声图像与其DCT变换系数
然而，当噪声方差本身很大时，噪声的能量和无噪声图像本身的能量级别相当，那么在变换域中无噪声图像部分次低频的变换系数反而有可能被噪声给掩盖掉。例如，图 3 是一幅在低照度环境下真实拍摄的图像，其实际上是用于叠加生成图 1 中的无噪声图像的上百张副本之一。可以看到，在空间域中，图像本身的信息已经极大地被噪声破坏，而在变换域中，我们也可以看到与图 2 人工噪声非常相似的杂乱的噪声，这也在一定程度上表明图像在低照度环境下的噪声确实以 AWGN为 主。另外，相比于图 2，图 3 变换域中的噪声幅度明显要大得多，普遍可以达到 80-100，所以其噪声的标准差应该大概在 40-50 左右。回顾图 1 无噪声图像的变换系数，其实有很大一部分的低频系数的幅度都在 100 以内，如果这时我们直接通过设置阈值的方法去除那些幅度低于 2~3 标准差的变换系数，就会把很多本身有用的信息给丢弃掉，最终我们得到的可能只是一幅过度平滑的图像，没有任何的细节，这无疑是不可接受的。况且，因为 AWGN 在变换域中仍然是 AWGN，所以不仅是高频的系数，我们所保留的低频系数本身也已经是被噪声污染了的。例如在图 3 中，最大的系数幅度（注意这里不包含最低频的 4 × 4 4times4 4×4 系数）只在 300-400 之间，而在图 1 中其应该在 400-500 的范围，所以即便我们只保留低频的信息还是会有比较大的失真的。
为此，类似于非局部均值（NLM）算法，一种被称之为协同滤波（Collaborative Filtering）的方法被提出。其原理也比较简单，即把那些相似的块叠加起来形成第 3 维，在前两维的变换基础上，再进行第 3 维的正交变换，之后就可以通过设置阈值的方法把那些低于一定幅度的变换系数去掉。那么，这种方法相比于前面的二维变换又有怎样的改进呢？其实，其基本的依据还是在于 AWGN 经过正交变换之后仍然是与原始噪声同分布的 AWGN，也就是说，无论我们进行多少次正交变换，噪声的强度都是不会发生改变的。然而，因为我们要处理的是相似的块，所以这些块本身无噪声的部分在第 3 维上是几乎一样的，也就是低频分量占据了主要部分，从而经过第 3 维的变换之后我们可以使得无噪声信号的能量更加集中，具体表现在属于无噪声信号的那部分系数幅度将变得更加的大。而因为噪声信号经过第 3 维的变换之后幅度没有变化，所以两者的差距将进一步地拉大，这时进行阈值置零操作就能够减少误杀的情况，而在那些没有被置零的系数中噪声信号的能量比例也进一步地缩小，从而更好地保留图像本身的高频信息。
图4 不同数量相似块对变换系数幅度的影响
为了便于理解，图 4 展示了在不同数量的相似块下有噪声图像的变换系数幅度的变化情况。注意，为了方便作图，这里只使用了图 1 中无噪声图像的某一列，其变换系数如左上角所示，并且直流系数已经被置零。为了模拟不同的相似块，这里只是简单重复地为该列加上标准差为 σ = 50 sigma = 50 σ=50 的高斯白噪声，然后把它们按列组合成一张二维图像，这时进行行方向的变换也就相当于我们前面所说的协同变换。因为经过行方向的协同变换后，变换系数的能量主要集中到第 1 列，所以图 4 只给出第 1 列的变换系数的绝对值，并且同样把直流分量置零。
可以看到，当只有单纯一个添加了噪声的块时，如图 4 右上角所示，由于噪声的方差很大，甚至达到了无噪声图像本身的量级，经过变换之后其系数和原始的相比除了直流分量已经完全无法分辨了，自然也不可能保留太多的细节。而随着相似块数量的不断增加，如图 4 下半部分，噪声的变换系数的量级并没有发生改变，始终保持在 100 以下，而其他原本属于无噪声图像的变换系数则随着能量的集中而不断地增大，从而与噪声的变换系数分离开来，这样我们就可以方便地抑制属于噪声信号的那部分能量，而不会过多地影响原本属于无噪声信号的那一部分，更好地保留图像的细节。
这时，你可能会想，如果我们直接把这些相似块的变换系数加起来求平均不是更好吗？其实，对于线性变换来说，变换域的相加和原始空间域的相加是等价的。我们已经说过，由于不同的相似块并不是完全一样的，所以直接把它们叠加起来只会让图像更加糟糕。图 4 只是展示了一种理想的情况，即所有相似块对应的无噪声图像都是一样的，所以左下角和右下角的变换系数才会和左上角那么相似。另外，这里也存在一个矛盾，即越多的相似块能够让低频分量的能量更加集中，但也会导致块与块之间的差异点越来越多，从而引入一部分的较高频的信息，但这部分系数的能量通常又不足以与噪声系数的能量区分开来，所以这部分信息也很有可能和噪声一起被抹去，经过反变换之后丢失更多的细节。
相对于 NLM 算法来说，协同滤波同样需要寻找足够的相似块。但是，NLM 是以每个像素点为操作单位的，也就是需要为每个像素找到足够多的相似块，然后把这些相似块的中心像素进行加权叠加，从而把该像素的噪声方差降低到足够小。所以，虽然其时间复杂度很高，但能够比较好地保留图像的细节。相反，协同滤波可以一次性处理所有块的所有像素，这样速度无疑要快很多，但是在细节保留方面则要逊色不少，因为该算法本身并没有把噪声的方差缩小，而无论我们怎样把低频分量的能量集中，总会有一些细节部分的的能量会被淹没在噪声当中，无法恢复。
另外，为了寻找相似块，我们总不能把整幅图像作为一个参考块，所以协同滤波本身也只是处理了图像的一小部分区域。为了保证图像中所有像素都能被覆盖到，最保险的方法是将图像拆分为规整的矩形块，然后使用栅格扫描的方式，对每一个区域都进行一次相似块匹配以及协同滤波。但这时又有一个问题，因为每个区域的相似块都是在整幅图像上寻找的，所以不同区域的相似块总可能包含同样的像素，也就是同一个像素通常会被进行多次的协同滤波。所以，一个很直观的改进方法就是把这些不同的协同滤波的结果进行整合，从而对当前像素的无噪声值做出更恰当的估计，相应地在使得图像更加纯净的同时能够更好地保留自身的细节。这就是下面我们所要说的 BM3D 降噪算法的基本思路。实际上，为了达到更好的降噪效果，通常来说这些拆分的区域本身也应该是有重叠的，这样可以进一步提高某个像素被多次滤波的可能，这种方法称之为过完备（Overcompleteness），并且被包括 BM3D 等多种算法所采用。。
3 BM3D算法原理
在第 2 节中，我们其实已经对 BM3D 算法的基本理论依据做了详细的描述，其主要可分为三个步骤，即首先对每个参考块进行相似块匹配（Block-Matching）并分别得到一个三维的组合，然后对其进行协同变换和滤波（3D-Transform），最后对各个参考块对应组合的滤波结果进行整合（Aggregation），从而得到最终的降噪结果。在此之上，为了进一步改善图像的质量，BM3D 实际进行了两次降噪，即将以上三个步骤再重复了一遍，但是具体的块匹配标准、滤波方式以及整合权重等等会有一些区别，其基本流程如图 5 所示。
图5 BM3D算法基本流程
在 Step1 中，协同滤波主要通过设置硬阈值（Hard-thresholding）的方法来去除噪声的能量，这也就是我们在第 2 节中所介绍的方法。这种方法相对来说比较简单，但是一方面可能误杀掉部分细节的信息，另一方面也无法去除噪声在低频系数中的能量。另外，在最初的有噪声图像中，块匹配的误差其实是一个随机变量，根据式 (2.4) 可知其方差大概正比于噪声方差的平方，当噪声方差较大时块匹配的质量并不是很理想，也就可能引入更多块与块之间不同的细节信息，并且被硬阈值滤波给误杀掉。所以，Step1 通常只是作为一种基础估计（Basic Estimate），但由于过完备以及块整合的缘故，其结果也是也是可以接受的。
作为改善手段，Step2 使用了 Step1 中基础估计的结果来进行相似块的匹配，因为噪声已经极大地被消除了，所以我们认为这时的匹配结果是更加准确的。注意，这里并不是直接复用 Step1 中已经找到的相似块的位置，而是利用基础估计图像重新搜寻一遍，所以 Step2 和 Step1 所需的时间其实是差不多的。但是，不同于 Step1 使用硬阈值的方法来进行协同滤波，对基础估计图像再进行一次硬阈值滤波其实是没有多大用处的，因为这样只会让得到的图像更加模糊。当然，这时我们可能会想使用在基础估计图像上找到的相似块位置对原始有噪声图像进行硬阈值滤波，但这种方法始终不能避免前面 Step1 所说的缺点，即误杀细节，以及无法去除低频系数的噪声。实际上，Step2 采用了经验维纳收缩（Emprical Wiener Shrinkage）的方法，其根据基础估计图像上的协同变换系数的功率谱以及噪声的强度，对原始有噪声图像同样位置 3D 块的协同变换系数进行收缩，然后对收缩后的系数进行反变换以及整合，直接得到最终所要的降噪结果。所以在这个过程中，我们实际处理的还是输入的原始有噪声图像，而基础估计只是作为一种辅助，用于更准确地确定原始图像中相似块的位置，以及相应的系数收缩比例。具体的原理将会在后面进行介绍。
为了方便阅读原论文，这里的所用到符号会与其保持一致。对于一幅添加了 AWGN 的图像，其可表示为
z ( x ) = y ( x ) η ( x ) , η ( x ) ∼ N ( 0 , σ 2 ) , x ∈ X ⊂ Z 2 . (3.1) z(x) = y(x) eta (x),;eta (x) sim N(0,;{sigma ^2}),;x in X subset {mathbb{Z}^2}.tag{3.1} z(x)=y(x) η(x),η(x)∼N(0,σ2),x∈X⊂Z2.(3.1)
注意，这里所说的图像都是 2D 单通道图像，对于彩色图像的降噪将在后面进行介绍。我们用 Z x Z_x Zx 来代表有噪声图像 z z z 上的一个左上角坐标为 x x x 的大小为 N 1 × N 1 N_1times N_1 N1×N1 的正方形块，而 Z S {bf{Z}}_S ZS 为一个 3D 块，所包含的每个 2D 块的左上角坐标由集合 S S S 来定义。为了表示每个 3D 块中的参考块，我们会在相应的坐标下添加 R R R (Reference) 的标识，即 x R x_R xR。同样，对于无噪声图像 y y y 上的一个块，我们也可以使用 Y x Y_x Yx 来表示。因为 Step1 和 Step2 的步骤基本一致，所以它们所用到的一些参数的意义也是一样的，为了区分两者，我们会分别加上 “ ht text{ht} ht” 和 “ wie text{wie} wie” 的标识。对于 Step1 的基础估计图像，我们用 y ^ baisc {hat y^{ {text{baisc}}}} y^baisc来表示，Step2 所得的最终估计图像则为 y ^ final {hat y^{ {text{final}}}} y^final。
3.1 Step1
接下来我们将对 BM3D 算法的 Step1 进行详细的介绍。首先有必要再次说明参考块的具体意义，不同于 NLM 算法会把每个像素的邻域都作为一个参考块，为了降低计算的复杂度，一般在算法实现中，我们会从图像的第一个 N 1 × N 1 N_1times N_1 N1×N1 块开始，采用栅格扫描的方式，每次平移几列（行），例如用 N s t e p < N 1 N_{step} lt N_1 Nstep<N1 来表示，一般为 3~6，直至铺满整幅图像。因为我们必须保证所有的像素都至少参与一次协同滤波，所以所有的块都会被选择一次作为参考块，而参考块的重叠也是过完备（Overcompleteness）的一种体现。在第 2 节中我们已经提到，对于无噪声图像 y y y 上某个大小为 N 1 ht × N 1 ht N_1^{ {text{ht}}} times N_1^{ {text{ht}}} N1ht×N1ht 的参考块 Y x R {Y_{ {x_R}}} YxR 和其他位置同样大小的块 Y x Y_x Yx，记其匹配误差或者块距离为
d ideal ( Y x R , Y x ) = ∥ Y x R − Y x ∥ 2 2 ( N 1 ht ) 2 . (3.2) {d^{ {text{ideal}}}}left( { {Y_{ {x_R}}},{Y_x}} right) = frac{ {left| { {Y_{ {x_R}}} – {Y_x}} right|_2^2}}{ { { {left( {N_1^{ {text{ht}}}} right)}^2}}}.tag{3.2} dideal(YxR,Yx)=(N1ht)2∥YxR−Yx∥22.(3.2)
如果这两个块是没有重叠的，也就是每个像素都是独立的，那么在有噪声图像 z z z 上对应位置的两个块的的匹配误差的数学期望以及方差分别为
E [ d noisy ( Z x R , Z x ) ] = d ideal ( Y x R , Y x ) 2 σ 2 . (3.3) Eleft[ { {d^{ {text{noisy}}}}left( { {Z_{ {x_R}}},{Z_x}} right)} right] = {d^{ {text{ideal}}}}left( { {Y_{ {x_R}}},{Y_x}} right) 2{sigma ^2}.tag{3.3} E[dnoisy(ZxR,Zx)]=dideal(YxR,Yx) 2σ2.(3.3)
D [ d noisy ( Z x R , Z x ) ] = 8 σ 4 ( N 1 ht ) 2 8 σ 2 d ideal ( Y x R , Y x ) ( N 1 ht ) 2 . (3.4) Dleft[ { {d^{ {text{noisy}}}}left( { {Z_{ {x_R}}},{Z_x}} right)} right] = frac{ {8{sigma ^4}}}{ { { {left( {N_1^{ {text{ht}}}} right)}^2}}} frac{ {8{sigma ^2}{d^{ {text{ideal}}}}left( { {Y_{ {x_R}}},{Y_x}} right)}}{ { { {left( {N_1^{ {text{ht}}}} right)}^2}}}.tag{3.4} D[dnoisy(ZxR,Zx)]=(N1ht)28σ4 (N1ht)28σ2dideal(YxR,Yx).(3.4)
所以，当噪声方差很大或者块比较小的时候，直接在有噪声图像 z z z 上寻找相似块是很不准确的，有可能无噪声图像 y y y 上两个差距很大的块在有噪声图像 z z z 上由于方差很大的原因差距显得比较小，从而被判定为相似，造成错误的匹配。
为了解决这个问题，论文提出先对有噪声图像 z z z 上的两个块进行可分离的二维正交变换，然后将那些幅度小于一定阈值的系数置零，这样两个块的匹配误差就可表示为这些系数的均方误差，即
d ( Z x R , Z x ) = ∥ Υ ′ ( T 2 D ht ( Z x R ) ) − Υ ′ ( T 2 D ht ( Z x ) ) ∥ 2 2 ( N 1 ht ) 2 . (3.5) dleft( { {Z_{ {x_R}}},{Z_x}} right) = frac{ {left| {Upsilon ‘left( {mathcal{T}_{2D}^{ {text{ht}}}left( { {Z_{ {x_R}}}} right)} right) – Upsilon ‘left( {mathcal{T}_{2D}^{ {text{ht}}}left( { {Z_x}} right)} right)} right|_2^2}}{ { { {left( {N_1^{ {text{ht}}}} right)}^2}}}.tag{3.5} d(ZxR,Zx)=(N1ht)2∥∥Υ′(T2Dht(ZxR))−Υ′(T2Dht(Zx))∥∥22.(3.5)
其中， Υ ′ Upsilon ‘ Υ′ 为硬阈值操作，相应的阈值为 λ 2 D σ {lambda _{2D}}sigma λ2Dσ，一般为两个标准差。 T 2 D ht mathcal{T}_{2D}^{ {text{ht}}} T2Dht 为相应的可分离二维正交变换。我们前面已经有提到，AWGN 经过正交变换之后还是 AWGN，而无噪声图像由于其在变换域中的稀疏特性，一般只有在少数低频的位置上有较大的值。所以经过二维正交变换之后，通过设置硬阈值的方法，我们可以去除绝大部分噪声的能量。虽然也有一部分图像本身的能量被去除，但其影响基本可以忽略，因此这时我们可以使得块匹配更加准确。注意，我们不需要对系数进行反变换再在空间域上计算两个块的匹配误差，因为对于正交变换以及其反变换来说，变换域和空间域的距离其实是一样的。
实际上，并不是所有的情况都需要先对变换系数进行硬阈值操作再计算块距离，论文中也指出一般在噪声的标准差 σ > 40 sigma > 40 σ>40 时才需要设置一个非零的阈值，对于阈值为零的情况，式 (3.5) 所得的块距离与直接在空间域计算是等价的。另外，我们可以先计算所有可能的候选块的 2D 变换系数并缓存起来，因为相邻的参考块的需要搜索的候选块大多数是一样的，所以这部分结果也可以进行复用。当然，我们不需要把整幅图像的候选块的变换都存储起来，因为全搜索不可能应用于整幅图像，而是只有某个较小的搜索窗口。当搜索窗口移出某个候选块时，我们就可以删除其变换系数以节省存储空间了。另一方面，因为我们已经计算得到所有相似块的 2D 变换系数，如果协同变换采用的是可分离的 3D 正交变换，那么这些 2D 变换的结果就可以进行复用，我们只需再进行第 3 维的 1D 变换即可，从而进一步降低其复杂度。
回到 Step1 的块匹配，根据式 (3.5) 计算得到当前参考块与其他所有块（实际上只有某个邻域内的块会被考虑）的匹配误差后，我们只保留那些误差小于一定阈值的块，并得到相应的坐标集合
S x R ht = { x ∈ X : d ( Z x R , Z x ) ⩽ τ match ht } . (3.6) S_{ {x_R}}^{ {text{ht}}} = left{ {x in X:dleft( { {Z_{ {x_R}}},{Z_x}} right) leqslant tau _{ {text{match}}}^{ {text{ht}}}} right}.tag{3.6} SxRht={ x∈X:d(ZxR,Zx)⩽τmatchht}.(3.6)
其中 τ match ht tau _{ {text{match}}}^{ {text{ht}}} τmatchht 为判定两个块相似的最大误差阈值，一般根据经验来设置。很明显当前参考块肯定会被判定为相似块，因为其匹配误差为 0，所以坐标集合中至少有 1 个元素。这时我们把所有的相似块堆叠起来，就可以得到一个形状为 N 1 ht × N 1 ht × ∣ S x R ht ∣ N_1^{ {text{ht}}} times N_1^{ {text{ht}}} times left| {S_{ {x_R}}^{ {text{ht}}}} right| N1ht×N1ht×∣∣SxRht∣∣ 的 3D 数组 Z S x R ht { {mathbf{Z}}_{S_{ {x_R}}^{ {text{ht}}}}} ZSxRht，其中 ∣ S x R ht ∣ left| {S_{ {x_R}}^{ {text{ht}}}} right| ∣∣SxRht∣∣ 代表该集合元素的个数。注意，通常来说该 3D 数组中块的顺序并不是很重要，但有时候因为第 3 维的正交变换可能需要输入的长度为 2 的幂，例如 Haar 变换，这时我们只保留集合中那些匹配误差最小的块以满足长度的要求。
得到参考块对应的 3D 数组后，我们就可以进行 3D 的协同变换和滤波，形式化地可表示为
Y ^ S x R ht ht = T 3 D ht − 1 ( Υ ( T 3 D ht ( Z S x R ht ) ) ) . (3.7) widehat {mathbf{Y}}_{S_{ {x_R}}^{ {text{ht}}}}^{ {text{ht}}} = mathcal{T}{_{3D}^{ {text{ht}}-1}}left( {Upsilon left( {mathcal{T}_{3D}^{ {text{ht}}}left( { { {mathbf{Z}}_{S_{ {x_R}}^{ {text{ht}}}}}} right)} right)} right).tag{3.7} Y SxRhtht=T3Dht−1(Υ(T3Dht(ZSxRht))).(3.7)
其中 T 3 D ht mathcal{T}_{3D}^{ {text{ht}}} T3Dht 为相应的 3D 正交变换，一般来说为了降低计算的复杂度我们会选用可分离的变换，因为前面在计算块匹配误差时二维的变换结果已经存在了，我们完全可以对其进行复用。所以， T 3 D ht mathcal{T}_{3D}^{ {text{ht}}} T3Dht 可以表示为前两维的 T 2 D ht mathcal{T}_{2D}^{ {text{ht}}} T2Dht 和第 3 维的 T 1 D ht mathcal{T}_{1D}^{ {text{ht}}} T1Dht 的组合。具体的变换类型没有太多的限制，一般前两维的变换会选用傅里叶变换类型的如 DCT 等，可以更好地捕捉图像块中周期性的信息，而第 3 维的变换则会选用小波变换类型的如 Haar 等，可以通过尺度缩放等更好地捕捉块与块之间的局部相似性，但实验数据也表明其性能差别并不是很大。 Υ Upsilon Υ 为另一个硬阈值操作，相应的阈值为 λ 3 D σ {lambda _{3D}}sigma λ3Dσ，通常为 2～3 个标准差。 Y ^ S x R ht ht widehat {mathbf{Y}}_{S_{ {x_R}}^{ {text{ht}}}}^{ {text{ht}}} Y SxRhtht 包含了 ∣ S x R ht ∣ left| {S_{ {x_R}}^{ {text{ht}}}} right| ∣∣SxRht∣∣ 个块 Y ^ x ht, x R widehat Y_x^{ {text{ht,}};{x_R}} Y xht,xR，表示经过式 (3.7) 的协同滤波所得的结果，对应的参考块的坐标为 x R x_R xR。注意，因为不同的参考块可能包含同样的相似块，所以一个块会对应多次的协同滤波结果，我们有必要通过其对应的参考块坐标进行区分，这些不同的协同滤波结果将会在后面的操作中进行整合，从而为每个像素获得 Step1 的基础估计。
得到每个参考块以及其对应的3D组的协同滤波结果后，因为参考块本身就是有重叠的，而且不同的参考块也可能包含重叠的相似块，为了得到每个像素的基础估计，我们有必要对这些来自不同组合的协同滤波结果进行整合，例如赋予这些不同的组合以不同的权值。为了获得一个比较好的结果，我们当然是希望每个组合的协同滤波都能够完美地去除噪声，但这明显是不可能的，因为经过多次的正交变换以后，AWGN 还是会以 AWGN 的形式残留在那些没有被去除的低频系数中，没有被置零的系数越多，残留的噪声能量就越多。另外，如果一个组合内的所有块都足够相似，那么其变换系数也本应该是非常稀疏的。所以，我们可以根据协同滤波后剩余的非零系数的个数来进行权值的分配，假设该 3D 组合内所有的像素都是独立的，记硬阈值操作后剩下的系数个数为 N hard x R N_{ {text{hard}}}^{ {x_R}} NhardxR，则有
w x R ht = { 1 σ 2 N hard x R , if N hard x R ⩾ 1 1. otherwise (3.8) w_{ {x_R}}^{ {text{ht}}} = left{ {begin{array}{c} {tfrac{1}{ { {sigma ^2}N_{ {text{hard}}}^{ {x_R}}}},}&{ {text{if }}N_{ {text{hard}}}^{ {x_R}} geqslant 1} \ {1.}&{ {text{otherwise}}} end{array}} right.tag{3.8} wxRht={ σ2NhardxR1,1.if NhardxR⩾1otherwise(3.8)
实际上，一个 3D 组合内的相似块通常都会有重叠，也就是说所有像素之间并不会完全独立，但为了简单起见我们还是会以式 (3.8) 来进行权值的计算。注意，这里的权值对于该 3D 组合内所有的像素都是一样的，但为了减少边界效应，通常都会为每个块加上一个窗函数，从而赋予中心像素更高的权值。论文中使用了 Kaiser 窗，其一维定义为
K ( n ) = I 0 ( β 1 − 4 n 2 ( M − 1 ) 2 ) / I 0 ( β ) , − M − 1 2 ⩽ n ⩽ M − 1 2 . (3.9) K(n) = {I_0}left( {beta sqrt {1 – frac{ {4{n^2}}}{ { { {left( {M – 1} right)}^2}}}} } right)/{I_0}left( beta right),; – frac{ {M – 1}}{2} leqslant n leqslant frac{ {M – 1}}{2}.tag{3.9} K(n)=I0(β1−(M−1)24n2 )/I0(β),−2M−1⩽n⩽2M−1.(3.9)
其中 I 0 I_0 I0 为零阶修正的贝塞尔函数，一般在 MATLAB 或者 Python.numpy 中可以直接调用。 n n n 为相对于中心点的距离，所以对于偶数长度的离散窗其应该为 1/2 的奇数倍。二维的 Kaiser 窗只需把列向量的一维 Kaiser 窗乘以其转置即可，在 BM3D 中其大小应该为 N 1 × N 1 N_1times N_1 N1×N1。注意这个窗是对单个块而言的，同一个 3D 组的各个相似块使用同样的窗函数。因为同一个像素可能包含在不同的相似块中，通过添加 Kaiser 窗也可以使得属于不同块不同位置的同一个像素获得不同的权值。最后，通过对某个像素所属的所有参考块以及对应的相似块的加权叠加，我们就可以得到其基础估计，即
y ^ basic ( x ) = ∑ x R ∈ X ∑ x m ∈ S x R ht w x R ht Y ^ x m ht, x R ( x ) ∑ x R ∈ X ∑ x m ∈ S x R ht w x R ht X x m ( x ) , ∀ x ∈ X . (3.10) {hat y^{ {text{basic}}}}(x) = frac{ {sumlimits_{ {x_R} in X} {sumlimits_{ {x_m} in S_{ {x_R}}^{ {text{ht}}}} {w_{ {x_R}}^{ {text{ht}}}widehat Y_{ {x_m}}^{ {text{ht, }}{x_R}}(x)} } }}{ {sumlimits_{ {x_R} in X} {sumlimits_{ {x_m} in S_{ {x_R}}^{ {text{ht}}}} {w_{ {x_R}}^{ {text{ht}}}{mathcal{X}_{ {x_m}}}(x)} } }},;forall x in X.tag{3.10} y^basic(x)=xR∈X∑xm∈SxRht∑wxRhtXxm(x)xR∈X∑xm∈SxRht∑wxRhtY xmht, xR(x),∀x∈X.(3.10)
其中我们假设每个块都已经通过补零的方法填充到和原图一样的大小，并且通过 X x m : X → { 0 , 1 } {mathcal{X}_{ {x_m}}}:X to { 0,;1} Xxm:X→{ 0,1} 来判断某个像素 x x x 是否在块 x m x_m xm 上。从式中也可以看出，一个像素既可以属于不同的参考块，也可以属于同一个参考块的不同相似块，而同一个块也可以是不同的参考块的相似块。通过这种过完备的协同滤波方法，我们一方面可以尽可能地去除噪声的能量，同时也能较好地保持图像中的细节，并且可以保持较低的计算复杂度，因为这里主要是以块为操作单位的，相比于 NLM 算法可以减少很大一部分寻找相似块的时间。
3.2 Step2
其实，经过 Step1 后，我们得到的基础估计图像已经能够极大地去除噪声的能量。虽然前面也说了，简单的硬阈值操作并不能去除残留在低频系数上噪声，并且很多的细节信息也有可能被一同抹杀掉，不过经过反变换之后这部分噪声的能量会被分散稀释到各个像素中，另外通过整合不同的协同滤波的结果也可以在一定程度上弥补这两个内在的缺陷。这时我们可能会想，一些细节信息会被抹杀掉，其实一部分原因是 3D 组合内的块不够相似，从而该细节部分的能量不够集中而无法与噪声区分开来，而块与块不够相似的一部分原因是输入图像中的噪声过于强大，无法对块匹配误差进行准确的估计。因为通过 Step1 我们已经能够得到一张比较纯净的图像，那作为一种改善手段是不是可以通过在基础估计图像上寻找相似块，然后重复 Step1 的操作呢？这似乎可以形成一种迭代的过程。但实际上，这种方法的性能提升并不是很明显，因为我们始终无法避免硬阈值滤波的缺点。为此，不同于硬阈值滤波直接将小于阈值的系数置零而完全保留那些较大的系数，BM3D 引入了一种称为经验维纳收缩（Emprical Wiener Shrinkage）的方法，其实际上是根据基础估计图像上的变换系数的功率谱来对原始有噪声图像的变换系数进行收缩，从而能够更好地去除包括低频系数中的噪声能量，同时尽可能地保留原本属于无噪声图像的那部分高频系数。
Step2 与 Step1 的流程基本一致，首先也是为每一个参考块寻找足够的相似块，注意这时的参考块大小不一定与 Step1 相同。不同的地方在于，Step2 是直接在基础估计图像上寻找相似块的，因为这时候的噪声已经比较小了，我们不需要像式 (3.5) 那样先进行变换然后硬阈值去除部分幅度较小的系数，即
S x R wie = { x ∈ X : ∥ Y ^ x R basic − Y ^ x basic ∥ 2 2 ( N 1 wie ) 2 < τ match wie } . (3.11) S_{ {x_R}}^{ {text{wie}}} = left{ {x in X:frac{ {left| {widehat Y_{ {x_R}}^{ {text{basic}}} – widehat Y_x^{ {text{basic}}}} right|_2^2}}{ { { {left( {N_1^{ {text{wie}}}} right)}^2}}} < tau _{ {text{match}}}^{ {text{wie}}}} right}.tag{3.11} SxRwie=⎩⎪⎨⎪⎧x∈X:(N1wie)2∥∥∥Y xRbasic−Y xbasic∥∥∥22<τmatchwie⎭⎪⎬⎪⎫.(3.11)
注意，不同于 Step1，这里将产生两个 3D 组合，即 Y ^ S x R wie basic widehat {mathbf{Y}}_{S_{ {x_R}}^{ {text{wie}}}}^{ {text{basic}}} Y SxRwiebasic 和 Z S x R wie { {mathbf{Z}}_{S_{ {x_R}}^{ {text{wie}}}}} ZSxRwie，前者为基础估计图像上的相似块所组成的 3D 组合，后者则为原始有噪声图像上同样位置的块。这时，我们就可以介绍经验维纳收缩的原理了（详见链接3）。
维纳滤波是一种使得有噪声信号经过滤波后的输出与原始无噪声信号均方误差最小的线性滤波器，其讨论一般是在频域中的，一种常用的定义为
H ( u , v ) = H d ∗ ( u , v ) ∣ H d ( u , v ) ∣ 2 P n ( u , v ) P s ( u , v ) . (3.12) H(u,v) = frac{ {H_d^ * (u,v)}}{ { { {left| { {H_d}(u,v)} right|}^2} frac{ { {P_n}(u,v)}}{ { {P_s}(u,v)}}}}.tag{3.12} H(u,v)=∣Hd(u,v)∣2 Ps(u,v)Pn(u,v)Hd∗(u,v).(3.12)
其中 H ( u , v ) H(u,v) H(u,v) 为维纳滤波器的频率响应函数， H d ( u , v ) H_d(u,v) Hd(u,v) 为退化系统的频率响应函数。退化系统也就是原始无噪声信号在被噪声干扰之前经过的某些前处理，例如低通滤波器，乃至由于相机运动所导致的动态模糊都可以归为退化系统，但在本文中我们需要处理的主要是噪声，所以可认为其恒为 1。 P n ( u , v ) P_n(u,v) Pn(u,v) 和 P s ( u , v ) P_s(u,v) Ps(u,v) 分别代表噪声和原始无噪声信号的功率谱，简单地说就是其频谱幅值的平方。然而问题在于，我们不可能知道噪声的具体数值，因为这是一个随机变量。所以，为了实现上面所描述的维纳滤波器，我们只能通过一些经验来进行估计。这里的经验包括两个部分，一是基础估计图像，因为我们认为它已经将绝大部分噪声去除了，所以从经验上讲它就是无噪声图像，这样就可得到 P s ( u , v ) P_s(u,v) Ps(u,v)；二是噪声的方差，通过噪声的方差我们可以估计噪声功率谱的数学期望，即 P n ( u , v ) P_n(u,v) Pn(u,v)。
前面我们已经多次提到，AWGN 经过正交变换以后仍然是 AWGN，而正交变换满足线性叠加，所以有噪声信号在变换域中仍然可以表示为无噪声信号的变换系数加上 AWGN 的结果。具体地，我们用 T y T_y Ty 来表示无噪声图像的变换系数， T z T_z Tz 为相应有噪声图像的变换系数，那么有
T z ( x ) = T y ( x ) η ( x ) , η ( x ) ∼ N ( 0 , σ 2 ) . (3.13) {T_z}(x) = {T_y}(x) eta (x),;eta (x) sim mathcal{N}(0,;{sigma ^2}).tag{3.13} Tz(x)=Ty(x) η(x),η(x)∼N(0,σ2).(3.13)
那么可求得噪声以及有噪声图像的功率谱的数学期望分别为
E ( η 2 ( x ) ) = D ( η ( x ) ) E 2 ( η ( x ) ) = σ 2 , E ( T z 2 ( x ) ) = E [ ( T y ( x ) η ( x ) ) 2 ] = T y 2 ( x ) σ 2 . (3.14) begin{aligned} & Eleft( { {eta ^2}(x)} right) = Dleft( {eta (x)} right) E^2left( { {eta }(x)} right) = {sigma ^2}, \ & Eleft( {T_z^2(x)} right) = Eleft[ { { {left( { {T_y}(x) eta (x)} right)}^2}} right] = T_y^2(x) {sigma ^2}. \ end{aligned} tag{3.14} E(η2(x))=D(η(x)) E2(η(x))=σ2,E(Tz2(x))=E[(Ty(x) η(x))2]=Ty2(x) σ2.(3.14)
所以，噪声信号的功率谱的数学期望刚好为其方差，而有噪声图像的功率谱的数学期望则为无噪声图像的功率谱与噪声方差的和。虽然噪声既可能为正也可能为负，也就是某个变换系数的功率谱既可能增加也可能减少，但是从数学期望上来说还是在增加的，增加的量就是噪声的方差。因为我们已经反复提到了基础估计图像可以作为原始无噪声图像的经验估计，所以可以认为无噪声图像的频谱 T y T_y Ty 可用基础估计图像的频谱来代替。那么，一种直观的滤波方法可定义为
T z ′ ( x ) = T y ^ b a s i c 2 ( x ) T y ^ b a s i c 2 ( x ) σ 2 T z ( x ) . (3.15) {T’_z}(x) = frac{ {T_{ { {hat y}^{basic}}}^2(x)}}{ {T_{ { {hat y}^{basic}}}^2(x) {sigma ^2}}}{T_z}(x).tag{3.15} Tz′(x)=Ty^basic2(x) σ2Ty^basic2(x)Tz(x).(3.15)
也就是说我们可以根据基础估计图像变换系数的功率占比来对有噪声图像的变换系数进行调整。但我们也发现，实际上所有的变换系数的幅度都是在缩小的，然而这也是无奈之举，因为我们并不知道噪声真正的值。虽然有噪声图像的某些变换系数本身就因为噪声的原因变小了，经过式 (3.15) 的缩小后似乎问题进一步加剧，但是由于噪声是随机信号，从数学期望上讲这种处理是有效的，因为很多变换系上的噪声功率是明显大于其方差的。那么，另一个问题在于，为什么式 (3.15) 使用功率谱比例对变换系数本身进行缩小，而不需要对该比例进行开根号？并且这种方法与维纳滤波即式 (3.12) 又有什么联系？其实，式 (3.15) 与式 (3.12) 是一样的，将基础估计图像和噪声的功率谱的经验估计代入式 (3.12) 可得
H ( x ) = 1 1 P n ( x ) P s ( x ) = P s ( x ) P s ( x ) P n ( x ) = T y ^ b a s i c 2 ( x ) T y ^ b a s i c 2 ( x ) σ 2 . (3.16) H(x) = frac{1}{ {1 frac{ { {P_n}(x)}}{ { {P_s}(x)}}}} = frac{ { {P_s}(x)}}{ { {P_s}(x) {P_n}(x)}} = frac{ {T_{ { {hat y}^{basic}}}^2(x)}}{ {T_{ { {hat y}^{basic}}}^2(x) {sigma ^2}}}.tag{3.16} H(x)=1 Ps(x)Pn(x)1=Ps(x) Pn(x)Ps(x)=Ty^basic2(x) σ2Ty^basic2(x).(3.16)
因为维纳滤波是均方误差最优的线性滤波器，所以从数学期望上讲式 (3.15) 的滤波结果可以使得输出的图像与原始无噪声信号的均方误差最小。因为所有变换系数的幅度都是缩小的，这时的维纳滤波也称为经验维纳收缩。
回到 BM3D 算法中，分别得到基础估计图像以及原始有噪声图像上的同样位置的两个 3D 组合 Y ^ S x R wie basic widehat {mathbf{Y}}_{S_{ {x_R}}^{ {text{wie}}}}^{ {text{basic}}} Y SxRwiebasic 和 Z S x R wie { {mathbf{Z}}_{S_{ {x_R}}^{ {text{wie}}}}} ZSxRwie 后，根据式 (3.15)，可得
Y ^ S x R wie wie = T 3 D wie − 1 ( W S x R wie ∘ T 3 D wie ( Z S x R wie ) ) , w h e r e W S x R wie = ∣ T 3 D wie ( Y ^ S x R wie basic ) ∣ 2 ∣ T 3 D wie ( Y ^ S x R wie basic ) ∣ 2 σ 2 . (3.17) begin{gathered} widehat {mathbf{Y}}_{S_{ {x_R}}^{ {text{wie}}}}^{ {text{wie}}} = mathcal{T}{_{3D}^{ {text{wie}}-1}}left( { { {mathbf{W}}_{S_{ {x_R}}^{ {text{wie}}}}} circ mathcal{T}_{3D}^{ {text{wie}}}left( { { {mathbf{Z}}_{S_{ {x_R}}^{ {text{wie}}}}}} right)} right),quad \ where;{ {mathbf{W}}_{S_{ {x_R}}^{ {text{wie}}}}} = frac{ { { {left| {mathcal{T}_{3D}^{ {text{wie}}}left( {widehat {mathbf{Y}}_{S_{ {x_R}}^{ {text{wie}}}}^{ {text{basic}}}} right)} right|}^2}}}{ { { {left| {mathcal{T}_{3D}^{ {text{wie}}}left( {widehat {mathbf{Y}}_{S_{ {x_R}}^{ {text{wie}}}}^{ {text{basic}}}} right)} right|}^2} {sigma ^2}}}. \ end{gathered} tag{3.17} Y SxRwiewie=T3Dwie−1(WSxRwie∘T3Dwie(ZSxRwie)),whereWSxRwie=∣∣∣T3Dwie(Y SxRwiebasic)∣∣∣2 σ2∣∣∣T3Dwie(Y SxRwiebasic)∣∣∣2.(3.17)
其中 ∣ ⋅ ∣ left| cdot right| ∣⋅∣ 代表复数的模，不过这里一般使用的是实数变换。注意虽然式中使用了矩阵形式的维纳收缩系数，但所有的运算 element-wise 的。对于低频的分量，由于相似块使得图像原本的能量更加集中，而噪声的强度始终不变，所以其收缩系数接近 1；对于高频的分量，图像原本的能量非常小，所以这时的收缩系数接近 0；而对于那些中间频率的分量，由于图像本身存在一些较高频的信息，相似块之间也存在一定的差异，所以这部分的能量较低，并且与噪声的区别并不是很大，这时通过维纳收缩既可以抑制部分噪声的能量，也可以尽保留图像本身的细节信息。所以，相比于 Step1 的硬阈值滤波，Step2 能够抑制所有频段的噪声，同时也不会完全地抹杀掉高频的信息，从而获得更好的降噪效果。
类似于 Step1，对图像中每一个参考块都进行维纳滤波后，我们同样需要对这些结果进行整合。在 Step1 中，我们根据每个 3D 组合经过硬阈值滤波后剩下的非零系数个数来分配该组合滤波结果的权值。如果非零的系数越多，证明其稀疏性越差，即包含更多的高频信息，而残留的噪声能量也就越大，所以应该赋予更低的权值。同样，因为无噪声图像的变换系数中只有很少的一部分明显大于零，所以式 (3.17) 的收缩系数矩阵大部分也应该为零。假设 3D 组合内所有像素都是独立的，那么残留的噪声与收缩系数矩阵的二阶范数成正比，所以类似于式 (3.9)，定义该组合的权值为
w x R wie = σ − 2 ∥ W S x R wie ∥ 2 − 2 . (3.18) w_{ {x_R}}^{ {text{wie}}} = {sigma ^{ – 2}}left| { { {mathbf{W}}_{S_{ {x_R}}^{ {text{wie}}}}}} right|_2^{ – 2}.tag{3.18} wxRwie=σ−2∥∥∥WSxRwie∥∥∥2−2.(3.18)
同样地，我们可以添加一个 Kaiser 窗来使得块中心的像素具有更高的权重，从而减小边界效应。那么，类似与式 (3.10)，Step2 所得的最终估计图像可表示为
y ^ final ( x ) = ∑ x R ∈ X ∑ x m ∈ S x R wie w x R wie Y ^ x m wie, x R ( x ) ∑ x R ∈ X ∑ x m ∈ S x R wie w x R wie X x m ( x ) , ∀ x ∈ X . (3.19) {hat y^{ {text{final}}}}(x) = frac{ {sumlimits_{ {x_R} in X} {sumlimits_{ {x_m} in S_{ {x_R}}^{ {text{wie}}}} {w_{ {x_R}}^{ {text{wie}}}widehat Y_{ {x_m}}^{ {text{wie, }}{x_R}}(x)} } }}{ {sumlimits_{ {x_R} in X} {sumlimits_{ {x_m} in S_{ {x_R}}^{ {text{wie}}}} {w_{ {x_R}}^{ {text{wie}}}{mathcal{X}_{ {x_m}}}(x)} } }},;forall x in X.tag{3.19} y^final(x)=xR∈X∑xm∈SxRwie∑wxRwieXxm(x)xR∈X∑xm∈SxRwie∑wxRwieY xmwie, xR(x),∀x∈X.(3.19)
虽然 Step2 和 Step1 一样都是直接对原始的有噪声图像进行降噪的，但是 Step2 借助了 Step1 所得的基础估计图像，一方面能够更加准确地找到相似块的位置，另一方面可作为原始无噪声图像的经验估计进行均方误差最优的维纳滤波，实验表明这种方法确实能够提高接近 0.5~1dB 的 PSNR 质量（论文数据），可以说是比较有效的。但要注意的是，因为 Step2 和 Step1 的整个流程基本相似，所以其所需计算时间也是差不多的，至于是否值得为了以上提升质量而多花一倍的时间需要根据实际的需求进行讨论。得益于巧妙的过完备像素整合方法，BM3D 的降噪质量是要明显优于 NLM 算法的，在计算复杂度方面需要根据具体的配置而论，因为 NLM 需要为每个像素的邻域寻找相似块，而 BM3D 是以块为单位的，块与块之间存在一定的距离，所以在寻找相似块上花费的时间要节省不少。
3.3 高效的实现方法
前面我们主要介绍了 BM3D 的基本原理，可以看到其复杂度还是比较高的，所以为了保证其实用性，我们有必要设计一种比较高效的实现方法，其中有一些已经在前面的讨论中有所提及。总来说其主要包括以下几个部分，注意因为 Step1 和 Step2 的流程基本相同，所以除非特别指出，以下参数的意义对于两个 Step 是一样的，但具体的值可能不同：
1) 减少需要处理的参考块数量。实际上，根据前面的算法描述，其大部分的时间都花在相似块的匹配上。不同于 NLM 算法将图像每一个像素的邻域都作为一个参考块来寻找相似块，BM3D 算法从图像左上角第一个 N 1 × N 1 N_1times N_1 N1×N1 的块开始，采用栅格扫描的方法，每次水平平移 N s t e p N_{step} Nstep 个像素，直至到达图像右侧再回到左侧第一个块垂直平移 N s t e p N_{step} Nstep 个像素，直到铺满整幅图像。这样，每个块都会作为一次参考块以保证每个像素都至少参与一次协同滤波，图像中参考块的数量被缩减到 W × H / N s t e p 2 W times H/N_{step}^2 W×H/Nstep2。
2) 降低建立 3D 组合的复杂度。虽然前面已经减少了参考块的数量，但是进行全图全搜索的复杂度还是很高的，所以我们限制最大的搜索窗口为 N S × N S N_Stimes N_S NS×NS，其中心为参考块的左上角坐标 x R x_R xR。更进一步，我们还可以采用预测搜索（Predictive Search）的方法，即搜索窗口不再是矩形，如果当前位置的候选块属于相似块，那么其全搜索区域将增加一个大小为 N P R × N P R N_{PR}times N_{PR} NPR×NPR（ N P R < < N S N_{PR}<<N_S NPR<<NS）的窗口，其中心为当前候选块位置沿着参考块处理的方向（例如栅格扫描方向）平移 N s t e p N_{step} Nstep 个像素。注意，总的搜索区域还是在 N S × N S N_Stimes N_S NS×NS 的范围内。当然，预测搜索的方法可能会带来一定的性能下降，所以我们可以设置一个步长 N F S N_{FS} NFS，每 N F S N_{FS} NFS 个参考块中的最后一个参考块都会进行一次区域为 N S × N S N_Stimes N_S NS×NS 的全搜索，从而达到一个折中的目的。为了减少后续协同滤波的计算复杂度，我们还可以设置保留的最大相似块数量 N 2 N_2 N2，注意其一般是 2 的幂以满足一些变换的长度要求。
3) 降低正交变换的复杂度。所有的包括 Step1 和 Step2 所使用的 3D 协同变换都为可分离的正交变换，这样该 3D 变换就可分解为各个相似块独自的 2D 变换以及相似块之间的 1D 变换。正如前面所提到的，对于 Step1 每个参考块 N S × N S N_Stimes N_S NS×NS 的邻域内的所有候选块的 2D 变换都会先缓存起来，一方面可用于式 (3.5) 中的块匹配误差计算，也可以直接用于后面的 3D 协同变换。因为相邻参考块的搜索邻域有很大的重叠，所以上一个参考块的候选块的 2D 变换结果可以给下一个参考块复用，从而每个候选块只需要进行一次 2D 变换。对于 Step2 因为我们直接在空间域上计算块距离，所以没必要对每个候选块都进行变换，但当前参考块的相似块的 2D 变换结果也可给相邻的参考块复用，如果其刚好也是相邻参考块的相似块。
4) 高效的像素整合方法。对于式 (3.10) 和 (3.19) 的像素整合，我们没必要把所有参考块的协同滤波结果都存储起来，因为这样需要很大的存储空间。实际上，我们可以用两个与图像同样尺寸的矩阵来分别存储分子和分母，其中分子为协同滤波结果与加权系数的乘积的累积，分母为加权系数的累积，最后我们只需把两个矩阵按元素相除即可得到最终整合的图像。
5) 降低边界效应。某些 2D 变换（如 DCT 等）可能会导致边界效应的产生，如块效应等等。为了减小边界效应的影响，BM3D 算法引入了一个 Kaiser 窗，具体定义如式 (3.9) 所示，即赋予所得的 3D 协同滤波结果的每一个 2D 块中心像素更高的权重，而边界像素的权重则更低。
论文根据应用环境的不同给出了三种基本配置，包括快速配置，以及标准配置中的低噪声和高噪声配置，具体参数如表 1 所示。快速配置为了加快算法的运行速度，在参考块数量、搜索窗口、相似块数量等方面做了一些妥协，并且使用了预测搜索的方法以减少相似块匹配的时间，可以看到相比于标准配置能够提升 5~8 倍的速度，当然图像质量也会有所下降。标准配置不使用预测搜索，并且也增大了全搜索窗口的大小。对于噪声比较强的情况，配置也相应地增大了参考块的大小，同时使用了变换域中的距离来计算两个块之间的误差，从而提高相似块的匹配准确度。
表1 BM3D基本配置参数
为了测试不同的变换类型对最终图像质量的影响，论文也做了一些实验，结果如表 2 所示。注意，所有的测试都基于表 1 的标准配置参数，只有表 2 中相应列的变换类型被替换，例如 T 2 D ht mathcal{T}_{2D}^{ {text{ht}}} T2Dht 为 Haar 变换表示表 1 中 Step1 的 2D 变换类型被替换为 Haar，但 Step2 以及其他的变换类型不变。噪声的标准差为 σ = 25 sigma=25 σ=25 。另外，为了更好地展示 Step2 所带来的质量提升，表 2 中每个样例（Boats 和 Lena）的第 1 列为 Step1 所得的基础估计图像的 PSNR，而剩余两列都为 Step2 所得最终估计图像的 PSNR。具体的变换类型原理这里不详细介绍。可以看到，使用不同的变换类型所带来的质量变化其实并不是很大，对于 Step2 所得的最终输出大多数都在 0.1dB 的变化范围内，除了 DC-only，也就是在第 3 维上进行简单的求和叠加。但是，DC rand 表明，直接在第 3 维上求和叠加，然后加上一些随机数，同样可以达到使用正交变换的效果，这也表明这些相似块的相似性确实能够得到比较好的保证，但依然会存在一些细节上差异，使得简单的叠加并不能很好地恢复图像。除此以外，对比第 1 列和其他两列，我们也可以看到 Step2 的维纳滤波确实可以带来明显的质量提升，在时间允许的情况下是十分值得的。
表2 不同的变换类型对结果的影响
3.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 块内不同的帧进行补偿。除此以外，还有一些比如自适应的块大小以及变换类型等等的改进，总之，这些扩展的基本原理都是相似的，即寻找相似块，然后进行协同滤波，最后把同一个像素的不同滤波结果整合。这里不再详述。
代码实现
这里是本人关于BM3D算法的C 实现，https://github.com/ZoengMingWong/BM3D_step1_cpp，目前只做了 Step1，不过效果也是不错的。里面主要采用了OOP的编程风格，整体阅读和理解也会相对轻松。
发布者：全栈程序员栈长，转载请注明出处：https://javaforall.cn/133964.html原文链接：https://javaforall.cn
