corrct.solvers

Solvers for the tomographic reconstruction problem.

@author: Nicola VIGANÒ, Computational Imaging group, CWI, The Netherlands, and ESRF - The European Synchrotron, Grenoble, France

Module Contents

Classes

SolutionInfo

Reconstruction info.

Solver

Initialize the base solver class.

FBP

Implementation of the Filtered Back-Projection (FBP) algorithm.

SART

Solver class implementing the Simultaneous Algebraic Reconstruction Technique (SART) algorithm.

MLEM

Initialize the MLEM solver class.

SIRT

Initialize the SIRT solver class.

PDHG

Initialize the PDHG solver class.

Data

eps

NDArrayFloat

API

corrct.solvers.eps

None

corrct.solvers.NDArrayFloat

None

class corrct.solvers.SolutionInfo(method: str, max_iterations: int, tolerance: Union[float, numpy.floating, None], residual0: float = np.inf, residual0_cv: float = np.inf)[source]

Reconstruction info.

Initialization

method: str

None

iterations: int

None

max_iterations: int

None

residual0: Union[float, numpy.floating]

None

residual0_cv: Union[float, numpy.floating]

None

residuals: corrct.solvers.NDArrayFloat

None

residuals_cv: corrct.solvers.NDArrayFloat

None

tolerance: Union[float, numpy.floating, None]

None

property residuals_rel: corrct.solvers.NDArrayFloat
property residuals_cv_rel: corrct.solvers.NDArrayFloat
class corrct.solvers.Solver(verbose: bool = False, leave_progress: bool = True, relaxation: float = 1.0, tolerance: Optional[float] = None, data_term: Union[str, corrct.data_terms.DataFidelityBase] = 'l2', data_term_test: Union[str, corrct.data_terms.DataFidelityBase, None] = None)[source]

Bases: abc.ABC

Initialize the base solver class.

Parameters

verbose : bool, optional Turn on verbose output. The default is False. tolerance : Optional[float], optional Tolerance on the data residual for computing when to stop iterations. The default is None. relaxation : float, optional The relaxation length. The default is 1.0. data_term : Union[str, data_terms.DataFidelityBase], optional Data fidelity term for computing the data residual. The default is “l2”. data_term_test : Optional[data_terms.DataFidelityBase], optional The data fidelity to be used for the test set. If None, it will use the same as for the rest of the data. The default is None.

Initialization

info() str[source]

Return the solver info.

Returns

str Solver info string.

upper() str[source]

Return the upper case name of the solver.

Returns

str Upper case string name of the solver.

lower() str[source]

Return the lower case name of the solver.

Returns

str Lower case string name of the solver.

abstract __call__(A: corrct.operators.BaseTransform, b: corrct.solvers.NDArrayFloat, *args: Any, **kwds: Any) tuple[corrct.solvers.NDArrayFloat, corrct.solvers.SolutionInfo][source]

Execute the reconstruction of the data.

Parameters

A : operators.BaseTransform The projection operator. b : NDArrayFloat The data to be reconstructed.

Returns

Tuple[NDArrayFloat, SolutionInfo] The reconstruction and related information.

static _initialize_data_fidelity_function(data_term: Union[str, corrct.data_terms.DataFidelityBase]) corrct.data_terms.DataFidelityBase[source]
static _initialize_regularizer(regularizer: Union[corrct.regularizers.BaseRegularizer, None, collections.abc.Sequence[corrct.regularizers.BaseRegularizer]]) collections.abc.Sequence[corrct.regularizers.BaseRegularizer][source]
static _initialize_b_masks(b: corrct.solvers.NDArrayFloat, b_mask: Optional[corrct.solvers.NDArrayFloat], b_test_mask: Optional[corrct.solvers.NDArrayFloat]) tuple[Optional[corrct.solvers.NDArrayFloat], Optional[corrct.solvers.NDArrayFloat]][source]
class corrct.solvers.FBP(verbose: bool = False, leave_progress: bool = False, regularizer: Union[collections.abc.Sequence[corrct.regularizers.BaseRegularizer], corrct.regularizers.BaseRegularizer, None] = None, data_term: Union[str, corrct.data_terms.DataFidelityBase] = 'l2', fbp_filter: Union[str, corrct.solvers.NDArrayFloat, corrct.filters.Filter] = 'ramp', pad_mode: str = 'constant')[source]

Bases: corrct.solvers.Solver

Implementation of the Filtered Back-Projection (FBP) algorithm.

Initialization

Initialize the Filtered Back-Projection (FBP) algorithm.

Parameters

