NLEM算法图像去噪源码程序

function  imgDenoised  =   NLEM(imgNoisy, h, P, S)

[m, n] = size(imgNoisy);

N  = 2*P + 1;
h2 = h * h;

kNN = ceil ( (2*S + 1)^2 / 2 ); % top 50% of the neighbors used 

imgPad  = padarray(imgNoisy, [P P], 'symmetric');

imgDenoised = zeros(m,n);

fprintf('Looping over pixels....\n');
for i = 1 : m 
    ii = i + P;
    for j = 1 : n
        jj = j + P;
        
        patchRef = imgPad(ii - P : ii + P, jj - P : jj + P);  
        
        pmin = max(i - S, 1) + P;
        qmin = max(j - S, 1) + P;
        pmax = min(i + S, m) + P;
        qmax = min(j + S, n) + P;
        
        s1 = pmax - pmin + 1;
        s2 = qmax - qmin + 1;
        
        w  =  zeros(s1, s2);
        neighborPatches =  zeros(s1, s2, N^2);
        patchCurrent    =  zeros(2*P + 1, 2*P + 1);
        patch           =  zeros(N, N);  %#ok<*NASGU>
        
        u = 0;
        for p = pmin : pmax
            u = u + 1;
            v = 0;
            for q = qmin : qmax
                v = v + 1;
                patchCurrent = imgPad(p - P : p + P, q - P : q + P);
                d2 = sum(sum( (patchRef - patchCurrent) ...
                    .* (patchRef - patchCurrent) ));
                w(u, v)  =  exp(- d2 / h2);
                patch  =  imgPad(p - P : p + P, q - P : q + P);
                neighborPatches(u, v, :) = reshape(patch, [N^2  1]);
            end
        end
        
        % sort weights
        w  = reshape(w, [s1*s2  1]);
        [w , I] = sort(w, 'descend');
        
        L = min(kNN, s1*s2);  
        
        f       =  reshape(neighborPatches, [s1*s2   N^2])';
        median  =  findEuclideanMedian(f(: , I(1 : L)), w(1 : L));
        patch_est  =  reshape(median, [N  N]);
        imgDenoised(i, j) =  patch_est(P + 1, P + 1);
    end
end