【图像去躁】基于BM3D实现图像去噪matlab源码

一、 总览

BM3D是2007年TIP的文章，题目是Image denoising by sparse 3D transform-domain collaborative ltering，论文、项目的地址是http\://www.cs.tut.fi/\~foi/GCF-BM3D/，提供matlab代码；http://www.ipol.im/pub/art/2012/l-bm3d/为C++的实现，这篇论文里面有对BM3D实现的更为详细的说明。

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

---

找相似块，得到块集合

---

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

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

二、 算法说明

1. 流程图

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

• S1 Basic estimate:

  • S11 Block-wise estimate

  • S111 Grouping 为目标块找到相似块，块集合无顺序限制；

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

  • S12 Aggregation 加权平均相似块(S112处理后)，叠加后得到basic estimate后的目标块；

• S2 Final estimate:

  • S21 Block-wise estimate

  • S211 Grouping 使用S111步的块集合，以及已由S1处理后的图像重新计算块集合

  • S212 Collaborative Wiener filtering 域变换后，使用维纳滤波

  • S22 Aggregation

2. S111 Grouping 找相似块

本文使用L2距离判断，采用a normalized 2D linear transform与hard-thresholding对块距离进行预处理操作，公式如下：\ \ 因为含有噪声的图像，计算L2距离，噪声使两个相似块的差异变大，可能找到错误的相似块。\ 其中，

• x为像素点，X为图像
• 目标块， 搜索块
• 是S1中的块大小
• 硬阈值操作，阈值设为
• 是归一化后的二维线性变换

根据距离就可以找到相似块集合，由式(5)所示，\ \ 其中，

• 确定是否相似的超参数
• 为 的相似块集合

3. S112 Collaborative hard-thresholding 相似块降噪

此处使用归一化的3D线性变换(normalized 3D linear transform)来降低相似块的噪声，然后使用反变换 得到S112步处理后的相似块 ，如下式(6)所示，\ \ 其中，

• 为归一化的3D线性变换， 为其反变换
• 硬阈值操作，阈值设为 
•  为处理后的相似块集合

a) 为什么要这么做

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

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

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

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

4. S12 Aggregation 加权生成目标块

同样是为了降低相似块的信息冗余与自身噪声，所以我认为与S112的目的一致。加权平均后得到目标块的像素值，如式(12)(10)所示：\

其中，

• 为权重
•  为公式(6)硬阈值操作后的非0系数的个数
• ，在实际搜索相似块时，需要对原图padding，由padding得到的相似块不纳入加权平均中

5. S211 Grouping

由S1初步处理后的图像，重新计算L2距离，得到相似块集合，此处不再使用硬阈值处理，公式如式(7)所示：

其中，

• 为S2步的判断是否相似的超参数
• 为相似块集合

6. S212 Collaborative Wiener filtering

由S211步得到的相似块集合 ，经过某个域变换后，使用维纳收缩系数(Wiener shrinkage coefcients) 加权，后经过反变换，得到块集合 。如式(9)、(8)所示。

\

7. S22 Aggregation

S12已经说明过了，S22与S12近似，公式如式(11)所示，而最终结果，与式(9)的计算步骤一致。

三、 结论

1. 加速

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

2. 难点

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

最终的去噪结果如下图所示，可以发现，噪声被很好地去除了，图像边缘保留完整，图像纹理得到了很好的还原。\ \ 但这些都是自行添加高斯白噪声产生的测试图像，若实际运用在图像降噪中，原始图像不会有这么多噪声，因此就不需要BM3D两步去噪。那么可以将BM3D的两步拆开，采用前步的硬阈值、2D变换寻找相似块、1D变换升至3D域再加权平均，或后步直接使用维纳滤波，或许就已经有很好的效果了。

下一步会总结BM4D、VBM4D的东西。

四、 C-BM3D

针对彩色图，本文将图像的RGB色彩空间转换为YUV色彩空间，因为YUV的Y分量拥有较其余分量更多的图像信息(边缘、材质、纹理等)，并且拥有更高的SNR(信噪比)，而U、V分量拥有更多的低频信息。\ 因此对于C-BM3D，本文使用Y分量搜寻相似块，U、V分量使用Y分量的相似块位置信息。\ \ 可以仿照他转换色彩空间、从Y分量搜索相似块的方法，来对彩色图像进行降噪处理。

五、 VBM3D

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

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

html  

``` clear all;clc ;

pauseTime = 1;

filePaths = 'cameraman256.png';

noise_leval = [10,15,20,25,30,35,40,45,50,55,60,65,70];

for i = 1:length(noise_leval)

PSNRs = [];

SSIMs = [];

sigma = noise_leval(i);

for j = 1:length(filePaths)

y = imread(filePaths);

if length(size(y)) > 2

y = rgb2gray(y);

end

y = im2double(y);

z = y + (sigma/255)*randn(size(y));

% 生成噪声图像

[PSNR,SSIM,y_est] = BM3D(y, z, sigma, 'np', 0);

PSNRs(j) = PSNR;

SSIMs(j) = SSIM;

imshow(cat(2,im2uint8(y),im2uint8(z),im2uint8(y_est)));

title([num2str(sigma),num2str(PSNR,'%2.2f'),'dB',' ',num2str(SSIMs(j),'%2.4f')])

drawnow;

pause(pauseTime)

end

disp(['sigma:', sigma,' psnr:', mean(PSNRs),' ssim:', mean(SSIMs)]);

end ```