% Construct tomographic measurement matrices
% using Matlab's radon.m and a brute-force approach.
%
% Samuli Siltanen Jan 2015

% Choose the size of the image (attenuation function)
N = 16;

% Choose the angles between 0 and 180 degrees
T  = 20;
th = ([0:(T-1)]/T)*180;

% Find out the size of the matrix A
sinog     = radon(zeros(N),th);
[row,col] = size(sinog);

% Initialize matrix
A = zeros(row*col,N^2);

% Construct A column by column
for jjj=1:N^2
    tmp      = zeros(N);
    tmp(jjj) = 1;
    sinog    = radon(tmp,th);
    A(:,jjj) = sinog(:);
    
    % Monitor the run
    if mod(jjj,100)==0
        disp([jjj,N^2])
    end
    
    % Take a look
    %     figure(1)
    %     clf
    %     subplot(1,2,1)
    %     imagesc(tmp)
    %     axis equal
    %     axis off
    %     colormap gray
    %     subplot(1,2,2)
    %     imagesc(sinog)
    %     axis equal
    %     axis off
    %     colormap gray
    %     drawnow
    
end

% Compute the singular value decomposition
[U,D,V] = svd(A);

% Save the matrix to file
savecommand = ['save data/XRmatrix_',num2str(N),'_',num2str(T),' A N T th U D V'];
eval(savecommand)
