通俗语言说BM3D

大家好，又见面了，我是你们的朋友全栈君。

随着友商某以摄像著称的旗舰机型的发布，其SOC中ISP5.0采用的所谓单反级降噪算法BM3D一下火热起来，本文试图用尽量通俗易懂的语言从算法原理的角度揭开BM3D算法的神秘面纱。

本文结构如下：
1.前言
2.硬阈值滤波原理介绍
3.维纳滤波原理介绍
4.BM3D原理详述
5.空间域降噪原理
6.再谈3D降噪（3DNR）

1.前言
图像去噪是计算机视觉前处理中很重要的一个环节，对于手机camera来讲，去噪的好坏直接影响最终图像的质量，图像去噪算法经历了传统的空间域去噪，基于傅立叶变换/离散余弦变换的频率域滤波降噪，基于变分法及模拟热对流的偏微分方程降噪方法，小波/多尺度几何变换（超小波）的多尺度降噪，以及现在非常火热的基于CNN（深度卷积神经网络）的AI降噪等阶段的演进。 BM3D降噪是芬兰的坦佩雷理工大学（Pampere University of Technology）的Kosadin、Alessandro、Vladimir、Karen等人2007年提出的基于传统方法的图像降噪算法，该方法的去噪性能目前是非AI图像降噪中去噪效果最好的，无愧于state-of-art denoising performance的称号。BM3D的基本思想来自于自然图像中本身有很多相似的重复结构的观察结果，采用图像块匹配的方式对这些相似的结构进行收集聚合，然后对其正交变换，得到它们的一个稀疏表示，充分利用稀疏性和结构相似性，进行滤波处理，BM3D的去噪能够充分保留图像的结构和细节，得到很好的信噪比。

在这里插入图片描述

自然图像自相似块结构

由于在BM3D降噪算法的两个很重要的步骤中分别用到了变换域硬阈值滤波和维纳滤波，因此在介绍BM3D之前，在本文的第2部分及第三部分分别简要介绍硬阈值滤波的基本原理和维纳滤波的基本原理。

2.硬阈值滤波原理介绍
硬阈值滤波是著名信号处理调和分析领域专家Donoho 1995年提出的在小波域对白噪声进行去噪的方法，它的基本假设是白噪声在小波的各个尺度中均匀分布，但是相对于主要信号的系数来说很小，于是可以通过一个硬阈值来将其区分开来，小于阈值的系数，将其置零，大于阈值的系数保持不变，通过这样的方法可以达到对信号进行去噪的目的，其基本流程如下图所示：

其具体过程描述如下：
（1）将信号从时间域或者空间域通过正交变换变换到变换域
（2）对变换域系数进行硬阈值滤波
（3）将硬阈值滤波后的变换域系数通过正交反变换还原到原来的时间域或者空间域得到去噪后的信号
这里的硬阈值滤波就是对变换域系数做一个简单的判断，如果其系数绝对值大于给定的阈值，那么保持系数不变，如果系数绝对值小于给点阈值，那么将系数置零。

这里的信号可以是一维信号，也可以是二维三维甚至是高维信号，其中本文要谈的图像就是二维信号，这里的正交变换可以是任意针对信号处理设计的正交变换，常用的傅立叶变换，离散余弦变换，小波变换，多尺度几何分析（超小波）等均在正交变换的范围内。

为了更形象理解上述过程，这里分别给出一维及二维信号的硬阈值滤波过程

close all
x = 0:0.02:3*pi;
y = 5 * cos(1 * x) + 8 * cos(3 * x)  + 7 * cos(5 * x) + 6 * cos(7 * x) + 5 * cos(10 * x);

y1 = y + 2.5 * randn(1, length(x));

freq_f1 = dct(y1);
freq_f  = dct(y);

hard_thresh = 50;

freq_signal_index = abs(freq_f1) > hard_thresh;
freq_f1_signal = freq_f1 .* freq_signal_index;

y1_hat = idct(freq_f1_signal);

y_hard = 50 * ones(1, length(x));


figure
subplot(211)
plot(x, y)
title('Original Signal')
subplot(212)
plot(x, y1)
title('Dirty Signal')

figure
subplot(511)
plot(x, y)
title('Original Signal')
subplot(512)
plot(x, y1)
title('Dirty Signal')
subplot(513)
plot(freq_f1)
title('DCT Coef of Dirty Signal')
hold on
plot(y_hard, 'r')
subplot(514)
plot(freq_f1_signal)
title('Dct Coef after Hard Shrinkage')
subplot(515)
plot(x, y1_hat)
title('Dirty Signal after Hard Shrinkage')

