%%%%%% stempo_fbp_example.m %%%%%%
%
% Example code for reconstructing the STEMPO phantom data using FBP
% (filtered back-projection).
% 
%%%%%%
% 
% Requirements:
% ASTRA Toolbox
% https://www.astra-toolbox.com/
% Recommended v1.9, problems with v2.0
%
% HelTomo Toolbox
% https://github.com/Diagonalizable/HelTomo
% v2.0
%
%%%%%%
%
% Created 8.9.2022 - Last edited 19.9.2022
% Tommi Heikkilä
% University of Helsinki

% Clear workspace
clear all
close all

%% Load data
dataType = 'static';
binning = 16;

switch dataType
    case 'static'
        % Static scan for comparison and computing ground truth
        load(sprintf('stempo_static_2d_b%d.mat',binning))
    case 'cont360'
        % Continuous 360 projection scan with 1 degree angle interval
        load(sprintf('stempo_cont360_2d_b%d.mat',binning))
    case 'seq8x45'
        % Sequence of 8x45 projection scans with 8 degree angle interval
        load(sprintf('stempo_seq8x45_2d_b%d.mat',binning))
end

%% Look at the sinogram

figure;
imagesc(CtData.sinogram)
title('Sinogram')
xlabel('Detector bins')
ylabel('Projection angles')
% Note the 2D sinogram orientation determined by ASTRA

%% Filtered back-projection

% Reconstruction resolution (based on binning, other options are possible)
Nx = 2240/binning;
Ny = 2240/binning;

% Use CUDA (if available)
useCUDA = true;

    % We need to just pass the CtData structure and desired resolution.
    % Available projection angles are already stored in
    % CtData.parameters.angles
if useCUDA
    fprintf('Using CUDA \n')
    fbp = tomorecon_2d_fan_fbp_astra_cuda(CtData, Nx, Ny);
else
    try % CPU implementation of fanbeam FBP has issues with ASTRA 2.0
        fbp = tomorecon_2d_fan_fbp_astra(CtData, Nx, Ny);
    catch e
        fprintf(' Error: %s \n', e.message);
        fbp = nan;
    end
end

figure;
imagesc(fbp)
title(sprintf('FBP reconstruction using %d angles', CtData.parameters.numberImages))
axis equal
axis off
% Note the reconstruction orientation determined by ASTRA

%% Sparse angle FBP

% Sparse subset of projection angles
anglesSparse = 2:6:359;

% This returns a new CtData structure with the correct parameters
CtDataSparse = subsample_sinogram(CtData, anglesSparse);

if useCUDA
    fprintf('Using CUDA \n')
    fbpSparse = tomorecon_2d_fan_fbp_astra_cuda(CtDataSparse, Nx, Ny);
else
    try % CPU implementation of fanbeam FBP has issues with ASTRA 2.0
        fbpSparse = tomorecon_2d_fan_fbp_astra(CtDataSparse, Nx, Ny);
    catch e
        fprintf(' Error: %s \n', e.message);
        fbpSparse = nan;
    end
end

figure;
imagesc(fbpSparse)
title(sprintf('Sparse angle FBP reconstruction using %d angles', CtDataSparse.parameters.numberImages))
axis equal
axis off
% Note the reconstruction orientation determined by ASTRA
