Features

This page contains a list of features and modules implemented in the current CASToR version.

CASToR v3.2

Main features:

• Iterative optimization-based reconstruction algorithms.
• PET (including Time of Flight), SPECT (parallel-hole and convergent-beam collimators) and CT geometries.
• Reconstruction of multi-frame and respiratory/cardiac gated acquisitions.
• Image-based motion management for gated acquisition and timestamp-based motion.
• Projector plug-in class system, including several methods (see below)
• Optimization algorithm plug-in class system, including several methods (see below) and a mechanism for optimization algorithms making use of sub iterations (e.g ADMM).
• Dynamic model class, dedicated to kinetic modelling / parametric image reconstruction (see below).
  Including an utility tool to apply the implemented dynamic model on a set of dynamic image (post-reconstruction kinetic analysis).
• An image-based deformation class, dedicated to image-based transformation for motion correction (see below).
• Plug-in class systems for image convolution and spatial regularization, and general image processing.
• A fully-bayesian reconstruction algorithm (RCP-GS)
• Two-level parallel CPU implementation using MPI (multi-computers) and openMP (multi-threads) libraries.
• Interfile I/O format for images.
• Utility tool to convert GATE Monte Carlo simulated data in ROOT format into the CASToR datafile format.
  SIMIND SPECT simulation platform also includes its own conversion program for CASToR.
• Easy and straightforward compilation procedure, as well as already compiled binaries for Unix and Windows 64-bits systems and Windows 32-bits systems.
• Utility tools for the exploration of datafile and geometry files.
• Toolkit for the computation of direct normalization factors using a normalization phantom.
• Event by event TOF resolution kernel for list-mode PET.
• User customizable fields in the Datafile

This version contains the following implemented modules.
More implementation details can be displayed when using the command-line options -help-projm, -help-opti, -help-pnlt and -help-conv with castor-recon for projectors, optimizers, penalties and image-based convolvers respectively.

Projectors

• Original Siddon:
  Line projector that computes the exact path length of a line through the voxels. Implemented from the following paper:
  RL Siddon. "Fast calculation of the exact radiological path for a three-dimensional CT array". Med Phys. 1985 Mar-Apr;12(2):252-5.
• Accelerated Siddon:
  Incremental version of the original algorithm proposed by R. L. Siddon. Implemented from the following paper:
  Filip Jacobs, et al. "A Fast Algorithm to Calculate the Exact Radiological Path Through a Pixel Or Voxel Space". Journal of Computing and Information Technology 6(1) · December 1998
• Joseph (not compatible with SPECT attenuation):
  Line projector that uses linear interpolation between pixels. Implemented from the following paper:
  Joseph PM. "An Improved Algorithm for Reprojecting Rays through Pixel Images". IEEE Trans Med Imaging. 1982;1(3):192-6.
• Distance-driven (not compatible with SPECT attenuation):
  Line projector based on computations of overlap between a pair of detection elements and voxels. Implemented from the following paper:
  B. De Man and S. Basu. "Distance-driven projection and backprojection in three dimensions", Phys. Med. Biol., vol. 49, pp. 2463-75, 2004.

Optimizers (reconstruction algorithms)

• MLEM:
  This optimizer is the standard MLEM (Maximum Likelihood Expectation Maximization).
  L. A. Shepp and Y.Vardi, "Maximum likelihood reconstruction for emission tomography" IEEE Trans. Med. Imaging, vol. MI-1, pp. 113–122, 1982.
  If subsets are used, it naturally becomes the OSEM optimizer. With transmission data (CT), the log-converted pre-corrected data are used as in J. Nuyts et al: "Iterative reconstruction for helical CT: a simulation study", Phys. Med. Biol., vol. 43, pp. 729-737, 1998.
• NEGML (only compatible with histogram, emission data):
  NEGML algorithm from K. Van Slambrouck et al, IEEE TMI, Jan 2015, vol. 34, pp. 126-136."
  Subsets can be used. This implementation only considers the psi parameter, but not the alpha image design parameter which is supposed to be 1 for all voxels. It implements the equation 17 from the paper.
