% Here we construct a tomographic X-ray measurement model using a sparse
% matrix A. (This is not very effective numerically, but it is used only
% for demonstration purposes.)
%
% Our specific example cases are N=16,32,50,64 and 128. The computed
% matrices are saved to a subfolder called data.
%
% Remark: this routine is computationally demanding with N greater than 64.
%
% Samuli Siltanen February 2012
% Create a subdirectory called 'data'. If it already exists, Matlab will
% show a warning. You don't need to care about the warning.
mkdir('.','data')
% Construct phantom. You can modify the resolution parameter N.
N = 16;
target = phantom('Modified Shepp-Logan',N);
% Choose measurement angles (given in degrees, not radians).
Nang = N;
angle0 = -90;
measang = angle0 + [0:(Nang-1)]/Nang*180;
% Initialize measurement matrix of size (M*P) x N^2, where M is the number of
% X-ray directions and P is the number of pixels that Matlab's Radon
% function gives.
P = length(radon(target,0));
M = length(measang);
A = sparse(M*P,N^2);
% Construct measurement matrix column by column. The trick is to construct
% targets with elements all 0 except for one element that equals 1.
for mmm = 1:M
for iii = 1:N^2
tmpvec = zeros(N^2,1);
tmpvec(iii) = 1;
A((mmm-1)*P+(1:P),iii) = radon(reshape(tmpvec,N,N),measang(mmm));
if mod(iii,100)==0
disp([mmm, M, iii, N^2])
end
end
end
% Test the result
Rtemp = radon(target,measang);
Rtemp = Rtemp(:);
Mtemp = A*target(:);
disp(['If this number is small, then the matrix A is OK: ', num2str(max(max(abs(Mtemp-Rtemp))))]);
% Save the result to file (with filename containing the resolution N)
savecommand = ['save data/RadonMatrix', num2str(N), ' A measang target N P Nang'];
eval(savecommand)