1 简介 BM3D算法借鉴了非局部均值（NL-Means）方法的非局部块匹配思想，在此基础上，将图像相似块堆叠成三维矩阵后进行协同滤波处理，再将处理结果聚合到原图像块的位置。BM3D算法由
1 简介
BM3D算法借鉴了非局部均值（NL-Means）方法的非局部块匹配思想，在此基础上，将图像相似块堆叠成三维矩阵后进行协同滤波处理，再将处理结果聚合到原图像块的位置。BM3D算法由两个大的步骤组成，即初步估计和最终估计阶段，每一阶段又分别包含三个部分：块匹配、协同滤波和聚合。
2 部分代码
%BM3D:ref:《[BM3D]Image denoising by sparse 3D transform-domain collaborative filtering》%相关符号主要借用自《An Analysis and Implementation of the BM3D Image Denoising Method》
clear;
clf;close all;
img1=imread('house.png');
% img1=imread('house2.png');
% img1=imread('peppers256.png');
% img1=imread('Cameraman256.png');
% img1=imread('Lena512.png');
% img1=imread('barbara.png');
% img1=imread('boat.png');
% img1=imread('man.png');
% img1=imread('boat.png');
% img1=imread('couple.png');
img1=double(img1);
[row,col]=size(img1);
if col>row
img=img1(:,:,1);
col=row;
else
img=img1;
end
% 噪声大小
sigma=25;
img_noise=img+sigma*randn(size(img));
figure;imshow(img_noise,[]);
psnr_noise=20*log10(255/sqrt(mean((img_noise(:)-img(:)).^2)))
%变换的方法 选择模式
tran_mode=0;%'fft'
% tran_mode=1;%'dct'
% tran_mode=2;%'dwt'
%tran_mode=3;%'db1' 还未实现
%用 strcmp还是比 直接比较数字 要慢。
tic
% first step
kHard=8;%块大小
pHard=4;%块移动间隔
lambda_distHard=0;%求相似的距离时，变换后，收缩的阈值
nHard=40;%搜索窗口大小
NHard=28;%最多相似块个数
tauHard=5000;%最大的相似距离for fft
beta=2;
% tauHard=50000;% for dct
if(tran_mode==0) %fft
lambda2d=400;
lambda1d=500;
lambda2d_wie=50;
lambda1d_wie=500;
% lambda2d=400;
% lambda1d=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
%kaiser 窗，实际代码中可能没有用到。
kaiser_win=kaiser(kHard,1)*kaiser(kHard,1)';
%图像分块，并且做变换，为找相似块做准备
[block,tran_block,block2row_idx,block2col_idx]=im2block(img_noise,kHard,pHard,lambda_distHard,0);
% block number row:行方向上的 块的个数
bn_r=floor((row-kHard)/pHard)+1;
bn_c=floor((col-kHard)/pHard)+1;
%基础估计的图像
img_basic_sum=zeros(row,col);
img_basic_weight=zeros(row,col);
%显示处理进度
is_disp_process=0;
process_step_total=bn_r/10;
process_step_cnt=0;
%基础估计
fprintf('first step\n');
%对每个块遍历
for i=1:1:bn_r
if((is_disp_process) &&(i>process_step_total))
process_step_total=process_step_total+bn_r/10;
process_step_cnt=process_step_cnt+1;
fprintf(' process:%d/10\n',process_step_cnt)
end
for j=1:1:bn_c
[sim_blk,sim_num,sim_blk_idx]=search_similar_block(i,j,block,tran_block,kHard,floor(nHard/pHard),...
bn_r,bn_c,tauHard,NHard);
%toc
% test fine_similar block
num_sim=size(sim_blk_idx,3);
%tic
%做3D变换，并且用阈值收缩
tran3d_blk_shrink=transform_3d(sim_blk,tran_mode,lambda2d,lambda1d);
%toc
NHard_P=nnz(tran3d_blk_shrink);%non_zero_num
if(NHard_P >1)
wHard_P=1/NHard_P;
else
wHard_P=1;
end
Wwin2D= kaiser(kHard, beta) * kaiser(kHard, beta)'; % Kaiser window used in the hard-thresholding part
wWien_P=Wwin2D*wHard_P;
%wHard_P=wHard_P*kaiser_win;% 要不要？
%tic
% 3D逆变换
blk_est=inv_transform_3d(tran3d_blk_shrink,tran_mode);
%pause
%可能是复数，虚部为0或接近于0，所以只取实部
%max(abs(imag(blk_est(:))))
blk_est=real(blk_est);
%toc
%tic
for k=1:sim_num
idx=sim_blk_idx(k);
ir=block2row_idx(idx);
jr=block2col_idx(idx);
%实在算不清了。利用提前保存好 块索引和左上角坐标的 对应关系
%ir=floor((idx-1)*pHard/col)+1;
%jr=(idx-1)*pHard-(i-1)*col+1;
img_basic_sum(ir:ir+kHard-1,jr:jr+kHard-1)=...
img_basic_sum(ir:ir+kHard-1,jr:jr+kHard-1)+wHard_P*blk_est(:,:,k);
img_basic_weight(ir:ir+kHard-1,jr:jr+kHard-1)=...
img_basic_weight(ir:ir+kHard-1,jr:jr+kHard-1)+wHard_P;
end
%toc
%pause
end
end
fprintf('first step end\n');
img_basic=img_basic_sum./img_basic_weight;
figure;
imshow(img_basic,[]);
psnr=20*log10(255/sqrt(mean((img_basic(:)-img(:)).^2)))
% second step
% 保持一致？
kWien=kHard;
pWien=pHard;
lambda_distWien=lambda_distHard;
nWien=nHard;%搜索窗口大小
NWien=NHard;%最多相似块个数
tauWien=tauHard;
sigma2=sigma*sigma;
[block_basic,tran_block_basic,block2row_idx_basic,block2col_idx_basic]=im2block(img_basic,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);
img_wien_weight=zeros(row,col);
process_step_total=bn_r/10;
process_step_cnt=0;
fprintf('second step\n');
for i=1:1:bn_r
if((is_disp_process) &&(i>process_step_total))
process_step_total=process_step_total+bn_r/10;
process_step_cnt=process_step_cnt+1;
fprintf(' process:%d/10\n',process_step_cnt)
end
for j=1:1:bn_c
[sim_blk_basic,sim_num,sim_blk_basic_idx]=search_similar_block(i,j,block_basic,tran_block_basic,kWien,floor(nWien/pWien),...
bn_r,bn_c,tauWien,NWien);
%对基础块进行3D变换，求得omega_P.
tran3d_blk_basic=transform_3d(sim_blk_basic,tran_mode,lambda2d_wie,lambda1d_wie);
omega_P=(tran3d_blk_basic.^2)./((tran3d_blk_basic.^2)+sigma2);
%用 基础块得到的相似块的索引，来找到 噪声块的相似块，并进行3D变换
tran3d_blk=transform_3d(block(:,:,sim_blk_basic_idx),tran_mode,lambda2d_wie,lambda1d_wie);
blk_est=inv_transform_3d(omega_P.*tran3d_blk,tran_mode);
%可能是复数，虚部为0或接近于0，所以只取实部
%max(abs(imag(blk_est(:))))
blk_est=real(blk_est);
NWien_P=nnz(omega_P); %IPOL文中8式中矩阵范数？应该如何求？？
if(NWien_P >1)
wWien_P=1/(NWien_P);
else
wWien_P=1;
end
% wWien_P=wWien_P/sigma2;
for k=1:sim_num
idx=sim_blk_basic_idx(k);
ir=block2row_idx_basic(idx);
jr=block2col_idx_basic(idx);
img_wien_sum(ir:ir+kWien-1,jr:jr+kWien-1)=...
img_wien_sum(ir:ir+kWien-1,jr:jr+kWien-1)+wWien_P*blk_est(:,:,k);
img_wien_weight(ir:ir+kWien-1,jr:jr+kWien-1)=...
img_wien_weight(ir:ir+kWien-1,jr:jr+kWien-1)+wWien_P;
end
end
end
fprintf('second step end\n');
img_wien=img_wien_sum./img_wien_weight;
figure;
imshow(img_wien,[]);
psnr=20*log10(255/sqrt(mean((img_wien(:)-img(:)).^2)))
toc
3 仿真结果
4 参考文献
[1]郑建明, and 王志伟. "BM3D的声呐图像去噪算法参数分析.".
博主简介：擅长智能优化算法、神经网络预测、信号处理、元胞自动机、图像处理、路径规划、无人机等多种领域的Matlab仿真，相关matlab代码问题可私信交流。
部分理论引用网络文献，若有侵权联系博主删除。
