% 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)
% w OMSORT wyznaczone wcześniej miejsca punktowych pomiarów, losow dajemy 0

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

qmf = MakeONFilter('Coiflet',2);%('Haar'); %('Coiflet',2); %było 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  - nie wykorzystujemy falkowej
% rekonstrukcji
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); %pomiary losowe


yq = y;
% quantization
%q_step = 0.2;
if q_step ~= 0,
    y = y/q_step; 
    ky = round(y);
else
    ky = y;
end

if recon == 1, %nie wykorzystujemy
    if q_step ~= 0,
        y= ky*q_step;
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%% 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));
    disp(sprintf('s0 cgsolve = %5.2f', psnr(X,reshape(s0, n, n)))); %dodatek
    % 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 - %nie wykorzystujemy
    %y = ky-sn; y = y*q_step;
    if q_step ~= 0,
        y = ky*q_step - ysn;
    else
        y = ky - ysn;
        ky = y; %****
    end
    sp = Phit(y);
    Xp = reshape(sp, n, n);
    round(Xp);
    Xp(Xp>255) = 255;
    Xp(Xp<0) = 0;
else
    Xp = 0;
end