%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    Compressive Sensing Reconstruction Example (2-d imaging example)    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Coded by Kun Qiu (kqiu@iastate.edu)
%Last updated (Dec. 29, 2009)



clear all
clc
close all

path(path, './subfunctions');
path(path,'./Algorithms');


RadialNum=30;                           %number of radial lines

%Generate phantom image
Img2D=phantom(128);
Img1D=Img2D(:);
[my,mx]=size(Img2D);
m=my*mx;                                    %number of underlying signal length

% number of radial lines in the Fourier domain
[M,Mh,mh,mhi]=LineMask(RadialNum,mx);
OMEGA=mhi;
smpl_pattern=fftshift(M);
%Sampling operator Phi
Phi=@(z) A_fhp(z,OMEGA);
Phit=@(z) At_fhp(z,OMEGA,mx);
%The sparsifying basis (Haar wavelet)
wav=daubcqf(2);
dwt_L=6;
Psi=@(z) midwt(z,wav,dwt_L);
Psit=@(z) mdwt(z,wav,dwt_L);
%Effective sensing operator H
H=@(z) H_idwt1d(Phi,z,wav,dwt_L,my,mx);
Ht=@(z) Ht_dwt1d(Phit,z,wav,dwt_L,my,mx);

% taking measurements
y=Phi(Img1D);
N=length(y);                                %number of measurements
Nratio=N/m;

phantom_coeff=Psit(Img2D);
r_init=length(find(abs(phantom_coeff(:))>0));

%solve by Back Projection
s_BackProj=Ht(y);
s_BackProj=reshape(s_BackProj,[my mx]);
Img2D_BackProj=Psi(s_BackProj);
PSNR_BackProj=psnr(Img2D,Img2D_BackProj);

%solve by IHTls
tic;
[s_IHT,Count_IHT]=hard_l0_Mterm_mod(y,H,m,r_init,1e-14,'P_trans',Ht,'step_size',1);
t_IHT=toc;
s_IHT=reshape(s_IHT,[my mx]);
Img2D_IHT=Psi(s_IHT);
PSNR_IHT=psnr(Img2D,Img2D_IHT);

%solve by NIHTls
tic;
[s_NIHT,Count_NIHT]=hard_l0_Mterm_mod(y,H,m,r_init,1e-14,'P_trans',Ht);
t_NIHT=toc;
s_NIHT=reshape(s_NIHT,[my mx]);
Img2D_NIHT=Psi(s_NIHT);
PSNR_NIHT=psnr(Img2D,Img2D_NIHT);

%solve by GPSR_BB
tic;
tau=0.001*max(abs(Ht(y)));
[s_GPSR,s_GPSRdb,objective,times,debias_start,mses]=GPSR_BB(y,H,tau,'AT',Ht,'Debias',1,'ToleranceA',1e-5,'Verbose',0);
if isempty(s_GPSRdb)
    s_GPSRdb=s_GPSR;
end
t_GPSR=toc;
s_GPSRdb=reshape(s_GPSRdb,[my mx]);
Img2D_GPSRdb=Psi(s_GPSRdb);
PSNR_GPSRdb=psnr(Img2D,Img2D_GPSRdb);


%solve by ExCoVls
tic;
[s_ExCoV,sigma2_hat,A_index_ExCoV,Count_ExCoV]=ExCoVls(H,Ht,y,'Visibility',1);
t_ExCoV=toc;
s_ExCoV=reshape(s_ExCoV,[my mx]);
Img2D_ExCoV=Psi(s_ExCoV);
PSNR_ExCoV=psnr(Img2D,Img2D_ExCoV);

%Plotting
figure
subplot(2,3,1)
imagesc(Img2D)
colormap(gray)
caxis([0 1]);
axis off
title('Original image');
subplot(2,3,2)
imagesc(smpl_pattern)
colormap(gray)
caxis([0 1]);
axis off
title(['Sampling Pattern']);
subplot(2,3,3)
imagesc(Img2D_BackProj)
colormap(gray)
caxis([0 1]);
axis off
title(['Back projectiom recovery (PSNR=',num2str(PSNR_BackProj),')']);
subplot(2,3,4)
imagesc(Img2D_NIHT)
colormap(gray)
caxis([0 1]);
axis off
title(['NIHT recovery (PSNR=',num2str(PSNR_NIHT),')']);
subplot(2,3,5)
imagesc(Img2D_GPSRdb)
colormap(gray)
caxis([0 1]);
axis off
title(['GPSR-DB recovery (PSNR=',num2str(PSNR_GPSRdb),')']);
subplot(2,3,6)
imagesc(Img2D_ExCoV)
colormap(gray)
caxis([0 1]);
axis off
title(['ExCoV recovery (PSNR=',num2str(PSNR_ExCoV),')']);


save(['Phantom_',num2str(RadialNum),'Radials.mat']);

