非局部均值滤波matlab代码理解

以下代码转载自：http://blog.csdn.net/frankgoogle/article/details/52209901

本文在以上博主的分析下进一步理解，主要是对高斯加权距离的理解。非局部均值滤波的代码如下：

function  [output]=NLmeansfilter(input,t,f,h)                

[m n]=size(input); %t:搜索框半径；f:相似框半径；h:滤波频率百分比

output=zeros(m,n);

input2 = padarray(input,[f f],'symmetric');              

 %将边缘对称折叠上去  

 % f :加宽的宽度值

 kernel = make_kernel(f);                              %计算得到一个高斯核，用于后续的计算

 kernel = kernel / sum(sum(kernel));             %sum：对矩阵k的每一列求和，k表示矩阵k的总和

 h=h*h;

 for i=1:m

 for j=1:n

         i1 = i+ f;

         j1 = j+ f;      

         W1= input2(i1-f:i1+f , j1-f:j1+f);     

         wmax=0;

         average=0;

         sweight=0;

         rmin = max(i1-t,f+1);      %确定相似框的边长，其实可以在代入数据的时候就设置搜索框和边界框边长的大小关系，就不用求最大最小值了

         rmax = min(i1+t,m+f);      %这段主要是控制不超出索引值

         smin = max(j1-t,f+1);

         smax = min(j1+t,n+f);

         for r=rmin:1:rmax

         for s=smin:1:smax                 

                if(r==i1 && s==j1)

                    continue;

                end;      

                W2= input2(r-f:r+f , s-f:s+f);                

                d = sum(sum(kernel.*(W1-W2).*(W1-W2)));

                w=exp(-d/h);                 

                if w>wmax      

                   wmax=w;                   

                end

                sweight = sweight + w;

                average = average + w*input2(r,s);                                  

         end

         end

        average = average + wmax*input2(i1,j1);

        sweight = sweight + wmax;       

        if sweight > 0

            output(i,j) = average / sweight;

        else

            output(i,j) = input(i,j);

        end                

 end

 end    

 function [kernel] = make_kernel(f)            %这是一个高斯核   

kernel=zeros(2*f+1,2*f+1);   

for d=1:f    

  value= 1 / (2*d+1)^2 ;    

  for i=-d:d

  for j=-d:d

    kernel(f+1-i,f+1-j)= kernel(f+1-i,f+1-j) + value ;

  end

  end

end

kernel = kernel ./ f;   

在运行过程中输入代码：

load('IMG_Normal002.mat')

imag=Img(:,:,50);

imshow(imag,[0 255])

fima=NLmeansfilter(imag,10,3,10);

4000



也就是对上图进行非局部均值滤波（这是一张脑部的MRI图像，大小为448×448）

其中，搜索框半径设置为10，相似框半径设置为3，滤波频率设置为10.

那么，程序运行大致过程如下。先扩充原图像为545×454（各边以相似框半径扩充）；然后以像素点p（4，4）为中心，构成一个7×7的相似框W1，设置搜索框范围为(行：4-14，列：4-14），再以点p1（4，5）为中心，构成同样大小的相似框W2。然后就这两个相似框的高斯加权距离d = sum(sum(kernel.*(W1-W2).*(W1-W2)));

这里为什么叫高斯加权距离呢？是因为(W1-W2).*(W1-W2）计算出两个相似框对应像素差值的平方W3，在用高斯核进行卷积，相当于把W3中的各像素点按高斯模板赋予权重（离中心点近的权重大些，远的就权重小些），然后再求和得到高斯加权距离。而普通的欧氏距离就是d = sum(sum((W1-W2).*(W1-W2)));

同理，循环进行，最后得到结果如下



至于h我现在还不知道该怎么理解，如果大神知道，望解惑，谢谢！