宝塔服务器面板,一键全能部署及管理,送你10850元礼包,点我领取

1. BM3D模型简介

BM3D模型是一个两阶段图像去噪方法,主要包含两个步骤:

1) 在噪声图像上,利用局部区域搜索相似块,并进行堆叠,在变换域DCT域、FFT域)利用硬阈值去噪方法对堆叠的图像块进行去噪,获得堆叠相似块的估计值,最后,根据均值权重进行聚合;

2) 通过步骤1) 获取初步估计的图像,在初步估计的图像上进行相似块的聚合; 然后,利用维纳协同滤波进行图像去噪,从而,获取最后的去噪结果

图像去噪序列——BM3D图像去噪模型实现-风君子博客图像去噪序列——BM3D图像去噪模型实现-风君子博客


2. 模型实现代码参考网络实现):

% BM3D_Color_Demo 
% BM3D 在彩色图像上去噪
% Author: HSW
% Date: 2018-05-06 
%

clc; 
close all; 
clear all; 

img_org = imread'timg.png'); 

figure1); 
imshowimg_org); 
title'原图像'); 

% 加噪声
sigma = 25; 
img_noise = doubleimg_org)+sigma * randnsizeimg_org));

figure; 
imshowimg_noise / 255, []); 
title'噪声图像'); 

img_denoise = BM3D_Colorimg_noise, 0, sigma, 0, 1); 

figure; 
imshowimg_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 = rgb2grayimg_org); 

figure1); 
imshowimg_gray); 
title'原图像'); 

% 加噪声
sigma = 25; 
img_noise = doubleimg_gray)+sigma * randnsizeimg_gray));

figure; 
imshowimg_noise / 255, []); 
title'噪声图像'); 

img_denoise = BM3D_Grayimg_noise, 0, sigma, 1); 

figure; 
imshowimg_denoise / 255, []); 
title'去噪图像'); 

