function f = FNNMM2(g, h, alpha_ep, beta)
% function u = FNNMM2(g, h, alpha_ep, beta, m)
%
% Alternating minimization algorithm based on fitting to Df (Algorithm 2) 
% to solve the following nonconvex nonsmooth minimization problem for image
% restoration and reconstruction :
%
%   min || Hf - g ||^2 + beta * psi(f) + omega * || Df - u ||^2
%        + beta * alpha_ep * sum || u_i || 
%  
% Proposed in the paper "Fast nonconvex nonsmooth minimization methods for
% image restoration and reconstruction" 
%
% Input:    g -- observed blurry and noisy image
%           h -- convolution kernel
%           alpha_ep -- positive parameter
%           beta -- regularization parameter
%           
% Output:   u -- restored image
%           or message of failure.
%
% This program is based on ADM2TVL2.m from FTVd package provided by Rice 
% University
%
%   Copyright 2010 April
%   Mila Nikolova,  Centre de Mathematiques et de Leurs Applications
%   Michael Ng,     CMIV and Dept. Math, Hong Kong Baptist Unviersity
%   ChiPan Tam,     CMIV and Dept. Math, Hong Kong Baptist Unviersity



[m n] = size(g);

% Set omega
omega = 0.05; 
omega_increasing_rate = 1.8;
mu = 2 * omega / (beta * alpha_ep);


% Transform variables into frequency domain
h_full = zeros(m); 
[mh nh] = size(h);
h_full(1:mh,1:nh) = h;  
h_full = circshift(h_full, -round([(nh-1)/2 (mh-1)/2]));
eigHT = conj(fft2(h_full));
eigHtH = abs(eigHT).^2;
HTg = real(ifft2(eigHT.* fft2(g)));
fft2HTg = eigHT.* fft2(g);
eigDtD=abs(psf2otf([1,-1],[m m])).^2 + abs(psf2otf([1;-1],[m m])).^2;


% Initialization of variables
n = 10; ep = 0:1/n:1; f = g;
tol = 1e-4; rho = 1; maxit = 1000; 
[D,Dt] = defDDt;
Lam1 = zeros(m);
Lam2 = Lam1;
gamma = 1.618;
[D1f,D2f] = D(f);


% Algorithm 2
for k = 1:n+1

    j = 1; relerr = tol + 1; 

    while j < maxit && relerr > tol 
        
        % Shrinkage step in ref. [30]
        Z1 = D1f + Lam1/mu;
        Z2 = D2f + Lam2/mu;
        V = Z1.^2 + Z2.^2;
        V = sqrt(V);
        V(V==0) = 1;
        V = max(V - 1/mu, 0)./V;
        Y1 = Z1.*V;
        Y2 = Z2.*V;

        f_old = f;
        
        if ep(k) == 0

            % Computation of f according to ref. [30]
            f = HTg + omega*(Dt(Y1,Y2)- Dt(Lam1,Lam2)/mu);
            f = fft2(f)./(eigHtH + omega*eigDtD);
            f = real(ifft2(f));

        else

             % Quasi-Newton Method
             fft2f = fft2(f);
             gradJ = 2*eigHtH.*fft2f  - 2*fft2HTg  + 2*omega*eigDtD.*fft2f ...
                 + fft2(beta*(phiprim(sqrt(D1f.^2+D2f.^2),alpha_ep, ep(k)).*-divTV_fft(eigDtD,fft2f,D1f,D2f)) ...
                 - 2*omega*(Dt(Y1,Y2) - Dt(Lam1,Lam2)/mu));
             df = ifft2(gradJ ./ (2*eigHtH + 2*omega*eigDtD)); df = -real(df);
             f = f + rho*df;

        end

        % Compute relative error
        relerr = norm(f(:) - f_old(:),'fro')/norm(f(:),'fro');
        j = j + 1;

        % Update lambda
        [D1f,D2f] = D(f);
        Lam1 = Lam1 - gamma*mu*(Y1 - D1f);
        Lam2 = Lam2 - gamma*mu*(Y2 - D2f);
        
    end

    % Update omega
    omega = omega * omega_increasing_rate;
    mu = 2 * omega / (beta * alpha_ep);
    
    % Error message
    if j >= maxit
        fprintf('\n\nAlgorithm 2: Fail to converge after %d iterations\n', maxit);
        return
    end
    
end



%----------------------- Nested functions -----------------------------------
     function [D,Dt] = defDDt
        % defines finite difference operator D
        % and its transpose operator

        D = @(U) ForwardD(U);
        Dt = @(X,Y) Dive(X,Y);

        function [Dux,Duy] = ForwardD(U)
            % Forward finite difference operator
            Dux = [diff(U,1,2), U(:,1) - U(:,end)];
            Duy = [diff(U,1,1); U(1,:) - U(end,:)];
        end

        function DtXY = Dive(X,Y)
            % Transpose of the forward finite difference operator
            DtXY = [X(:,end) - X(:, 1), -diff(X,1,2)];
            DtXY = DtXY + [Y(end,:) - Y(1, :); -diff(Y,1,1)];
        end
     end



    function y = phiprim(x, alpha_ep, ep)
        y = -ep*alpha_ep^2*(2+ep*alpha_ep*x) .*abs(x) ./ ((1+ep*alpha_ep*x).^2);
    end



    function y = divTV_fft(eigsDtD, fft2f, D1f, D2f)
        TVf = sqrt(D1f.^2+D2f.^2);
        y = (TVf>0) .* -real(ifft2(eigsDtD.*fft2f)) ./ (TVf+(TVf==0));
    end


end