【图像去噪】基于BM3D算法实现图像去噪含Matlab源码

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 stepkHard=8;%块大小pHard=4;%块移动间隔lambda_distHard=0;%求相似的距离时，变换后，收缩的阈值nHard=40;%搜索窗口大小NHard=28;%最多相似块个数tauHard=5000;%最大的相似距离for fftbeta=2;% tauHard=50000;% for dctif(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  endendfprintf('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  endendfprintf('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代码问题可私信交流。

部分理论引用网络文献，若有侵权联系博主删除。