#!/bin/bash

##
##
##  This is the Unix script to run the PET example for dynamic reconstruction incorporating:
##  - A dynamic model (Patlak)
##  - Timestamp-based motion correction (using rigid motion correction)
##  Dataset is list-mode generated from analytical projection of Zubal digital phantom with attenuation and noise
##   Zubal, I.G., Harrell, C.R, Smith, E.O, Rattner, Z., Gindi, G. and Hoffer, P.B.,
##   Computerized three-dimensional segmented human anatomy, Medical Physics, 21(2):299-302, 1994.
##  Activity in voxels represent typical FDG time activity curves.
##  The dynamic protocol starts at 2100 seconds, to reconstruct 5 frames of 300s each
##  5 motions are simulated at 2280s, 2460s, 3060s, 3180s and 3480s
##
##  This script will run the example and generate:
##  - 2 parametric images: "example_pet_dmodel-rmotion_reco5D_it2_par*" corresponding to Patlak slope and intercept
##  - 5 reconstructed images (each frame): "example_pet_dmodel-rmotion_reco5D_it*" 
##  - 5 sensitivity images (each frame) "example_pet_dmodel-rmotion_reco5D_sensitivity_fr*"
##
##

###########################
# Some checks
###########################

# Set the CASToR reconstruction program
recon="castor-recon"
# Test the existency of the CASToR reconstruction program in the PATH
type ${recon} > /dev/null
if [ $? != 0 ]
then
  echo "***** In order to run the benchmark script, please add the CASToR binary folder into your PATH !"
  exit 1
fi

###########################
# Command-line options
###########################

# To get help about the command-line options, run the program without argument or with '-h',
# '-help' or '--help' options. Then specific help about the different parts of the project
# can be prompted using specific help options.

# General verbose level (the verbose level relative to each specific part of the CASToR project
# can be customize separately if desired using options '-vb-XXXX'; to get specific help about
# these options, run the main program with the option '-help-misc')
verbose="-vb 2"

# The data file header (here, a pet list-mode data file). You can have a look at this ascii
# header to see the different fields used by the reconstruction program.
datafile="-df example_pet_dmodel-rmotion_df.cdh"

# Attenuation image, for atn correction (cm-1)
atn="-atn mumap_sebal_5mm.hdr"

# The output file base name. All files saved on the disk by the program will use this base name
# with an appended suffix and an extension. As an alternative, you can use the option '-dout'
# instead of '-fout'; this will create a folder of the provided name, and all files will be
# saved inside this folder. Note that the CASToR programs create a '.log' file which logs
# everything that is prompted on the screen.
output="-fout example_pet_dmodel-rmotion_reco5D"

# This is an option to specify that we want to save the image of the last iteration only. If not
# specified, by default, all iterations will be saved. You can alse use the '-oit' option to
# specifiy the series of iterations that you want to save. Run the reconstruction program with
# the '-help-out' option to get details about all options specific to the output settings.
last_it="-oit -1"

# Number of iterations (2) and subsets (28). Using this option, you can also define sequences
# of iterations with different number of subsets. For example '-it 2:28,1:14,10:1' will perform
# 2 iterations using 28 subsets, followed by 1 iteration using 14 subsets and followed by 10
# iterations using 1 subset.
iteration="-it 2:16"

# Number of voxels of the reconstructed image (X,Y,Z)
vox_number="-dim 128,128,64"

# Size of the field-of-view in mm (X,Y,Z). Alternatively, you can use the '-vox' option to specify
# the size of the voxels.
vox_size="-vox 2.2,2.2,2.8"

# The reconstruction algorithm. Here, we use an iterative optimization algorithm. If you want to
# get specific help on how to use this option or other algorithms, run the program with the
# '-help-algo' option. If you want to get the full list of all implemented optimization algorithms,
# run the program with the '-help-opti' option. From that list, you will get the description of
# how to parameterize each optimization algorithm. Here, for the MLEM algorithm, by not specifying
# any parameter, the default configuration file is used to parameterize the algorithm.
optimizer="-opti MLEM"

# The projection algorithm. Here the incremental Siddon projector is used. To get help about how to use
# this option, run the program with the '-help-proj' option and to get a list of all implemented
# projectors, use the '-help-projm' option.
projector="-proj joseph"

# Temporal splitting of the dynamic data into different images.
# The starting time of each frame must be entered. By default the end of the frame will be the beginning
# of the next frame, unless a duration is entered after a ':'. For the last frame, a duration must be provided.
# Here he choose to split the dataset in five 300s frames, starting from 2100 to 3300 seconds
# after the beginning of acquisition. 5 images will be reconstructed.
framing="-frm 2100,2400,2700,3000,3300:300"

# Dynamic model to be applied inside reconstruction between the frames. 
# The model type must be entered followed by its parameters 
# (provided either with set of parameters, or a configuration file containing these parameters) 
# Here the Patlak model is used with the nested optimization method (see Patlak help for more info)
dmodel="-dynamic-model Patlak:dmodel_patlak_nested.cf"
#dmodel="-dynamic-model Patlak,23682.79,25228.74,26636.99,27923.61,29101.16,5.4,4.91,4.48,4.1,3.75"

# Setup the motion correction, deformation model, and its parameters. 
# Different option can be used to define the type of motion:
# -rm: respiratory (cyclic) motion, for gated dataset
# -cm: cardial (cyclic) motion, for gated dataset
# -rcm: respiratory and cardiac (cyclic) motion, for dual-gated dataset
# -im: involuntary patient motion (chronological). This type of motion involves deformation at specific timestamps.
# We use here a model for rigid transformation, requiring 6 parameters (3 translations and 3 rotations)
# The parameters are provided from a configuration file (rigid_transformation_par.txt)
motion="-im deformationRigid:rigid_transformation_par.txt"

# This option defines the splitting of the dynamic dataset.
# A file must be provided. The information it contains depends on the type of dataset and motion to correct: 
# - physiological motion: number of events in each gates of a respiratory/cardiac gated dataset
# - involuntary motion:   timestamp for each transformation to perform. This is usually associated with head/body motion. 
# The file provided here contains the number of transformation, and the timestamps for theses transformation.
# This option is mandatory if the -im motion is enabled.
gating="-g motion_timestamp.txt"

# Parallel computation using the OpenMP library. If CASToR was not compiled using OpenMP, a warning will
# be displayed if this option is used, specifying that only one thread will be used. Here we specify '0'
# in order to let the computer choose the number of threads with respect to the available ressources.
# One can also manually specify the number of threads to be used. To get details about the computation
# settings, run the program with the '-help-comp' option.
thread="-th 0"

# We want to mask the extrem slices when saving the image. To do that, we use the '-slice-out' option
# in which we give the number of slices that we want to mask at both sides of the FOV.
out_mask_axial="-slice-out 3"

# We want to mask the edge of the image in the transaxial view. The number correspond to the percentage
# of the eliptical transaxial FOV to be kept while saving the image.
out_fov="-fov-out 90."

# This option merge output dynamic images (reconstructed or parametric) in one dynamic file instead of 3D separate files 
merge=""
#merge="-omd"


###########################
# Launch the reconstruction
###########################

# Launch the reconstruction
echo "=============================================================================================="
echo "Reconstruction is going on. Should take from one to several minutes depending on the hardware."
echo "=============================================================================================="
${recon} ${verbose} ${datafile} ${atn} ${output} ${last_it} ${iteration} ${vox_number} ${vox_size} ${optimizer} ${projector} ${framing} ${dmodel} ${motion} ${gating} ${thread} ${out_mask_axial} ${out_fov}

# Finished
echo ""
exit 0

