function iml = admmmat(b,K,D,mu,mask)
    % iml - latent image
    % b - blurred image
    % K - kernal matrix
    % mu - weight of gradient sparsity prior
    % mask - for removing saturated region from data term
    [h,w,~] = size(b);
    b = b(:);
    % masking out the saturated part -- 26-Oct-2014
    b = mask*b;
    K = mask*K;
    
    %% to solve:  1/2||b-kl||_2^2 + mu||Dl||_H
    % AL: 1/(2mu)||b-kl||_2^2 + ||s||_1 + lambda'(s-Dl) + rho/2||s-Dl||_2^2
    % <==> 1/(2mu)||b-kl||_2^2 + ||s||_1 + rho/2||s-Dl+lambda/rho||_2^2
    l = b;
    s = zeros(2*h*w,1); % slack
    lambda = s; % lagrangian
    rho = 5; % penalty
    KtK = K'*K;
    DtD = D'*D;
    Ktb = K'*b;
    LHS = 1/mu*KtK + rho*DtD;
    for i = 1:20
        lold = l;
        RHS = 1/mu*Ktb + rho*D'*(s+lambda/rho);
        % l = ( LHS ) \ ( 1/mu*Ktb + rho*D'*(s+lambda/rho) ); % too slow (10x compared to cg) and require huge memory
        l = cg(l, LHS, RHS);
        % l = pcg(LHS,RHS);%,1e-3,10);
        s = ProxH( D*l-lambda/rho, rho, 0 );
        
        lambda = lambda + rho*(s-D*l);
        
        if norm(l-lold)/norm(l) < 1e-4
            break;
        end
        
    end
    
    iml = reshape(l,h,w);


end

function x = cg(B, A, b)
    k = B(:);				           % Current iterate

    r = A*k - b;
    p = -r;                        % Current conjugate
    iter = 0;

    maxit = 20;
    tol = 1e-1;
    % ---------- Begin iteration ------------------------
    while ((norm(r) > tol) && (iter < maxit))
        rr = r' * r;
        Ap = A * p;
        pAp = p' * Ap;
        alpha = rr / (pAp+0.001);    % Compute the step length

        k = k + alpha * p;       % Update x

        r = r + alpha * Ap;   % Update residual 

        rr_new = r' * r;
        rho = rr_new / (rr+0.001);  % Compute beta

        p = rho * p - r;        % Update cg vector 
        iter = iter +1;
    end
    x = k;
end

function ret = ProxH(u,s,eps)
    tmp = s*u./(1 + eps * s);
    ret = u - 1/s * (s*u./(1 + eps * s)) ./ max(1, sqrt(tmp.^2));
end





    
    
