pytomography.algorithms.preconditioned_gradient_ascent#

This module consists of preconditioned gradient ascent (PGA) algorithms: these algorithms are both statistical (since they depend on a likelihood function dependent on the imaging system) and iterative. Common clinical reconstruction algorithms, such as OSEM, correspond to a subclass of PGA algorithms. PGA algorithms are characterized by the update rule \(f^{n+1} = f^{n} + C^{n}(f^{n}) \left[\nabla_{f} L(g^n|f^{n}) - \beta \nabla_{f} V(f^{n}) \right]\) where \(L(g^n|f^{n})\) is the likelihood function, \(V(f^{n})\) is the prior function, \(C^{n}(f^{n})\) is the preconditioner, and \(\beta\) is a scalar used to scale the prior function.

Module 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]\).

LinearPreconditionedGradientAscentAlgorithm

Implementation of a special case of PreconditionedGradientAscentAlgorithm whereby \(C^{n}(f^n) = D^{n} f^{n}\)

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]\)

RBIEM

Implementation of the rescaled block iterative expectation maximum algorithm

RBIMAP

Implementation of the rescaled block iterative maximum a posteriori algorithm

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\).

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})\).

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.
class pytomography.algorithms.preconditioned_gradient_ascent.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.preconditioned_gradient_ascent.LinearPreconditionedGradientAscentAlgorithm(likelihood, prior=None, object_initial=None, addition_after_iteration=0, **kwargs)[source]#

Bases: PreconditionedGradientAscentAlgorithm

Implementation of a special case of PreconditionedGradientAscentAlgorithm whereby \(C^{n}(f^n) = D^{n} f^{n}\)

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.

abstract _linear_preconditioner_factor(n_iter, n_subset)[source]#

Implementation of object independent scaling factor \(D^{n}\) in \(C^{n}(f^{n}) = D^{n} f^{n}\)

Parameters:

  • n_iter (int) – iteration number

  • n_subset (int) – subset number

Raises:

NotImplementedError –

.

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

Computes the preconditioner \(C^{n}(f^n) = D^{n} \text{diag}\left(f^{n}\right)\) using the associated _linear_preconditioner_factor method.

Parameters:

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

  • n_iter (int) – Iteration \(n\)

  • n_subset (int) – Subset \(m\)

Returns:

Preconditioner factor

Return type:

torch.Tensor

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

Estimates the uncertainty of the sum of voxels in a reconstructed image. Calling this method requires a masked region mask as well as an instance of DataStorageCallback that has been used in a reconstruction algorithm: this data storage contains the estimated object and associated forward projection at each subiteration number.

Parameters:

  • mask (torch.Tensor) – Masked region of the reconstructed object: a boolean Tensor.

  • data_storage_callback (Callback) – Callback that has been used in a reconstruction algorithm.

  • 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.

  • post_recon_filter (Transform | None) –

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]

_compute_Q(input, data_storage_callback, n)[source]#

Computes the operation of \(Q\) on an input object; this is a helper function for compute_uncertainty. For more details, see the uncertainty paper.

Parameters:

  • input (torch.Tensor) – Object on which Q operates

  • data_storage_callback (Callback) – Data storage callback containing all objects and forward projections at each subiteration

  • n (int) – Subiteration number

Returns:

Resulting output object from the operation of \(Q\) on the input object

Return type:

torch.Tensor

_compute_B(input, data_storage_callback, n, include_additive_term=False)[source]#

Computes the operation of \(B\) on an input object; this is a helper function for compute_uncertainty. For more details, see the uncertainty paper.

Parameters:

  • input (torch.Tensor) – Object on which B operates

  • data_storage_callback (Callback) – Data storage callback containing all objects and forward projections at each subiteration

  • n (int) – Subiteration number

  • include_additive_term (bool, optional) – Whether or not to include uncertainty estimation for the additive term. Defaults to False.

Returns:

Resulting output projections from the operation of \(B\) on the input object

Return type:

torch.Tensor

class pytomography.algorithms.preconditioned_gradient_ascent.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.preconditioned_gradient_ascent.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.preconditioned_gradient_ascent.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.preconditioned_gradient_ascent.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.preconditioned_gradient_ascent.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.preconditioned_gradient_ascent.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.preconditioned_gradient_ascent.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.preconditioned_gradient_ascent.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.preconditioned_gradient_ascent.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]