% Makes CS sensing and reconstruction considering the DWT as the sparse domain
%
% Inputs:
%	M: Number of measurements
%	q_step: quantization step
%	epsilon: optimization parameter
%	ImgName: image file name 
%
% Writen by: Adriana Schulz, Visgraf - IMPA aschulz@impa.br
% Last Update: January 2009

function [Xp,ky,OM_CS,yq] = cs_l1wavart(M, q_step, epsilon,X,recon,losow,OMSORT,ysn)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%  INITIALIZATION  %%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%addpath ../Measurements
%addpath ../Optimization
%addpath ../Utils
%addpath ../Data

%clear all;

% IMAGE COMPRESSED
x = X(:);
n = size(X,1);
N = n*n;

qmf = MakeONFilter('Coiflet',2);
%S = FWT2_PO(X,3,qmf);
%s = S(:);

% for repeatable experiments
if losow == 1,
    load RANDOM_STATES
    rand('state', rand_state);
    randn('state', randn_state);

    % NOISELET MEASUREMENTS -- for random projection, avoid mean
    q = randperm(N)';
    OM_CS = q(1:M);
else
    OM_CS = OMSORT;
end

% Making measurement matrix Phi
Phi = @(z) A_noiselet(z, OM_CS);
Phit = @(z) At_noiselet(z, OM_CS, N);

% Making measurement matrix Theta
Theta = @(z) A_noiseletWLOrth(z, OM_CS, N, qmf );
Thetat = @(z) At_noiseletWLOrth(z, OM_CS, N, qmf );


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%% SENSING THE IMAGE %%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

y = Phi(x);  % y = Theta(s);


yq = y;
% quantization
%q_step = 0.2;
y = y/q_step; 
ky = round(y); 

if recon == 1,

    y= ky*q_step;

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%% CS RECONSTRUCTION %%%%%%%%%%% 
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

    % CS Reconstruction

    % min l2 for cs
    % PPt = @(z) Theta(Thetat(z));
    % s0 = Thetat(cgsolve(PPt, y, 1e-8, 200));
    PPt = @(z) Phi(Phit(z));
    s0 = Phit(cgsolve(PPt, y, 1e-8, 200));

    % parameters for optimization code
    lbtol = 918;
    mu = 5;
    cgtol = 1e-8;
    cgmaxiter = 800;

    % cs recovery

    %epsilon = 1e-3*norm(y);
    %epsilon = 1.5;
    t0 = clock;
    %sp = l1qc_logbarrier_noShow(s0, Theta, Thetat, y, epsilon, lbtol, mu, cgtol, cgmaxiter);
    sp = tvqc_logbarrier_noShow(s0, Phi, Phit, y, epsilon, lbtol, mu, cgtol, cgmaxiter);%x0, Phi, Phit, y, epsilon, lbtol, mu, cgtol, cgmaxiter);
    t1 = clock;

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%% DISPLAYING RESULT %%%%%%%%%%% 
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    Sp =reshape(sp, n, n);
    Xp = Sp;% IWT2_PO(Sp,3,qmf);

    round(Xp);
    Xp(Xp>255) = 255;
    Xp(Xp<0) = 0;
elseif recon == 2, %rekonstrukcja odwrotną noiselet
    %y = ky-sn; y = y*q_step;    
    y = ky*q_step - ysn;
    sp = Phit(y);
    Xp = reshape(sp, n, n);
    round(Xp);
    Xp(Xp>255) = 255;
    Xp(Xp<0) = 0;
else
    Xp = 0;
end
% % % disp(sprintf('PSNR = %5.8f', psnr(X,Xp)));
% % % 
% % % pol = zeros(1,N);
% % % pol(OM_CS) = 1;
% % % ent_pol = entropy(pol);
% % % disp(sprintf('Entropy (pomiary losowe, polozenia): %d %d',M*entropy(ky)/N,ent_pol));
% % % disp(sprintf('Elapsed time = %8.2f seconds', etime(t1,t0)));
% % % disp(sprintf('Epsilon = %5.2f Meas Number = %d, qerror = %5.2f',epsilon,M,q_step));
%disp(sprintf('Elapsed time = %8.2f seconds', etime(t1,t0)));

%filename = ['l1DCT_noS_' ImgName '_M' int2str(M) '_Ps' int2str(q_step*100) '_Ep'  int2str(epsilon*1000)  '.mat'];
%save (filename, 'snr1', 'Xp', 'M', 'q_step', 'y', 'epsilon', 'ImgName',
%'q', 's0');