%% Tomography SVD limited data.
%
% We explore the singular values and singular vectors of the 2D Radon transform on a disk and their dependence on the data.
%  We will use AIRTools II https://github.com/jakobsj/AIRToolsII [2] and from https://github.com/jakobsj/FA18  VectorToDiskDomain.m and MatrixToDiskDomain.m etc, code for our book [1]

addpath('FA18')
addpath('AirToolsII')

%
%
% Set up parallel beam data as before

% Set the parameters for the test problem.
N = 50;           % The image is N-times-N..
theta = 0:2:178;  % Projection angles.
p = 75;           % No. of parallel rays.
% Create the test problem.
[A,b,x] = paralleltomo(N,theta,p);

% Now convert to circular image
%
%

xc= VectorToDiskDomain(x);
Ac=MatrixToDiskDomain(A);

%
% Normalized singular values on a log scale

ss=svd(full(Ac)); loglog(ss/ss(1));title('Singular values on a log log scale');

%
% Now lets see how big a matrix we can SVD to see if the singular values seem to be converging to the analytical formula. How big you can go will depend on your computer!

% Set the parameters for the test problem.
N = 100;           % The image is N-times-N..
theta = 0:1:178;  % Projection angles.
p = 150;           % No. of parallel rays.
% Create the test problem.
[A,b,x] = paralleltomo(N,theta,p);
Ac=MatrixToDiskDomain(A);
hold on
ss1=svd(full(Ac)); loglog(ss1/ss1(1)); 

%
% Now we consider limited angle data. Try 30 degrees
%

 hold off;
 As=Ac(1:(30*p),:);
ss2=svd(full(As));
loglog(ss2(1:3000)/ss2(1));title('Limited angle log log scale')
semilogy(ss2(1:3000)/ss2(1));title('Limited angle log lin scale')

% We see that with limited angle data the inverse problem is severely ill-posed. More about this in [1].
%
% Going back to the fulldata problem lets look at the singular functions (singular vectors).

% Set the parameters for the test problem.
N = 50;           % The image is N-times-N..
theta = 0:2:178;  % Projection  angles.
p = 75;           % No. of parallel rays.
% Create the test problem.
[A,b,x] = paralleltomo(N,theta,p);
xc= VectorToDiskDomain(x);
Ac=MatrixToDiskDomain(A);
[U,S,V]=svd(full(Ac));

%
%
% The V singular vectors correspond to circular images

hold off
tl = tiledlayout(4,4,TileSpacing='tight');
for i=1:16
    nexttile(tl);imagesc(VectorToRectDomain(V(:,i),N));colormap gray,axis image off
end
title(tl,'Right singular functions')

% These correspond to Zernike disk functions as detailed in [1].
%
% The left singular functions area combination of Fourier and Chebychev in sinogram space.

figure
tl1 = tiledlayout(4,4,TileSpacing='tight');
for i=1:16
    nexttile(tl1);imagesc(reshape(U(:,i),p,length(theta)));colormap gray,axis image off
end
title(tl1,'Left singular functions')

%% References
% [1] Hansen, Per Christian, Jakob Jorgensen, and William RB Lionheart, eds. Computed Tomography: Algorithms, Insight, and Just Enough Theory. Society for Industrial and Applied Mathematics, 2021. https://bookstore.siam.org/FA18
% [2] P. C. Hansen and J. S. Jorgensen, AIR Tools II: Algebraic Iterative Reconstruction Methods, Improved Implementation, Numerical Algorithms, 79 (2018), pp.  107-137. doi:10.1007/s11075-017-0430-x.