%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% This script is used to obtain the results in the paper
%%%% "Composite Splitting Algorithms"
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 clear; close all;
addpath('Rice_Wavelet_Toolbox_2.4');
rand('state',0); randn('state',0);
% F = double(imread('pics/Heart.png'));
F = double(imread('MRI_Coronal_Brain.jpg'));
% F = double(imread('pics/3DMR_Chest.jpg'));
% F = double(imread('pics/3DMR_Renal_Arteries.jpg'));
% F=F/255;

F=imresize(F,[256 256]);
F=im256(F); 
% figure; imshow(F, [])
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Preparing the data and operators
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for ii=1:1
    printf('ii=%d\n', ii);
[m n] = size(F); N = m*n;f = F;d = 4.4;   %%%% choose 3.3, 4, 4.4, 
[OMEGA] = RandMask_rect(double(m/d),double(n/d),m,n);
% Number of projections
k = length(OMEGA);
[mask] = RandMask_InverseTransfer(OMEGA,m,n); 
% partial FFT tranform
R = @(x) A_fhp_rect(x, OMEGA, m, n);
RT = @(x) At_fhp_rect(x, OMEGA, m, n);
R2 = @(x) A_fhp_rect2D(x, OMEGA, m, n);
RT2 = @(x) At_fhp_rect2D(x, OMEGA, m, n);
% define the function handles that compute 
% the products by W (inverse DWT) and W' (DWT)
wav = daubcqf(2,'min');W = @(x) midwt(x,wav);WT = @(x) mdwt(x,wav);
AO = @(x) R(W(x)); AOT = @(x) WT(RT(x)); 
A = A_operator(@(x) AO(x), @(x) AOT(x));
Phi = A_operator(@(x) WT(x), @(x) W(x));    % notice Phi=WT
RR = A_operator(@(x) R2(x), @(x) RT2(x));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sigma = 0.01; noise = sigma*randn(2*k+1,1); b = R(f) + noise;
data = zeros(m,n);data(1,1) = b(1); KK = length(b);
data(OMEGA) = b(2:(KK+1)/2) + i*b((KK+3)/2:KK);
im_poc = real(ifft2(data))*sqrt(m*n);
alpha = 1e-3; beta = 3.5e-2; 
xs = WT(f);input.f=f; input.xs=xs;
input.n1=m;input.n2=n;
input.alpha=alpha;input.beta=beta;
input.A=A;input.Phi=Phi;input.b=b;
input.maxitr=25000;input.maxiitr= 25000;
input.stopCriterion=2; 
input.tolA=1e-8; input.no=60; %%%%  Change no for different iteration numbers
input.L=1;input.num=1; 
input.l=min(f(:)); input.u=max(f(:));
input2=input; input2.A=RR;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% TVCMRI  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('calling the function TVCMRI.....\n');
input.f=f(:);input.xs=xs(:);
out = TVCMRI(input);

out2=out; x = out.x; im2=Phi'*x; im2=reshape(im2, [m, n]);
fprintf(1,'Iter_time = %.2fsec, Iteration Num = %d, snr = %f \n',out.iter_time, input.no, out.snr(end));
fprintf(1,'k = %d, rec_err = %10.6e, samp.ratio = %f, func.value = %f\n\n',k, norm(x(:)-xs(:))/norm(xs(:)),k/(m*n),out.funv(end));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% RecPF  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('calling the function RecPF.....\n');
opts.maxItr = input.no; opts.gamma = 1.6;opts.beta = 10;
opts.relchg_tol = 5e-4; opts.real_sol = true; 
opts.normalize = false; 
FB = fft2(F)/sqrt(m*n);
bb=FB(OMEGA)+1*randn(length(OMEGA), 1); 
pick=false(m,n);pick(OMEGA)=true;

[U,out] = RecPF(m,n,input.alpha,input.beta,pick,bb,2,opts,WT,W,range(F(:)),F);

out3=out; im3 = U; x=Phi*U(:);
fprintf(1,'Iter_time = %.2fsec, Iteration Num = %d, snr = %f \n',out.xtime(end), out.iter, out.snr(end));
fprintf(1,'k = %d, rec_err = %10.6e, samp.ratio = %f, func.value = %f\n\n', k, norm(x(:)-xs(:))/norm(xs(:)),k/(m*n),out.funv(end));


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% FCSAy 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('calling the function TV_FCSAy.....\n');

out = TVS_FCSAy(input2);

out4=out;im4 = out.y; x=Phi*im4; im4=reshape(im4, [m, n]);
fprintf(1,'Iter_time = %.2fsec, Iteration Num = %d, snr = %f \n',out.xtime(end), input.no, out.snr(end));
fprintf(1,'k = %d, rec_err = %10.6e, samp.ratio = %f, func.value = %f\n\n',...
    k, norm(x(:)-xs(:))/norm(xs(:)),k/(m*n),out.funv(end));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% FCSAm  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('calling the function NLTV_FCSA.....\n');

% profile on
out = NLTVS_FCSA(input2);
% profile viewer

