%% Tomography with Total Variation
% In this tutorial we generate a parallel beam Radon transform tomography problem (isin Air Tools II) with noisy data and use a primal dual solver from flexBox https://github.com/HendrikMuenster/flexBox by Hendrik Dirks [1]. For more details of optimization methods sutable for TV regularized inverse problem solution see [1, Chapter 13].

addpath('flexBox','AirToolsII')

%
% Set up standard paralel beam problem

% Set the parameters for the test problem.
% Set to be close to the exampleRadon.m in flexBox.
N = 128;           % The image is N-times-N..
theta = 1:3:180;;  % Projection angles.
p = 185;           % No. of parallel rays.
% Create the test problem.
[A,b,x] = paralleltomo(N,theta,p);
tt=tiledlayout(1,2);
nexttile(tt); imagesc(reshape(x,N,N))
colormap gray,axis image off
title('Original image');
nexttile(tt); imagesc(reshape(b,p,length(theta)));
axis image off
title('Original sinogram');

%
%
% Add Gaussian pseudo random noise to the sinogram

bn = b + randn(size(b))*0.001*norm(b);
figure, imagesc(reshape(bn,p, length(theta))); axis  off;
title('Sinogram with noise')

%
% Now set up the flexbox optimization problem as in exampleRadon.m

main = flexBox;
main.params.tryCPP = 0; %change, if C++ module is compiled
%add primal var u
numberU = main.addPrimalVar([N,N]);% Image is N by N

%add data-fidelity: 1/2\|Au-f\|_2^2
main.addTerm(L2dataTermOperator(1,A,bn),numberU);
%add regularizer
%main.addTerm(L2identity(1,[N,N]),numberU); %Tikhonov-regularization
%main.addTerm(L2gradient(2,[N,N]),numberU); %Smoothness
main.addTerm(L1gradientIso(0.08,[N,N]),numberU); %TV-regularization
%box constraint ensures the result to stay in [0,1]
main.addTerm(boxConstraint(0,1,[N,N]),numberU);

%run minimization algorithm
tic;main.runAlgorithm;toc

%get result
result = main.getPrimal(numberU);

%
% Display results

imagesc(result);axis image;colormap(gray);title('TV Reconstructed Image')
axis off image

%
% Compare with Matlab (image processing toolbox's) iradon.

imagesc(iradon(reshape(b,p,length(theta)),theta));axis off image;
title('Reconstruction with Matlab iradon')

%% References
% [1] Dirks, Hendrik. "A Flexible Primal-Dual Toolbox", arXiv e-prints (2016): arXiv-1603.
% [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.
% [3] 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