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(