clear all

% Save figures: 1/0 (default 0).
savfig = 0;

% Phantom 
load phantom;
%x = x(1:2:end,1:2:end); % Subsample for a faster demo. Comment if undesired.
n = length(x);		% Image width.

% Measurement (sensing) operator: Fourier samples along radial lines.
dict = 'Fourier';
L = 22;  			% Number of radial lines in the Fourier domain.
[M,Omega] = LineMask(L,n);
m = length(Omega);  		% Number of measurements m.
tightFrame = 1;
be = 1/tightFrame;

% Observed noisy data.
sigma = 0.01;
y0 = FastMeasure2D(x(:), dict, Omega);
SNR = 20*log10(std(y0)/sigma);
y = y0 + sigma*randn(size(y0));

% Solve Lasso in l1-penalized form with FISTA. 
% lambda chosen such that ||r|| <= epsilon.
options.mu    = be;		% Descent step-size for FB.
options.method = 'fista';
options.niter = 100;
options.verb  = 0;
lambda = 2E-3;
K  = @(x)FastMeasure2D(x, dict, Omega);
KS = @(x)FastMeasureAdjoint2D(x, n*n, dict, Omega);
ProxF = @(x,mu)reshape(perform_proxexact_anisotv2D(reshape(x,n,n),mu*lambda),n*n,1);
GradG = @(x)KS(K(x) - y);
tvnorm= @(x)sum(sum(sqrt(sum(grad(reshape(x,n,n)).^2,3))));
%options.report = @(x)norm(y-K(x))^2/2+lambda*tv(x);
tic;[xhatlassofista,Efista] = fb(zeros(n*n,1), ProxF, GradG, be, options);timelassofista=toc

% Solve BP with DR.
options.gamma = 0.1;    % Stepsize parameter for Douglas-Rachford iteration.
options.tau   = 1;    % No relaxation.
options.niter = 100;
options.verb  = 0;
ProxF = @(x,ga)reshape(perform_proxexact_anisotv2D(reshape(x,n,n),ga),n*n,1);
ProxG = @(x,ga)x + KS(y0 - K(x))/tightFrame;
tic;xhatBPDR = real(dr(zeros(n*n,1),ProxF,ProxG,options));timeBPDR=toc

% Solve BPDN with DR.
options.gamma = 0.1;    % Stepsize parameter for Douglas-Rachford iteration.
options.tau   = 1;    % No relaxation.
options.niter = 100;
options.verb  = 0;
epsilon = n*sigma/3;	      
%epsilon = n*sigma;
ProxF = @(x,ga)reshape(perform_proxexact_anisotv2D(reshape(x,n,n),ga),n*n,1);
ProxG = prox_F_tightframe(@(x,ga)compute_l2ball_projection(x,epsilon),K,KS,tightFrame,y);
tic;xhatBPDNDR = real(dr(zeros(n*n,1),ProxF,ProxG,options));timeBPDNDR=toc

% Solve Dantzig-Selector with PD.
options.mu    = .99*be;      % Stepsize on the dual.
options.nu    = .99*be;      % Stepsize on the primal.
options.tau   = 1;	     % No relaxation.
options.niter = 100;
options.verb  = 0;
delta = lambda;       % Same parameter as Lasso.      
ProxG = @(x,mu)reshape(perform_proxexact_anisotv2D(reshape(x,n,n),mu),n*n,1);
z = KS(y);
ProxF = @(x,nu)z+compute_linfball_projection(x-z,delta);
ProxFS= prox_conjugate(ProxF);
K  = @(x)FastMeasureAdjoint2D(FastMeasure2D(x, dict, Omega), n*n, dict, Omega);
KS = K;
GradH = @(x)0;
%options.report = @(x,u)tv(x)+sum(z.*u)+delta*norm(u,1);
tic;xhatDSPD = real(primaldual(zeros(n*n,1), K,  KS, ProxFS, ProxG, GradH, options));timeDSPD=toc


fprintf('%s\n','*'*ones(1,90));
fprintf('%40sSummary\n',' ');
fprintf('%s\n','*'*ones(1,90));
fprintf('%-10s%-15s\t%-15s\t%-15s\t%-15s\t%-15s\n',' ','Original','LASSO-FISTA','BP-DR','BPDN-DR','DS-PD');
fprintf('%-10s%-15g\t%-15g\t%-15g\t%-15g\t%-15g\n','||x||_0:',length(find(abs(grad(reshape(x,n,n))))),			...
							      length(find(abs(grad(reshape(xhatlassofista,n,n))))),	...
							      length(find(abs(grad(reshape(xhatBPDR,n,n))))),		...
							      length(find(abs(grad(reshape(xhatBPDNDR,n,n))))),		...
							      length(find(abs(grad(reshape(xhatDSPD,n,n))))));

fprintf('%-10s%-15g\t%-15g\t%-15g\t%-15g\t%-15g\n','||x||_TV:',tvnorm(x),		...
							       tvnorm(xhatlassofista), 	...
							       tvnorm(xhatBPDR), 	...
							       tvnorm(xhatBPDNDR),	...
							       tvnorm(xhatDSPD));

fprintf('%-25s\t%-15g\t%-15g\t%-15g\t%-15g\n','CPU (s):', timelassofista, timeBPDR, timeBPDNDR, timeDSPD);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure
subplot(231)
yft = zeros(n,n);
yft(find(M)) = y;
imagesc(log(1+abs(fftshift(yft))));axis image;axis off
title(sprintf('Fourier measurements SNR=%g dB',SNR));

subplot(232)
imagesc(reshape(x,n,n));axis image;axis off
title(sprintf('Original image'));

subplot(233)
imagesc(reshape(xhatlassofista,n,n));axis image;axis off
title(sprintf('LASSO-FISTA PSNR=%g dB',psnr(x(:),xhatlassofista(:))));

subplot(234)
imagesc(reshape(xhatBPDR,n,n));axis image;axis off
title(sprintf('BP-DR (noiseless) PSNR=%g dB',psnr(x(:),xhatBPDR(:))));

subplot(235)
imagesc(reshape(xhatBPDNDR,n,n));axis image;axis off
title(sprintf('BPDN-DR PSNR=%g dB',psnr(x(:),xhatBPDNDR(:))));

subplot(236)
imagesc(reshape(xhatDSPD,n,n));axis image;axis off
title(sprintf('DS-PD PSNR=%g dB',psnr(x(:),xhatDSPD(:))));

colormap('gray')

if savfig
  saveas(gcf,'2D/Datasets/testsCSFourierTVSheppLogan2D.fig','fig');
  saveas(gcf,'2D/Datasets/testsCSFourierTVSheppLogan2D.eps','epsc');
end