%%%%%% stempo_fdk_example.m %%%%%%
%
% Example code for reconstructing the STEMPO phantom data using FDK
% (Feldkamp-Davis-Kress) algorithm. FDK is the industry standard analytic
% reconstruction method for (spatially) 3D data.
% 
%%%%%%
% 
% Requirements:
% ASTRA Toolbox
% https://www.astra-toolbox.com/
% Recommended v1.9 or higher
%
% HelTomo Toolbox
% https://github.com/Diagonalizable/HelTomo
% v2.0
%
% CUDA capable graphics card
% See https://developer.nvidia.com/cuda-gpus#compute or try Matlab command:
% 'gpuDevice'
%
%%%%%%
%
% Created 15.9.2022 - Last edited 15.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_3d_b%d.mat',binning))
    case 'cont360'
        % Continuous 360 projection scan with 1 degree angle interval
        load(sprintf('stempo_cont360_3d_b%d.mat',binning))
    case 'seq8x45'
        % Sequence of 8x45 projection scans with 8 degree angle interval
        load(sprintf('stempo_seq8x45_3d_b%d.mat',binning))
end

%% Look at the projection images

P = CtData.parameters.numberImages;
figure(1);
for i = 1:P
    imagesc(squeeze(CtData.sinogram(:,i,:))')
    title(sprintf('Projection image from angle: %.1f deg.', CtData.parameters.angles(i)))
    colormap gray
    axis off
    drawnow
    pause(0.02)
end

%% Full angle FDK

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

% We need to just pass the CtData structure and desired resolution.
% Available projection angles are already stored in CtData.parameters.angles
fdk = tomorecon_3d_cone_fdk_astra_cuda(CtData, Nx, Ny, Nz);

%% Look at the reconstructed volume

VSconfig = defineVSconfig(0.11);
figure(2);
% Cut the volume in half
VS = volshow(fdk(1:Nx/2,:,:), VSconfig);
pause(1)

% Remove unwanted parts for better visualization
fdkCut = fdk(:,:,round(0.2*Nz):round(0.8*Nz));
VSconfig = defineVSconfig(0.15);
figure(3);
VS = volshow(fdkCut, VSconfig);

%% Helper function
function VSconfig = defineVSconfig(isoValue)
    % Helper function for defining reasonable config values for the volshow
    % function
    % Most of these are the default values given by the Volume Viewer app.
    % Threshold value 'isoValue' for the isosurface must be given by user.
    VSconfig.CameraPosition = [10/3 10/3 2];
    VSconfig.CameraUpVector = [0 0 1];
    VSconfig.CameraTarget = [0 0 0];
    VSconfig.CameraViewAngle = 15;
    VSconfig.BackgroundColor = [1 1 1];
    VSconfig.Renderer = 'Isosurface';
    VSconfig.Colormap = repmat(linspace(0,1,256)',1,3);
    VSconfig.Alphamap = linspace(0,1,256)';
    VSconfig.Lighting = true;
    VSconfig.IsosurfaceColor = [0.5 0.5 0.5];
    VSconfig.Isovalue = isoValue;
end