function img_denoise = BM3D_Colorimg_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] = sizeimg_noise);
img_trans = rgb2otherimg_noise, color_mode);
% First Step 参数
kHard           = 8;          % 块大小
pHard           = 4;          % 块移动间隔
lambda_distHard = 0;          % 求相似的距离时,变换后,收缩的阈值
nHard           = 40;         % 搜索窗口大小
NHard           = 28;         % 最多相似块个数
tauHard         = 5000;       % 最大的相似距离for fft
% kaiser窗口的参数,实际上并没有特别大的影响
beta=2;
Wwin2D = kaiserkHard, beta) * kaiserkHard, 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] = im2blockimg_trans:,:,1), kHard, pHard, lambda_distHard, 0);
[block_ch2, tran_block_ch2, block2row_idx_ch2, block2col_idx_ch2] = im2blockimg_trans:,:,2), kHard, pHard, lambda_distHard, 0); 
[block_ch3, tran_block_ch3, block2row_idx_ch3, block2col_idx_ch3] = im2blockimg_trans:,:,3), kHard, pHard, lambda_distHard, 0); 
%bn_r和bn_c为行和列上的图像块个数
bn_r = floorrow - kHard) / pHard) + 1;
bn_c = floorcol - kHard) / pHard) + 1;
%基础估计的图像
img_basic_sum = zerosrow, col, 3);
img_basic_weight = zerosrow, col, 3);
%对每个块遍历
for i=1:bn_r
for j=1:bn_c
% 利用亮度通道进行相似块搜索
[sim_blk_ch1, sim_num, sim_blk_idx] = search_similar_blocki, j, block_ch1, tran_block_ch1, floornHard/pHard), bn_r, bn_c, tauHard, NHard);
% 进行亮度通道处理
% 协同滤波: 公式2)
tran3d_blk_shrink_ch1 = transform_3dsim_blk_ch1, tran_mode, lambda2d, lambda1d);
tran3d_blk_shrink_ch2 = transform_3dblock_ch2:,:,sim_blk_idx), tran_mode, lambda2d, lambda1d); 
tran3d_blk_shrink_ch3 = transform_3dblock_ch3:,:,sim_blk_idx), tran_mode, lambda2d, lambda1d); 
% 聚合: 公式3)中的说明
NHard_P_ch1 = nnztran3d_blk_shrink_ch1);
NHard_P_ch2 = nnztran3d_blk_shrink_ch2); 
NHard_P_ch3 = nnztran3d_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_3dtran3d_blk_shrink_ch1,tran_mode);
blk_est_ch1 = realblk_est_ch1);
blk_est_ch2 = inv_transform_3dtran3d_blk_shrink_ch2, tran_mode); 
blk_est_ch2 = realblk_est_ch2); 
blk_est_ch3 = inv_transform_3dtran3d_blk_shrink_ch3, tran_mode); 
blk_est_ch3 = realblk_est_ch3); 
% 公式3): 对亮度通道,即第1个通道
for k=1:sim_num
idx = sim_blk_idxk);
ir = block2row_idx_ch1idx);
jr = block2col_idx_ch1idx);
img_basic_sumir:ir+kHard-1, jr:jr+kHard-1, 1) = img_basic_sumir:ir+kHard-1, jr:jr+kHard-1, 1) + wHard_P_ch1 * blk_est_ch1:, :, k);
img_basic_weightir:ir+kHard-1, jr:jr+kHard-1, 1) = img_basic_weightir:ir+kHard-1, jr:jr+kHard-1, 1) + wHard_P_ch1;
img_basic_sumir:ir+kHard-1, jr:jr+kHard-1, 2) = img_basic_sumir:ir+kHard-1, jr:jr+kHard-1, 2) + wHard_P_ch2 * blk_est_ch2:, :, k);
img_basic_weightir:ir+kHard-1, jr:jr+kHard-1, 2) = img_basic_weightir:ir+kHard-1, jr:jr+kHard-1, 2) + wHard_P_ch2;
img_basic_sumir:ir+kHard-1, jr:jr+kHard-1, 3) = img_basic_sumir:ir+kHard-1, jr:jr+kHard-1, 3) + wHard_P_ch3 * blk_est_ch3:, :, k);
img_basic_weightir:ir+kHard-1, jr:jr+kHard-1, 3) = img_basic_weightir: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 = other2rgbimg_basic, color_mode); 
imshowimg_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] = im2blockimg_basic:, :, 1), kWien, pWien, lambda_distWien, 0);
[block_basic_ch2,tran_block_basic_ch2,block2row_idx_basic_ch3,block2col_idx_basic_ch2] = im2blockimg_basic:, :, 2), kWien, pWien, lambda_distWien, 0);
[block_basic_ch3,tran_block_basic_ch3,block2row_idx_basic_ch3,block2col_idx_basic_ch3] = im2blockimg_basic:, :, 3), kWien, pWien, lambda_distWien, 0);
bn_r = floorrow - kWien) / pWien) + 1;
bn_c = floorcol - kWien) / pWien) + 1;
img_wien_sum = zerosrow, col, 3);
img_wien_weight = zerosrow, 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_blocki, j, block_basic_ch1, tran_block_basic_ch1, floornWien/pWien), bn_r, bn_c, tauWien, NWien);
% 公式6)
tran3d_blk_basic_ch1 = transform_3dsim_blk_basic_ch1, tran_mode, lambda2d_wie, lambda1d_wie);
tran3d_blk_basic_ch2 = transform_3dblock_basic_ch2:, :, sim_blk_basic_idx), tran_mode, lambda2d_wie, lambda1d_wie); 
tran3d_blk_basic_ch3 = transform_3dblock_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_3dblock_ch1:, :, sim_blk_basic_idx), tran_mode, lambda2d_wie, lambda1d_wie);
tran3d_blk_ch2 = transform_3dblock_ch2:, :, sim_blk_basic_idx), tran_mode, lambda2d_wie, lambda1d_wie); 
tran3d_blk_ch3 = transform_3dblock_ch3:, :, sim_blk_basic_idx), tran_mode, lambda2d_wie, lambda1d_wie); 
blk_est_ch1 = inv_transform_3domega_P_ch1 .* tran3d_blk_ch1, tran_mode);
blk_est_ch2 = inv_transform_3domega_P_ch2 .* tran3d_blk_ch2, tran_mode); 
blk_est_ch3 = inv_transform_3domega_P_ch3 .* tran3d_blk_ch3, tran_mode); 
blk_est_ch1 = realblk_est_ch1);
blk_est_ch2 = realblk_est_ch2); 
blk_est_ch3 = realblk_est_ch3); 
NWien_P_ch1 = nnzomega_P_ch1);
NWien_P_ch2 = nnzomega_P_ch2); 
NWien_P_ch3 = nnzomega_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_idxk);
ir=block2row_idx_basic_ch1idx);
jr=block2col_idx_basic_ch1idx);
img_wien_sumir:ir+kWien-1, jr:jr+kWien-1, 1) = img_wien_sumir:ir+kWien-1, jr:jr+kWien-1, 1) + wWien_P_ch1 * blk_est_ch1:, :, k);
img_wien_weightir:ir+kWien-1, jr:jr+kWien-1, 1) = img_wien_weightir:ir+kWien-1, jr:jr+kWien-1, 1) + wWien_P_ch1;
img_wien_sumir:ir+kWien-1, jr:jr+kWien-1, 2) = img_wien_sumir:ir+kWien-1, jr:jr+kWien-1, 2) + wWien_P_ch2 * blk_est_ch2:, :, k);
img_wien_weightir:ir+kWien-1, jr:jr+kWien-1, 2) = img_wien_weightir:ir+kWien-1, jr:jr+kWien-1, 2) + wWien_P_ch2;
img_wien_sumir:ir+kWien-1, jr:jr+kWien-1, 3) = img_wien_sumir:ir+kWien-1, jr:jr+kWien-1, 3) + wWien_P_ch3 * blk_est_ch3:, :, k);
img_wien_weightir:ir+kWien-1, jr:jr+kWien-1, 3) = img_wien_weightir:ir+kWien-1, jr:jr+kWien-1, 3) + wWien_P_ch3;
end
end
end
img_other = img_wien_sum ./ img_wien_weight; 
img_denoise = other2rgbimg_other, color_mode);
fprintf'BM3D: Second Stage End\n');

