BM3D

This doc introduces the process of bm3d algorithm in detail.

BM3D Algorithm (centered)

1 总览

BM3D是2007年TIP的文章，题目是Image denoising by sparse 3D transform-domain collaborative filtering，论文、项目的地址是https://webpages.tuni.fi/foi/GCF-BM3D/，提供python、matlab代码；An Analysis and Implementation of the BM3D Image Denoising Method为C++实现，这篇论文里面有对BM3D实现的更为详细的说明。

处理灰度图的BM3D以及它的变体CBM3D（彩色图）、VBM3D（时域）是图像去噪领域公认的去噪效果（PSNR）最好的，而BM4D、VBM4D等也都是沿袭BM3D的基于块处理（block-wise estimate）的思想，但其计算时间复杂度极大，或许只能用于离线处理（offline），当然后续有文章进行优化（代码、算法），这里就不再提及。

BM3D算是NLM(non-local mean)的升级版本，因为它主要用到了非局部块匹配的思想，首先找相似块，不同于传统NLM使用L2距离，它用了硬阈值线性变换降低了L2距离的复杂度；找到相似块后，NLM是做一个均值处理，而BM3D则是将相似块域转换，提出Collaborative filtering降低相似块自身含有的噪声（NLM做均值，引入了相似块的噪声），并在aggregation处对相似块加权处理，得到降噪后的目标块。

BM3D是一种基于增强变换域中稀疏表示的新型图像去噪策略。通过将图像的类似2D片段分组为3D数据组来实现稀疏性的增强。协同过滤是为处理这些3D组而开发的一种特殊程序。它包括三个连续的步骤：组的3D变换，变换谱的收缩和逆3D变换。借此，我们获得由一组联合滤波的2D片段组成的组的3D估计。由于分组块之间的相似性，变换可以实现真实信号的高度稀疏表示，从而可以通过收缩很好地分离噪声。通过这种方式，协同过滤甚至可以显示分组片段共享的嵌套细节，同时保留每个片段的基本独特特征。

BM3D的复杂之处在于，首先找相似块，其次是采用了两次block-wise estimate，复杂度相比NLM翻了1倍，且含有域变换操作。

2 算法说明

2.1 算法架构

我们将首先描述如何处理灰度图像。对彩色图像的扩展将在后面进行解释。在下面所有的工作中，我们工作在高斯白噪声的情况下，其中方差用σ2表示。该算法可分为两个主要步骤：1。第一步是在协同滤波过程中使用硬阈值法估计去噪图像。这一步中的参数用硬指数表示；2。第二步既基于原始的噪声图像，又基于第一步得到的基本估计。它使用维纳滤波。因此，第二步用指数维纳来表示。该算法将总结在接下来的两个图1和图2中。图3显示了补丁，一个以参考补丁P为中心的搜索窗口，并说明了导致多个估计的补丁重叠。

2.2 流程

分为两步，Basic estimate与Final estimate，两步除Collaborative filtering不同外，其余子步骤近似。其中：

• S1 Basic estimate:

  • S11 Grouping 为目标块找到$N^{hard}$个相似块；

  • S12 Collaborative hard-thresholding 使用域变换方法，得到“变换后的相似块集合”，采用协同（由相似块共同确定）硬阈值策略“弱化”相似块的噪声，后反变换回原始块域；

  • S13 Aggregation 加权平均相似块，叠加后得到basic estimate后的目标块；

• S2 Final estimate:

  • S21 Grouping 使用S1处理后的图像重新计算块集合，并以相同的顺序相乘噪声图像集合

  • S22 Collaborative Wiener filtering 对两个集合都应用三维变换。利用基本估计的能谱作为真实（导频）能谱，对噪声能谱进行维纳滤波。通过对过滤后的值应用逆三维变换来生成所有分组块的估计值，并将块的估计值返回到其原始位置。

  • S23 Aggregation 加权平均相似块

协同过滤(Collaborative filtering)是BM3D分组和过滤过程的名称。实现四个步骤：1）找到与给定图像块相似的图像块，并将其分组为三维块2）三维线性变换；3）变换谱系数收缩；4）三维反变换。因此，这个3D滤波器可以同时过滤掉3D块中的所有2D图像补丁。

