% Truncated SVD reconstruction for tomography. 
% Note: WITH inverse crime
%
% Samuli Siltanen Jan 2015

% Choose noise level
noiselevel = 0.1;

% Choose the size of the image (attenuation function)
N = 32;

% Choose the angles between 0 and 180 degrees
T  = 50;

% Load precomputed stuff
loadcommand = ['load data/XRmatrix_',num2str(N),'_',num2str(T),' A N T th U D V'];
eval(loadcommand)
svals = diag(D);

% Simulate tomographic data WITH inverse crim
ph   = phantom(N);
mIC  = A*ph(:);
mnIC = mIC + noiselevel*max(abs(mIC))*randn(size(mIC));

% Take a look at the matrices
figure(2)
clf
subplot(1,2,1)
spy(A)
axis equal
subplot(1,2,2)
semilogy([1:min(size(A))],svals,'r')
xlim([1 min(size(A))])

% Compute truncated SVD reconstruction
% Number of singular values to use
r_alpha = 200; 
% Build the matrix Dplus
Dplus = zeros(size(D.'));
Dplus(1:r_alpha,1:r_alpha) = diag(1./svals(1:r_alpha));
% Compute reconstructions
recIC    = V*Dplus*(U.')*mIC;
recIC    = reshape(recIC,N,N);
recnIC   = V*Dplus*(U.')*mnIC;
recnIC   = reshape(recnIC,N,N);


% Take a look
figure(1)
clf
imagesc([ph,recIC,recnIC])
colormap gray
axis equal
axis off
