NLM图像去噪算法

关于NLM图像去噪算法，效果可以，但运行效率实在太慢了，能够较好的去除噪声，但同时也损失了很多细节。记录一下，希望以后有方法能够改善其速度和细节保留的效果。

function im_rec=nlmeans_filt2D(im,sigmas,ksize,ssize,noise_std)
%%% 2D NL-means Filter. 
% im：待滤波的输入图像（灰度图）
% sigmas：高斯核的标准差：5
% ksize：滤波器的大小（以像素为单位）7
% ssize：搜索窗口的大小（以像素为单位）21
% noise_std：噪声的标准差
%==============计算滤波器和搜索窗口的一半大小===================
half_ksize=floor(ksize/2);
half_ssize=floor(ssize/2);
%%% 初始化一个扩展了搜索窗口大小的零矩阵 xm，以处理边界情况.
%%% 并将原始图像 xn 复制到 xm 的中心位置。
[M,N]=size(im);
im_plus=zeros(M+ssize-1,N+ssize-1);
im_plus(half_ssize+1:M+half_ssize,half_ssize+1:N+half_ssize)=im;
im_plus(1:half_ssize,:)=im_plus(ssize:-1:half_ssize+2,:);
im_plus(M+half_ssize+1:M+ssize-1,:)=im_plus(M+half_ssize-1:-1:M,:);
im_plus(:,1:half_ssize)=im_plus(:,ssize:-1:half_ssize+2);
im_plus(:,N+half_ssize+1:N+ssize-1)=im_plus(:,N+half_ssize-1:-1:N);
%%% 生成二维高斯核. 
gauss_win=gauss_ker2D(sigmas,ksize);
% gauss_win=ones(ksize,ksize);
%%% NL-means Filter Implementation.
%%%%%%%%%%%%%%%%计算滤波器的参数%%%%%%%%%%%%%%%%
% 较小的 filt_h 值意味着滤波器的作用范围更小，只影响邻近像素的值，可能会产生较锐利的图像细节。
% 较大的 filt_h 值意味着滤波器的作用范围更大，影响到更远处的像素值，可能导致图像边缘变得模糊。
filt_h=0.65*noise_std;

[M,N]=size(im_plus);
%%% 遍历图像中的每个像素，但不包括边界像素
%%% 在当前像素位置获取局部块和搜索窗口；计算每个搜索窗口内的权重；
%%% 用归一化权重加权求和得到滤波后的像素值，并将其存储在结果图像 im_rec 中。
for ii=half_ssize+1:M-half_ssize
   for jj=half_ssize+1:N-half_ssize
      xtemp=im_plus(ii-half_ksize:ii+half_ksize,jj-half_ksize:jj+half_ksize);
      search_win=im_plus(ii-half_ssize:ii+half_ssize,jj-half_ssize:jj+half_ssize);
      for kr=1:(ssize-ksize+1)
         for kc=1:(ssize-ksize+1)
            euclid_dist=(xtemp-search_win(kr:kr+ksize-1,kc:kc+ksize-1)).^2;
            wt_dist=gauss_win.*euclid_dist;
            sq_dist=sum(sum((wt_dist)))/(ksize^2);
            weight(kr,kc)=exp(-max(sq_dist-(2*noise_std^2),0)/filt_h^2);
         end
      end
      sum_wt=sum(sum(weight));
      weightn=weight/sum_wt;
      sum_pix=sum(sum(search_win(half_ksize+1:ssize-half_ksize,half_ksize+1:ssize-half_ksize).*weightn));
      im_rec(ii-half_ssize,jj-half_ssize)=sum_pix;
   end
end
end
function h=gauss_ker2D(sigmas,ksize)
%%% 2D Gaussian Kernel Generation.
% for x=-floor(ksize/2):floor(ksize/2)
%    for y=-floor(ksize/2):floor(ksize/2)
%       h1(x+floor(ksize/2)+1,y+floor(ksize/2)+1)=exp(-(x^2+y^2)/(2*(sigmas)^2));
%    end
% end
%% 2nd Method.
x=-floor(ksize/2):floor(ksize/2);
x2=repmat(x.^2,ksize,1);
x3=x2'+x2;
h=exp(-(x3/(2*(sigmas)^2)));
end

% 获取噪声的标准差noise_std
function noise_std=mean_noise(x)
% x:输入图像
    wname='db8'; %% db8 sym8 db16 coif5 bior6.8
    sorh='s'; % s or h or t -> trimmed
    if(size(x,3)==3)
        x=rgb2gray(x);
    end
    %%% Gaussian Noise addition.
    mean_val=0;
    noise_std=20; %% 10,20,30,40,50.
    sizeA = size(x);
    randn('seed',212096); %%% Results depend on 'seed' of the random noise.
    xn = double(x) + (noise_std*randn(sizeA)) + mean_val;
    xn = max(0,min(xn,255));

    %%% PSNR Computation of Noisy Image.
%     xn_mse=sum(sum((double(x)-double(xn)).^2))/(M*N);
%     xn_psnr=10*log10(255^2./xn_mse);
    %%% Noise Level Estimation using Robust Median Estimator.
    if(isequal(wname,'DCHWT'))
        dchw=dchwtf2(xn,1);
        tt1=dchw{1}(:)';
    else
        % ca 是逼近系数（Approximation coefficients），表示信号的低频部分。
        % ch 是水平细节系数（Horizontal detail coefficients），表示信号在水平方向上的高频细节部分。
        % cv 是垂直细节系数（Vertical detail coefficients），表示信号在垂直方向上的高频细节部分。
        % cd 是对角线细节系数（Diagonal detail coefficients），表示信号在对角线方向上的高频细节部分。
        [ca,ch,cv,cd]=dwt2(xn,wname); %二维离散小波变换
        tt1=cd(:)';
    end
    median_hh2=median(abs(tt1)); %% HH1->Subband containing finest level diagonal details.
    std_dev2=(median_hh2/0.6745);
 end