function img_denoise = BM3D_Grayimg_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] = sizeimg_noise);
% First Step 参数
kHard           = 8;          % 块大小
pHard           = 4;          % 块移动间隔
lambda_distHard = 0;          % 求相似的距离时,变换后,收缩的阈值
nHard           = 40;         % 搜索窗口大小
NHard           = 28;         % 最多相似块个数
tauHard         = 5000;       % 最大的相似距离for fft
% kaiser窗口的参数,实际上并没有特别大的影响
beta=2;
Wwin2D = kaiserkHard, beta) * kaiserkHard, beta)';
% Second Step参数
kWien           = kHard;
pWien           = pHard;
lambda_distWien = lambda_distHard;
nWien           = nHard;
NWien           = NHard;
tauWien         = tauHard;
sigma2          = sigma*sigma;
iftran_mode==0)        %fft
lambda2d=400;
lambda1d=500;
lambda2d_wie=50;
lambda1d_wie=500;
elseiftran_mode == 1)  %dct
lambda2d=50;
lambda1d=80;
lambda2d_wie=20;
lambda1d_wie=60;
elseiftran_mode == 2)  %dwt
lambda2d=50;
lambda1d=80;
lambda2d_wie=20;
lambda1d_wie=60;
end
%block为原始图像块, tran_block为FFT变换且硬阈值截断后的频域系数频域, 计算距离的时候采用的是变换块)
[block,tran_block,block2row_idx,block2col_idx]=im2blockimg_noise,kHard,pHard,lambda_distHard,0);
%bn_r和bn_c为行和列上的图像块个数
bn_r=floorrow-kHard)/pHard)+1;
bn_c=floorcol-kHard)/pHard)+1;
%基础估计的图像
img_basic_sum=zerosrow,col);
img_basic_weight=zerosrow,col);
%basic处理
fprintf'BM3D: First Stage Start...\n');
%对每个块遍历
for i=1:bn_r
for j=1:bn_c
[sim_blk,sim_num,sim_blk_idx]=search_similar_blocki,j,block,tran_block,floornHard/pHard),bn_r,bn_c,tauHard,NHard);
% 协同滤波: 公式2)
tran3d_blk_shrink=transform_3dsim_blk,tran_mode,lambda2d,lambda1d);
% 聚合: 公式3)中的说明
NHard_P=nnztran3d_blk_shrink);
ifNHard_P >1)
wHard_P=1/NHard_P;
else
wHard_P=1;
end
blk_est =inv_transform_3dtran3d_blk_shrink,tran_mode);
blk_est=realblk_est);
% 公式3)
for k=1:sim_num
idx=sim_blk_idxk);
ir=block2row_idxidx);
jr=block2col_idxidx);
img_basic_sumir:ir+kHard-1,jr:jr+kHard-1) = img_basic_sumir:ir+kHard-1,jr:jr+kHard-1) + wHard_P*blk_est:,:,k);
img_basic_weightir:ir+kHard-1,jr:jr+kHard-1) = img_basic_weightir:ir+kHard-1,jr:jr+kHard-1) + wHard_P;
end
end
end
fprintf'BM3D: First Stage End...\n');
img_basic=img_basic_sum./img_basic_weight;
if isDisplay
figure;
imshowimg_basic,[]);
title'BM3D:Fist Stage Result');
end
[block_basic,tran_block_basic,block2row_idx_basic,block2col_idx_basic] = im2blockimg_basic,kWien,pWien,lambda_distWien,0);
bn_r=floorrow-kWien)/pWien)+1;
bn_c=floorcol-kWien)/pWien)+1;
img_wien_sum=zerosrow,col);
img_wien_weight=zerosrow,col);
fprintf'BM3D: Second Stage Start...\n');
for i=1:1:bn_r
for j=1:1:bn_c
% 公式5)
[sim_blk_basic,sim_num,sim_blk_basic_idx] = search_similar_blocki,j,block_basic,tran_block_basic,floornWien/pWien),bn_r,bn_c,tauWien,NWien);
% 公式6)
tran3d_blk_basic = transform_3dsim_blk_basic,tran_mode,lambda2d_wie,lambda1d_wie);
omega_P=tran3d_blk_basic.^2)./tran3d_blk_basic.^2)+sigma2);
% 公式7)
tran3d_blk = transform_3dblock:,:,sim_blk_basic_idx),tran_mode,lambda2d_wie,lambda1d_wie);
blk_est=inv_transform_3domega_P.*tran3d_blk,tran_mode);
blk_est=realblk_est);
NWien_P=nnzomega_P);
ifNWien_P >1)
wWien_P=1/NWien_P);
else
wWien_P=1;
end
% 公式8)
for k=1:sim_num
idx=sim_blk_basic_idxk);
ir=block2row_idx_basicidx);
jr=block2col_idx_basicidx);
img_wien_sumir:ir+kWien-1,jr:jr+kWien-1) = img_wien_sumir:ir+kWien-1,jr:jr+kWien-1) + wWien_P*blk_est:,:,k);
img_wien_weightir:ir+kWien-1,jr:jr+kWien-1) = img_wien_weightir:ir+kWien-1,jr:jr+kWien-1) + wWien_P;
end
end
end
fprintf'BM3D: Second Stage End\n');
img_denoise = img_wien_sum./img_wien_weight;

