% 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)

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

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

%clear all;

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

%reconstructed measurements
cy = ry.*quantyz;%q_step
if bzaokrag,
    y = ky;
else
    y = ky*q_step;
end
y = [cy;y];
 
 
 Phi = @(z) A_lpnlet(z, n, OM, OM_CS);
 Phit = @(z) At_lpnlet(z, n, OM, OM_CS);


%min l2 for cs ********
PPt = @(z) Phi(Phit(z));
x0 = Phit(cgsolve(PPt, y, 1e-8, 200)); %1e-8,200
% 
xp = x0; %kombinowany
disp(sprintf('x0 cgsolve = %5.2f', psnr(X,reshape(x0, n, n)))); %dodatek

% 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, Phi, Phit, y, epsilon, lbtol, mu, cgtol, cgmaxiter);
xp = tvqc_logbarrier(x0, Phi, Phit, y, epsilon, lbtol, mu, cgtol, cgmaxiter);
%xp = l1qc_logbarrier_noShow(x0, Phi, Phit, y, 0.9, lbtol, mu, cgtol, cgmaxiter);
t1 = clock;

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