function [output,params] = NLmeans(in, search, similarity, lambda, lambda_adaptive)

if((search(1) < similarity(1)) || (search(2) < similarity(2)))
    error('similarity should not be bigger than search. Terminating...');
end

[ts1, ts2] = size(search);
if((ts1 ~= 1) && (ts2 ~= 2))
    error('search should be a row vector with two entries.');
end

[ts1, ts2] = size(similarity);
if((ts1 ~= 1) && (ts2 ~= 2))
    error('similarity should be a row vector with two entries.');
end

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

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

% Building string with parameters
params =          'Method:     NL-Means\n';
params = [params, 'search:     ', num2str(search), '\n'];
params = [params, 'similarity: ', num2str(similarity), '\n'];
params = [params, 'lambda:     ', num2str(lambda), '\n'];
params = [params, ' |adaptive: ', lambda_adaptive, '\n'];

search_y = search(1);
search_x = search(2);

sim_y = similarity(1);
sim_x = similarity(2);

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

lambda2 = 1 / (lambda * lambda);

for temp_i=1:m
    for temp_j=1:n
        % Jump to actual entries in the matrix
        i = temp_i + sim_x;
        j = temp_j + sim_y;
        
        sum_weights = 0;
        sum_values = 0;
        
        N1 = input(i-sim_x:i+sim_x, j-sim_y:j+sim_y);
        
        amin = max(i-search_x,sim_x+1);
        amax = min(i+search_x,sim_x+m);
        bmin = max(j-search_y,sim_y+1);
        bmax = min(j+search_y,sim_y+n);

        exps = zeros(amax-amin+1,bmax-bmin+1);
        
        for a=amin:1:amax
            for b=bmin:1:bmax
                if(a==i && b==j)
                    continue;
                end
                
                N2 = input(a-sim_x:a+sim_x, b-sim_y:b+sim_y);
                
                distance = sum(sum((N1-N2).*(N1-N2)));
                
                c = a+1-amin;
                d = b+1-bmin;
                
                if(strcmp(lambda_adaptive,'y'))
                    templambda = compute_lambda(lambda,N2);
                    lambda2 = 1 / (templambda * templambda);
                end
                
                exps(c,d) = exp(-distance * lambda2);
            end
        end
        
        Z = sum(exps(:));
        
        if (Z == 0)
            wm = exps;
        else
            wm = exps ./ Z;
        end
        
        for a=amin:1:amax
            for b=bmin:1:bmax
                c = a+1-amin;
                d = b+1-bmin;
                
                w = wm(c,d);
                sum_weights = sum_weights + w;
                sum_values = sum_values + w*input(a,b);
            end
        end
        
        % Doing the inpainting step
        if(sum_weights > 0)
            output(temp_i,temp_j) = sum_values / sum_weights;
        else
            output(temp_i,temp_j) = in(temp_i,temp_j);
        end
    end
end

end

function [output] = compute_lambda(lambda,B)
output = lambda * sqrt(std(B(:)));
end