function [block,transform_block,block2row_idx,block2col_idx] =im2blockimg,k,p,lambda2D,delta)
% 实现图像分块
% Inputs:
%        k: 块大小
%        p: 块移动步长
%        lambda_2D: 收缩阈值
%        delta: 收缩阈值
%  Outputs:
%        block: 返回的块
%        transform_block: 变换后的块
%        block2row_idx: 块索引与图像块的左上角行坐标对应关系
%        block2col_idx: 块索引与图像块的左上角列坐标对应关系
%
[row,col] = sizeimg);
% 频域去噪中的硬阈值,实际上原文中,对于噪声方差小于40时thres = 0, 具体见公式1)的说明第2点即距离计算)
thres = lambda2D*delta*sqrt2*logrow*col));
% r_num 和 c_num分别表示行和列上可以采集的块的数目
r_num = floorrow-k)/p)+1;
c_num = floorcol-k)/p)+1;
block = zerosk,k,r_num*c_num);
block2row_idx = [];
block2col_idx = [];
cnt = 1;
for i = 0:r_num-1
rs = 1+i*p;
for j = 0:c_num-1
cs = 1+j*p;
block:,:,cnt) = imgrs:rs+k-1,cs:cs+k-1);
block2row_idxcnt) = rs;
block2col_idxcnt) = cs;
tr_b = fft2block:,:,cnt));
idx = findabstr_b)<thres);
tr_bidx) = 0;
transform_block:,:,cnt) = tr_b;
cnt = cnt+1;
end
end
end