2.21 S11 基础估计 - Grouping 相似块分组

在有噪声的图像中对于每个目标块通过块匹配来实现分组，在附近寻找最多$N^{hard}$个相似的图块，为了避免噪点的影响，将图块经过2D变换后再用欧氏距离衡量相似程度。按距离从小到大排序后取最多前$N^{hard}$个且两个块的距离小于固定阈值，这些块被认为相似并堆叠成一个三维数组。

具体定义及公式：

用P表示大小为 $k^{hard}\times k^{hard}$ 的参照块，以P为中心的 $n^{hard}\times n^{hard}$ 邻域中搜索原始的噪声图像。相似块的集合被简单地定义为 \(\mathcal{P}(P)=\{ Q:d(P,Q)\le \tau^{hard}_{match} \}\) 其中，

• $\tau^{hard}_{match}$ 为d的最大距离阈值，两个块被认为相似的程度；该参数的选择是基于对理想差值的可接受值的确定性推测，主要忽略了信号的噪声分量。

• $d(P,Q)=\frac{

  \gamma’(\tau_{2D}^{hard}(P))-\gamma’(\tau_{2D}^{hard}(Q)

  _2^2}{(k^{hard})^2}$ ，$

  \cdot

  _2$是L2norm。d是块之间的归一化二次距离，用以描述当前块和参考块的相似度；为了避免相对大的方差或小的块导致的错误集合，使用粗预滤波来测量块距离，这种预滤波是通过在两个块上应用归一化的二维线性变换，然后对得到的系数进行硬阈值化。

• $\gamma’$ 是一个阈值为$\lambda^{hard}{2D}\sigma$的硬阈值算子。对于$\sigma≤40$，$\lambda^{hard}{2D}\sigma=0$，它简单地将绝对值低于阈值$\lambda^{hard}_{2D}\sigma$的补丁系数设为零；
• $\sigma^2$ 是零均值高斯噪声的方差。

通过在两个块上应用归一化的二维线性变换，然后对得到的系数进行硬阈值化，

2.22 S12 基础估计 - 3D协同硬阈值滤波

在三维变换域中，通过硬阈值法实现了协同滤波。通过硬阈值化进行有效的噪声衰减，然后进行反变换，产生一个块级估计的三维阵列。

具体定义及公式：

先对每个三维矩阵中的二维的块进行二维变换（DCT变换等）后再在矩阵的第三个维度进行一维变换（通常为Hadamard Transform）。变换后采用硬阈值的方式将小于超参数 $λ^{hard}_{3D}\sigma$ 的成分置为0。同时统计非零成分的数量作为后续权重的参考。最后通过在第三维的一维反变换和二维反变换得到处理后的图像块。这一过程的公式如下： \(\mathbb{P}(P)^{hard}=\tau_{3D}^{hard^{-1}}(\gamma(\tau_{3D}^{hard}(\mathbb{P}(P))))\) 其中 $\gamma$ 为硬阈值操作：

$$ \begin{split} \gamma(x)= \left { \begin{array}{ll} 0 & if\space |x|\le \lambda^{hard}_{3D}\sigma\

x                                 & othersise \end{array} \right. \end{split} $$ $\tau_{3D}^{hard}$ 为归一化的3D线性变换（normalized 3D linear transform），降低相似块的噪声

为了实际目的，P的三维变换$\tau_{3D}^{hard}$是由两个转换：2d转换表示$\tau^{hard}{2D}$应用在集合$\mathcal{P}(P)$的每个块，和1维转换表示$\tau^{hard}{1D}$应用三维。

传统方法，如NLM，由空域得到近似块，然后对近似块的每个像素一一对应去平均，作为目标块每个像素的值。但是，上述策略对于如下场景并不合适：

i. 某些相似块拥有的噪声更小，相比其它相似块，该块的“权重”应更大，而不是简单取平均

ii. 相似块图像信息冗余，从空域上看，两个有重叠区域的相似块，简单平均会造成目标块信息重复。

因此采用“Collaborative filtering by shrinkage in transform domain”的方式，能够加强相似块的稀疏性，同时降低相似块的噪声。

2.23 S13 基础估计 - 聚合

通过之前的滤波，图像中的大部分噪声被滤除，接下来打算将滤波后的图像块恢复到原来的位置。因为组之间的图像块之间存在重叠相交的情况，即过完备。如何复原对去噪的结果也有很大的影响。因此在存在相交的位置，采用权值平均的方法来进行聚合复原，即对于同一个像素点的不同估计值赋予它们一个权值，取加权平均结果作为最终的结果。

具体定义及公式：

当协同过滤完成时，得到每个块的估计值，然后得到每个像素的可变数量的估计值。这些估计值保存在一个缓冲区中： $$ \forall{Q}\in P,\forall x \in Q,\begin{split} { \begin{array}{ll} v(x)=v(x)+w_P^{hard}u_{Q,P}^{hard} (x) \

\delta(x)=\delta(x)+w_P^{hard}  \end{array}

\end{split} $$ 其中，

• v(resp.$\delta$) 是分子(分母)在分组步骤结束时获得的图像的基本估计的一部分；
• $u_{Q,P}^{hard}(x)$ 为在参考patch P的协同滤波过程中得到的属于patch Q的像素x值的估计值；
• $w_P^{hard}=\begin{split}{\begin{array}{ll}(\sigma^2N_P^{hard})^{-1}&if \space N_P^{hard}\ge1\1&otherwise\end{array}\end{split}$ , 噪声更大的块级估计应该被授予更小的权重。
• $N_P^{hard}$ 为硬阈值后三维块中非零数的数量：

为了降低相似块的信息冗余与自身噪声。加权平均后得到目标块的像素值（在实际搜索相似块时，需要对原图padding，由padding得到的相似块不纳入加权平均中）。

使用Kaiser窗口（参数$\beta$）作为聚合中的权重的一部分，以减少当使用某些二维变换（例如，二维DCT、二维DFT或周期小波）时可能出现的边界效应。

2.24 S21 最终估计 - 块匹配

经过第一阶段的基本估计后已经极大地消除了噪点，得到一幅基本估计滤波后的图像$u^{basic}$。对基础估计图的块进行类似步骤S11的分组，得到新的三维数组$\mathbb{P}^{basic}(P)$，然后在原始含噪图像的相同位置按照基本滤波后的分组，将原始含噪图像也进行相同的分组$\mathbb{P}(P)$。

基础估计图的相似块的集合被定义为 \(\mathcal{P}^{basic}(P)=\{Q:d(P,Q)<\tau^{wien}_{match}\}\) 其中，

• $\tau^{wien}_{match}$ 为d的最大距离阈值

• $d(P,Q)=\frac{

  P-Q

  _2^2}{(k^{wien})^2}$

即形成了两个新的3D组：

• $\mathbb{P}^{basic}(P)$ 通过从基本估计$u^{basic}$中将块堆叠在一起
• $\mathbb{P}(P)$ 通过从原始的噪声图像中以相同的顺序将块堆叠在一起。

为了优化，在两个三维组中保留最多 $N^{wien}$ 个块。

2.25 S22 最终估计 - 3D协同维纳滤波

将 $\mathbb{P}^{basic}(P)$ 作为指导图像块，当作未经过噪声污染的原始图像，$\mathbb{P}(P)$ 当作要去噪的图像，对两组分别进行3D变换，然后进行维纳滤波，维纳滤波后，对$\mathbb{P}(P)$组进行3D反变换到空间域图像块，并将图像块的图像块放回其原始位置。

也就是，对含基础估计3D数组的第三维，即图块叠起来后，每个图块同一个位置的像素点构成的数组，进行DCT变换，利用如下公式得到系数。两个三维矩阵都进行二维和一维变换，用维纳滤波（Wiener Filtering）将噪声图形成的三维矩阵进行系数收缩，该系数$w_P$通过基础估计的三维矩阵的值以及噪声强度得出。这一过程用如下公式表示： \(\mathbb{P}^{wien}(P)=\tau_{3D}^{wien^{-1}}(w_P\cdot\tau_{3D}^{wien}(\mathbb{P}(P)))\)
\(w_P(\xi)=\frac{|\tau_{3D}^{wien}(\mathbb{P}^{basic}(P))(\xi)|^2}{|\tau_{3D}^{wien}(\mathbb{P}^{basic}(P))(\xi)|^2+\sigma^2}\)

有效的协同过滤可以通过在变换域中的收缩来实现。假设已经形成了类似信号片段的d + 1维组，则协同收缩包括以下步骤。

• 将d + 1维线性变换应用于该组。
• 收缩（例如通过软阈值和硬阈值或维纳滤波）变换系数以衰减噪声。
• 反转线性变换以生成所有分组片段的估计值

这些组的特征在于：

• 片段内相关性出现在每个分组片段的像素之间——自然图像的独特性
• 片段间相关性出现在不同片段的相应像素之间——分组片段之间的相似性

3D变换可以利用两种相关性，从而在组中生成真实信号的稀疏表示。 这种稀疏性使收缩非常有效地衰减了噪声，同时保留了信号的特征。

2.26 S23 最终估计 - 聚合

与第一步中一样，将系数与含噪3D图块相乘放回原处，最后做加权即可得到最终估计图。相对于基础估计图，最终估计图还原了更多原图的细节。

$$ \forall{Q}\in P,\forall x \in Q,\begin{split} { \begin{array}{ll} v(x)=v(x)+w_P^{wien}u_{Q,P}^{wien} (x) \

\delta(x)=\delta(x)+w_P^{wien}  \end{array}

\end{split} $$ 其中，

• v(resp.$\delta$) 是分子(分母)在上一步结束时得到的图像的最终估计的一部分；
• $u_{Q,P}^{wien}(x)$ 为在参考块P的协同滤波后得到的属于块Q的像素x处的估计值；

• $w_P^{wien}=\sigma^{-2}

  w_P

  ^{-2}_2$

采用$k^{wien}\times k^{wien}$的Kaiser窗口来减少边界效应。

则最终估计为 \(u^{final}(x)=\frac{\sum_Pw^{wien}_P\sum_{Q\in\mathcal{P}(P)}\mathcal{X}_{Q}(x)u_{Q,P}^{wien}(x)}{\sum_Pw^{wien}_P\sum_{Q\in\mathcal{P}(P)}\mathcal{X}_{Q}(x)}\)

2.3 参数

噪声没那么严重时，λ2d=0，窗口大小N1=8，严重时为λ2d=2，N1=12

3 结论

仔细分析BM3D的两个阶段滤波细节可知，第一个阶段对图像进行块匹配，3D协同滤波得到的图像，并没有直接贡献到最终的图像结果里，而是滤波后的结果噪声相对原始噪声更少，块匹配的效果更准，从而用来指导第二阶段的块匹配和3D协同滤波中的维纳滤波系统函数的权重计算，得到权重后真正参与到协同滤波的像素点还是原始的含噪图像，3D维纳滤波后反变换回到空间域后再对图像块进行加权，得到最终去噪结果，因此真正的去噪过程是第二阶段，去噪涉及到的技术有变换域的维纳滤波和空间域的线性加权。

3.1加速

在实际操作中，为加快BM3D的计算速度，在寻找相似块的步骤后，得到的块实际上已经进行了2D变换处理，然后再加上一个1D变换（文中使用1D-Haar离散小波变换），成为3D变换，使用2D+1D的变换方法替代直接3D变换。

1)在进行块匹配过程中，不是滑窗逐像素进行匹配，而是以一定步长进行匹配，可以减少计算量

2)降低建立3D组合的复杂度。构成3D组的时候，限制最大的匹配块数目 搜索匹配块时，限制最大的搜索区域

3)块匹配时，使用predictive search（具体可能要看下代码如何实现）

4)匹配时已经得到了2D变换，构成3D组时，再做1D变换就得到3D变换，不需要直接做3D变换

5)两个步骤中，进行hard_thresholding和Wiener_filtering时顺便把权重计算好，存在buffer里，最后Aggregation时就可以直接算

