% 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