function [blk_est]=inv_transform_3dblk_tran3d,tran_mode)
% 3D 逆变换
% Inputs:
%       blk_tran3d: 在频域中,硬阈值滤波的图像块
%       tran_mode: 变换方法
% Outputs:
%       blk_est:
%
global blk_tran1d_s;
global blk_2d_s;
[m,n,blk_num]=sizeblk_tran3d);
blk_invtran1d=zerosm,n,blk_num);
blk_est=zerosm,n,blk_num);
iftran_mode==0)    %fft
for i=1:1:m
for j=1:1:n
blk_invtran1di,j,:)=ifftblk_tran3di,j,:));
end
end
for i=1:1:blk_num
blk_est:,:,i)=ifft2blk_invtran1d:,:,i));
end
elseiftran_mode==1)  %dct
for i=1:1:m
for j=1:1:n
blk_invtran1di,j,:)=idctblk_tran3di,j,:));
end
end
for i=1:1:blk_num
blk_est:,:,i)=idct2blk_invtran1d:,:,i));
end
elseiftran_mode==2)    %dwt
blk_num=lengthblk_2d_s);
blk_c=waverec2blk_tran3d,blk_tran1d_s,'haar');
blk_est=[];
for i=1:1:blk_num
blk_est:,:,i)=waverec2blk_c:,i),blk_2d_s{i},'Bior1.5');
end
else
error'tran_mode error');
end
end

function img_trans = other2rgbimg_in, color_mode)
% 将RGB颜色空间转为其他颜色空间
% Inputs:
%        img_in: RGB颜色空间图像
%        color_mode: 彩色图像去噪时采用的颜色空间, 默认值为0, color_mode: = 0, YUV; = 1, YCbCr; = 2, OPP
% Outputs:
%        img_trans: 其他颜色空间
%
% Author: HSW
% Date: 2018-05-06
img_trans = zerossizeimg_in));
[row, col, dims] = sizeimg_in);
if color_mode == 0
color_tran = [0.30, 0.59, 0.11; -0.15, -0.29, 0.44; 0.61, -0.51, -0.10];
color_tran_inv = invcolor_tran);
for i = 1:row
for j = 1:col
other = [img_ini, j, 1); img_ini, j, 2); img_ini, j, 3)];
img_transi, j, :) = color_tran_inv * other;
end
end
elseif color_mode == 1
color_tran = [0.30, 0.59, 0.11; -0.17, -0.33, 0.50; 0.50, -0.42, -0.08];
color_tran_inv = invcolor_tran);
for i = 1:row
for j = 1:col
other = [img_ini, j, 1); img_ini, j, 2); img_ini, j, 3)];
img_transi, j, :) = color_tran_inv * other;
end
end
elseif color_mode == 2
color_tran = [1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0; 1.0 / 2.0, 0, -1.0 / 2.0; 1.0 / 4.0, -1.0 / 2.0, 1.0 / 4.0];
color_tran_inv = invcolor_tran);
for i = 1:row
for j = 1:col
other = [img_ini, j, 1); img_ini, j, 2); img_ini, j, 3)];
img_transi, j, :) = color_tran_inv * other;
end
end
end
end

