% 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 [XX] = cs_l1dctrecon(ry, ky, quantyz, q_step, epsilon, X, OM, OM_CS,bzaokrag)
% quantyz i  q_step są ustawione na 0 - nieważna opcja w tym momencie
qmf = MakeONFilter('Coiflet',2);%('Haar');%
x = X(:);
n = size(X,1);
N = n*n;

%reconstructed measurements
if q_step ~= 0
    cy = ry.*quantyz;%q_step
    if bzaokrag,
        y = ky;
    else
        y = ky*q_step;
    end
else %bez kwantyzacji pomiarów
    cy = ry;
    y = ky;
end
if ry ~= 0, %są pomiary kosinusowe, więc do rekonstrukcji będziemy używać łacznych procedur rekonstrukcji pomiarów kosinusowych i punktowych
    y = [cy;y];
    Phi = @(z) A_lpnlet(z, n, OM, OM_CS);
    Phit = @(z) At_lpnlet(z, n, OM, OM_CS);
else %tylko punktowe pomiary są rekonstruowane
    % Making measurement matrix Phi
    Phi = @(z) A_noiselet(z, OM_CS);
    Phit = @(z) At_noiselet(z, OM_CS, N);
end
 
 
%  Phi = @(z) A_lpnlet(z, n, OM, OM_CS);
%  Phit = @(z) At_lpnlet(z, n, OM, OM_CS);
 
%Making measurement matrix Theta - nieużywane
%Theta = @(z) A_noiseletWLOrth(z, OM_CS, N, qmf );
%Thetat = @(z) At_noiseletWLOrth(z, OM_CS, N, qmf );
%Theta = @(z) A_noiseletWLOrth(z, [OM;OM_CS], N, qmf );
%Thetat = @(z) At_noiseletWLOrth(z,[OM;OM_CS], N, qmf );
Theta = @(z) A_noiseletWLOrthart(z, OM, OM_CS, n, qmf );
Thetat = @(z) At_noiseletWLOrthart(z, OM, OM_CS, n, qmf );



%min l2 for cs ******** wstępna rekonstrukcja z min normy l2
%PPt = @(z) Phi(Phit(z));
%x0 = Phit(cgsolve(PPt, y, 1e-8, 200)); %1e-8,200
 PPt = @(z) Theta(Thetat(z));
 x0 = Thetat(cgsolve(PPt, y, 1e-8, 200));
% 
xp = x0; %kombinowany
S = FWT2_PO(reshape(x,n,n),3,qmf);
s = S(:);
disp(sprintf('x0 cgsolve = %5.2f', psnr(s,x0))); %dodatek niesitotny - jakość wstępnej, czasowo bardzo korzystnej rekonstrukcji, ale często nienajlepszej
%disp(sprintf('x0 cgsolve = %5.2f', psnr(X,reshape(x0, n, n)))); %dodatek niesitotny - jakość wstępnej, czasowo bardzo korzystnej rekonstrukcji, ale często nienajlepszej

% parameters for optimization code
lbtol = 918;     % 1e-3*tv(X);
mu = 5;
cgtol = 1e-8;
cgmaxiter = 800;


% cs recovery ******
%epsilon = 1e-3*norm(y); %3
t0 = clock;
xp = l1qc_logbarrier_noShow(x0, Theta, Thetat, y, 0.9, lbtol, mu, cgtol, cgmaxiter);
%xp = tvqc_logbarrier_noShow(x0, Theta, Thetat, y, epsilon, lbtol, mu, cgtol, cgmaxiter); %próba z falkami
% % Sp =reshape(sp, n, n);
% % XX = IWT2_PO(Sp,3,qmf);

%xp = tvqc_logbarrier_noShow(x0, Phi, Phit, y, epsilon, lbtol, mu, cgtol, cgmaxiter); %zasadnicza procedura rekonstrukcji total variance, lepiej się sprawdza
%xp = l1qc_logbarrier_noShow(x0, Phi, Phit, y, 0.9, lbtol, mu, cgtol, cgmaxiter);
XX = reshape(xp, n, n);

XX(XX>255) = 255;
XX(XX<0) = 0;
% 
t1 = clock;
disp(sprintf('Elapsed time = %8.2f seconds', etime(t1,t0)));