out5=out;im5 = out.y; x=Phi*im5; im5=reshape(im5, [m, n]);
fprintf(1,'Iter_time = %.2fsec, Iteration Num = %d, snr = %f \n',out.xtime(end), input.no, out.snr(end));
fprintf(1,'k = %d, rec_err = %10.6e, samp.ratio = %f, func.value = %f\n\n',...
    k, norm(x(:)-xs(:))/norm(xs(:)),k/(m*n),out.funv(end));


SNR(:,ii,2)=out2.snr; 
SNR(:,ii,3)=out3.snr; 
SNR(:,ii,4)=out4.snr; 
SNR(:,ii,5)=out5.snr;

Time(:,ii,2)=out2.xtime; 
Time(:,ii,3)=out3.xtime; 
Time(:,ii,4)=out4.xtime; 
Time(:,ii,5)=out5.xtime;

Object(:,ii,2)=out2.funv; 
Object(:,ii,3)=out3.funv; 
Object(:,ii,4)=out4.funv; 
Object(:,ii,5)=out5.funv;

Rate(:,ii,2)=out2.xrate; 
Rate(:,ii,3)=out3.xrate; 
Rate(:,ii,4)=out4.xrate; 
Rate(:,ii,5)=out5.xrate;

end
SNR_M(:,:)=mean(SNR, 2); SNR_S(:,:)=std(SNR, 0, 2);
Time_M(:,:)=mean(Time, 2); Time_S(:,:)=std(Time, 0, 2);
Object_M(:,:)=mean(Object, 2); Object_S(:,:)=std(Object, 0, 2);
Rate_M(:,:)=mean(Rate, 2); Rate_S(:,:)=std(Rate, 0, 2);

no1=2; no2=0;
num=size(SNR_M,1);
start=2;
% 
% figure;subplot(2,1,1);imshow([f f-im2 f-im5 f-im7],[]);
% % imshow([f abs(f-im2) abs(f-im5) abs(f-im7)],[]);
% subplot(2,1,2);imshow([f im2 im5 im7],[]);


figure; hold on;
plot(SNR_M(start:num,2), 'gd-', 'linewidth', no1,'MarkerSize',1);
plot(SNR_M(start:num,3), 'mo--', 'linewidth', no1,'MarkerSize',1);
plot(SNR_M(start:num,4), 'bh-', 'linewidth', no1,'MarkerSize',1);
plot(SNR_M(start:num,5), 'r-', 'linewidth', no1,'MarkerSize',1);

xlabel('Iterations');ylabel('SNR');box on;
% legend('TVCMRI', 'RecPF ', 'CSA', 'FCSA','ASA', 'FASA',  'Location','SouthEast');
legend('TVCMRI', 'RecPF ', 'FCSA', 'FCSANL', 'Location','SouthEast');

textobj = findobj('type', 'text');
set(textobj, 'fontsize', 16);
h_xlabel = get(gca,'XLabel');
set(h_xlabel,'FontSize',16); 
h_xlabel = get(gca,'YLabel');
set(h_xlabel,'FontSize',16); 


figure; hold on;
plot(Time_M(start:num,2), SNR_M(start:num,2), 'g-', 'linewidth', no1,'MarkerSize',1);
plot(Time_M(start:num,3), SNR_M(start:num,3), 'mo--', 'linewidth', no1,'MarkerSize',1);
plot(Time_M(start:num,4), SNR_M(start:num,4), 'bh-.', 'linewidth', no1+0,'MarkerSize',1);
plot(Time_M(start:num,5), SNR_M(start:num,5), 'r-', 'linewidth', no1,'MarkerSize',1);

xlabel('CPU Time (s)');ylabel('SNR');box on;
% legend('TVCMRI', 'RecPF ', 'CSA', 'FCSA','ASA', 'FASA',  'Location','SouthEast');
legend('TVCMRI', 'RecPF ', 'FCSA', 'FCSANL','Location','SouthEast');
textobj = findobj('type', 'text');
set(textobj, 'fontsize', 16);
h_xlabel = get(gca,'XLabel');
set(h_xlabel,'FontSize',16); 
h_xlabel = get(gca,'YLabel');
set(h_xlabel,'FontSize',16); 
% 


% step=2;
% figure; hold on;
% % semilogy([start+step:num], Object_M(start+step:num,2), 'g-.', 'linewidth', no1);
% % plot(Object_M(start+step:num,3), 'm--', 'linewidth', no1);
% % semilogy([start+step:num], Object_M(start+step:num,4), 'b-.', 'linewidth', no1);
% semilogy([start+step:num], Object_M(start+step:num,5), 'k-', 'linewidth', no1+1);
% % semilogy([start+step:num], Object_M(start+step:num,6), 'm-', 'linewidth', no1);
% semilogy([start+step:num], Object_M(start+step:num,7), 'r-', 'linewidth', no1);
% 
% xlabel('Iterations'); ylabel('Objective');box on;
% % legend('TVCMRI', 'RecPF ', 'CSA', 'FCSA', 'ASA', 'FASA', 'Location','SouthEast');
% legend('TVCMRI', 'FCSA', 'FCSANL', 'Location','SouthEast');
% 
% 
% textobj = findobj('type', 'text');
% set(textobj, 'fontsize', 16);
% h_xlabel = get(gca,'XLabel');
% set(h_xlabel,'FontSize',16); 
% h_xlabel = get(gca,'YLabel');
% set(h_xlabel,'FontSize',16); 