function [result] = NLmeansR(in, search, similarity, lambda)

% Determine the size of the input
[m, n] = size(in);

% Allocating memory for the result
result = zeros(m,n);

% Dummies
input = padarray(in, [similarity similarity], 'symmetric');

lambda2 = 1 / (lambda * lambda);

kernel = make_kernel(similarity);
kernel = kernel / sum(sum(kernel));

for temp_i=1:m
    for temp_j=1:n
        % Jump to actual entries in the matrix
        i = temp_i + similarity;
        j = temp_j + similarity;
        
        sum_weights = 0;
        sum_values = 0;
        wmax = 0;
        
        N1 = double(input(i-similarity:i+similarity, j-similarity:j+similarity));
        
        smin = max(i-search,similarity+1);
        smax = min(i+search,similarity+m);
        tmin = max(j-search,similarity+1);
        tmax = min(j+search,similarity+n);
        
        for s=smin:1:smax
            for t=tmin:1:tmax
                if(s==i && t==j)
                    continue;
                end
                
                N2 = double(input(s-similarity:s+similarity, t-similarity:t+similarity));
                
                distance = sum(sum(kernel.*(N1-N2).*(N1-N2)));
                
                w = exp(-distance * lambda2);
                
                if w>wmax
                    wmax = w;
                end
                
                sum_weights = sum_weights + w;
                sum_values = sum_values + w*input(s,t);
            end
        end
        
        % Stabalizing
        sum_weights = sum_weights + wmax;
        sum_values = sum_values + wmax*input(i,j);

        % Doing the inpainting step
        if(sum_weights > 0)
            result(temp_i,temp_j) = sum_values / sum_weights;
        else
            result(temp_i,temp_j) = in(temp_i,temp_j);
        end
    end
end

end

function [kernel] = make_kernel(similarity)

% Allocating memory
kernel = zeros(2*similarity+1,2*similarity+1);

for d=1:similarity
    value = 1 / (2*d+1)^2;
    for i=-d:d
        for j=-d:d
            kernel(similarity+1-i,similarity+1-j) = kernel(similarity+1-i,similarity+1-j) + value;
        end
    end
end

kernel = double(kernel ./ similarity);

end