% Try truncated SVD for tomography and save the results. Note that in this
% routine we are committing "inverse crime," meaning that we use the same
% model (matrix A) both for simulating data and in the reconstruction.
%
% Note that you have to run first the routine tomo01_RadonMatrix_comp.m
% with the particular values of j and Nang you are using here.
%
% Samuli Siltanen Feb 2017

% Choose the size of the unknown. The image has size (2^j)x(2^j);
% let's denote n = (2^j)*(2^j).
j = 5;

% Choose the angles for tomographic projections
Nang   = 33; 

% Load precomputed results of the routine tomo01_RadonMatrix_comp.m
% with the particular values of j and Nang above
loadcommand = ['load data/RadonMatrix_', num2str(j), '_', num2str(Nang), ' A j M n Nang angles U D V'];
eval(loadcommand)

% Evaluate the phantom
f = SqPhantom(M);

% Use radon.m to create a sinogram
% Note: we are committing "inverse crime!"
m = A*f(:);

% Use TSVD to reconstruct
svals   = diag(D);
% Truncated SVD:
r_alpha = 20;
% Naive inversion:
%r_alpha = length(svals);
Dp_alpha = zeros(size(A));
for iii = 1:r_alpha
    Dp_alpha(iii,iii) = 1/svals(iii);
end
f0 = V*(Dp_alpha.')*(U.')*m;

% Noisy data option 
%f0 = V*(Dp_alpha.')*(U.')*(m+10*randn(size(m)));

% Save results to disc
recon = reshape(f0,M,M);
savecommand = ['save data/TSVD_', num2str(j), '_', num2str(Nang), ' j M n Nang angles r_alpha recon f'];
eval(savecommand)

% Plot the results
tomo02_firstTSVD_plot