Iterative Reconstruction

This page describes the iterative reconstruction algorithm in LEAP. These algorithms can be applied to any CT geometry. Iterative reconstruction algorithms may outperform Filtered Backprojection (FBP) reconstruction when the reconstruction problem is severly ill-posed such as very noisey data, not enough projections, or projections that only cover a limited angular range. Sometimes iterative reconstruction can underperform FBP due to overfitting when the line integral model is not an accurate model. Inaccurate models can be the result of incorrect CT geometry specification, strong beam hardening or scatter effects, etc.

We shall separate the following iterative reconstruction algorithms into the following categories

  • Algebraic: These algorithms use algebraic means to solve the reconstruction problem. The time per iteration is fast, but the algorithms are pretty basic.

    • SIRT

    • SART

  • Statistical (Transmission): These algorithm model the statistics of transmission CT data and are very good at reconstructing noisey data.

    • RWLS

    • MLTR

  • Statistical (Emission): These algorithm model the statistics of emission CT data and are very good at reconstructing noisey data.

    • MLEM

    • OSEM

  • Special Purpose: These algorithms are good at solving reconstruction problems where one has a small number of projections and/or the projections only cover a limited angular range.

    • RDLS

    • ASDPOCS

We shall use the following notation in this section

  • \(g\): projection (attenuation) data

  • \(f\): reconstruction volume

  • \(P\): forward projection operator

  • \(P^T\): backprojection operator (adjoint of forward projection)

  • \(\boldsymbol{1}\): vector of all ones

  • \(R(\cdot)\): regularization functional (e.g., TV)

For more information on how to specify regularization functions in LEAP, see the Filter Sequence description.

leapctype.tomographicModels.SIRT(self, g, f, numIter, mask=None)

Simultaneous Iterative Reconstruction Technique reconstruction

This is the same algorithm as a SART reconstruction with one subset. The SIRT algorithm is performed using the following update equation

\[\begin{eqnarray} f_{n+1} &:=& f_n + \frac{0.9}{P^T 1} P^T\left[ \frac{Pf_n - g}{P1} \right] \end{eqnarray}\]

The CT geometry parameters and the CT volume parameters must be set prior to running this function.

Parameters:

  • g (C contiguous float32 numpy or torch array) – projection data

  • f (C contiguous float32 numpy or torch array) – volume data

  • numIter (int) – number of iterations

  • mask (C contiguous float32 numpy or torch array) – projection data to mask out bad data, etc. (zero values indicate projection data pixels not to use)

Returns:

f, the same as the input with the same name

leapctype.tomographicModels.SART(self, g, f, numIter, numSubsets=1, mask=None, nonnegativityConstraint=True)

Simultaneous Algebraic Reconstruction Technique reconstruction

The SART algorithm is performed using two nested loops. The inner loop is performed like the SIRT algorithm (see SIRT algorithm documentation for a description of the algorithm), but the successive updates of the reconstructed volume are done with a subset of the projection angles. Once every subset of is complete, the process starts over again. Using these ordered subsets reduces the time it takes for this algorithm to converge.

If one wishes to combine this algorithm with regularization, (e.g., TV), please see ASDPOCS.

The CT geometry parameters and the CT volume parameters must be set prior to running this function.

Parameters:

  • g (C contiguous float32 numpy or torch array) – projection data

  • f (C contiguous float32 numpy or torch array) – volume data

  • numIter (int) – number of iterations

  • numSubsets (int) – number of subsets

  • mask (C contiguous float32 numpy or torch array) – projection data to mask out bad data, etc. (zero values indicate projection data pixels not to use)

  • nonnegativityConstraint (bool) – if true constrains values of reconstruction to be nonnegative

Returns:

f, the same as the input with the same name

leapctype.tomographicModels.ASDPOCS(self, g, f, numIter, numSubsets, numTV, filters=None, mask=None, nonnegativityConstraint=True)

