% Example computations related to limited-angle X-ray tomography.
% The following routines must be precomputed:
% Xraytomo_01_phantom_comp.m and Xraytomo_06_limitedmatrix_comp.m
%
% Samuli Siltanen March 2010

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

% Take a look at the signal
figure(1)
clf
imagesc(x)
colormap gray
axis('square')
axis('off')
colorbar
title('Signal x: cross section of I-beam')
drawnow
pause

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

% Take a look at matrix A
figure(2)
clf
spy(A)
axis square
title('Nonzero elements of matrix A')
pause

% Construct ideal measurement m 
m  = A*x(:);
m  = reshape(m,P,length(measang));

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Demonstration of unsuccessful naive inversion

% Take a look at the reconstruction
figure(3)
clf
subplot(1,2,1)
imagesc(x)
colormap gray
axis off
title('Signal x')
axis square
subplot(1,2,2)
imagesc(naive_recon)
colormap gray
axis off
title('Naive least squares reconstruction from ideal data')
axis square
