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