% This is an example code for the open access Gel phantom data for dynamic
% tomography. It performs a simple Tikhonov regularization on some of the 
% time frames of the dataset using the conjugate gradient method.
%
% This routine requires the ASTRA Toolbox, Spot operator Toolbox and
% HelTomo Toolbox to work properly! A CUDA compatible graphics card (GPU)
% is also required!
%
% Tommi Heikkilä, University of Helsinki - 5.3.2020
%
%%%     %%%     %%%     %%% PROGRAM STARTS HERE %%%     %%%     %%%     %%% 

% Load the data structure from the file.
load GelPhantomData_b4.mat

% Choose desired time frames
frames = [1,7,14];

% Define reconstruction size
N = 256;

% Regularization parameter
alpha = 0.1;

%% Case 1: 120 projections
% Desired measurement angles in degrees
angles = [0:3:357];

% Downsample the data from the first measurement
CtData = subsample_sinogram(GelPhantomData_b4(1), angles);

% Create the ASTRA Spot forward operator
A = create_ct_operator_2d_fan_astra_cuda(CtData, N, N);

% Normalize the operator
normA = normest(A);
A = A/normA;

% Define the forward model as a function
fun = @(x) A.'*(A*x(:))+alpha*x(:);

% Preallocate the reconstruction
recn120 = zeros(N,N,length(frames));

% Go through the time frames
for i = 1:length(frames)
    % Pick the desired data and downsample it
    CtData = subsample_sinogram(GelPhantomData_b4(frames(i)), angles);
    % Normalize sinogram
    m = CtData.sinogram / normA;
    % Back-project the data
    b = A.'*m(:);
    % Perform conjugate gradient algorithm
    recn = pcg(fun, b);     % This is a column vector!
    % Save the reconstructed time frame
    recn120(:,:,i) = reshape(recn, N, N);
end

%% Case 2: 20 projections
% Measurement angles in degrees
angles = [0:18:342];

% Downsample the data from the first measurement
CtData = subsample_sinogram(GelPhantomData_b4(1), angles);

% Create the ASTRA Spot forward operator
A = create_ct_operator_2d_fan_astra(CtData, N, N);

% Normalize the operator
normA = normest(A);
A = A/normA;

% Define the forward model as a function
fun = @(x) A.'*(A*x(:))+alpha*x(:);

% Preallocate the reconstruction
recn20 = zeros(N,N,length(frames));
% Go through the time frames
for i = 1:length(frames)
    % Pick the desired data and downsample it
    CtData = subsample_sinogram(GelPhantomData_b4(frames(i)), angles);
    % Normalize sinogram
    m = CtData.sinogram / normA;
    % Back-project the data
    b = A.'*m(:);
    % Perform conjugate gradient algorithm
    recn = pcg(fun, b);     % This is a column vector!
    % Save the reconstructed time frame
    recn20(:,:,i) = reshape(recn, N, N);
end

%% Case 3: FBP from original projections
% Preallocate the reconstruction
recnFBP = zeros(N,N,length(frames));
% Go through the time frames
for i = 1:length(frames)
    % Pick the desired data and reconstruct it
    recnFBP(:,:,i) = tomorecon_2d_fan_fbp_astra_cuda(GelPhantomData_b4(frames(i)),N, N);
end

%% Look at the results and the corresponding sinograms
% Normalize the reconstructions to interval [0,1] for viewing
recn120 = recn120/max(recn120(:));
recn20 = recn20/max(recn20(:));
recnFBP = recnFBP/max(recnFBP(:));

for i = 1:length(frames)
    figure(i);
    hold on
    % 120 projection sinogram
    subplot(2,2,1)
    imagesc(GelPhantomData_b4(frames(i)).sinogram(1:3:358,:))
    colormap gray
    axis square
    axis off
    title({'Sinogram,'; ['t=', num2str(frames(i)), ', 120 projections']})
    
    % 120 projection reconstruction
    subplot(2,2,2)
    imagesc(recn120(:,:,i),[0,1])
    colormap gray
    axis square
    axis off
    title({'Tikhonov regularization,'; ['t=', num2str(frames(i)), ', 120 projections']})
    
    % 20 projection sinogram
    subplot(2,2,3)
    imagesc(GelPhantomData_b4(frames(i)).sinogram(1:18:343,:))
    colormap gray
    axis square
    axis off
    title({'Sinogram,'; ['t=', num2str(frames(i)), ', 20 projections']})
    
    % 20 projection reconstruction
    subplot(2,2,4)
    imagesc(recn20(:,:,i),[0,1])
    colormap gray
    axis square
    axis off
    title({'Tikhonov regularization,'; ['t=', num2str(frames(i)), ', 20 projections']})
    
    hold off
end

% FBP reconstructions
figure(length(frames)+1);
for i = 1:length(frames)
    hold on
    subplot(1,length(frames),i)
    imagesc(recnFBP(:,:,i),[0,1])
    colormap gray
    axis square
    axis off
    title({'Filtered back-projection,'; ['t=', num2str(frames(i)), ', 360 projections']})
    hold off
end
