非局部均值图像去噪的原始算法

自己实现的非局部均值图像去噪的原始算法，对理解算法原理有一定帮助

function G = nlm（F win adj h）
% input: F:    the original image;
%          win: the size of window region
%          adj: the size of adjacent region
%          h:    a degree of filtering
% output G:  the output image
% BY Liu Fuxing

[mn]=size（F）;
G=zeros（mn）;

for i=1:m
    i
    for j=1:n
        % （ij） is the center point
        
        % set the size of adjacent region
        ax1=floor（adj/2）;
        ax2=floor（adj/2）;
        ay1=floor（adj/2）;
        ay2=floor（adj/2）;
        
        if j<=floor（adj/2）
            ay1=j-1;
        elseif j+floor（adj/2）>n
            ay2=n-j;
        end
        
        if i<=floor（adj/2）
            ax1=i-1;
        elseif i+floor（adj/2）>m
            ax2=m-i;
        end
        
        % set the size of window region
        wx1=floor（win/2）;
        wx2=floor（win/2）;
        wy1=floor（win/2）;
        wy2=floor（win/2）;
        
        if j<=floor（win/2）
            wy1=j-1;
        elseif j+floor（win/2）>n
            wy2=n-j;
        end
        
        if i<=floor（win/2）
            wx1=i-1;
        elseif i+floor（win/2）>m
            wx2=m-i;
        end
        
        k=0;
        weigh=0*ones（1win^2）;
        % choose an adjacent region
        for ii=i-wx1+ax1:i+wx2-ax2
            for jj=j-wy1+ay1:j+wy2-ay2
                % （iijj） is the center of certain adjacent region
                if ii-ax1>0 & ii+ax2<=m & jj-ay1>0 & jj+ay2<=n
                    w=0;
                    % calculate each element
                    for iii=ii-ax1:ii+ax2
                        for jjj=jj-ay1:jj+ay2
                            % （iiijjj） is a pixel in a certain
                            % adjacent region
                            if iii>0 & iii<=m & jjj>0 & jjj<=n
                                dist=（F（iiijjj）-F（iii-ii+ijjj-jj+j））^2;
                                w=w+exp（-dist/（h^2））;
                            end
                        end
                    end
                    k=k+1;
                    weigh（k）=w;
                end
            end
        end
        
        % calculate the weight
        total=0;
        for t=1:k
            total=total+weigh（k）;
        end
        weigh=weigh/total;
        
        % reconstruct
        t=0;
        for ii=i-wx1+ax1:i+wx2-ax2
            for jj=j-wy1+ay1:j+wy2-ay2
                if ii-ax1>0 & ii+ax2<=m & jj-ay1>0 & jj+ay2<=n
                    t=t+1;
                    G（ij）=G（ij）+F（iijj）*weigh（t）;
                end
            end
        end
        
    end
end

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       2652  2013-04-23 17:14  nlm.m

----------- ---------  ---------- -----  ----

                 2652                    1