% 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      = 32;
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

A =full(A);

% 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)

