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/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/contentcsg/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 快速版'});