pytomography.algorithms#

This module contains all the available reconstruction algorithms in PyTomography.

Submodules#

Package Contents#

Classes#

PreconditionedGradientAscentAlgorithm

Generic class for preconditioned gradient ascent algorithms: i.e. those that have the form \(f^{n+1} = f^{n} + C^{n}(f^{n}) \left[\nabla_{f} L(g^n|f^{n}) - \beta \nabla_{f} V(f^{n}) \right]\).

OSEM

Implementation of the ordered subset expectation maximum algorithm \(f^{n+1} = f^{n} + \frac{f^n}{H_n^T} \nabla_{f} L(g^n|f^{n})\).

OSMAPOSL

Implementation of the ordered subset maximum a posteriori one step late algorithm \(f^{n+1} = f^{n} + \frac{f^n}{H_n^T+\nabla_f V(f^n)} \left[ \nabla_{f} L(g^n|f^{n}) - \nabla_f V(f^n) \right]\)

BSREM

Implementation of the block sequential regularized expectation maximum algorithm \(f^{n+1} = f^{n} + \frac{\alpha(n)}{\omega_n H^T 1} \left[\nabla_{f} L(g^n|f^{n}) - \nabla_f V(f^n) \right]\)

KEM

Implementation of the ordered subset expectation maximum algorithm \(\alpha^{n+1} = \alpha^{n} + \frac{\alpha^n}{\tilde{H}_n^T} \nabla_{f} L(g^n|\alpha^{n})\) and where the final predicted object is \(f^n = K \hat{\alpha}^{n}\). The system matrix \(\tilde{H}\) includes the kernel transform \(K\).

RBIEM

Implementation of the rescaled block iterative expectation maximum algorithm

RBIMAP

Implementation of the rescaled block iterative maximum a posteriori algorithm

SART

Implementation of the SART algorithm. This algorithm takes as input the system matrix and projections (as opposed to a likelihood). This is an implementation of equation 3 of https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8506772/

PGAAMultiBedSPECT

Assistant class for performing reconstruction on multi-bed SPECT data. This class is a wrapper around a reconstruction algorithm that is called for each bed and then the results are stitched together.

MLEM

Implementation of the maximum likelihood expectation maximum algorithm \(f^{n+1} = f^{n} + \frac{f^n}{H^T} \nabla_{f} L(g|f^{n})\).

FilteredBackProjection

Implementation of filtered back projection reconstruction \(\hat{f} = \frac{\pi}{N_{\text{proj}}} \mathcal{R}^{-1}\mathcal{F}^{-1}\Pi\mathcal{F} g\) where \(N_{\text{proj}}\) is the number of projections, \(\mathcal{R}\) is the 3D radon transform, \(\mathcal{F}\) is the 2D Fourier transform (applied to each projection seperately), and \(\Pi\) is the filter applied in Fourier space, which is by default the ramp filter.

DIPRecon

