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/2for 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_simii = 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_simii = 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;endend
endbasic = 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/2for 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_simii = 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_simii = 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;endend
endfinal = 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 快速版'});