6)使用Kaiser window来减少边界效应

每个像素的复杂度（计算量很大)

3.2难点

文中提到的2D变换与各种超参数，并没有一个确定值，对于真实视频去噪，使用的2D变换与超参数可能与文章实验设置不同，因此需要进行微调，也就是需要大量实验的积累。

最终的去噪结果可以发现，噪声被很好地去除了，图像边缘保留完整，图像纹理得到了很好的还原。

4 CBM3D

处理彩色图像时，将图像的RGB色彩空间转换为YUV色彩空间，Y通道进行匹配，匹配的结果直接在UV分量上使用。因为YUV的Y分量拥有较其余分量更多的图像信息（边缘、材质、纹理等），并且拥有更高的SNR（信噪比），而U、V分量拥有更多的低频信息，人眼对这些色度的变化相对不会那么敏感。因此对于C-BM3D，使用Y分量搜寻相似块，U、V分量使用Y分量的相似块位置信息，后续滤波和聚合等过程则分通道进行。 \(A_{YCbCr}=[\begin{matrix}0.30 & 0.59 &0.11\\-0.17&-0.33&0.50\\0.5&-0.42&-0.08\end{matrix}]\\ A_{opp}=[\begin{matrix}\frac{1}{3} & \frac{1}{3} &\frac{1}{3}\\\frac{1}{\sqrt{6}}&0&\frac{-1}{\sqrt{6}}\\\frac{1}{3\sqrt{2}}&\frac{-\sqrt{2}}{3}&\frac{1}{3\sqrt{2}}\end{matrix}]\)

