% Build a matrix for parallel beam tomographic measurement.
% This is a simple method based on the use of Matlab's radon.m
% routine. There surely exist more effective ways to do this.
%
% 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;
n = 2^(2*j);
% Choose the angles for tomographic projections
Nang = 33; % odd number is preferred
angles = [0:(Nang-1)]*360/Nang;
% Let's find out the eventual size of matrix A. Namely, the size is
% k x n, where n is as above and k is the number of X-rays. The number of
% X-rays depends on Matlab's inner workings in routine radon.m.
tmp = radon(zeros(2^j,2^j),angles);
k = length(tmp(:));
% Initialize matrix A as sparse matrix
A = sparse(k,n);
% Construct matrix A column by column. This is a computationally
% wasteful method, but we'll just do it.
for iii = 1:n
% Construct iii'th unit vector
unitvec = zeros(2^j,2^j);
unitvec(iii) = 1;
% Apply radon.m
tmp = radon(unitvec,angles);
A(:,iii) = sparse(tmp(:));
% Take a look
figure(1)
clf
subplot(1,3,1)
imagesc(unitvec)
axis equal
title(['Unit vector ',num2str(iii)],'fontsize',16)
subplot(1,3,2)
imagesc(tmp)
axis equal
title('Current sinogram','fontsize',16)
subplot(1,3,3)
spy(A)
drawnow
title('Nonzero elements in A','fontsize',16)
pause
end