% Makes CS sensing and reconstruction considering TV minimization
%
% 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 [snr1, y] = cs_tv(M, q_step, epsilon, ImgName)

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

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

X = double(imread(ImgName));
x = X(:);

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

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

q = randperm(N)';
OM_CS = q(1:M);


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

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

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


% quantization
y = y/q_step; 
y = round(y); 
y= y*q_step;
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%% CS RECONSTRUCTION %%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

 % CS Recontruction
 
 % min l2 for cs
 PPt = @(z) Phi(Phit(z));
 x0 = Phit(cgsolve(PPt, y, 1e-8, 200));
 
 % parameters for optimization code
 lbtol = 918; 
 mu = 5;
 cgtol = 1e-8;
 cgmaxiter = 800;
 
 % cs recovery
 t0 = clock;
 xp = tvqc_logbarrier_noShow(x0, Phi, Phit, y, epsilon, lbtol, mu, cgtol, cgmaxiter);
 t1 = clock;
 

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

Xp =reshape(xp, n, n);

snr1 = psnr(X,Xp);

disp(sprintf('PSNR = %5.8f', psnr(X,Xp)));

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