BM3D 图像降噪算法与 Python 实现

Image denoising with block-matching and 3D filtering. SPID 2006

https://www.spiedigitallibrary.org/conference-proceedings-of-spie/6064/606414/Image-denoising-with-block-matching-and-3D-filtering/10.1117/12.643267.short

是啥

BM3D 是一种基于分块匹配和三维协同滤波的图像降噪算法，与之前讲的 NLM 算法有以下两点不同之处

降噪时算法操作的单位是图像块，而不是 NLM 的逐像素处理，速度比 NLM 快
对于相似块，块内先进行正交变换，然后叠加在一起变成三维矩阵后再块间进行正交变换，文中称为三维滤波（也称协同滤波**）**

BM3D 的优缺点

优点：

在保持较高降噪效果的同时，算法的计算速度相对较快，适合大规模图像处理。
BM3D算法采用分块匹配策略，可以处理非线性噪声，适用范围广。
算法具有较强的自适应性和鲁棒性，可以有效处理多种不同类型的图像噪声

缺点：

BM3D算法的参数调节相对较为复杂，需要一定的经验和技巧才能得到最佳的降噪效果。
算法对于高斯噪声的降噪效果较好，但对于其他噪声类型的降噪效果可能不如其他算法。
算法在处理低对比度图像时可能会出现失真和伪影等问题，需要进行一定的处理和优化。

怎么做

算法流程

BM3D 算法流程分为两步，这两步的流程是一致的，在一些计算上会有差异

分块，把图像分成一个个大小相同的块，每个块大小为 NxN
匹配，对于每个块，通过相似性度量找到最相似的 K 个块
聚合，对上一步获得的块进行排序，通过硬阈值（第二步为 Wiener Filtering）的方式获得 M 个块
协同滤波，将聚合获得的块组合成三维矩阵，先块内正交后块间正交，最后卡阈值获得降噪块
重构，将所有的块拼起来，获得降噪后的图片

大体和 NLM 流程相似，主要区别是协同滤波，它的优点是可以更好的保存图像细节

协同滤波

由上面流程所述，对于相似块，协同滤波块内先进行正交变换，然后叠加在一起变成三维矩阵后再块间进行正交变换。

在这里，正交变换是指将图像信号从时域转换成频域，如离散余弦变换（DCT）和傅里叶变换，如下图所示

如何生成右边频域图的 python 代码如下

import cv2
import numpy as np
import matplotlib.pyplot as plt# Set up a figure twice as tall as it is wide
fig = plt.figure()
fig.suptitle('DCT Visualization', fontsize=16)# Second subplot
ax = fig.add_subplot(1, 1, 1, projection='3d')image = cv2.imread("assets/a0016.jpg", 0)
# resize 到 100 方便后面可视化
image = cv2.resize(image, (100, 100))
# 加入噪音
# size = image.shape
# noise = np.random.normal(0, 0.001 ** 0.5, size)
# image = image + 1 * noiseimg_dct = np.float32(image)
dct = cv2.dct(img_dct)
# 方便可视化，这里除于100，并把大于10的都置零
dct_out = abs(dct) / 100
dct_out[dct_out > 8] = 0h, w = dct_out.shapeX = np.arange(0, w, 1)
Y = np.arange(0, h, 1)
X, Y = np.meshgrid(X, Y)surf = ax.plot_surface(X, Y, dct_out, rstride=1, cstride=1,linewidth=0, antialiased=False)
ax.set_zlim(-4, 8)plt.show()

一般来说噪音都是高斯随机分布，会遍布整个区域，包括高频和低频，而非噪音一般在低频区域，通过将小于某个值置零的方法（也就是硬阈值）可以去除掉噪音。关于时域转频域的相关原理后续单独展开说

但若噪音方与低频非噪音相近，硬阈值则会抹掉图像细节，因此对块间也进行了正交变换，再次提取相似块相同位置的低频信息

为什么块间正交变换后就可以保留更多的图像细节呢？

我们假设找到的相似块，他们图像细节抛除噪音后是一致的，那么我们在横向计算正交变换的时候，理论上非噪音的部分还是会聚集在低频区域，噪音还是会遍布整个范围，因此可以再次使用硬阈值

Python 代码实现

下面代码只实现了第一步