图中的红线即为给定的阈值，小于红线的系数全部改成零，然后反变换回来，可以将脏信号恢复到原始信号

从上到下的信号处理过程如下图所示

从图片可以看出，经过简单的硬阈值滤波，可以很好的过滤掉噪声信号，当然，这里的滤波结果太完美，实际上在噪声很大的时候，会残留一些噪声信号，但是信号的整体信噪比还是会有很大的提升的。

对于二维信号，这里主要指的是数字图像信号，其原理也类似，下面也给出一个直观的例子，

close all
x = linspace(0, 3*pi, 128);
y = linspace(0, 3*pi, 128);
[X, Y] = meshgrid(x, y);
Z = 6 + cos(X).*cos(Y) + 8 * cos(3*X).*cos(3*Y) + 9 * cos(5*X).*cos(5*Y) + 7 * cos(7*X).*cos(7*Y);
Z1 = 10.5 * randn(128, 128) + Z;

freq_Z = dct2(Z);
freq_Z1 = dct2(Z1);

hard_thresh = 45;
freq_Z1_mask = abs(freq_Z1) > hard_thresh;
freq_Z1_hat = freq_Z1 .* freq_Z1_mask;
Z1_hat = idct2(freq_Z1_hat);


%surf(Z1)

Z_img = uint8((Z+30)/60 * 255);
Z1_img = uint8((Z1+30)/60 * 255);
Z1_hat_img = uint8((Z1_hat + 30)/60 * 255);

imshow(Z_img)
title('Original 2D Signal')
figure, imshow(Z1_img)
title('2D Signal with Noise')
figure, imshow(Z1_hat_img)
title('2D Signal After Hard Shrinkage')
figure, mesh(freq_Z)
title('2D DCT Coef of Original 2D Signal')
figure, mesh(freq_Z1)
title('2D DCT Coef of Noisy 2D Signal')
figure, mesh(freq_Z1_hat)
title('2D DCT Coef after Hard Shrinkage')

figure
subplot(151)
imshow(Z_img)
title('Original 2D Signal')
subplot(152)
imshow(Z1_img)
title('2D Signal with Noise')
subplot(153)
surf(abs(freq_Z1))
title('2D DCT Coef of Noisy 2D Signal')
subplot(154)
surf(abs(freq_Z1_hat))
title('2D DCT Coef after Hard Shrinkage')
subplot(155)
imshow(Z1_hat_img)
title('2D Signal After Hard Shrinkage')

3.维纳滤波原理介绍
维纳滤波器是采用统计的方法对平稳信号进行滤波的一种方法，其基本思想是设计一个滤波器使得信号经过滤波器后的输出信号和原始信号误差在统计意义上最小，其流程如下图所示：

在信号的采集或者传输的过程中，信道里存在噪声，比如拍照过程，传感器中本身存在热噪声，同时在光子计数过程中也存在着泊松噪声，这样我们的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$ ，于是给出在正交变换域中的维纳滤波系统函数公式：
$\frac{T_x^2}{T_x^2+\sigma^2}$
其中 $T_x$ 为信号经过正交变换后的系数，比如傅立叶变换系数，离散余弦变换系数，小波变换系数等
在工程中，我们是没法拿到原始信号的，因此也无法拿到原始信号的正交变换系数放到上述维纳滤波器中进行滤波，最简单的方法便是将噪声信号当作原始信号进行正交变换带入维纳滤波公式得到滤波结果，但是这样的效果通常不怎么理想，下面给出这样的一组实验结果

从图像效果上来看，直接拿含噪信号当作原始信号带入到维纳滤波系统函数中进行维纳滤波，效果并不太理想，从2D-DCT系数上来看也是如此。因此对上述维纳滤波进行改进，对含噪信号进行硬阈值滤波，得到的信号当作原始信号带入到维纳滤波系统函数中，然后对为经过任何操作的含噪信号进行维纳滤波，这样得到的结果更好，下面给出这样的一组实验：

clear
clc
close all
x = linspace(0, 3*pi, 128);
y = linspace(0, 3*pi, 128);
[X, Y] = meshgrid(x, y);
Z = 6 + cos(X).*cos(Y) + 8 * cos(3*X).*cos(3*Y) + 9 * cos(5*X).*cos(5*Y) + 7 * cos(7*X).*cos(7*Y);
Z1 = 10.5 * randn(128, 128) + Z;

freq_Z = dct2(Z);
freq_Z1 = dct2(Z1);

sigma = 20.5^2;
hard_thresh = 35;