• AML (AB-EMML) (only compatible with histogram, emission data):
  This optimizer is the AML algorithm derived from the AB-EMML of C. Byrne, "Iterative algorithms for deblurring and deconvolution with constraints" Inverse Problems, 1998, vol. 14, pp. 1455-67
  The bound B is taken to infinity, so only the bound A can be parameterized. This bound must be quantitative (same unit as the reconstructed image). It is provided as a single value and thus assuming a uniform bound. Subsets can be used.
  With a negative or null bound, this algorithm implements the equation 6 of A. Rahmim et al, Phys. Med. Biol., 2012, vol. 57, pp. 733-55.
  If a positive bound is provided, then we suppose that the bound A is taken to minus infinity. In that case, this algorithm implements the equation 22 of K. Van Slambrouck et al, IEEE TMI, Jan 2015, vol. 34, pp. 126-136.
• Landweber (only compatible with histogram data):
  This optimizer implements the standard Landweber algorithm for least-squares optimization. With transmission data, it uses the log-converted model to derive the update. The relaxation parameter is not automatically set, so it often requires some trials and errors to find an optimal setting. This algorithm is particularly slow to converge.
  L. Landweber "An Iteration Formula for Fredholm Integral Equations of the First Kind" American Journal of Mathematics, Vol. 73, No. 3 (Jul., 1951), pp. 615-624
• MLTR (only compatible with histogram, transmission data):
  This optimizer is a version of the MLTR algorithm implemented from the equation 16 of the paper from K. Van Slambrouck and J. Nuyts: Reconstruction scheme for accelerated maximum lihelihood reconstruction: the patchwork structure, IEEE Trans. Nucl. Sci., vol. 61, pp. 173-81, 2014.
  An additional empiric relaxation factor has been added onto the additive update. Its value onto the first and last updates can be parameterized. Its value for all the intermediate updates is computed linearly from these first and last provided values. Subsets can be used.
• One-Step-Late:
  This optimizer is the algorithm proposed by P. J. Green, IEEE TMI, Mar 1990, vol. 9, pp. 84-93.
  Subsets can be used as for OSEM. It accepts penalty terms that have a derivative order of at least one. Without penalty, it is stricly equivalent to the MLEM algorithm (thus can used with transmission data, see MLEM description).
• BSREM (only compatible with histogram, emission data):
  This optimizer is the Block Sequential Regularized Expectation Maximization (BSREM) algorithm from S. Ahn and J. Fessler, IEEE TMI, May 2003, vol. 22, pp. 613-626. Its abbreviated name in this paper is BSREM-II.
  This algorithm is the only one to have proven convergence using subsets. Its implementation is entirely based on the reference paper. It may have numerical problems when a full field-of-view is used, because of the sharp sensitivity loss at the edges of the cylindrical field-of-view. As it is simply based on the gradient, penalty terms must have a derivative order of at least one. Without penalty, it reduces to OSEM but where the sensitivity is not dependent on the current subset. This is a requirement of the algorithm, explaining why it starts by computing the global sensitivity before going through iterations.
• Penalized Preconditioned Gradient ML (only compatible with histogram, emission data):
  This optimizer is the Penalized Preconditioned Gradient algorithm from J. Nuyts et al, IEEE TNS, Feb 2002, vol. 49, pp. 56-60.
  As usually described by its inventor, it is a heuristic but effective gradient ascent algorithm for penalized maximum-likelihood reconstruction. It addresses the shortcoming of One-Step-Late when large penalty strengths can create numerical problems. Penalty terms must have a derivative order of at least two. Subsets can be used as for OSEM. Without penalty, it is equivalent to the gradient ascent form of the MLEM algorithm.
• Modified EM for Penalized ML (only compatible with emission data):
  This optimizer is based on the algorithm from A. De Pierro, IEEE TMI, vol. 14, pp. 132-137, 1995.
  This algorithm uses optimization transfer techniques to derive an exact and convergent algorithm for maximum likelihood reconstruction including a Markov Random Field penalty with different potential functions. The algorithm is convergent and is numerically robust to high penalty strength. It is stricly equivalent to MLEM without penalty, but can be unstable with extremly low penalty strength. Currently, it only implements the quadratic penalty. Subsets can be used as for OSEM, without proof of convergence though.
• ADMM:
  This optimizer is the ADMM with non-negativity constraint on projection space. It was proposed by H. Lim, Y. K. Dewaraja, and J. A. Fessler, Physics in Medicine & Biology, 2018. It is for now developed without subsets and implemented using the same equations as in the original paper, also with one inner iteration.
  The algorithm can be used without penalty, with quadratic penalty or with Markov Random Field (MRF) penalty with quadratic potential.
