% Example computations related to X-ray tomography.
% Demonstration of solution by truncated SVD.
%
% The following routines must be run before this file:
% Xraytomo_01_phantom_comp.m, Xraytomo_02_matrix_comp.m and
% Xraytomo_04_SVD_comp.m.
%
% Samuli Siltanen March 2010

% Load the measured object (signal x) from disc
load data/tomo_phantom phantom N
x = phantom;

% Load measurement matrix A from disc
load data/RadonMatrix A measang
[Nrow,tmp] = size(A);
P          = Nrow/length(measang); 

% Construct noisy measurements mn
m          = A*x(:); 
m          = reshape(m,P,length(measang));
noiselevel = 0.1; 
mn         = m + noiselevel*randn(size(m));

% Take a look at the measurement
figure(1)
clf
subplot(1,2,1)
imagesc(x)
colormap gray
title('Signal x')
axis square
subplot(1,2,2)
imagesc(mn)
colormap gray
title('Noisy measurement')
axis square
drawnow
pause

% Load singular value decomposition of A
load data/RadonMatrix_SVD U D V

% Take a look at matrix A and its SVD
figure(2)
clf
subplot(1,2,1)
spy(A)
axis square
title('Nonzero elements of matrix A')
subplot(1,2,2)
semilogy(diag(D),'r','linewidth',3)
axis([1 length(diag(D)) 0 max(diag(D))])
set(gca,'xtick',[1 round(length(diag(D))/2) length(diag(D))])
axis square
title(['Singular values of matrix A'])
drawnow
pause

% Reconstruct using truncated SVD (Ns singular values included)
Ns               = 600;
[row,col]        = size(D.');
Dplus            = sparse(row,col);
svals            = diag(D);
Dplus(1:Ns,1:Ns) = diag(1./svals(1:Ns));
recn             = V*Dplus*U.'*mn(:);
recn             = reshape(recn,N,N);

% Take a look at the reconstruction
figure(3)
clf
subplot(1,2,1)
imagesc(x,[0,1])
colormap gray
title('Signal x')
axis square
axis off
subplot(1,2,2)
imagesc(recn,[0,1])
colormap gray
title(['Truncated SVD, Ns = ', num2str(Ns)])
axis square
axis off

% Take a look at singular vectors (columns of matrix V)
figure(4)
subplot(2,3,1)
imagesc(reshape(V(:,1),N,N))
title('First singular vector')
colormap gray
axis square
axis off
subplot(2,3,2)
imagesc(reshape(V(:,2),N,N))
title('Second singular vector')
colormap gray
axis square
axis off
subplot(2,3,3)
imagesc(reshape(V(:,10),N,N))
title('Tenth singular vector')
colormap gray
axis square
axis off
subplot(2,3,4)
imagesc(reshape(V(:,200),N,N))
title('Singular vector number 200')
colormap gray
axis square
axis off
subplot(2,3,5)
imagesc(reshape(V(:,end-10),N,N))
title('Tenth-to-last singular vector')
colormap gray
axis square
axis off
subplot(2,3,6)
imagesc(reshape(V(:,end),N,N))
title('Last singular vector')
colormap gray
axis square
axis off

figure(3)