BM3D 图像降噪算法与 Python 实现

Edvard 的视界深度之旅

已于 2023-05-31 09:13:46 修改

阅读量3.5k

点赞数 5

分类专栏： # 图像降噪图形图像算法文章标签：算法论文阅读计算机视觉图像处理 python

于 2023-05-16 11:26:59 首次发布

本文链接：https://blog.csdn.net/edvard_hua/article/details/130701681

版权

图形图像算法同时被 2 个专栏收录

9 篇文章

订阅专栏

图像降噪

4 篇文章

订阅专栏

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 * noise

img_dct = np.float32(image)
dct = cv2.dct(img_dct)
# 方便可视化，这里除于100，并把大于10的都置零
dct_out = abs(dct) / 100
dct_out[dct_out > 8] = 0

h, w = dct_out.shape

X = 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 math


def 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_str
    print(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, :])  # 二维kaiser

diffT = 100  # 允许纳入数组的最大diff
coefT = 0.0005

diffArray = 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 - searchsize

newh = math.floor((size[0] - newh - blocksize) / blockstep) * blockstep + newh
neww = math.floor((size[1] - neww - blocksize) / blockstep) * blockstep + neww

imnoise = 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 = 1
        diffArray[0] = 0
        # 不算自己
        for sy in range(systart, syend, blockstep):
            for sx in range(sxstart, sxend, blockstep):
                if sy == y and sx == x:
                    continue

                diff = np.sum(np.abs(
                    imnoise[y: y + blocksize, x: x + blocksize] - imnoise[sy: sy + blocksize, sx: sx + blocksize]))

                if diff > diffT:
                    continue
                # 塞入

                changeid = 0
                if hasboxnum < blockmaxnum - 1:
                    changeid = hasboxnum
                    hasboxnum = hasboxnum + 1
                else:
                    # 排序
                    for difid in range(1, blockmaxnum - 1):

                        if diff < diffArray[difid]:
                            changeid = difid

                if changeid != 0:
                    similarBoxArray[:, :, changeid] = imnoise[sy: sy + blocksize, sx: sx + blocksize]
                    diffArray[changeid] = diff

        # 开始做dct
        for difid in range(1, hasboxnum):
            similarBoxArray[:, :, difid] = cv2.dct(similarBoxArray[:, :, difid])

        # 开始做1d,阈值操作，计算非零个数
        notzeronum = 0
        for y1d in range(0, blocksize):
            for x1d in range(0, blocksize):
                temp3ddct = cv2.dct(similarBoxArray[y1d, x1d, :])
                zeroidx = np.abs(temp3ddct) < coefT

                temp3ddct[zeroidx] = 0
                notzeronum = notzeronum + temp3ddct[temp3ddct != 0].size
                similarBoxArray[y1d, x1d, :] = cv2.idct(temp3ddct)[:, 0]

        if notzeronum < 1:
            notzeronum = 1

        # 最终恢复图像所用到的分子分母
        for difid in range(1, hasboxnum):
            # 求权重：kaiser * w
            weight = kai / ((sigma ** 2) * notzeronum)
            imidct = cv2.idct(similarBoxArray[:, :, difid])

            # 分子，分母
            imnum[y: y + blocksize, x: x + blocksize] += imidct * weight
            imden[y: y + blocksize, x: x + blocksize] += weight

        x = x + blocksize

    y = y + blocksize

# 完成去噪
imout = imnum / imden

cv2.imshow('in', imnoise)
cv2.waitKey(0)
cv2.imshow('out', imout)
cv2.waitKey(0)