function img_trans = rgb2otherimg_in, color_mode)
% 将RGB颜色空间转为其他颜色空间
% Inputs:
%        img_in: RGB颜色空间图像
%        color_mode: 彩色图像去噪时采用的颜色空间, 默认值为0, color_mode: = 0, YUV; = 1, YCbCr; = 2, OPP
% Outputs:
%        img_trans: 其他颜色空间
%
% Author: HSW
% Date: 2018-05-06
img_trans = zerossizeimg_in));
[row, col, dims] = sizeimg_in);
if color_mode == 0
color_tran = [0.30, 0.59, 0.11; -0.15, -0.29, 0.44; 0.61, -0.51, -0.10];
for i = 1:row
for j = 1:col
rgb = [img_ini, j, 1); img_ini, j, 2); img_ini, j, 3)]; 
img_transi, j, :) = color_tran * rgb)';
end
end
elseif color_mode == 1
color_tran = [0.30, 0.59, 0.11; -0.17, -0.33, 0.50; 0.50, -0.42, -0.08];
for i = 1:row
for j = 1:col
rgb = [img_ini, j, 1); img_ini, j, 2); img_ini, j, 3)]; 
img_transi, j, :) = color_tran * rgb)';
end
end
elseif color_mode == 2
color_tran = [1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0; 1.0 / 2.0, 0, -1.0 / 2.0; 1.0 / 4.0, -1.0 / 2.0, 1.0 / 4.0];
for i = 1:row
for j = 1:col
rgb = [img_ini, j, 1); img_ini, j, 2); img_ini, j, 3)]; 
img_transi, j, :) = color_tran * rgb)';
end
end
end
end

function [sim_blk,sim_num,sim_blk_idx]=search_similar_blockik,jk,block,tran_block,np,bn_r,bn_c,tau,max_sim_num)
% 搜索相似块
% Inputs:
%       ik, jk: 待搜索相似块的索引
%       block: 图像块集合
%       tran_block: 图像块FFT硬阈值过滤后的FFT系数
%       k: 图像块大小
%       np: floornHard / pHard), 其中nHard表示图像的搜索区域大小, pHard表示块的移动步长
%       bn_r, bn_c: 图像总的行/列可以采集图像块的数目
%       tau: 图像块相似性判断阈值,见公式1)
%       max_sim_num: 最多保留相似块的数目
% Ouputs:
%       sim_blk:
%       sim_num:
%       sim_blk_idx:
%
% 搜索窗口的左上角,右下角的块索引
in_s = maxik-floornp/2),1);
jn_s = maxjk-floornp/2),1);
in_e = minik+floornp/2),bn_r);
jn_e = minjk+floornp/2),bn_c);
% 当前参考块
ref_blk = tran_block:,:,ik-1)*bn_c+jk));
ii = in_s:1:in_e;
jj = jn_s:1:jn_e;
[II,JJ] = meshgridii,jj);
IDX = II-1)*bn_c+JJ;
blk_idx=IDX:);
% 收缩范围内的全部图像块
cur_blk=tran_block:,:,blk_idx);
cnt=sizecur_blk,3);
ref_blk_mat=repmatref_blk,[1,1,cnt]);
delta_blk=cur_blk-ref_blk_mat;
dist=sumsumdelta_blk.*delta_blk,1),2);
[dist_sort,dist_idx]=sortdist);
% 最大相似块是真实相似块和目标参数相似块的最小值
max_num=mincnt,max_sim_num);
ifdist_sortmax_num)<tau)
sim_num=max_num;
else
sim_num=sumdist_sort1:max_num)<tau);
end
cnt_idx=dist_idx1:sim_num);
sim_blk_idx=blk_idxcnt_idx);
sim_blk=block:,:,sim_blk_idx);
end

