function [ R ] = deconv_alm( B, h , b, K)
% Reference:
% http://ieeexplore.ieee.org/iel5/97/26750/01193031.pdf?arnumber=1193031#search=%22%22Adaptive%20Landweber%20Method%22%22

% Constants
n = length(h);
M = size(B,1);
N = size(B,2);


% Build blur operator
H = conv_matrix(h, B);

% Loop over all color channels
for c = 1:size(B,3)
    % Reshape pictures matrixes to stacked vector
    g = reshape(B(:,:,c), N * M,1);

    % (1) Initilalize
    % Attention: Due to Matlab restrictions iteration indexes are shifted +1
    f{1} = zeros(size(g,1),1);
    %b = 0.1;
    %K = 10000;

    % (2) First iteration
    f{2} = f{1} + b * H'*(g - H*f{1});
    % (3) Second iteration
    f{3} = f{2} + b * H'*(g - H*f{2});

    % (4)
    if norm(gradient(f{2})) <= 0
        a{3} = norm(gradient(f{3})) / norm(gradient(f{2}));
    else
        a{3} = b;
    end

    % (5)
    for k = 4 : K+1
        k
        % (a)
        f{k} = f{k-1} + max(a{k-1},b) * H'*(g - H*f{k-1});
        % (b)
        if norm(gradient(f{k-1})) <= 0
            a{k} = norm(gradient(f{k})) / norm(gradient(f{k-1}));
        else
            a{k} = b;
        end
    end


    % Convert image back to Matlab notation
    R(:,:,c) = reshape(f{k}, size(B,1), size(B,2));
end