Adaptive Steepest Descent-Projection onto Convex Subsets reconstruction

This algorithm combines SART with regularization (e.g., TV; see "filters" argument). See SART and SIRT documentation for more information. This algorithm solves the following optimization problem

\[\begin{eqnarray} \text{minimize } R(f) \text{ subject to } \|Pf - g\|_2^2 < \varepsilon \end{eqnarray}\]

The CT geometry parameters and the CT volume parameters must be set prior to running this function. This function actually implements the iTV reconstruction method which is a slight varition to ASDPOCS which we find works slightly better.

Here is the reference Ritschl, Ludwig, and Marc Kachelriess. “Improved total variation regularized image reconstruction (iTV) applied to clinical CT data.” In Medical Imaging 2011: Physics of Medical Imaging, vol. 7961, pp. 786-798. SPIE, 2011.

Parameters:

  • g (C contiguous float32 numpy or torch array) – projection data

  • f (C contiguous float32 numpy or torch array) – volume data

  • numIter (int) – number of iterations

  • numSubsets (int) – number of subsets

  • numTV (int) – number of TV diffusion steps, larger numbers perform stronger regularization/ smoothing

  • filters (filterSequence object) – list of regularization filters

  • mask (C contiguous float32 numpy or torch array) – projection data to mask out bad data, etc. (zero values indicate projection data pixels not to use)

  • nonnegativityConstraint (bool) – if true constrains values of reconstruction to be nonnegative

Returns:

f, the same as the input with the same name

leapctype.tomographicModels.LS(self, g, f, numIter, preconditioner=None, nonnegativityConstraint=True)

Least Squares reconstruction

This function minimizes the Least Squares cost function using Preconditioned Conjugate Gradient. The optional preconditioner is the Separable Quadratic Surrogate for the Hessian of the cost function which is given by (P*P1)^-1, where 1 is a volume of all ones, P is forward projection, and P* is backprojection. The Least Squares cost function is given by the following

\[\begin{eqnarray} C_{LS}(f) &:=& \frac{1}{2} \| Pf - g \|^2 \end{eqnarray}\]

The CT geometry parameters and the CT volume parameters must be set prior to running this function.

Parameters:

  • g (C contiguous float32 numpy or torch array) – projection data

  • f (C contiguous float32 numpy or torch array) – volume data

  • numIter (int) – number of iterations

  • preconditioner (string) – specifies the preconditioner as ‘SQS’, ‘RAMP’, or ‘SARR’

  • nonnegativityConstraint (bool) – if true constrains values of reconstruction to be nonnegative

Returns:

f, the same as the input with the same name

leapctype.tomographicModels.WLS(self, g, f, numIter, W=None, preconditioner=None, nonnegativityConstraint=True)

Weighted Least Squares reconstruction

This function minimizes the Weighted Least Squares cost function using Preconditioned Conjugate Gradient. The optional preconditioner is the Separable Quadratic Surrogate for the Hessian of the cost function which is given by (P*WP1)^-1, where 1 is a volume of all ones, W are the weights, P is forward projection, and P* is backprojection. The Weighted Least Squares cost function is given by the following

\[\begin{eqnarray} C_{WLS}(f) &:=& \frac{1}{2} (Pf - g)^T W (Pf - g) \end{eqnarray}\]

The CT geometry parameters and the CT volume parameters must be set prior to running this function.

Parameters:

  • g (C contiguous float32 numpy or torch array) – projection data

  • f (C contiguous float32 numpy or torch array) – volume data

  • numIter (int) – number of iterations

  • W (C contiguous float32 numpy array) – weights, should be the same size as g, if not given, W=exp(-g); can also be used to mask out bad data

  • preconditioner (string) – specifies the preconditioner as ‘SQS’, ‘RAMP’, or ‘SARR’

  • nonnegativityConstraint (bool) – if true constrains values of reconstruction to be nonnegative

Returns:

f, the same as the input with the same name

leapctype.tomographicModels.RLS(self, g, f, numIter, filters=None, preconditioner=None, nonnegativityConstraint=True)

Regularized Least Squares reconstruction

This function minimizes the Regularized Least Squares cost function using Preconditioned Conjugate Gradient. The optional preconditioner is the Separable Quadratic Surrogate for the Hessian of the cost function which is given by (P*P1)^-1, where 1 is a volume of all ones, P is forward projection, and P* is backprojection. The Regularized Least Squares cost function is given by the following

\[\begin{eqnarray} C_{RLS}(f) &:=& \frac{1}{2} \| Pf - g \|^2 + R(f) \end{eqnarray}\]

The CT geometry parameters and the CT volume parameters must be set prior to running this function.

Parameters:

  • g (C contiguous float32 numpy or torch array) – projection data

  • f (C contiguous float32 numpy or torch array) – volume data

  • numIter (int) – number of iterations

  • filters (filterSequence object) – list of differentiable regularization filters

  • preconditioner (string) – specifies the preconditioner as ‘SQS’, ‘RAMP’, or ‘SARR’

  • nonnegativityConstraint (bool) – if true constrains values of reconstruction to be nonnegative

Returns:

f, the same as the input with the same name

leapctype.tomographicModels.RWLS(self, g, f, numIter, filters=None, W=None, preconditioner=None, nonnegativityConstraint=True)

Regularized Weighted Least Squares reconstruction

This function minimizes the Regularized Weighted Least Squares cost function using Preconditioned Conjugate Gradient. The optional preconditioner is the Separable Quadratic Surrogate for the Hessian of the cost function which is given by (P*WP1)^-1, where 1 is a volume of all ones, W are the weights, P is forward projection, and P* is backprojection. The Regularized Weighted Least Squares cost function is given by the following

\[\begin{eqnarray} C_{RWLS}(f) &:=& \frac{1}{2} (Pf - g)^T W (Pf - g) + R(f) \end{eqnarray}\]

The CT geometry parameters and the CT volume parameters must be set prior to running this function.

Parameters:

  • g (C contiguous float32 numpy or torch array) – projection data

  • f (C contiguous float32 numpy or torch array) – volume data

  • numIter (int) – number of iterations

  • filters (filterSequence object) – list of differentiable regularization filters

  • W (C contiguous float32 numpy array) – weights, should be the same size as g, if not given, W:=exp(-g); can also be used to mask out bad data

  • preconditioner (string) – specifies the preconditioner as ‘SQS’, ‘RAMP’, or ‘SARR’

  • nonnegativityConstraint (bool) – if true constrains values of reconstruction to be nonnegative

Returns:

f, the same as the input with the same name

leapctype.tomographicModels.DLS(self, g, f, numIter, preconditionerFWHM=1.0, nonnegativityConstraint=False, dimDeriv=2)

Derivative Least Squares reconstruction

See documentation for RDLS because this is the same algorithm without the regularization.

leapctype.tomographicModels.RDLS(self, g, f, numIter, filters=None, preconditionerFWHM=1.0, nonnegativityConstraint=False, dimDeriv=1)

Regularized Derivative Least Squares reconstruction

This function minimizes the Regularized Derivative Least Squares cost function using Preconditioned Conjugate Gradient. The optional preconditioner is a 2D blurring for each z-slice. The Regularized Weighted Least Squares cost function is given by the following

\[\begin{eqnarray} C_{RDLS}(f) &:=& \frac{1}{2} (Pf - g)^T \Delta (Pf - g) + R(f) \end{eqnarray}\]

where \(\Delta\) is the Laplacian operator. The CT geometry parameters and the CT volume parameters must be set prior to running this function.