function [val]=thres_shrinkdata,thres)
% 进行阈值截断: 即 datai) < thres ? datai) = 0 : datai) = datai)
% Inputs:
%       data: 阈值截断前的数据
%       thres: 阈值
% Ouputs:
%       val: 阈值截断后的数据
% 
val=data;
idx=findabsdata)<thres);
validx)=0;
end

function blk_tran3d = transform_3dblk_3d,tran_mode,lambda2d,lambda1d)
% 进行3D变换,即Collaborative Filtering: 在图像块内进行2D变换,在图像块间进行1D变换
% 公式2)
% Inputs:
%        blk_3d:
%        tran_mode:
% Ouputs:
%
global blk_tran1d_s;
global blk_2d_s;
[m,n,blk_num]=sizeblk_3d);
%变换不同时,可能需要修改??
blk_2d_shrink=zerosm,n,blk_num);
blk_1d_shrink=zerosm,n,blk_num);
iftran_mode==0)    %fft
for i=1:1:blk_num
blk_tran2d = fft2blk_3d:,:,i));
blk_2d_shrink:,:,i) = thres_shrinkblk_tran2d,lambda2d);
end
for i=1:1:m
for j=1:1:n
blk_tran1d = fftblk_2d_shrinki,j,:));
blk_1d_shrinki,j,:) = thres_shrinkblk_tran1d,lambda1d);
end
end
blk_tran3d=blk_1d_shrink;
elseiftran_mode==1)  %dct
for i=1:1:blk_num
blk_tran2d=dct2blk_3d:,:,i));
blk_2d_shrink:,:,i)=thres_shrinkblk_tran2d,lambda2d);
end
for i=1:1:m
for j=1:1:n
blk_tran1d=dctblk_2d_shrinki,j,:));
blk_1d_shrinki,j,:)=thres_shrinkblk_tran1d,lambda1d);
end
end
blk_tran3d=blk_1d_shrink;
elseiftran_mode==2)    %dwt
blk_2d_s={};
blk_2d_shrink=[];%zeros)
for i=1:1:blk_num
[blk_tran2d_c,blk_tran2d_s]=wavedec2blk_3d:,:,i),2,'Bior1.5');
blk_2d_shrink:,i)=thres_shrinkblk_tran2d_c,lambda2d);
blk_2d_s{i}=blk_tran2d_s;
end
%这里应该用 wavedec.因为是对1维??
[blk_tran1d_c,blk_tran1d_s]=wavedec2blk_2d_shrink,1,'haar');
blk_tran3d=thres_shrinkblk_tran1d_c,lambda1d);
%   elseifstrcmptran_mode,'db1')) %还未实现
%       blk_2d_s={};
%       blk_2d_shrink=[];%zeros)
%       for i=1:1:blk_num
%           [blk_tran2d_cA,blk_tran2d_cH,blk_tran2d_cV,blk_tran2d_cD]=...
%               dwt2blk_3d:,:,i),'db1');
%           blk_2d_shrink:,i)=thres_shrinkblk_tran2d_c,lambda2d);
%           blk_2d_s{i}=blk_tran2d_s;
%       end
%       [blk_tran1d_c,blk_tran1d_s]=wavedec2blk_2d_shrink,1,'haar');
%       blk_tran3d=thres_shrinkblk_tran1d_c,lambda1d);
else
error'tran_mode error');
end
end

3. 模型效果:

3.1 灰度图像

图像去噪序列——BM3D图像去噪模型实现-风君子博客图像去噪序列——BM3D图像去噪模型实现-风君子博客

图像去噪序列——BM3D图像去噪模型实现-风君子博客图像去噪序列——BM3D图像去噪模型实现-风君子博客

3.2 彩色图像

图像去噪序列——BM3D图像去噪模型实现-风君子博客图像去噪序列——BM3D图像去噪模型实现-风君子博客

图像去噪序列——BM3D图像去噪模型实现-风君子博客图像去噪序列——BM3D图像去噪模型实现-风君子博客