• APPGML:
  This optimizer is the A-Penalized Preconditioned Gradient algorithm (A-PPGML) from M. Millardet et al, IEEE TRPMS, 2021, vol. 6, pp. 629-640. It is based on the same idea as AML, which is allowing for negative values in the image by shifting an algorithm enforcing positive values. APPGML can be seen as a shifted version of PPGML with a lower bound A. This algorithm allows for negative image values in case the provided bound is also negative. The shift is applied uniformly in the whole image unless a shift image is provided. This shift image allows to customize the shift value applied in each voxel. The algorithm can also use penalized reconstruction because it is based on PPGML.
• MLMUMAP:
  This optimizer is a modified MLEM algorithm where the image voxel values are attenuation factors and the data are taken as the logarithm of ACF values associated to PET histogram data. It is numerically implemented in the multiplicative form (as opposed to the gradient form). It truncates ACF lower than 1 to 1 to satisfy the positivity constraint. Note that although physicists are used to attenuation values in 1/cm, the reconstructed image units are 1/mm as the whole castor code is using mm. Subsets can be used.
   

Penalties (that can be associated with some optimizers)

• MRF (Markov Random Field):
  The Markov Random Field (MRF) penalty is implemented for several types of neighborhood, potential functions, similarity and proximity factors. Neighborhood shape can be defined as nearest neighbors or a sphere. Proximity factors can be uniform or based on the distance between neighbors. Similarity factors can be uniform or based on the asymmetrical Bowsher's method. The potential function can be quadratic, based on relative differences (Nuyts et al, IEEE Trans. Nucl. Sci., vol. 49, pp. 56-60, 2002) or based on many different published functions.
• MRP (Median Root Prior):
  The Median Root Prior (MRP) is implemented for several types of neighborhood as for MRF. The current implementation is the simplest one, based on the following reference: S. Alenius and U Ruotsalainen, "Bayesian image reconstruction for emission tomography based on the median root prior", European Journal of Nuclear Medicine, vol. 24, pp. 258-265, 1997.

Dynamic Models

• LinearModel (Base class for linear models):
  - Several level of application of the models (dynamic frames, respiratory/cardiac gates).
  - Several optimisation methods (nested EM, NNLS, direct (within system matrix) ).
  - Interpolation of arterial input curve interpolation on framing protocols.
• Patlak (Patlak Model):
  - Automatic Patlak basis function computation from arterial input curve.
  - Several optimisation methods (nested EM, NNLS, direct (within system matrix), linear regression ).
• Spectral (Spectral Model):
  - Parameterization of spectral function coefficients (number, rate, etc..).
  - Interpolation of arterial input function on framing protocol.
• _1TCM (1 Tissue Compartment Model):
  - Several optimisation methods (NNLS, LS with ridge-regression).
  - Several integration methods (Weighed parabola overlapping, trapezoid).

Image-based Deformation Models

• deformationRigid (Rigid deformation model ):
  Transformation performed through vectors containing 3 translation and 3 rotation parameters

Image-based convolution

An image convolver is called through the -conv option. The provided argument describes the convolver to be used, its options, and when to include it within the algorithm. The 'when' parameter is an argument describing when to include the convolver within the algorithm. It is a list of keywords separated by commas. The following keywords can be used:
 -forward  (include convolver into the forward model; a convolution of the current estimate is forward-projected)
 -backward (include convolver into the backward model; a convolution of the correction terms is used for the update)
 -post     (apply convolver before saving the image; the convolved image is not put back as the estimate for the next update)
 -psf      (include both 'forward' and 'backward'; the standard image-based PSF modelling)
 -sieve    (include both 'psf' and 'post'; the standard method of sieve)
 -intra    (apply convolver to the updated image and use it as the current estimate for the next update)

Currently, only a non-variant Gaussian convolver is implemented:

• Stationary Gaussian:
  This convolver is based on a classic stationary Gaussian kernel. It can be anisotropic so that the transaxial and axial FWHM can be different. One also has to choose the number of sigmas of the Gaussian distribution that will be included in the kernel. The following options can be used (in this particular order when provided as a list, separated by commas):
  -transaxial FWHM: to set the transaxial FWHM (in mm)
  -axial FWHM: to set the axial FWHM (in mm)
  -number of sigmas: to set the number of sigmas included in the kernel (recommendations: at least 3. and maximum 5.)