Parameters:

  • g (C contiguous float32 numpy or torch array) – projection data

  • f (C contiguous float32 numpy or torch array) – volume data

  • numIter (int) – number of iterations

  • filters (filterSequence object) – list of differentiable regularization filters

  • preconditionerFWHM (float) – specifies the FWHM of the blur preconditioner

  • nonnegativityConstraint (bool) – whether to apply a nonnegativity constraint

  • dimDeriv (int) – number of dimensions (1 or 2) to apply the Laplacian derivative

Returns:

f, the same as the input with the same name

leapctype.tomographicModels.MLTR(self, g, f, numIter, numSubsets=1, filters=None, mask=None)

Maximum Likelihood Transmission reconstruction

This function maximizes the Maximum Likelihood function of CT transmission data which assumes a Poisson noise model. This algorithm best models the noise for very low transmission/ low count rate data. The MLTR cost function is given by the following

\[\begin{eqnarray} C_{MLTR}(f) &:=& \left< -t\log\left(e^{-Pf}\right) + e^{-Pf} , 1 \right> \end{eqnarray}\]

where \(t = e^{-g}\) is the transmission data and "1" is a vector of all ones. The inner product notation is just use for simplicity, but all it is really doing is performing a sum over all the elements. The CT geometry parameters and the CT volume parameters must be set prior to running this function.

Parameters:

  • g (C contiguous float32 numpy or torch array) – projection data

  • f (C contiguous float32 numpy or torch array) – volume data

  • numIter (int) – number of iterations

  • numSubsets (int) – number of subsets (reduces the time it takes for this algorithm to converge)

  • filters (filterSequence object) – list of differentiable regularization filters

  • mask (C contiguous float32 numpy or torch array) – projection data to mask out bad data, etc. (zero values indicate projection data pixels not to use)

Returns:

f, the same as the input with the same name

leapctype.tomographicModels.MLEM(self, g, f, numIter, filters=None, mask=None)

Maximum Likelihood-Expectation Maximization reconstruction

This algorithm performs reconstruction with the following update equation

\[\begin{eqnarray} f_{n+1} &:=& \frac{f_n}{P^T 1 + R'(f_n)} P^T\left[ \frac{g}{Pf_n} \right] \end{eqnarray}\]

where R’(f) is the gradient of the regularization term(s).

The CT geometry parameters and the CT volume parameters must be set prior to running this function. This reconstruction algorithms assumes the projection data, g, is Poisson distributed which is the correct model for SPECT data. CT projection data is not Poisson distributed because of the application of the -log

Parameters:

  • g (C contiguous float32 numpy or torch array) – projection data

  • f (C contiguous float32 numpy or torch array) – volume data

  • numIter (int) – number of iterations

  • filters (filterSequence object) – list of differentiable regularization filters

  • mask (C contiguous float32 numpy or torch array) – projection data to mask out bad data, etc.

Returns:

f, the same as the input with the same name

leapctype.tomographicModels.OSEM(self, g, f, numIter, numSubsets=1, filters=None, mask=None)

Ordered Subsets-Expectation Maximization reconstruction

The OSEM algorithm is performed using two nested loops. The inner loop is performed like the MLEM algorithm (see MLEM algorithm documentation for a description of the algorithm), but the successive updates of the reconstructed volume are done with a subset of the projection angles. Once every subset of is complete, the process starts over again. Using these ordered subsets reduces the time it takes for this algorithm to converge.

The CT geometry parameters and the CT volume parameters must be set prior to running this function. This reconstruction algorithms assumes the projection data, g, is Poisson distributed which is the correct model for SPECT data. CT projection data is not Poisson distributed because of the application of the -log

Parameters:

  • g (C contiguous float32 numpy or torch array) – projection data

  • f (C contiguous float32 numpy or torch array) – volume data

  • numIter (int) – number of iterations

  • numSubsets (int) – number of subsets

  • filters (filterSequence object) – list of differentiable regularization filters

  • mask (C contiguous float32 numpy or torch array) – projection data to mask out bad data, etc.

Returns:

f, the same as the input with the same name