BM3D 图像降噪快速算法的 MATLAB 实现

1. 快速 BM3D 算法流程（概述）

步骤	操作	加速技巧
① 分组	块匹配 + 堆叠	FFT 互相关
② 协同滤波	3D 变换 + 硬阈值	FFT 沿第三维
③ 聚合	加权平均	稀疏矩阵累加

---

2. 核心函数（单文件版）

保存为 bm3d_fast.m 即可调用：

function [img_denoised] = bm3d_fast(img_noisy, sigma)
% 快速 BM3D 图像降噪（纯 MATLAB，FFT 加速）
% 输入：img_noisy 灰度图 0~255
%       sigma     噪声标准差
% 输出：img_denoised 同尺寸

img_noisy = double(img_noisy);
[H,W] = size(img_noisy);

%% 参数（与作者论文一致）
block_size = 8;            % 块尺寸
step       = 3;            % 滑动步长
max_similar = 16;          % 最大相似块数
tau_hard   = 2.7*sigma;    % 硬阈值系数
tau_wien   = 2.5*sigma;    % Wiener 阈值

%% Step1：基础估计（硬阈值）
basic = bm3d_step1(img_noisy, sigma, block_size, step, max_similar, tau_hard);

%% Step2：最终估计（Wiener 协同滤波）
img_denoised = bm3d_step2(img_noisy, basic, sigma, block_size, step, max_similar, tau_wien);

img_denoised = uint8(img_denoised);
end

---

3. Step1：基础估计（硬阈值）

function basic = bm3d_step1(img, sigma, bs, step, max_sim, tau)
[H,W] = size(img);
basic = zeros(H,W);  weight = zeros(H,W);

% 预计算 FFT 加速互相关
img_fft = fft2(img);

for i = bs/2+1 : step : H-bs/2
    for j = bs/2+1 : step : W-bs/2
        % 当前块
        block = img(i-bs/2:i+bs/2-1, j-bs/2:j+bs/2-1);
        block_fft = fft2(block);
        
        % 快速块匹配（FFT 互相关）
        corr = real(ifft2(block_fft .* conj(img_fft)));
        [vals, idx] = sort(corr(:),'descend');
        idx = idx(1:max_sim);  % 最相似块
        [di,dj] = ind2sub([H,W], idx);
        
        % 堆叠 3D 组
        group = zeros(bs,bs,max_sim);
        for k = 1:max_sim
            ii = di(k)-bs/2; jj = dj(k)-bs/2;
            group(:,:,k) = img(ii+1:ii+bs, jj+1:jj+bs);
        end
        
        % 3D 变换（FFT 沿第三维）+ 硬阈值
        group_fft = fft(group,[],3);
        group_fft(abs(group_fft) < tau*sigma) = 0;
        group_est = real(ifft(group_fft,[],3));
        
        % 加权聚合
        w = 1/(sigma^2 * max_sim + 1e-6);
        for k = 1:max_sim
            ii = di(k)-bs/2; jj = dj(k)-bs/2;
            basic(ii+1:ii+bs, jj+1:jj+bs) = basic(ii+1:ii+bs, jj+1:jj+bs) + w * group_est(:,:,k);
            weight(ii+1:ii+bs, jj+1:jj+bs) = weight(ii+1:ii+bs, jj+1:jj+bs) + w;
        end
    end
end

basic = basic ./ (weight + 1e-6);
end

---

4. Step2：最终估计（Wiener 协同滤波）

function final = bm3d_step2(img, basic, sigma, bs, step, max_sim, tau)
[H,W] = size(img);
final = zeros(H,W);  weight = zeros(H,W);

basic_fft = fft2(basic);

for i = bs/2+1 : step : H-bs/2
    for j = bs/2+1 : step : W-bs/2
        % 基础图块匹配
        block_b = basic(i-bs/2:i+bs/2-1, j-bs/2:j+bs/2-1);
        block_b_fft = fft2(block_b);
        corr = real(ifft2(block_b_fft .* conj(basic_fft)));
        [vals, idx] = sort(corr(:),'descend');
        idx = idx(1:max_sim);
        [di,dj] = ind2sub([H,W], idx);
        
        % 两组：原图 + 基础图
        group_noisy = zeros(bs,bs,max_sim);
        group_basic = zeros(bs,bs,max_sim);
        for k = 1:max_sim
            ii = di(k)-bs/2; jj = dj(k)-bs/2;
            group_noisy(:,:,k) = img(ii+1:ii+bs, jj+1:jj+bs);
            group_basic(:,:,k) = basic(ii+1:ii+bs, jj+1:jj+bs);
        end
        
        % 3D FFT + Wiener 系数
        fft_noisy = fft(group_noisy,[],3);
        fft_basic = fft(group_basic,[],3);
        power = abs(fft_basic).^2;
        wiener_coef = power ./ (power + sigma^2 + 1e-6);
        fft_est = fft_noisy .* wiener_coef;
        group_est = real(ifft(fft_est,[],3));
        
        % 加权聚合
        w = 1/(sigma^2 + 1e-6);
        for k = 1:max_sim
            ii = di(k)-bs/2; jj = dj(k)-bs/2;
            final(ii+1:ii+bs, jj+1:jj+bs) = final(ii+1:ii+bs, jj+1:jj+bs) + w * group_est(:,:,k);
            weight(ii+1:ii+bs, jj+1:jj+bs) = weight(ii+1:ii+bs, jj+1:jj+bs) + w;
        end
    end
end

final = final ./ (weight + 1e-6);
end

参考代码 BM3D图像降噪快速算法 www.youwenfan.com/contentcnl/53384.html

5. 运行

img = imread('lena_gray.png');
img_noisy = imnoise(img,'gaussian',0,0.01);   % σ=10
img_denoised = bm3d_fast(img_noisy, 10);

figure; montage({img, img_noisy, img_denoised});
title({'原图','噪声图','BM3D 快速版'});