function [output,params] = NLmeansFast(in, search, similarity, lambda, lambda_adaptive, eta1, eta2)

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 (fast)\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'];
params = [params, 'eta1:       ', num2str(eta1), '\n'];
params = [params, 'eta2:       ', num2str(eta2), '\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');
[mi, ni] = size(input);

lambda2 = 1 / (lambda * lambda);

[FX,FY] = gradient(in);
normgradient = sqrt(FX.*FX + FY.*FY);

sigma_nabla = compute_sigma(normgradient);

% compute mean values for gradients
mdir_x = zeros(m,n);
mdir_y = zeros(m,n);
atans = zeros(m,n);
theta  = zeros(mi*ni*(2*search_x+1)*(2*search_y+1),1);

dir_x = padarray(0.5*(input(3:mi,1:ni)-input(1:mi-2,1:ni)),[1 0]);
dir_y = padarray(0.5*(input(1:mi,3:ni)-input(1:mi,1:ni-2)),[0 1]);

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;
        
        temp_dir_x = dir_x(i-sim_x+1:i+sim_x-1, j-sim_y:j+sim_y);
        temp_dir_y = dir_y(i-sim_x:i+sim_x, j-sim_y+1:j+sim_y-1);
        
        mdir_x(temp_i,temp_j) = mean(temp_dir_x(:));
        mdir_y(temp_i,temp_j) = mean(temp_dir_y(:));
        
        atans(temp_i,temp_j) = atan2(mdir_y(temp_i,temp_j),mdir_x(temp_i,temp_j));
    end
end

x = 0;
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;
        
        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);
        
        for a=amin:1:amax
            for b=bmin:1:bmax
                if(a==i && b==j)
                    continue;
                end
                x = x+1;
                c = a - sim_x;
                d = b - sim_y;
                
                theta(x) = abs(atans(temp_i,temp_j) - atans(c,d));
                if(theta(x) > pi)
                    theta(x) = 2*pi - theta(x);
                end
            end
        end
    end
end
sigma_theta = compute_sigma(theta);

nbhnorm = zeros(m,n);
nbhmean = zeros(m,n);
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;
        
        nbh = input(i-sim_x:i+sim_x, j-sim_y:j+sim_y);
        
        nbhnorm(temp_i,temp_j) = double(norm(nbh));
        nbhmean(temp_i,temp_j) = double(mean(nbh(:)));
    end
end

x = 0;
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);
        
        % statement1: \| \nabla v(i) \| < \sigma_\nabla
        statement1 = (nbhnorm(temp_i,temp_j) < sigma_nabla);
        
        for a=amin:1:amax
            for b=bmin:1:bmax
                if(a==i && b==j)
                    continue;
                end
                c = a+1-amin;
                d = b+1-bmin;
                x = x+1;
                
                % statement2: \| \nabla v(j) \| < \sigma_\nabla
                statement2 = (nbhnorm(a-sim_x,b-sim_y) < sigma_nabla);
                
                % statement3: | \theta(i,j) | < \sigma_\theta
                statement3 = (theta(x) < sigma_theta);
                
                helper = nbhmean(temp_i,temp_j) / nbhmean(a-sim_x,b-sim_y);
                % statement4 & statement5: \eta_1 < \frac{\overline{v}(i)}{\overline{v}(j)} < \eta_2
                statement4 = (helper > eta1);
                statement5 = (helper < eta2);
                
                if((statement1 || statement2 || statement3) && (statement4 && statement5))
                    N2 = input(a-sim_x:a+sim_x, b-sim_y:b+sim_y);
                    
                    distance = sum(sum((N1-N2).*(N1-N2)));
                    
                    if(strcmp(lambda_adaptive,'y'))
                        templambda = compute_lambda(lambda,N2);
                        lambda2 = 1 / (templambda * templambda);
                    end
                
                    exps(c,d) = exp(-distance * lambda2);
                else
                    exps(c,d) = 0;
                end
            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


function [sigma] = compute_sigma(input)
temp = median(input(:));
sigma = 1.4826 * median(abs(input(:) - temp));
end