BM3Dpython代码

以下是使用 Python 实现 BM3D 算法的示例代码：

import numpy as np
from skimage.measure import compare_psnr
from skimage import io
import matplotlib.pyplot as plt

def block_matching(img, block_size, search_size):
    # img: input image
    # block_size: size of the blocks to match
    # search_size: size of the search window
    h, w = img.shape
    b = block_size // 2
    s = search_size // 2
    search_area = img[b:h-b, b:w-b]
    matches = np.zeros((h-2*b, w-2*b), dtype=np.int32)
    for i in range(-s, s+1):
        for j in range(-s, s+1):
            shifted = img[b+i:h-b+i, b+j:w-b+j]
            diff = (search_area - shifted) ** 2
            ssd = np.sum(diff, axis=(2,3))
            min_idx = np.argmin(ssd, axis=0)
            matches += min_idx
    return matches

def bm3d(img, sigma, block_size=8, search_size=21, window_size=35, threshold=2.7, hard_threshold=False):
    # img: input image
    # sigma: noise level
    # block_size: size of the blocks to match
    # search_size: size of the search window
    # window_size: size of the 3D transform window
    # threshold: threshold value for soft thresholding
    # hard_threshold: whether to use hard thresholding instead of soft thresholding
    h, w = img.shape
    b = block_size // 2
    s = search_size // 2
    img_pad = np.pad(img, b, mode='symmetric')
    matches = block_matching(img_pad, block_size, search_size)
    img_denoise = np.zeros((h, w))
    weights = np.zeros((h, w))
    for i in range(b, h-b):
        for j in range(b, w-b):
            block = img_pad[i-b:i+b+1, j-b:j+b+1]
            match = matches[i-b, j-b]
            similar_blocks = img_pad[match-b:match+b+1, j-b:j+b+1]
            similar_blocks = similar_blocks.reshape(block_size, block_size, -1)
            mean = np.mean(similar_blocks, axis=2, keepdims=True)
            similar_blocks -= mean
            cov = np.matmul(similar_blocks, similar_blocks.transpose((0,1,3,2))) / (block_size ** 2 - 1)
            eigvals, eigvecs = np.linalg.eigh(cov)
            eigvals = np.sqrt(np.maximum(eigvals, 0))
            transform = np.matmul(eigvecs, eigvals[:, :, None] * eigvecs.transpose((0,1,3,2)))
            transform = np.moveaxis(transform, [2,3], [0,1])
            weights_window = np.zeros((window_size, window_size))
            for k in range(-s, s+1):
                for l in range(-s, s+1):
                    ii = i + k
                    jj = j + l
                    if ii < b or ii >= h+b or jj < b or jj >= w+b:
                        continue
                    block2 = img_pad[ii-b:ii+b+1, jj-b:jj+b+1]
                    match2 = matches[ii-b, jj-b]
                    if match2 != match:
                        continue
                    similar_blocks2 = img_pad[match2-b:match2+b+1, jj-b:jj+b+1]
                    similar_blocks2 = similar_blocks2.reshape(block_size, block_size, -1)
                    mean2 = np.mean(similar_blocks2, axis=2, keepdims=True)
                    similar_blocks2 -= mean2
                    transform2 = np.matmul(eigvecs, eigvals[:, :, None] * eigvecs.transpose((0,1,3,2)))
                    transform2 = np.moveaxis(transform2, [2,3], [0,1])
                    diff = block - block2
                    dist = np.sum(diff ** 2) / (block_size ** 2)
                    weight = np.exp(-dist / (2 * sigma ** 2))
                    weights_window += weight
                    if hard_threshold:
                        mask = np.abs(transform - transform2) > threshold * sigma
                        transform = np.where(mask, transform, 0)
                    else:
                        transform = transform - threshold * sigma ** 2 / (sigma ** 2 + dist) * (transform - transform2)
            if np.sum(weights_window) == 0:
                continue
            weights_window = weights_window / np.sum(weights_window)
            weights[i-b:i+b+1, j-b:j+b+1] += weights_window
            transform = np.sum(transform * weights_window[..., None, None], axis=(0,1))
            transform = transform / np.sqrt(np.sum(transform ** 2))
            block_denoise = np.matmul(transform, block.reshape(-1))
            block_denoise = np.clip(block_denoise, 0, 255)
            img_denoise[i-b:i+b+1, j-b:j+b+1] += block_denoise.reshape(block_size, block_size) * weights_window
    img_denoise /= weights
    img_denoise = np.clip(img_denoise, 0, 255).astype(np.uint8)
    return img_denoise

这个代码实现了 BM3D 算法的主要步骤，包括块匹配、3D 变换、软阈值滤波等。您可以通过调整函数参数来控制算法的性能。