 function [Kernel,zh] = NLM(y,R,h,ksz)
% Build NLM K and W matrices for the selected patches
% y: input image
% R: pach size
[M,N] = size(y);

% 
wrad = R; % patch radius
        
% h2 = 5; % 1.5
wsz = 2*wrad+1; % Search window

% ksz = 7; % Patch sizes (default 7)
krad = (ksz-1)/2;

img = EdgeMirror(y, [R+krad,R+krad]);


% W = zeros(wsz^2,wsz^2);
Dinv = zeros(wsz^2); 
% K = W;

G = fspecial('gaussian',[ksz ksz],1.2);
G = G(:);
G = G./sum(G);

for i=1:M
    
    for j=1:N
        
    loc_p = [i+R+krad,j+R+krad];
    
    yp = img(loc_p(1)-R :loc_p(1)+R , loc_p(2)-R :loc_p(2)+R);
    
    rowMin = loc_p(1,1) - R - krad;
    rowMax = loc_p(1,1) + R + krad;
    colMin = loc_p(1,2) - R - krad;
    colMax = loc_p(1,2) + R + krad;
    z = img(rowMin:rowMax,colMin:colMax);
    Z = im2col(z,[ksz ksz],'sliding');
    Z = Z.*repmat(G,[1 wsz^2]); % Gaussian weighted patches
    
%     for cen = 1:wsz^2
        Zc = Z(:,(wsz^2-1)/2+1);
        Zc = repmat(Zc,[1 wsz^2]);
        Ker = exp(-sum((Zc - Z).^2)/h^2);
%         d2 =  mean((Zc - Z).^2);
%         Ker = exp(-max(d2-2*sigma^2,0)/h^2);

        zh(i,j) = (Ker/sum(Ker))*yp(:);
        Kernel(:,i,j) = Ker; 
    end
    
end

 end