winer_filter_ideal =  freq_Z.^2 ./ (freq_Z.^2 + sigma);
winer_filter = freq_Z1.^2 ./ (freq_Z1.^2 + sigma);
freq_Z1_hat_ideal = winer_filter_ideal .* freq_Z1;
freq_Z1_hat = winer_filter .* freq_Z1;
Z1_hat = dct2(freq_Z1_hat);
Z1_hat_ideal = dct2(freq_Z1_hat_ideal);

freq_Z1_mask = abs(freq_Z1) > hard_thresh;
freq_Z1_hard = freq_Z1 .* freq_Z1_mask;
freq_Z1_hat_hard = freq_Z1_hard.^2 ./ (freq_Z1_hard.^2 + sigma) .* freq_Z1;
Z1_hat_hard = dct2(freq_Z1_hat_hard);


Z_img = uint8((Z+30)/60 * 255);
Z1_img = uint8((Z1+30)/60 * 255);
Z1_hat_img = uint8((Z1_hat + 30)/60 * 255);
Z1_hat_img_ideal = uint8((Z1_hat_ideal + 30)/60 * 255);
Z1_hat_img_hard = uint8((Z1_hat_hard + 30)/60 * 255);

figure, imshow(Z_img);
title('Original 2D Signal')
figure, imshow(Z1_img);
title('Noisy 2D Signal')
figure, imshow(Z1_hat_img);
title('2D Signal After Wiener Filtering')
figure, imshow(Z1_hat_img_ideal)
title('2D Signal After Hard Wiener Filtering')
figure, mesh(freq_Z);
title('2D DCT Coef of Noisy Signal')
figure, mesh(freq_Z1);
title('2D DCT Coef of Original Signal')
figure, mesh(freq_Z1_hat);
title('2D DCT Coef After Wiener Filtering')
figure, mesh(freq_Z1_hat_ideal)
title('2D DCT Coef After Hard Wiener Filtering')

figure, mesh(freq_Z1_hat_hard)
figure, imshow(Z1_hat_img_hard)

这次可以看到实验结果有了很明显的提升，无论从图像效果，还是2D-DCT频谱系数上来看均是如此。

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

4.BM3D原理详述
BM3D全称Block Matching and 3D collaborative filtering，基于块匹配的3D协同滤波，算法的中心思想是充分利用自然图像中丰富的自相似结构来进行图像降噪。

BM3D降噪算法过程分为两个阶段，第一个阶段称为基本估计滤波，第二个阶段称为最终估计滤波，其流程如下图所示：

第一阶段：
a)块匹配估计：i)块匹配分组 ii)3D协同硬阈值滤波
b)聚合加权滤波
第二阶段：
a)块匹配估计：i)块匹配分组 ii)3D协同维纳滤波
b)聚合加权滤波

下面按照流程里的数据流向，详细说明BM3D滤波过程

1)第一阶段块匹配：按照划窗的方式在图像中搜索与给定图像块的相似图像块，如果和图像的相似度高，那么将其加入到该分组中，如果相似度不高，舍弃，继续寻找下一个相似的图像块，需要特殊说明的是，图像的相似度与图像的差异成反相关关系，在相似度比较之前先对图像块做一个硬阈值收缩，去掉因为高频噪声带来的影响，块匹配的动态过程如下图所示

2)第一阶段3D协同硬阈值滤波：将上一个过程中得到的图像块进行3D的正交变换，在变换域中进行硬阈值滤波，然后反变换到空间域中，得到3D协同硬阈值滤波后的图像块，然后将滤波后的图像块分别放回其图像原始位置

这里的3D正交变换充分利用了图像块内部的相关性和图像块之间的相关性，在上图中第一二维空间组成的图像块内部，一个2D的正交变换揭示了图像内部的相关性，所有图像块经过第一个2D正交变换后，对应位置的像素点沿着第三维串起来，在第三维方向进行1D的正交变换，这个正交变换揭示了图像块之间的相关性，充分利用了同一组图像块相似来进行系数表示，通过这样的先2D正交变换再1D正交变换构成的3D正交变换硬阈值滤波，可以很好的抑制图像的噪声，同时保留了该组图像块中各个图像块自身的独特的细节。

3)第一阶段聚合加权滤波：
经过第2)中的3D协同滤波后，将图像块放回到其原始的位置，但是由于搜索相似图像块的时候，存在图像块重叠的情况，如上图中红色区域中的几个图像块，彼此之间有重叠区域，那么重叠区域里像素点会有不止一个去噪后的结果，到底取哪一个呢，处理这个问题的方法就叫做聚合加权滤波，即对于同一个像素点的不同估计值赋予它们一个权值，取加权平均结果作为最终的结果。

