clear all

% Titles: 'book' (subfigures), or 'toolbox' (meaningful titles).
titles = 'toolbox';

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

% Image 
x = double(imread('Mondrian.tif'));
%x = x(1:2:end,1:2:end);
n = length(x);  % Image size n x n.
x = x(:) - mean(x(:));

% Dictionary (here DWT).
qmf=MakeONFilter('Haar');
dict='PO2';pars1=0;pars2=qmf;pars3=0;
p = SizeOfDict2(n,dict,pars1,pars2,pars3);
tightFrame = p/(n*n); 	% The tight frame constant is p/(n*n) x constant of the sensing matrix (e.g. 1 for Fourier or RST).
be = 1/tightFrame;

% Measurement (sensing) operator.
dictCS = 'Fourier';
L = 100;  % number of radial lines in the Fourier domain.
[M,Omega] = LineMask(L,n);
m = length(Omega);  	% Number of measurements m.

% Observed noisy data.
SNR = 30;
y0 = FastMeasure2D(x, dictCS, Omega);
sigma = std(y0)*10^(-SNR/20);
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 = 500;
options.verb  = 0;
lambda = sigma;
K  = @(x)FastMeasure2D(reshape(FastS2(x,n,dict,pars1,pars2,pars3),n*n,1), dictCS, Omega);
KS = @(x)FastA2(reshape(FastMeasureAdjoint2D(x,n*n,dictCS,Omega),n,n),dict,pars1,pars2,pars3);
ProxF = @(x,mu)SoftThresh(x,mu*lambda);
GradG = @(x)KS(K(x) - y);
%options.report = @(x)norm(y-K(x))^2/2+lambda*norm(x,1);
tic;xhatlassofista = real(fb(zeros(p,1), ProxF, GradG, be, options));timelassofista=toc

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

% Solve BPDN with DR.
options.gamma = 1;    % Stepsize parameter for Douglas-Rachford iteration.
options.tau   = 1;    % No relaxation.
options.niter = 500;
options.verb  = 0;
%epsilon = n*sigma*sqrt(1+2*sqrt(2/(n*n)));	      
epsilon = n*sigma/2;
ProxF = @(x,ga)SoftThresh(x,ga);
ProxG = prox_F_tightframe(@(x,ga)compute_l2ball_projection(x,epsilon),K,KS,tightFrame,y);
tic;xhatBPDNDR real(dr(zeros(p,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 = 500;
options.verb  = 0;
delta = lambda;       % Same parameter as Lasso.      
ProxG = @(x,mu)SoftThresh(x,mu);
z = KS(y);
ProxF = @(x,nu)z+compute_linfball_projection(x-z,delta);
ProxFS= prox_conjugate(ProxF);
KK  = @(x)KS(K(x));
KKS = KK;
GradH = @(x)0;
%options.report = @(x,u)norm(x,1)+sum(z.*u)+delta*norm(u,1);
tic;xhatDSPD = real(primaldual(zeros(p,1), KK,  KKS, ProxFS, ProxG, GradH, options));timeDSPD=toc

% Compute recovered images.
NAME = NthList(dict, 1);
PAR1 = NthList(pars1, 1);
PAR2 = NthList(pars2, 1);
PAR3 = NthList(pars3, 1);
C = xhatlassofista;
yhatlassofista = eval(['Fast' NAME, 'Synthesis(C(:), n, PAR1, PAR2, PAR3)']);
C = xhatBPDR;
yhatBPDR = eval(['Fast' NAME, 'Synthesis(C(:), n, PAR1, PAR2, PAR3)']);
C = xhatBPDNDR;
yhatBPDNDR = eval(['Fast' NAME, 'Synthesis(C(:), n, PAR1, PAR2, PAR3)']);
C = xhatDSPD;
yhatDSPD = eval(['Fast' NAME, 'Synthesis(C(:), n, PAR1, PAR2, PAR3)']);

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

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

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure
subplot(231)
imagesc(reshape(x,n,n));axis image;axis off
if strcmp(titles,'book'),	h=title(['(a)']);set(h,'fontsize',18);
else				title(sprintf('Original image'));
end

subplot(232)
xft = fftshift(fft2(reshape(x,n,n)));
zft = xft.*fftshift(M);
imagesc(log(1+abs(zft)));axis image;axis off
if strcmp(titles,'book'),	h=title(['(b)']);set(h,'fontsize',18);
else				title(sprintf('Projections %s m/n=%.2f SNR=%g dB',dictCS,m/(n*n),SNR));
end

subplot(233); 
imagesc(yhatlassofista);axis image;axis off
if strcmp(titles,'book'),	h=title(['(c)']);set(h,'fontsize',18);
else				title(sprintf('LASSO-FISTA PSNR=%g dB',psnr(x(:),yhatlassofista(:))));
end

subplot(234); 
imagesc(yhatBPDR);axis image;axis off
if strcmp(titles,'book'),	h=title(['(d)']);set(h,'fontsize',18);
else				title(sprintf('BP-DR (noiseless) PSNR=%g dB',psnr(x(:),yhatBPDR(:))));
end

subplot(235); 
imagesc(yhatBPDNDR);axis image;axis off
if strcmp(titles,'book'),	h=title(['(e)']);set(h,'fontsize',18);
else				title(sprintf('BPDN-DR PSNR=%g dB',psnr(x(:),yhatBPDNDR(:))));
end

subplot(236); 
imagesc(yhatDSPD);axis image;axis off
if strcmp(titles,'book'),	h=title(['(f)']);set(h,'fontsize',18);
else				title(sprintf('DS-PD PSNR=%g dB',psnr(x(:),yhatDSPD(:))));
end

colormap('gray')

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