verbose : bool, optional Turn on verbose output. The default is False. leave_progress: bool, optional Leave the progress bar after the computation is finished. The default is True. regularizer : Sequence[regularizers.BaseRegularizer] | regularizers.BaseRegularizer | None, optional NOT USED, only exposed for compatibility reasons. data_term : Union[str, data_terms.DataFidelityBase], optional NOT USED, only exposed for compatibility reasons. fbp_filter : Union[str, NDArrayFloat], optional FBP filter to use. Either a string from scikit-image’s list of iradon filters, or an array. The default is “ramp”. pad_mode: str, optional The padding mode to use for the linear convolution. The default is “constant”.

info() str[source]

Return the solver info.

Returns

str Solver info string.

__call__(A: corrct.operators.BaseTransform, b: corrct.solvers.NDArrayFloat, iterations: int = 0, x0: Optional[corrct.solvers.NDArrayFloat] = None, lower_limit: Union[float, corrct.solvers.NDArrayFloat, None] = None, upper_limit: Union[float, corrct.solvers.NDArrayFloat, None] = None, x_mask: Optional[corrct.solvers.NDArrayFloat] = None, b_mask: Optional[corrct.solvers.NDArrayFloat] = None) tuple[corrct.solvers.NDArrayFloat, corrct.solvers.SolutionInfo][source]

Reconstruct the data, using the FBP algorithm.

Parameters

A : BaseTransform Projection operator. b : NDArrayFloat Data to reconstruct. iterations : int Number of iterations. x0 : Optional[NDArrayFloat], optional Initial solution. The default is None. lower_limit : Union[float, NDArrayFloat], optional Lower clipping value. The default is None. upper_limit : Union[float, NDArrayFloat], optional Upper clipping value. The default is None. x_mask : Optional[NDArrayFloat], optional Solution mask. The default is None. b_mask : Optional[NDArrayFloat], optional Data mask. The default is None.

Raises

ValueError In case the data is 1D.

Returns

Tuple[NDArrayFloat, SolutionInfo] The reconstruction, and None.

class corrct.solvers.SART(verbose: bool = False, leave_progress: bool = True, relaxation: float = 1.0, tolerance: Optional[float] = None, data_term: Union[str, corrct.data_terms.DataFidelityBase] = 'l2', data_term_test: Union[str, corrct.data_terms.DataFidelityBase, None] = None)[source]

Bases: corrct.solvers.Solver

Solver class implementing the Simultaneous Algebraic Reconstruction Technique (SART) algorithm.

Initialization

compute_residual(A: Callable, b: corrct.solvers.NDArrayFloat, x: corrct.solvers.NDArrayFloat, A_num_rows: int, b_mask: Optional[corrct.solvers.NDArrayFloat]) corrct.solvers.NDArrayFloat[source]

Compute the solution residual.

Parameters

A : Callable The forward projector. b : NDArrayFloat The detector data. x : NDArrayFloat The current solution A_num_rows : int The number of projections. b_mask : Optional[NDArrayFloat] The mask to apply

Returns

NDArrayFloat The residual.

__call__(A: Union[Callable[[numpy.typing.NDArray, int], numpy.typing.NDArray], corrct.projectors.ProjectorUncorrected], b: corrct.solvers.NDArrayFloat, iterations: int, A_num_rows: Optional[int] = None, At: Optional[Callable] = None, x0: Optional[corrct.solvers.NDArrayFloat] = None, lower_limit: Union[float, corrct.solvers.NDArrayFloat, None] = None, upper_limit: Union[float, corrct.solvers.NDArrayFloat, None] = None, x_mask: Optional[corrct.solvers.NDArrayFloat] = None, b_mask: Optional[corrct.solvers.NDArrayFloat] = None) tuple[corrct.solvers.NDArrayFloat, corrct.solvers.SolutionInfo][source]

Reconstruct the data, using the SART algorithm.

Parameters

A : Union[Callable, BaseTransform] Projection operator. b : NDArrayFloat Data to reconstruct. iterations : int Number of iterations. A_num_rows : int Number of projections. x0 : Optional[NDArrayFloat], optional Initial solution. The default is None. At : Callable, optional The back-projection operator. This is only needed if the projection operator does not have an adjoint. The default is None. lower_limit : Union[float, NDArrayFloat], optional Lower clipping value. The default is None. upper_limit : Union[float, NDArrayFloat], optional Upper clipping value. The default is None. x_mask : Optional[NDArrayFloat], optional Solution mask. The default is None. b_mask : Optional[NDArrayFloat], optional Data mask. The default is None.

Returns

Tuple[NDArrayFloat, SolutionInfo] The reconstruction, and the residuals.

