% Configuration.
distance_source_origin = 300; % [mm]
distance_origin_detector = 100; % [mm]
detector_pixel_size = 1.05; % [mm]
detector_rows = 200; % Vertical size of detector [pixels].
detector_cols = 200; % Horizontal size of detector [pixels].
num_of_projections = 180;
angles = linspace2(0, 2 *pi, num_of_projections); %linspace2 from ASTRA toolbox

% Create phantom.
phantom = zeros(detector_rows, detector_cols, detector_cols);
hb = 110; % Height of beam [pixels].
wb = 40; % Width of beam [pixels].
hc = 100; % Height of cavity in beam [pixels].
wc = 30; % Width of cavity in beam [pixels].
phantom(detector_rows/2 - hb/2 : detector_rows/2 + hb/2,...
detector_cols/2 - wb/2 : detector_cols/2 + wb/2,...
detector_cols/2 - wb/2 : detector_cols/2 + wb/2) = 1;
phantom(detector_rows/2 - hc/2 : detector_rows/2 + hc/2,...
detector_cols/2 - wc/2 : detector_cols/2 + wc/2,...
detector_cols/2 - wc/2 : detector_cols/2 + wc/2) = 0;
phantom(detector_rows/2 - 5 : detector_rows/2 + 5,...
detector_cols/2 + wc/2 : detector_cols/2 + wb/2,...
detector_cols/2 - 5 : detector_cols/2 + 5) = 0;


vol_geom = astra_create_vol_geom(detector_rows, detector_cols, detector_cols);
phantom_id = astra_mex_data3d('create', '-vol', vol_geom, phantom);

% Create projections. With increasing angles, the projection are such that the
% object is rotated clockwise. Slice zero is at the top of the object. The
% projection from angle zero looks upwards from the bottom of the slice.
proj_geom = astra_create_proj_geom('cone', 1, 1, detector_rows, detector_cols, angles,...
(distance_source_origin + distance_origin_detector)/detector_pixel_size, 0);

[~, projections] = astra_create_sino3d_cuda(phantom_id, proj_geom, vol_geom);
projections = projections./max(projections(:));

%% Reconstruction:
projections_id = astra_mex_data3d('create', '-proj3d', proj_geom, projections);
reconstruction_id = astra_mex_data3d('create', '-vol', vol_geom);

cfg = astra_struct('FDK_CUDA');
cfg.ReconstructionDataId = reconstruction_id;
cfg.ProjectionDataId = projections_id;

% Create the algorithm object from the configuration structure
algorithm_id = astra_mex_algorithm('create', cfg);
% Run iterations of the algorithm
astra_mex_algorithm('iterate', algorithm_id , 1);
% Get the result
reconstruction = astra_mex_data3d('get', reconstruction_id);

% Cleanup.
astra_mex_data3d('clear');
astra_mex_projector('clear');
astra_mex_matrix('clear');
astra_mex_algorithm('clear');

slice1 = squeeze(reconstruction(fix(end/2),:,:));
slice2 = squeeze(reconstruction(:,fix(end/2),:));

figure;imshow(slice1,[]);
figure;imshow(slice2,[]);