NLM去噪算法代码

代码：

clear
clc

%导入图片
Lenaim = imread('lenano.gif'); %直接导入有噪声的灰色图片
% Lenayyy = rgb2gray(Lenaxxx);
Lenad = im2double(Lenaim);

%设置初始参数
a = 2; %设置高斯核函数的标准差
k = 0.35;
h = k*a; %滤波系数h与标准差a正相关

%similarity square neighborhood N = 7*7
SSNx = 7;
SSNy = 7;
%search window S = 21*21
SWx = 21;
SWy = 21;

[m,n] = size(Lenad); %获取噪声图片的长与宽

%计算高斯权重
Ga = zeros(SWx-SSNx+1,SWy-SSNy+1);
for i = 1:SWx-SSNx+1
    for j = 1:SWy-SSNy+1
        Ga(i,j) = (1/(2*pi*a^2))*exp((-1/(2*a^2))*((i-1-(SWx-SSNx)/2)^2+(j-1-(SWy-SSNy)/2)^2));
    end
end

%扩充边缘
Lenanosw = zeros((m+SWx-1),(n+SWy-1));
for i = 1:m
    for j = 1:n
        Lenanosw(i+(SWy-1)/2,j+(SWx-1)/2) = Lenad(i,j);
    end
end

%NLM去噪
NL = zeros(m,n);
for i = 1+(SWy-1)/2:m+(SWy-1)/2
    for j = 1+(SWx-1)/2:n+(SWx-1)/2
        %每一个i确定一个需要比较的SSNi
        Ni = zeros(SSNx,SSNy);
        for x = 1:SSNx
            for y = 1:SSNy
                Ni(x,y) = Lenanosw((i-(SSNy-1)/2+x-1),(j-(SSNx-1)/2+y-1));
            end
        end
        %每一个i确定一个比较范围SW
        SW = zeros(SWx,SWy);
        for x = 1:SWx
            for y = 1:SWy
                SW(x,y) = Lenanosw((i-(SWy-1)/2+x-1),(j-(SWx-1)/2+y-1));
            end
        end
        %在SW中进行比较
        Z = zeros(SWx-SSNx+1,SWy-SSNy+1);
        ZZ = 0;
        for x = 1+(SSNx-1)/2:(SWx-(SSNx-1))+(SSNx-1)/2
            for y = 1+(SSNy-1)/2:(SWy-(SSNy-1))+(SSNy-1)/2
                %在SW中选一个SSNj用来与SSNi进行比较，每一个xy确定一个一个SSNj
                Nj = zeros(SSNx,SSNy);
                for xx = 1:SSNx
                    for yy = 1:SSNy
                        Nj(xx,yy) = SW((x-(SSNy-1)/2+xx-1),(y-(SSNx-1)/2+yy-1));
                    end
                end
                %计算这个SSNj与SSNi之间的欧式距离
                NN = 0;
                for xx = 1:SSNx
                    for yy = 1:SSNy
                        NN = NN + (Ni(xx,yy)-Nj(xx,yy))^2;
                    end
                end
                %计算这个SSNj与SSNi之间的Zij与归一化常数Zi
                Z(x-(SSNx-1)/2,y-(SSNy-1)/2) = exp(((-1/h^2)*Ga(x-(SSNx-1)/2,y-(SSNy-1)/2))*NN); %加上高斯加权时求出来的权重
                %Z(x-(SSNx-1)/2,y-(SSNy-1)/2) = exp((-1/h^2)*NN); %没有加上高斯加权时求出来的权重
                ZZ = ZZ + Z(x-(SSNx-1)/2,y-(SSNy-1)/2);
            end
        end
        
        W = zeros(SWx-SSNx+1,SWy-SSNy+1);
        WW = 0;
        for xx = 1:SWx-SSNx+1
            for yy = 1:SWy-SSNy+1
                W(xx,yy) = (1/ZZ)*Z(xx,yy);
                WW = WW + W(xx,yy);
            end
        end
        for xx = 1:SWx-SSNx+1
            for yy = 1:SWy-SSNy+1
                NL(i-(SWy-1)/2,j-(SWx-1)/2) = NL(i-(SWy-1)/2,j-(SWx-1)/2) + W(xx,yy)*SW(xx+(SSNx-1)/2,yy+(SSNy-1)/2);
            end
        end
    end
end

Lena = mat2gray(NL);
imshow(Lena)