clear all

n=256;	% Image size.

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

%%%%%%%%%%%%%%% 2D synthetic image: texture+gaussians %%%%%%%%%%%%%%%%%%%%%%%%
% Make three patches of textured part.
ix = ((-n/2):(n/2-1))' * ones(1,n);
iy = ones(n,1) * ((-n/2):(n/2-1));
ix = ix./n; iy = iy./n;
imgdct=zeros(n,n);
imgdct(1:n/4,1:n/4)=cos(2*pi*ix(1:n/4,1:n/4)*32).*cos(2*pi*iy(1:n/4,1:n/4)*32);
imgdct(1:n/4,n-n/4+1:n)=cos(2*pi*ix(1:n/4,1:n/4)*32+2*pi*iy(1:n/4,1:n/4)*32);
imgdct(n/2-n/8+1:n/2+n/8,n/2-n/8+1:n/2+n/8)=cos(2*pi*ix(1:n/4,1:n/4)*32).*sin(2*pi*iy(1:n/4,1:n/4)*16);

% Make the gaussian part.
imgwav=exp(-160.*ix.^2- 160.*iy.^2 )+exp(-160.*(ix-0.25).^2- 160.*(iy-0.25).^2 )+exp(-640.*(ix-0.25).^2-640.*(iy+0.25).^2);

% Image = part1 + part2
%imgdct=imgdct(1:2:end,1:2:end);
%imgwav=imgwav(1:2:end,1:2:end);
y0 = imgdct+imgwav;
n = length(y0);
y0 = y0(:);

% Dictionary (here local DCT + UDWT).
qmf=MakeONFilter('Symmlet',6);
dict1='UDWT2';pars11=2;pars12=qmf;pars13=0;
dict2='RealFourier2';pars21=n/8;pars22=0;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 = 15;
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 = 100;
options.method = 'fista';
options.verb  = 0;
lambda = 1.4;
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);
options.report = @(x)norm(y-K(x))^2/2+lambda*norm(x,1);
tic;[xhatlassofista,Efista] = fb(zeros(p,1), ProxF, GradG, be, options);timelassofista=toc

% Solve BP with DR.
options.gamma = 0.5;	% Stepsize parameter for Douglas-Rachford iteration.
options.tau   = 1;	% No relaxation.
options.niter = 100;
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 = 0.5;	% Stepsize parameter for Douglas-Rachford iteration.
options.tau   = 1;	% No relaxation.
options.niter = 100;
options.verb  = 0;
%epsilon = n*sigma*sqrt(1+2*sqrt(2)/n);	      
epsilon = 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 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)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,EDS] = primaldual(zeros(p,1), KK,  KKS, ProxFS, ProxG, GradH, options);timeDSPD=toc

% Solve with MCA.
tic;yhatMCA=SolveMCA(reshape(y,n,n),dict,pars1,pars2,pars3,options.niter,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;
	C = xhatDSPD((ind+1):(ind+DimDomain));
	yhatDSPD(i,:,:) = eval(['Fast' NAME, 'Synthesis(C, n, PAR1, PAR2, PAR3)']);
	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\t%-15s\n',' ','LASSO-FISTA','BP-DR','BPDN-DR','DS-PD','MCA');
fprintf('%-10s\t%-15g\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))), 		...
							      length(find(abs(xhatMCA))));

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

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


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

subplot(6,3,2); 
imagesc(imgwav);axis image;axis off
title(sprintf('Original Gaussian part'));

subplot(6,3,3); 
imagesc(imgdct);axis image;axis off
title(sprintf('Original DCT part'));

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

subplot(6,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(6,3,10); 
imagesc(squeeze(sum(yhatBPDNDR,1)));axis image;axis off
title(sprintf('BPDN-DR PSNR=%g dB',psnr(y0,sum(yhatBPDNDR,1))));

subplot(6,3,13); 
imagesc(squeeze(sum(yhatDSPD,1)));axis image;axis off
title(sprintf('DS-PD PSNR=%g dB',psnr(y0,sum(yhatDSPD,1))));

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


% Display recovered DCT component.
subplot(6,3,5); 
imagesc(squeeze(yhatlassofista(1,:,:)));axis image;axis off
title(sprintf('LASSO-FISTA PSNR=%g dB',psnr(imgwav(:),squeeze(yhatlassofista(1,:,:)))));

subplot(6,3,8); 
imagesc(squeeze(yhatBPDR(1,:,:)));axis image;axis off
title(sprintf('BP-DR (noiseless) PSNR=%g dB',psnr(imgwav(:),squeeze(yhatBPDR(1,:,:)))));

subplot(6,3,11); 
imagesc(squeeze(yhatBPDNDR(1,:,:)));axis image;axis off
title(sprintf('BPDN-DR PSNR=%g dB',psnr(imgwav(:),squeeze(yhatBPDNDR(1,:,:)))));

subplot(6,3,14); 
imagesc(squeeze(yhatDSPD(1,:,:)));axis image;axis off
title(sprintf('DS-PD PSNR=%g dB',psnr(imgwav(:),squeeze(yhatDSPD(1,:,:)))));

subplot(6,3,17); 
imagesc(squeeze(yhatMCA(1,:,:)));axis image;axis off
title(sprintf('MCA PSNR=%g dB',psnr(imgwav(:),squeeze(yhatMCA(1,:,:)))));

% Display recovered Gaussian  component.
subplot(6,3,6); 
imagesc(squeeze(yhatlassofista(2,:,:)));axis image;axis off
title(sprintf('LASSO-FISTA PSNR=%g dB',psnr(imgdct(:),squeeze(yhatlassofista(2,:,:)))));

subplot(6,3,9); 
imagesc(squeeze(yhatBPDR(2,:,:)));axis image;axis off
title(sprintf('BP-DR (noiseless) PSNR=%g dB',psnr(imgdct(:),squeeze(yhatBPDR(2,:,:)))));

subplot(6,3,12); 
imagesc(squeeze(yhatBPDNDR(2,:,:)));axis image;axis off
title(sprintf('BPDN-DR PSNR=%g dB',psnr(imgdct(:),squeeze(yhatBPDNDR(2,:,:)))));

subplot(6,3,15); 
imagesc(squeeze(yhatDSPD(2,:,:)));axis image;axis off
title(sprintf('DS-PD PSNR=%g dB',psnr(imgdct(:),squeeze(yhatDSPD(2,:,:)))));

subplot(6,3,18); 
imagesc(squeeze(yhatMCA(2,:,:)));axis image;axis off
title(sprintf('MCA PSNR=%g dB',psnr(imgdct(:),squeeze(yhatMCA(2,:,:)))));

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