function [W_sym] = WM(loc_p,ksize,K_rec,R,kernel)

%%%%%%%%%

kr = (ksize - 1) / 2;


% K_rec = LARK(y,ksize,wsize);

%% Kernel Matrix
K_org = zeros((R*2+1)^2, (R*2+1)^2);
q = 1;
% ktmp = zeros(R*2+1, R*2+1);
for k = loc_p(2)-R : loc_p(2)+R
    for j = loc_p(1)-R : loc_p(1)+R
        ktmp = zeros(R*2+1, R*2+1);
        ktmp(R+1-kr:R+1+kr, R+1-kr:R+1+kr) = reshape(K_rec(:,j,k), ksize, ksize);
        ktmp = ishift2(ktmp, j-loc_p(1), k-loc_p(2)); % matrix shift with zero-padding
        K_org(q,:) = ktmp(:).';
        q = q + 1;
    end
end

%% Kernel Symmetrization
Ki = K_org; % original kernel matrix may be asymmetric, but positive and positive-definite

if kernel == 'LRK'
for k=1:2
    Ks = (Ki+Ki')/2;
    [V,D]=eig(Ks,'nobalance');
    Dpos = diag(max(diag(D),0));
    KsP=(V*Dpos*V');
    KsP= max(KsP,0); % symmetric K
    Ki = KsP;
end
end
%% Row-Stochastic (Asymmetric) Matrix
% W_asy = inv(diag(sum(Ki')))*Ki; % row-stochastic

%% Sinkhorn
A = Ki; % start from a positive symmetric kernel matrix
r0 = ones(length(A),1);
r = r0;
for k = 1:20
    c = 1./(A'*r);
    r = 1./(A*c);
end
CC = diag(c);
RR = diag(r);
W_sym = RR*A*CC; % doubly-stochastic




% 
% rec_vec = W_sym*patch(:);
% rec_patch = reshape (rec_vec,N,N);
% figure, imshow(rec_patch,[]);
% title('denoised patch');
% 
% figure, imshow(patch,[]);
% title('noisy patch');
end
