1. BM3D模型简介
BM3D模型是一个两阶段图像去噪方法，主要包含两个步骤:
(1) 在噪声图像上，利用局部区域搜索相似块，并进行堆叠，在变换域(DCT域、FFT域）利用硬阈值去噪方法对堆叠的图像块进行去噪，获得堆叠相似块的估计值，最后，根据均值权重进行聚合；
(2) 通过步骤(1) 获取初步估计的图像，在初步估计的图像上进行相似块的聚合； 然后，利用维纳协同滤波进行图像去噪，从而，获取最后的去噪结果
2. 模型实现(代码参考网络实现)：
% BM3D_Color_Demo % BM3D 在彩色图像上去噪 % Author: HSW % Date: 2018-05-06 % clc; close all; clear all; img_org = imread('timg.png'); figure(1); imshow(img_org); title('原图像'); % 加噪声 sigma = 25; img_noise = double(img_org)+sigma * randn(size(img_org)); figure; imshow(img_noise / 255, []); title('噪声图像'); img_denoise = BM3D_Color(img_noise, 0, sigma, 0, 1); figure; imshow(img_denoise / 255, []); title('去噪图像'); % BM3D_Gray_Demo % BM3D 在灰度图像上去噪 % Author: HSW % Date: 2018-05-06 % clc; close all; clear all; img_org = imread('timg.png'); img_gray = rgb2gray(img_org); figure(1); imshow(img_gray); title('原图像'); % 加噪声 sigma = 25; img_noise = double(img_gray)+sigma * randn(size(img_gray)); figure; imshow(img_noise / 255, []); title('噪声图像'); img_denoise = BM3D_Gray(img_noise, 0, sigma, 1); figure; imshow(img_denoise / 255, []); title('去噪图像'); function img_denoise = BM3D_Color(img_noise, tran_mode, sigma, color_mode, isDisplay) % BM3D实现去噪 % Inputs: % img_noise: 噪声图像 % tran_mode: 变换方法： 默认值为0, tran_mode: = 0, fft; = 1, dct; = 2, dwt, = 3, db1 % sigma: 噪声水平，默认值为10 % color_mode: 彩色图像去噪时采用的颜色空间, 默认值为0, color_mode: = 0, YUV; = 1, YCbCr; = 2, OPP % Ouputs: % img_out: 去噪图像 % 参考文献：An Analysis and Implementation of the BM3D Image Denoising Method % Inputs: % img_in: 噪声图像，必须为矩形方阵 % tran_mode: = 0, FFT; = 1, DCT; = 2, DWT, = 3, db1 % Outputs: % img_denoise: 去噪图像 % % if ~exist('isDisplay', 'var') isDisplay = 0; end if ~exist('color_mode', 'var') color_mode = 0; end if ~exist('sigma', 'var') sigma = 10; end if ~exist('tran_mode', 'var') tran_mode = 0; end [row, col, dims] = size(img_noise); img_trans = rgb2other(img_noise, color_mode); % First Step 参数 kHard = 8; % 块大小 pHard = 4; % 块移动间隔 lambda_distHard = 0; % 求相似的距离时，变换后，收缩的阈值 nHard = 40; % 搜索窗口大小 NHard = 28; % 最多相似块个数 tauHard = 5000; % 最大的相似距离for fft % kaiser窗口的参数，实际上并没有特别大的影响 beta=2; Wwin2D = kaiser(kHard, beta) * kaiser(kHard, beta)'; % Second Step参数 kWien = kHard; pWien = pHard; lambda_distWien = lambda_distHard; nWien = nHard; NWien = NHard; tauWien = tauHard; sigma2 = sigma*sigma; if tran_mode == 0 % FFT lambda2d=400; lambda1d=500; lambda2d_wie=50; lambda1d_wie=500; elseif tran_mode == 1 % DCT lambda2d=50; lambda1d=80; lambda2d_wie=20; lambda1d_wie=60; elseif tran_mode == 2 % DWT lambda2d=50; lambda1d=80; lambda2d_wie=20; lambda1d_wie=60; end fprintf('BM3D: First Stage Start...\n'); %block为原始图像块， tran_block为FFT变换且硬阈值截断后的频域系数(频域， 计算距离的时候采用的是变换块) [block_ch1, tran_block_ch1, block2row_idx_ch1, block2col_idx_ch1] = im2block(img_trans(:,:,1), kHard, pHard, lambda_distHard, 0); [block_ch2, tran_block_ch2, block2row_idx_ch2, block2col_idx_ch2] = im2block(img_trans(:,:,2), kHard, pHard, lambda_distHard, 0); [block_ch3, tran_block_ch3, block2row_idx_ch3, block2col_idx_ch3] = im2block(img_trans(:,:,3), kHard, pHard, lambda_distHard, 0); %bn_r和bn_c为行和列上的图像块个数 bn_r = floor((row - kHard) / pHard) + 1; bn_c = floor((col - kHard) / pHard) + 1; %基础估计的图像 img_basic_sum = zeros(row, col, 3); img_basic_weight = zeros(row, col, 3); %对每个块遍历 for i=1:bn_r for j=1:bn_c % 利用亮度通道进行相似块搜索 [sim_blk_ch1, sim_num, sim_blk_idx] = search_similar_block(i, j, block_ch1, tran_block_ch1, floor(nHard/pHard), bn_r, bn_c, tauHard, NHard); % 进行亮度通道处理 % 协同滤波： 公式(2) tran3d_blk_shrink_ch1 = transform_3d(sim_blk_ch1, tran_mode, lambda2d, lambda1d); tran3d_blk_shrink_ch2 = transform_3d(block_ch2(:,:,sim_blk_idx), tran_mode, lambda2d, lambda1d); tran3d_blk_shrink_ch3 = transform_3d(block_ch3(:,:,sim_blk_idx), tran_mode, lambda2d, lambda1d); % 聚合： 公式(3)中的说明 NHard_P_ch1 = nnz(tran3d_blk_shrink_ch1); NHard_P_ch2 = nnz(tran3d_blk_shrink_ch2); NHard_P_ch3 = nnz(tran3d_blk_shrink_ch3); if NHard_P_ch1 > 1 wHard_P_ch1 = 1 / NHard_P_ch1; else wHard_P_ch1 = 1; end if NHard_P_ch2 > 1 wHard_P_ch2 = 1 / NHard_P_ch2; else wHard_P_ch2 = 1; end if NHard_P_ch3 > 1 wHard_P_ch3 = 1 / NHard_P_ch3; else wHard_P_ch3 = 1; end blk_est_ch1 = inv_transform_3d(tran3d_blk_shrink_ch1,tran_mode); blk_est_ch1 = real(blk_est_ch1); blk_est_ch2 = inv_transform_3d(tran3d_blk_shrink_ch2, tran_mode); blk_est_ch2 = real(blk_est_ch2); blk_est_ch3 = inv_transform_3d(tran3d_blk_shrink_ch3, tran_mode); blk_est_ch3 = real(blk_est_ch3); % 公式(3): 对亮度通道，即第1个通道 for k=1:sim_num idx = sim_blk_idx(k); ir = block2row_idx_ch1(idx); jr = block2col_idx_ch1(idx); img_basic_sum(ir:ir+kHard-1, jr:jr+kHard-1, 1) = img_basic_sum(ir:ir+kHard-1, jr:jr+kHard-1, 1) + wHard_P_ch1 * blk_est_ch1(:, :, k); img_basic_weight(ir:ir+kHard-1, jr:jr+kHard-1, 1) = img_basic_weight(ir:ir+kHard-1, jr:jr+kHard-1, 1) + wHard_P_ch1; img_basic_sum(ir:ir+kHard-1, jr:jr+kHard-1, 2) = img_basic_sum(ir:ir+kHard-1, jr:jr+kHard-1, 2) + wHard_P_ch2 * blk_est_ch2(:, :, k); img_basic_weight(ir:ir+kHard-1, jr:jr+kHard-1, 2) = img_basic_weight(ir:ir+kHard-1, jr:jr+kHard-1, 2) + wHard_P_ch2; img_basic_sum(ir:ir+kHard-1, jr:jr+kHard-1, 3) = img_basic_sum(ir:ir+kHard-1, jr:jr+kHard-1, 3) + wHard_P_ch3 * blk_est_ch3(:, :, k); img_basic_weight(ir:ir+kHard-1, jr:jr+kHard-1, 3) = img_basic_weight(ir:ir+kHard-1, jr:jr+kHard-1, 3) + wHard_P_ch3; end end end img_basic = img_basic_sum ./ img_basic_weight; if isDisplay figure; img_rgb = other2rgb(img_basic, color_mode); imshow(img_rgb / 255.0 ,[]); title('BM3D:Fist Stage Result'); end fprintf('BM3D: First Stage End...\n'); fprintf('BM3D: Second Stage Start...\n'); [block_basic_ch1,tran_block_basic_ch1,block2row_idx_basic_ch1,block2col_idx_basic_ch1] = im2block(img_basic(:, :, 1), kWien, pWien, lambda_distWien, 0); [block_basic_ch2,tran_block_basic_ch2,block2row_idx_basic_ch3,block2col_idx_basic_ch2] = im2block(img_basic(:, :, 2), kWien, pWien, lambda_distWien, 0); [block_basic_ch3,tran_block_basic_ch3,block2row_idx_basic_ch3,block2col_idx_basic_ch3] = im2block(img_basic(:, :, 3), kWien, pWien, lambda_distWien, 0); bn_r = floor((row - kWien) / pWien) + 1; bn_c = floor((col - kWien) / pWien) + 1; img_wien_sum = zeros(row, col, 3); img_wien_weight = zeros(row, col, 3); for i=1:1:bn_r for j=1:1:bn_c % 公式(5), 利用亮度进行相似性搜索 [sim_blk_basic_ch1, sim_num, sim_blk_basic_idx] = search_similar_block(i, j, block_basic_ch1, tran_block_basic_ch1, floor(nWien/pWien), bn_r, bn_c, tauWien, NWien); % 公式(6) tran3d_blk_basic_ch1 = transform_3d(sim_blk_basic_ch1, tran_mode, lambda2d_wie, lambda1d_wie); tran3d_blk_basic_ch2 = transform_3d(block_basic_ch2(:, :, sim_blk_basic_idx), tran_mode, lambda2d_wie, lambda1d_wie); tran3d_blk_basic_ch3 = transform_3d(block_basic_ch3(:, :, sim_blk_basic_idx), tran_mode, lambda2d_wie, lambda1d_wie); omega_P_ch1 = (tran3d_blk_basic_ch1.^2) ./ ((tran3d_blk_basic_ch1.^2) + sigma2); omega_P_ch2 = (tran3d_blk_basic_ch2.^2) ./ ((tran3d_blk_basic_ch2.^2) + sigma2); omega_P_ch3 = (tran3d_blk_basic_ch3.^2) ./ ((tran3d_blk_basic_ch3.^2) + sigma2); % 公式(7) tran3d_blk_ch1 = transform_3d(block_ch1(:, :, sim_blk_basic_idx), tran_mode, lambda2d_wie, lambda1d_wie); tran3d_blk_ch2 = transform_3d(block_ch2(:, :, sim_blk_basic_idx), tran_mode, lambda2d_wie, lambda1d_wie); tran3d_blk_ch3 = transform_3d(block_ch3(:, :, sim_blk_basic_idx), tran_mode, lambda2d_wie, lambda1d_wie); blk_est_ch1 = inv_transform_3d(omega_P_ch1 .* tran3d_blk_ch1, tran_mode); blk_est_ch2 = inv_transform_3d(omega_P_ch2 .* tran3d_blk_ch2, tran_mode); blk_est_ch3 = inv_transform_3d(omega_P_ch3 .* tran3d_blk_ch3, tran_mode); blk_est_ch1 = real(blk_est_ch1); blk_est_ch2 = real(blk_est_ch2); blk_est_ch3 = real(blk_est_ch3); NWien_P_ch1 = nnz(omega_P_ch1); NWien_P_ch2 = nnz(omega_P_ch2); NWien_P_ch3 = nnz(omega_P_ch3); if NWien_P_ch1 > 1 wWien_P_ch1 = 1 / (NWien_P_ch1); else wWien_P_ch1 = 1; end if NWien_P_ch2 > 1 wWien_P_ch2 = 1/(NWien_P_ch2); else wWien_P_ch2 = 1; end if NWien_P_ch3 > 1 wWien_P_ch3 = 1 / (NWien_P_ch3); else wWien_P_ch3 = 1; end % 公式(8) for k=1:sim_num idx=sim_blk_basic_idx(k); ir=block2row_idx_basic_ch1(idx); jr=block2col_idx_basic_ch1(idx); img_wien_sum(ir:ir+kWien-1, jr:jr+kWien-1, 1) = img_wien_sum(ir:ir+kWien-1, jr:jr+kWien-1, 1) + wWien_P_ch1 * blk_est_ch1(:, :, k); img_wien_weight(ir:ir+kWien-1, jr:jr+kWien-1, 1) = img_wien_weight(ir:ir+kWien-1, jr:jr+kWien-1, 1) + wWien_P_ch1; img_wien_sum(ir:ir+kWien-1, jr:jr+kWien-1, 2) = img_wien_sum(ir:ir+kWien-1, jr:jr+kWien-1, 2) + wWien_P_ch2 * blk_est_ch2(:, :, k); img_wien_weight(ir:ir+kWien-1, jr:jr+kWien-1, 2) = img_wien_weight(ir:ir+kWien-1, jr:jr+kWien-1, 2) + wWien_P_ch2; img_wien_sum(ir:ir+kWien-1, jr:jr+kWien-1, 3) = img_wien_sum(ir:ir+kWien-1, jr:jr+kWien-1, 3) + wWien_P_ch3 * blk_est_ch3(:, :, k); img_wien_weight(ir:ir+kWien-1, jr:jr+kWien-1, 3) = img_wien_weight(ir:ir+kWien-1, jr:jr+kWien-1, 3) + wWien_P_ch3; end end end img_other = img_wien_sum ./ img_wien_weight; img_denoise = other2rgb(img_other, color_mode); fprintf('BM3D: Second Stage End\n'); function img_denoise = BM3D_Gray(img_noise, tran_mode, sigma, isDisplay) % 参考文献：An Analysis and Implementation of the BM3D Image Denoising Method % Inputs: % img_noise: 灰度噪声图像，必须为矩形方阵 % tran_mode: = 0, fft; = 1, dct; = 2, dwt, = 3, db1 % Outputs: % img_denoise: 去噪图像 % if ~exist('tran_mode', 'var') tran_mode = 0; end if ~exist('sigma', 'var') sigma = 10; end if ~exist('isDisplay', 'var') isDisplay = 0; end [row,col] = size(img_noise); % First Step 参数 kHard = 8; % 块大小 pHard = 4; % 块移动间隔 lambda_distHard = 0; % 求相似的距离时，变换后，收缩的阈值 nHard = 40; % 搜索窗口大小 NHard = 28; % 最多相似块个数 tauHard = 5000; % 最大的相似距离for fft % kaiser窗口的参数，实际上并没有特别大的影响 beta=2; Wwin2D = kaiser(kHard, beta) * kaiser(kHard, beta)'; % Second Step参数 kWien = kHard; pWien = pHard; lambda_distWien = lambda_distHard; nWien = nHard; NWien = NHard; tauWien = tauHard; sigma2 = sigma*sigma; if(tran_mode==0)