clear all

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

%%%%%%%%%%%%%%% 2D Barbara %%%%%%%%%%%%%%%%%%%%%%%%
y0 = double(imread('barbara_512x512.png'));
y0 = y0(end/2-end/4+1:end/2+end/4,end/2-end/4+1:end/2+end/4);
%y0 = y0(1:2:end,1:2:end);
n = length(y0);
y0 = y0(:);

% Dictionary (here local DCT + UDWT).
qmf=MakeONFilter('Symmlet',6);
dict1='UDWT2';pars11=2;pars12=qmf;pars13=0;
dict2='LDCT2iv';pars21=0;pars22=16;pars23=0;
dict=MakeList(dict1,dict2);
pars1=MakeList(pars11,pars21);
pars2=MakeList(pars12,pars22);
pars3=MakeList(pars13,pars23);

% Coefficient vector length.
p = SizeOfDict2(n,dict,pars1,pars2,pars3);
tightFrame = p/(n*n); % The dictionary is a tight frame with constant p/(n*n).
be = 1/tightFrame;

% Observed noisy data.
PSNR = 30;
sigma = max(abs(y0))*10^(-PSNR/20);
y = y0 + sigma*randn(size(y0));

% Solve Lasso in l1-penalized form with FISTA. 
% lambda chosen such that ||r|| <= epsilon.
options.niter = 200;
options.method = 'fista';
options.verb  = 0;
lambda = 5*sigma;
K  = @(x)reshape(FastS2(x,n,dict,pars1,pars2,pars3),n*n,1);
KS = @(x)FastA2(reshape(x,n,n),dict,pars1,pars2,pars3);
ProxF = @(x,mu)SoftThresh(x,mu*lambda);
GradG = @(x)KS(K(x) - y);
tic;[xhatlassofista,Efista] = fb(zeros(p,1), ProxF, GradG, be, options);timelassofista=toc

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

% Solve BPDN with DR.
options.gamma = 5;	% Stepsize parameter for Douglas-Rachford iteration.
options.tau   = 1;	% No relaxation.
options.niter = 200;
options.verb  = 0;
epsilon = 1.5*n*sigma;	      
ProxF = @(x,ga)SoftThresh(x,ga);
ProxG = prox_F_tightframe(@(x,ga)compute_l2ball_projection(x,epsilon),K,KS,tightFrame,y);
tic;[xhatBPDNDR,EBPDN] = dr(zeros(p,1),ProxF,ProxG,options);timeBPDNDR=toc


% Solve with MCA.
tic;yhatMCA=SolveMCA(reshape(y,n,n),dict,pars1,pars2,pars3,200,0,3*sigma,epsilon,0);timeMCA=toc

% Compute recovered signals and respective morphological components.
NumberOfDicts = LengthList(dict);
ind = 0;
for i = 1:NumberOfDicts,
	NAME = NthList(dict, i);
	PAR1 = NthList(pars1, i);
	PAR2 = NthList(pars2, i);
	PAR3 = NthList(pars3, i);
	DimDomain = SizeOfDict2(n, NAME, PAR1, PAR2, PAR3);
	C = xhatlassofista((ind+1):(ind+DimDomain));
	yhatlassofista(i,:,:) = eval(['Fast' NAME, 'Synthesis(C, n, PAR1, PAR2, PAR3)']);
	C = xhatBPDR((ind+1):(ind+DimDomain));
	yhatBPDR(i,:,:) = eval(['Fast' NAME, 'Synthesis(C, n, PAR1, PAR2, PAR3)']);
	C = xhatBPDNDR((ind+1):(ind+DimDomain));
	yhatBPDNDR(i,:,:) = eval(['Fast' NAME, 'Synthesis(C, n, PAR1, PAR2, PAR3)']);
	C = eval(['Fast' NAME, 'Analysis(squeeze(yhatMCA(i,:,:)), PAR1, PAR2, PAR3)']);
	xhatMCA((ind+1):(ind+DimDomain)) = C;
	ind = ind + DimDomain;
end

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','MCA');
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(xhatMCA))));

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

fprintf('%-10s\t%-15g\t%-15g\t%-15g\t%-15g\n','CPU (s):', timelassofista, timeBPDR, timeBPDNDR, timeMCA);fprintf('%s\n','*'*ones(1,90));


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure
% Display recovered images.
subplot(5,1,1); 
imagesc(reshape(y,n,n));axis image;axis off
title(sprintf('Noisy image PSNR_{in}=%g dB',PSNR));

subplot(5,3,4); 
imagesc(squeeze(sum(yhatlassofista,1)));axis image;axis off
title(sprintf('LASSO-FISTA PSNR=%g dB',psnr(y0,sum(yhatlassofista,1))));

subplot(5,3,7); 
imagesc(squeeze(sum(yhatBPDR,1)));axis image;axis off
title(sprintf('BP-DR (noiseless) PSNR=%g dB',psnr(y0,sum(yhatBPDR,1))));

subplot(5,3,10); 
imagesc(squeeze(sum(yhatBPDNDR,1)));axis image;axis off
title(sprintf('BPDN-DR PSNR=%g dB',psnr(y0,sum(yhatBPDNDR,1))));

subplot(5,3,13); 
imagesc(squeeze(sum(yhatMCA,1)));axis image;axis off
title(sprintf('MCA PSNR=%g dB',psnr(y0,sum(yhatMCA,1))));

% Display recovered cartoon component.
subplot(5,3,5); 
imagesc(squeeze(yhatlassofista(1,:,:)));axis image;axis off
title(sprintf('LASSO-FISTA Cartoon'));

subplot(5,3,8); 
imagesc(squeeze(yhatBPDR(1,:,:)));axis image;axis off
title(sprintf('BP-DR (noiseless) Cartoon'));

subplot(5,3,11); 
imagesc(squeeze(yhatBPDNDR(1,:,:)));axis image;axis off
title(sprintf('BPDN-DR Cartoon'));

subplot(5,3,14); 
imagesc(squeeze(yhatMCA(1,:,:)));axis image;axis off
title(sprintf('MCA Cartoon'));

% Display texture recovered component.
subplot(5,3,6); 
imagesc(squeeze(yhatlassofista(2,:,:)));axis image;axis off
title(sprintf('LASSO-FISTA Texture'));

subplot(5,3,9); 
imagesc(squeeze(yhatBPDR(2,:,:)));axis image;axis off
title(sprintf('BP-DR (noiseless) Texture'));

subplot(5,3,12); 
imagesc(squeeze(yhatBPDNDR(2,:,:)));axis image;axis off
title(sprintf('BPDN-DR Texture'));

subplot(5,3,15);
imagesc(squeeze(yhatMCA(2,:,:)));axis image;axis off
title(sprintf('MCA Texture'));

colormap('gray')

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