Implementation of the Deep Image Prior reconstruction technique (see https://ieeexplore.ieee.org/document/8581448). This reconstruction technique requires an instance of a user-defined prior_network that implements two functions: (i) a fit method that takes in an object (\(x\)) which the network f(z;\theta) is subsequently fit to, and (ii) a predict function that returns the current network prediction \(f(z;\theta)\). For more details, see the Deep Image Prior tutorial.
class pytomography.algorithms.PreconditionedGradientAscentAlgorithm(likelihood, prior=None, object_initial=None, addition_after_iteration=0, **kwargs)[source]#

Generic class for preconditioned gradient ascent algorithms: i.e. those that have the form \(f^{n+1} = f^{n} + C^{n}(f^{n}) \left[\nabla_{f} L(g^n|f^{n}) - \beta \nabla_{f} V(f^{n}) \right]\).

Parameters:

  • likelihood (Likelihood) – Likelihood class that facilitates computation of \(L(g^n|f^{n})\) and its associated derivatives.

  • prior (Prior, optional) – Prior class that faciliates the computation of function \(V(f)\) and its associated derivatives. If None, then no prior is used Defaults to None.

  • object_initial (torch.Tensor | None, optional) – Initial object for reconstruction algorithm. If None, then an object with 1 in every voxel is used. Defaults to None.

  • addition_after_iteration (float, optional) – Value to add to the object after each iteration. This prevents image voxels getting “locked” at values of 0 for certain algorithms. Defaults to 0.

_set_n_subsets(n_subsets)[source]#

Sets the number of subsets used in the reconstruction algorithm.

Parameters:

n_subsets (int) – Number of subsets

abstract _compute_preconditioner(object, n_iter, n_subset)[source]#

Computes the preconditioner factor \(C^{n}(f^{n})\). Must be implemented by any reconstruction algorithm that inherits from this generic class.

Parameters:

  • object (torch.Tensor) – Object \(f^n\)

  • n_iter (int) – Iteration number

  • n_subset (int) – Subset number

Raises:

NotImplementedError –

.

Return type:

None

_compute_callback(n_iter, n_subset)[source]#

Method for computing callbacks after each reconstruction iteration

Parameters:

  • n_iter (int) – Number of iterations

  • n_subset (int) – Number of subsets

__call__(n_iters, n_subsets=1, n_subset_specific=None, callback=None)[source]#

_summary_

Parameters:

  • Args –

  • n_iters (int) – Number of iterations

  • n_subsets (int) – Number of subsets

  • n_subset_specific (int) – Ignore all updates except for this subset.

  • callback (Callback, optional) – Callback function to be called after each subiteration. Defaults to None.

Returns:

Reconstructed object.

Return type:

torch.Tensor

class pytomography.algorithms.OSEM(likelihood, object_initial=None)[source]#

Bases: LinearPreconditionedGradientAscentAlgorithm

Implementation of the ordered subset expectation maximum algorithm \(f^{n+1} = f^{n} + \frac{f^n}{H_n^T} \nabla_{f} L(g^n|f^{n})\).

Parameters:

  • likelihood (Likelihood) – Likelihood function \(L\).

  • object_initial (torch.Tensor | None, optional) – Initial object for reconstruction algorithm. If None, then an object with 1 in every voxel is used. Defaults to None.

_linear_preconditioner_factor(n_iter, n_subset)[source]#

Computes the linear preconditioner factor \(D^n = 1/H_n^T 1\)

Parameters:

  • n_iter (int) – iteration number

  • n_subset (int) – subset number

Returns:

linear preconditioner factor

Return type:

torch.Tensor

class pytomography.algorithms.OSMAPOSL(likelihood, object_initial=None, prior=None)[source]#

Bases: PreconditionedGradientAscentAlgorithm

Implementation of the ordered subset maximum a posteriori one step late algorithm \(f^{n+1} = f^{n} + \frac{f^n}{H_n^T+\nabla_f V(f^n)} \left[ \nabla_{f} L(g^n|f^{n}) - \nabla_f V(f^n) \right]\)

Parameters:

  • likelihood (Likelihood) – Likelihood function \(L\).

  • object_initial (torch.Tensor | None, optional) – Initial object for reconstruction algorithm. If None, then an object with 1 in every voxel is used. Defaults to None.

  • prior (Prior, optional) – Prior class that faciliates the computation of function \(V(f)\) and its associated derivatives. If None, then no prior is used. Defaults to None.

_compute_preconditioner(object, n_iter, n_subset)[source]#

Computes the preconditioner factor \(C^n(f^n) = \frac{f^n}{H_n^T+\nabla_f V(f^n)}\)

Parameters:

  • object (torch.Tensor) – Object estimate \(f^n\)

  • n_iter (int) – iteration number

  • n_subset (int) – subset number

Returns:

preconditioner factor.

Return type:

torch.Tensor

class pytomography.algorithms.BSREM(likelihood, object_initial=None, prior=None, relaxation_sequence=lambda _: ..., addition_after_iteration=0.0001)[source]#

Bases: LinearPreconditionedGradientAscentAlgorithm

Implementation of the block sequential regularized expectation maximum algorithm \(f^{n+1} = f^{n} + \frac{\alpha(n)}{\omega_n H^T 1} \left[\nabla_{f} L(g^n|f^{n}) - \nabla_f V(f^n) \right]\)

Parameters:

  • likelihood (Likelihood) – likelihood function \(L\)

  • object_initial (torch.Tensor | None, optional) – Initial object for reconstruction algorithm. If None, then an object with 1 in every voxel is used. Defaults to None.

  • prior (Prior, optional) – Prior class that faciliates the computation of function \(V(f)\) and its associated derivatives. If None, then no prior is used. Defaults to None.

  • relaxation_sequence (Callable, optional) – Relxation sequence \(\alpha(n)\) used to scale future updates. Defaults to 1 for all \(n\). Note that when this function is provided, it takes the iteration number (not the subiteration) so that e.g. if 4 iterations and 8 subsets are used, it would call \(\alpha(4)\) for all 8 subiterations of the final iteration.

  • addition_after_iteration (float, optional) – Value to add to the object after each iteration. This prevents image voxels getting “locked” at values of 0. Defaults to 1e-4.

_linear_preconditioner_factor(n_iter, n_subset)[source]#

Computes the linear preconditioner factor \(D^n = 1/(\omega_n H^T 1)\) where \(\omega_n\) corresponds to the fraction of subsets at subiteration \(n\).

Parameters:

  • n_iter (int) – iteration number

  • n_subset (int) – subset number

Returns:

linear preconditioner factor

Return type:

torch.Tensor

class pytomography.algorithms.KEM(likelihood, object_initial=None)[source]#

Bases: OSEM

Implementation of the ordered subset expectation maximum algorithm \(\alpha^{n+1} = \alpha^{n} + \frac{\alpha^n}{\tilde{H}_n^T} \nabla_{f} L(g^n|\alpha^{n})\) and where the final predicted object is \(f^n = K \hat{\alpha}^{n}\). The system matrix \(\tilde{H}\) includes the kernel transform \(K\).

Parameters:

  • likelihood (Likelihood) – Likelihood function \(L\).

  • object_initial (torch.Tensor | None, optional) – Initial object for reconstruction algorithm. If None, then an object with 1 in every voxel is used. Defaults to None.

_compute_callback(n_iter, n_subset)[source]#

Method for computing callbacks after each reconstruction iteration. This is reimplemented for KEM because the callback needs to be called on \(f^n = K \hat{\alpha}^{n}\) as opposed to \(\hat{\alpha}^{n}\)

Parameters:

  • n_iter (int) – Number of iterations

  • n_subset (int) – Number of subsets

__call__(*args, **kwargs)[source]#

Reimplementation of the call method such that \(f^n = K \hat{\alpha}^{n}\) is returned as opposed to \(\hat{\alpha}^{n}\)

Returns:

reconstructed object

Return type:

torch.Tensor

class pytomography.algorithms.RBIEM(likelihood, object_initial=None, prior=None)[source]#

Bases: LinearPreconditionedGradientAscentAlgorithm

Implementation of the rescaled block iterative expectation maximum algorithm

Parameters:

  • likelihood (Likelihood) – Likelihood function \(L\).

  • object_initial (torch.Tensor | None, optional) – Initial object for reconstruction algorithm. If None, then an object with 1 in every voxel is used. Defaults to None.

  • prior (Prior, optional) – Prior class that faciliates the computation of function \(V(f)\) and its associated derivatives. If None, then no prior is used. Defaults to None.

_compute_preconditioner(object, n_iter, n_subset)[source]#

Computes the preconditioner factor \(C^n(f^n) = \frac{f^n}{H_n^T+\nabla_f V(f^n)}\)

Parameters:

  • object (torch.Tensor) – Object estimate \(f^n\)

  • n_iter (int) – iteration number

  • n_subset (int) – subset number

Returns:

preconditioner factor.

Return type:

torch.Tensor

class pytomography.algorithms.RBIMAP(likelihood, object_initial=None, prior=None)[source]#

Bases: PreconditionedGradientAscentAlgorithm

Implementation of the rescaled block iterative maximum a posteriori algorithm

Parameters:

  • likelihood (Likelihood) – Likelihood function \(L\).

  • object_initial (torch.Tensor | None, optional) – Initial object for reconstruction algorithm. If None, then an object with 1 in every voxel is used. Defaults to None.

  • prior (Prior, optional) – Prior class that faciliates the computation of function \(V(f)\) and its associated derivatives. If None, then no prior is used. Defaults to None.

_compute_preconditioner(object, n_iter, n_subset)[source]#

Computes the preconditioner factor \(C^n(f^n) = \frac{f^n}{H_n^T+\nabla_f V(f^n)}\)

Parameters:

  • object (torch.Tensor) – Object estimate \(f^n\)

  • n_iter (int) – iteration number

  • n_subset (int) – subset number

Returns:

preconditioner factor.

Return type:

torch.Tensor

class pytomography.algorithms.SART(system_matrix, projections, additive_term=None, object_initial=None, relaxation_sequence=lambda _: ...)[source]#

Bases: PreconditionedGradientAscentAlgorithm

Implementation of the SART algorithm. This algorithm takes as input the system matrix and projections (as opposed to a likelihood). This is an implementation of equation 3 of https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8506772/

Parameters:

  • system_matrix (SystemMatrix) – System matrix for the imaging system.

  • projections (torch.Tensor) – Projections for the imaging system.

  • additive_term (torch.Tensor | None, optional) – Additive term for the imaging system. If None, then no additive term is used. Defaults to None.

  • object_initial (torch.Tensor | None, optional) – Initial object for reconstruction algorithm. If None, then an object with 1 in every voxel is used. Defaults to None.

  • relaxation_sequence (collections.abc.Callable) –

_compute_preconditioner(object, n_iter, n_subset)[source]#

Computes the preconditioner factor \(C^n(f^n) = \frac{1}{H_n^T+\nabla_f V(f^n)}\)

Parameters:

  • object (torch.Tensor) – Object estimate \(f^n\)

  • n_iter (int) – iteration number

  • n_subset (int) – subset number

Returns:

preconditioner factor.

Return type:

torch.Tensor

class pytomography.algorithms.PGAAMultiBedSPECT(files_NM, reconstruction_algorithms)[source]#

Bases: PreconditionedGradientAscentAlgorithm

Assistant class for performing reconstruction on multi-bed SPECT data. This class is a wrapper around a reconstruction algorithm that is called for each bed and then the results are stitched together.

Parameters:

  • files_NM (Sequence[str]) – Sequence of SPECT raw data paths corresponding to each likelihood

  • reconstruction_algorithm (Algorithm) – Reconstruction algorithm used for reconstruction of each bed position

  • reconstruction_algorithms (collections.abc.Sequence[object]) –

__call__(n_iters, n_subsets, callback=None)[source]#

Perform reconstruction of each bed position for specified iteraitons and subsets, and return the stitched image

Parameters:

  • n_iters (int) – Number of iterations to perform reconstruction for.

  • n_subsets (int) – Number of subsets to perform reconstruction for.

  • callback (Callback | Sequence[Callback] | None, optional) – Callback function. If a single Callback is given, then the callback is computed for the stitched image. If a sequence of callbacks is given, then it must be the same length as the number of bed positions; each callback is called on the reconstruction for each bed position. If None, no Callback is used. Defaults to None.

Returns:

_description_

Return type:

torch.Tensor

_compute_callback(n_iter, n_subset)[source]#

Computes the callback at iteration n_iter and subset n_subset.

Parameters:

  • n_iter (int) – Iteration number

  • n_subset (int) – Subset index

_finalize_callback()[source]#

Finalizes callbacks after reconstruction. This method is called after the reconstruction algorithm has finished.

compute_uncertainty(mask, data_storage_callbacks, subiteration_number=None, return_pct=False, include_additive_term=False)[source]#

Estimates the uncertainty in a mask (should be same shape as the stitched image). Calling this method requires a sequence of DataStorageCallback instances that have been used in a reconstruction algorithm: these data storage contain required information for each bed position.

Parameters:

  • mask (torch.Tensor) – Masked region of the reconstructed object: a boolean Tensor. This mask should be the same shape as the stitched object.

  • data_storage_callbacks (Sequence[Callback]) – Sequence of data storage callbacks used in reconstruction corresponding to each bed position.

  • subiteration_number (int | None, optional) – Subiteration number to compute the uncertainty for. If None, then computes the uncertainty for the last iteration. Defaults to None.

  • return_pct (bool, optional) – If true, then additionally returns the percent uncertainty for the sum of counts. Defaults to False.

  • include_additive_term (bool) – Whether or not to include uncertainty contribution from the additive term. This requires the additive_term_variance_estimate as an argument to the initialized likelihood. Defaults to False.

Returns:

Absolute uncertainty in the sum of counts in the masked region (if return_pct is False) OR absolute uncertainty and relative uncertainty in percent (if return_pct is True)

Return type:

float | Sequence[float]

class pytomography.algorithms.MLEM(likelihood, object_initial=None)[source]#

Bases: OSEM

Implementation of the maximum likelihood expectation maximum algorithm \(f^{n+1} = f^{n} + \frac{f^n}{H^T} \nabla_{f} L(g|f^{n})\).

Parameters:

  • likelihood (Likelihood) – Likelihood function \(L\).

  • object_initial (torch.Tensor | None, optional) – Initial object for reconstruction algorithm. If None, then an object with 1 in every voxel is used. Defaults to None.

__call__(n_iters, callback=None)[source]#

_summary_

Parameters:

  • Args –

  • n_iters (int) – Number of iterations

  • n_subsets (int) – Number of subsets

  • n_subset_specific (int) – Ignore all updates except for this subset.

  • callback (Callback, optional) – Callback function to be called after each subiteration. Defaults to None.

Returns:

Reconstructed object.

Return type:

torch.Tensor

class pytomography.algorithms.FilteredBackProjection(projections, system_matrix, filter=RampFilter)[source]#

Implementation of filtered back projection reconstruction \(\hat{f} = \frac{\pi}{N_{\text{proj}}} \mathcal{R}^{-1}\mathcal{F}^{-1}\Pi\mathcal{F} g\) where \(N_{\text{proj}}\) is the number of projections, \(\mathcal{R}\) is the 3D radon transform, \(\mathcal{F}\) is the 2D Fourier transform (applied to each projection seperately), and \(\Pi\) is the filter applied in Fourier space, which is by default the ramp filter.

Parameters:

  • projections (torch.Tensor) – projection data \(g\) to be reconstructed

  • system_matrix (SystemMatrix) – system matrix for the imaging system. In FBP, phenomena such as attenuation and PSF should not be implemented in the system matrix

  • filter (Callable, optional) – Additional Fourier space filter (applied after Ramp Filter) used during reconstruction.

__call__(projections)[source]#

Applies reconstruction

Returns:

Reconstructed object prediction

Return type:

torch.tensor

class pytomography.algorithms.DIPRecon(likelihood, prior_network, rho=0.003)[source]#

Implementation of the Deep Image Prior reconstruction technique (see https://ieeexplore.ieee.org/document/8581448). This reconstruction technique requires an instance of a user-defined prior_network that implements two functions: (i) a fit method that takes in an object (\(x\)) which the network f(z;\theta) is subsequently fit to, and (ii) a predict function that returns the current network prediction \(f(z;\theta)\). For more details, see the Deep Image Prior tutorial.

Parameters:

  • likelihood (Likelihood) – Initialized likelihood function for the imaging system considered

  • prior_network (nn.Module) – User defined prior network that implements the neural network \(f(z;\theta)\) that predicts an object given a prior image \(z\). This network also implements a fit method that takes in an object and fits the network to the object (for a specified number of iterations: SubIt2 in the paper).

  • rho (float, optional) – Value of \(\rho\) used in the optimization procedure. Larger values of \(rho\) give larger weight to the neural network, while smaller values of \(rho\) give larger weight to the EM updates. Defaults to 1.

_compute_callback(n_iter, n_subset)[source]#

Method for computing callbacks after each reconstruction iteration

Parameters:

  • n_iter (int) – Number of iterations

  • n_subset (int) – Number of subsets

__call__(n_iters, subit1, n_subsets_osem=1, callback=None)[source]#

Implementation of Algorithm 1 in https://ieeexplore.ieee.org/document/8581448. This implementation gives the additional option to use ordered subsets. The quantity SubIt2 specified in the paper is controlled by the user-defined prior_network class.

Parameters:

  • n_iters (int) – Number of iterations (MaxIt in paper)

  • subit1 (int) – Number of OSEM iterations before retraining neural network (SubIt1 in paper)

  • n_subsets_osem (int, optional) – Number of subsets to use in OSEM reconstruction. Defaults to 1.

Returns:

Reconstructed image

Return type:

torch.Tensor