以下代码转载自：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);

也就是对上图进行非局部均值滤波（这是一张脑部的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我现在还不知道该怎么理解，如果大神知道，望解惑，谢谢！