% 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;
M = 2^j;
n = M*M;

% Choose the angles for tomographic projections
Nang   = 17; % 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(M,M),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
    % Monitor the run
    if mod(iii,round(n/10))==0
        disp([iii n])
    end
end

% Compute SVD
disp('Computing SVD')
tic
[U,D,V] = svd(full(A));
disp(['SVD computed in ',num2str(toc), ' seconds'])
svals = diag(D);

% Save the result to disc
%mkdir data
savecommand = ['save -v7.3 data/RadonMatrix_', num2str(j), '_', num2str(Nang), ' A j M n Nang angles U D V svals'];
eval(savecommand)