ISP调试工程师-面试问题-001-BM3D算法

BM3D是一种基于**分组稀疏表示**的图像降噪算法，该算法结合了**块匹配3D**（Block-Matching 3D，BM3D）和小波变换。BM3D 3D变换滤波的基本思路是将信号

，看作3D信号，将每一个2D图像看作一个3D信号的一个切片，然后在3D信号上进行**匹配和滤波**，最后将3D信号切片后形成的2D降噪图像拼起来，即得到最终的降噪结果。其中，匹配可以通过**块模式匹配**（block-matching）实现，滤波可以通过**非线性阈值处理（nonlinear thresholding）**和 3D变换域滤波（3D transform-based filtering）实现。

以下是基于python3.7实现对一张jpg图的BM 3D算法的代码。

import numpy as np
from scipy.fftpack import dct, idct
from skimage import data, img_as_float
import matplotlib.pyplot as plt

# 先加载一张测试图片
img = img_as_float(data.camera())

# 添加高斯噪声
sigma = 0.1
img_noise = img + sigma * np.random.randn(*img.shape)

# 首先对图像进行分块
block_size = 8
H, W = img_noise.shape[:2]
block_H, block_W = block_size, block_size
blocks = np.zeros((H//block_H, W//block_W, block_H, block_W))

# 分块循环
for i in range(0, H - block_H + 1, block_H):
    for j in range(0, W - block_W + 1, block_W):
        blocks[i//block_H, j//block_W] = img_noise[i:i+block_H, j:j+block_W]

# 经过分块后处理每块，遍历每个块
for i in range(blocks.shape[0]):
    for j in range(blocks.shape[1]):
        # 获取当前块
        block = blocks[i, j]

        # 第一步：去噪，首先使用2D DCT变换将块转换为DCT系数表示
        dct_block = dct(dct(block.T, norm='ortho').T, norm='ortho')
        
        # 第二步：去噪，对DCT系数进行软阈值（soft thresholding）去噪
        # 具体实现就是从DCT系数中减去软阈值（threshold）得到新的DCT系数
        T = 0.8 * sigma * np.sqrt(2*np.log(block.size))
        noisy_dct_block = np.copy(dct_block)
        for k in range(0, block_H):
            for l in range(0, block_W):
                if abs(dct_block[k, l]) < T:
                    noisy_dct_block[k, l] = 0
                else:
                    noisy_dct_block[k, l] -= np.sign(dct_block[k, l]) * T
        
        # 第三步：去噪，将DCT系数用IDCT（反DCT）变换恢复到空域
        idct_block = idct(idct(noisy_dct_block.T, norm='ortho').T, norm='ortho')

        # 储存恢复后的块
        blocks[i, j] = idct_block

# 再次循环各个块，将其放在完整图像上，完成去噪
denoised_img = np.zeros_like(img_noise)

for i in range(blocks.shape[0]):
    for j in range(blocks.shape[1]):
        denoised_img[i*block_H:(i+1)*block_H, j*block_W:(j+1)*block_W] = blocks[i, j]
        
# 绘制结果，分别显示原图、加噪图和去噪后的图像
fig, axs = plt.subplots(1, 3, figsize=(10, 3))
axs[0].imshow(img, cmap='gray')
axs[0].set_title('Original Image')
axs[1].imshow(img_noise, cmap='gray')
axs[1].set_title('Noisy Image (Sigma={:.2f})'.format(sigma))
axs[2].imshow(denoised_img, cmap='gray')
axs[2].set_title('BM3D Image (Sigma={:.2f})'.format(sigma))
plt.show()

得到的效果。