% Makes CS sensing and reconstruction considering the SVD 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 [snr1, y] = cd_l1dct(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;

[U S V] = svd(X);    %X = U*S*V',  S = U'*X*V
s = S(:);

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

% Making measurement matrix Theta
Theta = @(z) A_noiseletSVD(z, OM_CS, N, U, V);
Thetat = @(z) At_noiseletSVD(z, OM_CS, N,U, V);


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

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


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

% CS Reconstruction
 
% min l2 for cs
PPt = @(z) Theta(Thetat(z));
s0 = Thetat(cgsolve(PPt, y, 1e-8, 200));
 
% parameters for optimization code
lbtol = 918;
mu = 5;
cgtol = 1e-8;
cgmaxiter = 800;
 
% cs recovery
t0 = clock;
sp = l1qc_logbarrier_noShow(s0, Theta, Thetat, y, epsilon, lbtol, mu, cgtol, cgmaxiter);
t1 = clock;


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

Sp =reshape(sp, n, n);
Xp = U*Sp*V';

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