import cv2
import numpy as np
import mathdef process_bar(percent, start_str='', end_str='', total_length=0):bar = ''.join(["\033[31m%s\033[0m" % '   '] * int(percent * total_length)) + ''bar = '\r' + start_str + bar.ljust(total_length) + ' {:0>4.1f}%|'.format(percent * 100) + end_strprint(bar, end='', flush=True)image = cv2.imread('assets/a0016.jpg')
image = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)size = image.shape  # h w c# 生成噪声
image = np.array(image / 255, dtype=np.float32)
noise = np.random.normal(0, 0.001 ** 0.5, size)
imnoise = image + 1 * noise# cv2.imshow('result',imnoise)
# cv2.waitKey()sigma = 1
searchsize = 16  # 图像搜索区域半径
blocksize = 4  # 图像块尺寸
blockstep = 2  # 搜索步长
blockmaxnum = 8  # 相似块最大数量
searchblocksize = searchsize / blocksize  # 半径内块个数
kai = np.kaiser(blocksize, 5)
kai = np.dot(kai[:, None], kai[None, :])  # 二维kaiserdiffT = 100  # 允许纳入数组的最大diff
coefT = 0.0005diffArray = np.zeros(blockmaxnum)  # 相似块
similarBoxArray = np.zeros((blocksize, blocksize, blockmaxnum))
# jinbuqu csdn newland# 让图像尺寸符合倍数
newh = math.floor((size[0] - blocksize) / blockstep) * blockstep - searchsize
neww = math.floor((size[1] - blocksize) / blockstep) * blockstep - searchsizenewh = math.floor((size[0] - newh - blocksize) / blockstep) * blockstep + newh
neww = math.floor((size[1] - neww - blocksize) / blockstep) * blockstep + newwimnoise = imnoise[0:newh, 0:neww]# 初始化分子分母
imnum = np.zeros(imnoise.shape)
imden = np.zeros(imnoise.shape)# 将左上角作为块的坐标.每个块独立做一次去噪，只要关注一块的计算就行
for y in range(0, newh - blocksize, blockstep):  # 检查能不能完整走到底,范围对不对systart = max(0, y - searchsize)syend = min(newh - blocksize, y + searchsize - 1)process_bar(y / (newh - blocksize), start_str='进度', end_str="100", total_length=15)for x in range(0, neww - blocksize, blockstep):sxstart = max(0, x - searchsize)sxend = min(neww - blocksize, x + searchsize - 1)# 排序矩阵初始化similarBoxArray[:, :, 0] = imnoise[y: y + blocksize, x: x + blocksize]hasboxnum = 1diffArray[0] = 0# 不算自己for sy in range(systart, syend, blockstep):for sx in range(sxstart, sxend, blockstep):if sy == y and sx == x:continuediff = np.sum(np.abs(imnoise[y: y + blocksize, x: x + blocksize] - imnoise[sy: sy + blocksize, sx: sx + blocksize]))if diff > diffT:continue# 塞入changeid = 0if hasboxnum < blockmaxnum - 1:changeid = hasboxnumhasboxnum = hasboxnum + 1else:# 排序for difid in range(1, blockmaxnum - 1):if diff < diffArray[difid]:changeid = difidif changeid != 0:similarBoxArray[:, :, changeid] = imnoise[sy: sy + blocksize, sx: sx + blocksize]diffArray[changeid] = diff# 开始做dctfor difid in range(1, hasboxnum):similarBoxArray[:, :, difid] = cv2.dct(similarBoxArray[:, :, difid])# 开始做1d,阈值操作，计算非零个数notzeronum = 0for y1d in range(0, blocksize):for x1d in range(0, blocksize):temp3ddct = cv2.dct(similarBoxArray[y1d, x1d, :])zeroidx = np.abs(temp3ddct) < coefTtemp3ddct[zeroidx] = 0notzeronum = notzeronum + temp3ddct[temp3ddct != 0].sizesimilarBoxArray[y1d, x1d, :] = cv2.idct(temp3ddct)[:, 0]if notzeronum < 1:notzeronum = 1# 最终恢复图像所用到的分子分母for difid in range(1, hasboxnum):# 求权重：kaiser * wweight = kai / ((sigma ** 2) * notzeronum)imidct = cv2.idct(similarBoxArray[:, :, difid])# 分子，分母imnum[y: y + blocksize, x: x + blocksize] += imidct * weightimden[y: y + blocksize, x: x + blocksize] += weightx = x + blocksizey = y + blocksize# 完成去噪
imout = imnum / imdencv2.imshow('in', imnoise)
cv2.waitKey(0)
cv2.imshow('out', imout)
cv2.waitKey(0)

参考链接

下面这篇文章写的很详细，图文并茂
https://blog.csdn.net/qq_33552519/article/details/108632146

https://blog.csdn.net/sinat_29018995/article/details/123032100