4)第二阶段块匹配
经过第一阶段的基本估计滤波后得到一幅基本估计滤波后的图像，再对其进行块匹配估计，得到新的分组，然后在未经去噪的原始含噪图像的相同位置按照基本滤波后的分组，将原始含噪图像也进行相同的分组，具体的块匹配过程与第一阶段块匹配过程相似，这里便不再赘述，需要说明的是，这里得到两个块匹配分组，一个分组的图像块来自第一阶段基本估计滤波后的图像，我们称之为a组，另一个分组来自没有经过任何处理的含噪图像，我们称之为b组。
5)第二阶3D协同段维纳滤波
在得到第二阶段块匹配的两个分组后，我们将a组作为指导图像块，当作为经过噪声污染的原始图像，b组图像当作要去噪的图像，对两组分别进行3D正交变换，然后按照文本第三部分介绍的那样，进行维纳滤波，维纳滤波后，对b组进行3D反正交变换到空间域图像块，并将图像块的图像块放回其原始位置。
6)第二阶段聚合加权滤波
同第一阶段聚合加权滤波一样，这里对b组的图像数据进行聚合加权滤波，得到BM3D的最终结果。

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

通过简单的扩展，就可以将针对灰度图像滤波的BM3D扩展到针对彩色图像滤波的C-BM3D,这里的C代表颜色color。

至此BM3D降噪算法就介绍完毕了。

5.空间域降噪原理详述
在文章的最后，笔者想啰嗦几句，介绍一下空间域像素之间进行加权平均能够达到降噪目的的原因
（1）从空间域加权平均本质上是一个低通滤波器，权重如何分配决定率滤波其的滤波带块形状等性能指标
（2）从随机变量的角度来看，噪声信号通常是假设为独立同分布的随机变量，即：
假设有n个像素点参与加权，它们的噪声分量为 $x(i)服从N(0,\sigma^2), i.i.d, i=1,2,\cdots,n$ ，权重为 $\omega_1,\omega_2,\cdots,\omega_n$ ,且 $\omega_1+\omega_2+\cdots+\omega_n=1$ ，那么 $x=\omega_1x(1)+\omega_2x(2)+\cdots+\omega_nx(n)$ 由方差求和公式很容易得到
$\begin{aligned} var(x) &=\omega_1^2\sigma^2+\omega_2^2\sigma^2+\cdots+\omega_n^2\sigma^2 \\ &=(\omega_1^2 + \omega_2^2 + \cdots + \omega_n^2)\sigma^2 \end{aligned}$
从上式可知经过加权后噪声的方差变成原来的 $\omega_1^2+\omega_2^2+\cdots+\omega_n^2$ ，且有 $\omega_1+\omega_2+\cdots+\omega_n=1$ 约束条件，这样的一个优化问题很容易得出，在某个 $\omega=1$ 时，其他均为零，此时只有一个块或者像素进行加权，它的噪声方差最大，为 $\sigma^2$ ,通过拉格朗日乘子发很容易分析得出，在 $\omega_1=\omega_2=\cdots=\omega_n=\frac{1}{n}$ 时n个图像块或者像素点进行加权，噪声方差最小，变成 $\frac{1}{n}\sigma^2$ ,由这个数学过程分析，在采用空间域图像去噪时，需要尽可能搜索到相似的图像块或者像素点进行加权，这样彼此的权重接近，得到的加权结果噪声方差最小，图像效果更佳。

6.再谈3D降噪（3DNR）
从概率论的独立重复实验说起，比如我们统计扔一个骰子点数为1的概率是多少，我们可以让一个人扔一万次，然后统计结果中点数为1的个数，从而计算出扔骰子点数为1的概率（大数定律意义下），另一种方法是让1000个人，每人扔十次骰子，然后统计点数为1的概率；从上面两个实验中可以看出，我们要得到一个随机事件的很多样本，可以让这个事件采用实验的方法，让随机事件重复N次，得到统计规律，也可以构造很多m个相同的随机实验，然后每个随机试验重复N/m次，这样从整体上来看，总体上也是相同的随机事件重复了N，得到的统计规律和前面的结论一致。
回到我们要说的图像去噪，我们知道含噪图像信号是由一个确定不变的图像信号和一个噪声随机变量的和，那么我们观测到的图像信号本身也是一个随机变量，我们想通过加权的方法得到真实的信号，那么我们就需要这个随机变量的N个样本。这样的样本从哪里得到了，一种方法就是对图像进行N次观测，这种方法就对应了我们在工程上采用的多帧相同曝光图像进行融合的3D降噪算法，这种方法类似于一个人做了N次独立重复实验；另一种方法就是在本身里面去寻找相同的图像块，它的一种实现方法就是本文中提到的BM3D，这种方法近似类比于N个人，每个人做了一次独立重复实验。