5 VBM3D

对于视频去噪，一是可以将视频转化为单帧图像，然后使用图像去噪算法对单帧进行处理，然后融合还原成已去噪的视频；二是根据视频前后帧信息，某个像素点，前帧没有噪声，或噪声较少，那么就可以作为后帧的值，但视频中物体是运动的，如果按前后帧的同一位置的像素处理是不合理的，因此会引入运动补偿、跟踪的算法，对于实时处理来说，条件就有点苛刻。

VBM3D不含运动补偿，对中间帧的目标块搜索相似块，搜索对象是前后帧与中间帧，提出了predictive-search block-matching（PS-BM），用于前后帧的相似块搜索。具体而言，PS-BM，先以中间帧的目标块为中心、搜索半径NS的区域寻找相似块，然后在该块对应的前后帧的位置为中心、搜索半径NPR的区域寻找相似块，中间帧、前后帧的相似块构成块集合。其余步骤与BM3D无异。

6 Kaiser window

一种局部优化的具有较强能力的窗函数。

它是凯塞(Kaiser)采用修正的零阶贝塞尔函数实现的。已经证明，凯塞窗是十分接近最大能量比的窗函数.用它来设计有限单位脉冲响应(FIR)滤波器可以有效地调整窗函数频谱主瓣与旁瓣的相对关系，从而使设计的FIR滤波器能够同时满足通带波动、阻带衰耗和过渡带宽等多项设计指标的要求。在数字信号处理中，其定义如下 $$ \begin{split} w[n]=\left { \begin{array}{ll} \frac{I_0(\pi\alpha\sqrt{1-(\frac{2n}{N-1}-1)^2})}{I_0(\beta)} & 0\le n\le N-1\

0                                 & otherwise \end{array} \right. \end{split} $$ 其中：

• N 为序列的长度
• I0 是零阶的第一类修正贝索函数
• α 是任意非负实数，用来调整凯泽窗的外形。在频域上可以在主瓣（main-lobe）宽度及旁瓣（side lobe）大小中取拾，这是窗函数设计的重要考量因素。

若N为奇数，窗函数最大值会在 $w[(N-1)/2]=1$ ,若N为偶数，窗函数最大值会在 $w[N/2-1]=w[N/2]<1$ 。

凯塞窗函数的表示式为 $$ \begin{split} w(n)=\left { \begin{array}{ll} \frac{I_0(\alpha)}{I_0(\beta)} & |n|\le \frac{N-1}{2}\sigma\

0                                 & |n| > \frac{N-1}{2} \end{array} \right. \end{split} $$

7 维纳滤波原理

维纳滤波器是采用统计的方法对平稳信号进行滤波的一种方法，其基本思想是设计一个滤波器使得信号经过滤波器后的输出信号和原始信号误差在统计意义上最小，其流程可描述为信号+噪声产生含噪信号，经过滤波器后后输出新信号。

在信号的采集或者传输的过程中，信道里存在噪声，比如拍照过程，传感器中本身存在热噪声，同时在光子计数过程中也存在着泊松噪声，这样我们的ISP处理器得到的信号就是一个含噪信号，ISP去噪的过程就是为了从含有噪声的CMOS RAW数据中恢复出原始不含噪的照片。维纳滤波器想要达到的目的是经过滤波其得到的信号y和原始信号在统计意义上误差最小，即均方误差最小（MMSE）,可用如下公式表示：$min E(e^2)= min E((y-s)^2)$

图像信号通常是先采集然后处理的信号，因此可以用非因果的维纳滤波器，其系统函数在频率域的表达式为： $H(z)=\frac{P_{xy}(z)}{P_{ss}(z)+P_w(z)} $

这里 H(z) 为维纳滤波器系统函数的频域表达式， $P_{xy}(z)$ 为信号x与滤波后信号y的相关函数的功率谱， $P_{ss}(z)$ 为原始信号的功率谱， $P_w(z)$ 为噪声信号的功率谱。如果对统计信号处理不太了解，可以将维纳滤波器当作一个黑盒子，经过黑盒子后，信号便恢复到和原始信号误差统计意义上最小，也就达到我们去噪的目的了。

Donoho等人在1995年发表文章认为，==白噪声信号经过正交变换后，在每个正交基张成的子空间里仍然是白噪声==，而通常采用的正交变换基函数为归一化的基函数，其功率为单位功率，因此在正交变换域的某个基上的功率即为该系数的代数平方，基于此在正交变换域做维纳滤波的方法是点对点的系数收缩即可达到滤波的目的，而图像的噪声水平很容易采用统计的方法估计出来，我们记其方差为 $\sigma^2$，于是给出在正交变换域中的维纳滤波系统函数公式： $H = \frac{T_x^2}{T_x^2+\sigma^2}$

其中$T_x$为信号经过正交变换后的系数，比如傅立叶变换系数，离散余弦变换系数，小波变换系数等 在工程中，我们是没法拿到原始信号的，因此也无法拿到原始信号的正交变换系数放到上述维纳滤波器中进行滤波，最简单的方法便是将噪声信号当作原始信号进行正交变换带入维纳滤波公式得到滤波结果，但是这样的效果通常不怎么理想。

因此对上述维纳滤波进行改进，对含噪信号进行硬阈值滤波，得到的信号当作原始信号带入到维纳滤波系统函数中，然后对为经过任何操作的含噪信号进行维纳滤波，这样得到的结果更好。

有人会问，既然前面一部分已经指出，直接经过硬阈值滤波就可以得到很好的结果，为什么这里要做一次硬阈值滤波，然后放到维纳滤波器中进行维纳滤波呢？这是因为，在工程中想要找到区分信号和噪声的合适阈值很困难，很容易“误伤”小系数的真实信号，而维纳滤波里通过系统函数H可以分成以下三种情况：（1）在信号功率很小的时候噪声为主，那么H几乎等于零，H与含噪信号系数相乘，含噪信号的系数几乎等于零，很好抑制噪声信号；（2）在信号功率很大的时候信号为主，那么H几乎等于一，H与含噪信号系数相乘，信号几乎不变，信号得到了很好的保留；（3）在噪声功率和信号功率差不多的时候，这时候信号与噪声各占一半，H的数值约等于二分之一，H与含噪信号系数相乘，含噪系数近似恢复到原始信号系数；基于上面三种情况的分析，维纳滤波器可以减少信号的“误伤”，比直接进行硬阈值滤波会有更好的结果，另外维纳滤波是一种统计意义上误差最小的滤波方式。

Ref

Ref:

https://zhuanlan.zhihu.com/p/92973703

https://cloud.tencent.com/developer/article/2030568

https://cloud.tencent.com/developer/article/2091059

https://blog.csdn.net/qq_33552519/article/details/108632146

https://blog.csdn.net/zhognsc08/article/details/122770149

https://cloud.tencent.com/developer/article/2073115?from=article.detail.2030597

https://blog.csdn.net/haronchou/article/details/108867110 维纳滤波

An Analysis and Implementation of the BM3D Image Denoising Method

https://cloud.tencent.com/developer/article/1769748?from=article.detail.2030597 图像去噪综述

https://zhuanlan.zhihu.com/p/46399784 BM3D ：稀疏三维变换域协同过滤的图像去噪

原理&算法实现