class corrct.solvers.MLEM(verbose: bool = False, leave_progress: bool = True, tolerance: Optional[float] = None, regularizer: Union[collections.abc.Sequence[corrct.regularizers.BaseRegularizer], corrct.regularizers.BaseRegularizer, None] = None, data_term: Union[str, corrct.data_terms.DataFidelityBase] = 'kl', data_term_test: Union[str, corrct.data_terms.DataFidelityBase, None] = None)[source]

Bases: corrct.solvers.Solver

Initialize the MLEM solver class.

This class implements the Maximul Likelihood Expectation Maximization (MLEM) algorithm.

Parameters

verbose : bool, optional Turn on verbose output. The default is False. leave_progress: bool, optional Leave the progress bar after the computation is finished. The default is True. tolerance : Optional[float], optional Tolerance on the data residual for computing when to stop iterations. The default is None. regularizer : Sequence[regularizers.BaseRegularizer] | regularizers.BaseRegularizer | None, optional Regularizer to be used. The default is None. data_term : Union[str, data_terms.DataFidelityBase], optional Data fidelity term for computing the data residual. The default is “l2”. data_term_test : Optional[data_terms.DataFidelityBase], optional The data fidelity to be used for the test set. If None, it will use the same as for the rest of the data. The default is None.

Initialization

info() str[source]

Return the MLEM info.

Returns

str info string.

__call__(A: corrct.operators.BaseTransform, b: corrct.solvers.NDArrayFloat, iterations: int, x0: Optional[corrct.solvers.NDArrayFloat] = None, lower_limit: Union[float, corrct.solvers.NDArrayFloat, None] = None, upper_limit: Union[float, corrct.solvers.NDArrayFloat, None] = None, x_mask: Optional[corrct.solvers.NDArrayFloat] = None, b_mask: Optional[corrct.solvers.NDArrayFloat] = None, b_test_mask: Optional[corrct.solvers.NDArrayFloat] = None) tuple[corrct.solvers.NDArrayFloat, corrct.solvers.SolutionInfo][source]

Reconstruct the data, using the MLEM algorithm.

Parameters

A : BaseTransform Projection operator. b : NDArrayFloat Data to reconstruct. iterations : int Number of iterations. x0 : Optional[NDArrayFloat], optional Initial solution. The default is None. lower_limit : Union[float, NDArrayFloat], optional Lower clipping value. The default is None. upper_limit : Union[float, NDArrayFloat], optional Upper clipping value. The default is None. x_mask : Optional[NDArrayFloat], optional Solution mask. The default is None. b_mask : Optional[NDArrayFloat], optional Data mask. The default is None. b_test_mask : Optional[NDArrayFloat], optional Test data mask. The default is None.

Returns

Tuple[NDArrayFloat, SolutionInfo] The reconstruction, and the residuals.

class corrct.solvers.SIRT(verbose: bool = False, leave_progress: bool = True, relaxation: float = 1.95, tolerance: Optional[float] = None, regularizer: Union[collections.abc.Sequence[corrct.regularizers.BaseRegularizer], corrct.regularizers.BaseRegularizer, None] = None, data_term: Union[str, corrct.data_terms.DataFidelityBase] = 'l2', data_term_test: Union[str, corrct.data_terms.DataFidelityBase, None] = None)[source]

Bases: corrct.solvers.Solver

Initialize the SIRT solver class.

This class implements the Simultaneous Iterative Reconstruction Technique (SIRT) algorithm.

Parameters

verbose : bool, optional Turn on verbose output. The default is False. leave_progress: bool, optional Leave the progress bar after the computation is finished. The default is True. tolerance : Optional[float], optional Tolerance on the data residual for computing when to stop iterations. The default is None. relaxation : float, optional The relaxation length. The default is 1.95. regularizer : Sequence[regularizers.BaseRegularizer] | regularizers.BaseRegularizer | None, optional Regularizer to be used. The default is None. data_term : Union[str, data_terms.DataFidelityBase], optional Data fidelity term for computing the data residual. The default is “l2”. data_term_test : Optional[data_terms.DataFidelityBase], optional The data fidelity to be used for the test set. If None, it will use the same as for the rest of the data. The default is None.

Initialization

info() str[source]

Return the SIRT info.

Returns

str SIRT info string.

__call__(A: corrct.operators.BaseTransform, b: corrct.solvers.NDArrayFloat, iterations: int, x0: Optional[corrct.solvers.NDArrayFloat] = None, lower_limit: Union[float, corrct.solvers.NDArrayFloat, None] = None, upper_limit: Union[float, corrct.solvers.NDArrayFloat, None] = None, x_mask: Optional[corrct.solvers.NDArrayFloat] = None, b_mask: Optional[corrct.solvers.NDArrayFloat] = None, b_test_mask: Optional[corrct.solvers.NDArrayFloat] = None) tuple[corrct.solvers.NDArrayFloat, corrct.solvers.SolutionInfo][source]

