通俗语言说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）,可用如下公式表示：
$$minE(e^2)=minE((y-s{)}^2)$$
图像信号通常是先采集然后处理的信号，因此可以用非因果的维纳滤波器，其系统函数在频率域的表达式为：
$$H(z)=\frac{P_{xy}(z)}{P_{ss}(z)+P_w(z)}$$
这里$H(z)$为维纳滤波器系统函数的频域表达式， P x y ( z ) P_{xy}(z)