%% Tikhonov and generalized Tikhonov of the Radon transform in 2D
% For this tutorial you will need to download AIRTools II from https://github.com/jakobsj/AIRToolsII  and run airtoolsII_setup.m first. Airtools is a simple Matlab library to explore the algebraic approach to 2D radon transform (ie tomography) problems in Matlab. Described also in [2].
%
% Set up a parallel beam tomography problem with AIRTools II

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

%
% Show the original "phantom" image

imagesc(reshape(x,N,N)); colormap gray, axis image off
title('Test image')

%
% For this example with no noise backslash (least squares) solution works well

imagesc(reshape(A\b,N,N));axis image off

% Lets look at the singular values to see why that is by plotting the singular values on a log scale
%

s=svd(full(A)); loglog(s);

% As this is straight until it falls dramatically we conclude that this problem is mildly ill-posed.
%
% Now add some noise.

bn = b+ 0.001* randn(size(b))*norm(b);

% Plot with and without error

subplot(1,2,1);imagesc( reshape(b,p,length(theta)))
subplot(1,2,2);imagesc( reshape(bn,p,length(theta)))

%
% Reconstruct with backslash with no regularization

xnr=reshape(A\bn,N,N);figure,imagesc(reshape(A\bn,N,N));axis image off; colormap gray;
plot(xnr(N/2,:))

%
%

nA=norm(full(A));
TikFac=0.5*nA;
Aa =[A; TikFac*eye(N^2)];
bna=[bn;zeros(N^2,1)];
xnt = Aa\bna;
xntr= reshape(xnt,N,N);
imagesc(xntr); title('Simple Tikhonov regularization'); axis off
plot(xntr(N/2,:));title('Cross section through above')

%
% Make a derivative maltrix for the image. See [1,Chapter 12]

D1=diag(-ones(1,N))+ diag(ones(1,N-2),2);
D1(end,end)=0;D1=D1/N;
L =[kron(eye(N),D1);kron(D1,eye(N))];
L=sparse(L);
spy(L)

%
% Generalized Tikhonov using this matrix L penalizes ||Lx||^2 rather than just ||x||^2 so it constrains image values on adjacent pixels to be similar.

TikFac=0.5*nA;
Aa =[A; TikFac*L];
bna=[bn;zeros(2*N^2,1)];
xnt = Aa\bna;
xntr= reshape(xnt,N,N);
imagesc(xntr); title('Gen Tikhonov regularization');
plot(xntr(N/2,:));title('Cross section through above')

%%
% Suggestion:  1) Is it sensible to use the same regularization factor with I and with L? Should it be normalized differently? 2) was the choise of the regularization parameter here lucky? What hapens if you change it, or change the level of noise?  How well do the solutions fit the data (ie norm of the misfit to the data)?
%
%% 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.
%