function [kernel,zh] = BLT(y,R,hs,hr)
% Build Bilateral K and W matrices for the selected patches
% y: input image
% patches: parameters for the selected patches
[M,N] = size(y);

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

% R = patchs.R; % patch radius
% loc_p

% smoothing parameters
% hs = 20.0;
% hr = 20.0;

wsz = 2*R+1;
K = zeros(wsz^2,wsz^2);
W = K;
[x2, x1] = meshgrid(-R:R, -R:R);

    for j = 1 : M
        for k = 1 : N
            
            loc_p = [j+R,k+R];
    
            yp = img(loc_p(1)-R :loc_p(1)+R , loc_p(2)-R :loc_p(2)+R);

            % compute the radiometric and spatial distnces
            yd = yp - y(j,k);
            x1d = x1;
            x2d = x2;
            
            % compute the bilateral kernel weights
            bk = exp(-yd.^2 / (2*hr^2)) .* exp(-(x1d.^2 + x2d.^2) / (2*hs^2));
            
            % normalization
            ker = bk(:) / sum(bk(:));
            kernel(:,j,k) = ker;
            zh(j,k) =  ker'*yp(:);
        end
    end
    