Reconstruct the data, using the SIRT algorithm.

Parameters

A : BaseTransform Projection operator. b : NDArrayFloat Data to reconstruct. iterations : int Number of iterations. x0 : Optional[NDArrayFloat], optional Initial solution. The default is None. lower_limit : Union[float, NDArrayFloat], optional Lower clipping value. The default is None. upper_limit : Union[float, NDArrayFloat], optional Upper clipping value. The default is None. x_mask : Optional[NDArrayFloat], optional Solution mask. The default is None. b_mask : Optional[NDArrayFloat], optional Data mask. The default is None. b_test_mask : Optional[NDArrayFloat], optional Test data mask. The default is None.

Returns

Tuple[NDArrayFloat, SolutionInfo] The reconstruction, and the residuals.

class corrct.solvers.PDHG(verbose: bool = False, leave_progress: bool = True, tolerance: Optional[float] = None, relaxation: float = 0.95, regularizer: Union[collections.abc.Sequence[corrct.regularizers.BaseRegularizer], corrct.regularizers.BaseRegularizer, None] = None, data_term: Union[str, corrct.data_terms.DataFidelityBase] = 'l2', data_term_test: Union[str, corrct.data_terms.DataFidelityBase, None] = None)[source]

Bases: corrct.solvers.Solver

Initialize the PDHG solver class.

PDHG stands for primal-dual hybrid gradient algorithm from Chambolle and Pock.

Parameters

verbose : bool, optional Turn on verbose output. The default is False. leave_progress: bool, optional Leave the progress bar after the computation is finished. The default is True. tolerance : Optional[float], optional Tolerance on the data residual for computing when to stop iterations. The default is None. relaxation : float, optional The relaxation length. The default is 0.95. regularizer : Sequence[regularizers.BaseRegularizer] | regularizers.BaseRegularizer | None, optional Regularizer to be used. The default is None. data_term : Union[str, data_terms.DataFidelityBase], optional Data fidelity term for computing the data residual. The default is “l2”. data_term_test : Optional[data_terms.DataFidelityBase], optional The data fidelity to be used for the test set. If None, it will use the same as for the rest of the data. The default is None.

Initialization

info() str[source]

Return the PDHG info.

Returns

str PDHG info string.

static _initialize_data_fidelity_function(data_term: Union[str, corrct.data_terms.DataFidelityBase])[source]
power_method(A: corrct.operators.BaseTransform, b: corrct.solvers.NDArrayFloat, iterations: int = 5) tuple[numpy.floating, collections.abc.Sequence[int], numpy.typing.DTypeLike][source]

Compute the l2-norm of the operator A, with the power method.

Parameters

A : BaseTransform Operator whose l2-norm needs to be computed. b : NDArrayFloat The data vector. iterations : int, optional Number of power method iterations. The default is 5.

Returns

Tuple[float, Tuple[int], DTypeLike] The l2-norm of A, and the shape and type of the solution.

_get_data_sigma_tau_unpreconditioned(A: corrct.operators.BaseTransform, b: corrct.solvers.NDArrayFloat)[source]
__call__(A: corrct.operators.BaseTransform, b: corrct.solvers.NDArrayFloat, iterations: int, x0: Optional[corrct.solvers.NDArrayFloat] = None, lower_limit: Union[float, corrct.solvers.NDArrayFloat, None] = None, upper_limit: Union[float, corrct.solvers.NDArrayFloat, None] = None, x_mask: Optional[corrct.solvers.NDArrayFloat] = None, b_mask: Optional[corrct.solvers.NDArrayFloat] = None, b_test_mask: Optional[corrct.solvers.NDArrayFloat] = None, precondition: bool = True) tuple[corrct.solvers.NDArrayFloat, corrct.solvers.SolutionInfo][source]

Reconstruct the data, using the PDHG algorithm.

Parameters

A : BaseTransform Projection operator. b : NDArrayFloat Data to reconstruct. iterations : int Number of iterations. x0 : Optional[NDArrayFloat], optional Initial solution. The default is None. lower_limit : Union[float, NDArrayFloat], optional Lower clipping value. The default is None. upper_limit : Union[float, NDArrayFloat], optional Upper clipping value. The default is None. x_mask : Optional[NDArrayFloat], optional Solution mask. The default is None. b_mask : Optional[NDArrayFloat], optional Data mask. The default is None. b_test_mask : Optional[NDArrayFloat], optional Test data mask. The default is None. precondition : bool, optional Whether to use the preconditioned version of the algorithm. The default is True.

Returns

Tuple[NDArrayFloat, SolutionInfo] The reconstruction, and the residuals.