Pipelines¶

Reconstruction¶

mumott.pipelines.reconstruction.run_mitra(data_container, use_absorbances=True, use_sirt_weights=True, use_gpu=False, maxiter=20, ftol=None, **kwargs)[source]¶

Reconstruction pipeline for the Modular Iterative Tomographic Reconstruction Algorithm (MITRA). This is a versatile, configureable interface for tomographic reconstruction that allows for various optimizers, projectors, loss functions and regularizers to be supplied.

This is meant as a convenience interface for intermediate or advanced users to create customized reconstruction pipelines.

Parameters

data_container (DataContainer) – The DataContainer from loading the data set of interest.
use_absorbances (bool) – If True, the reconstruction will use the absorbances calculated from the diode, or absorbances provided via a keyword argument. If False, the data in data_container.data will be used.
use_sirt_weights (bool) – If True (default), SIRT or tensor SIRT weights will be computed for use in the reconstruction.
use_gpu (bool) – Whether to use GPU resources in computing the projections. Default is False, which means SAXSProjectorBilinear. If set to True, the method will use SAXSProjectorCUDA.
maxiter (int) – Maximum number of iterations for the gradient descent solution.
ftol (Optional[float]) – Tolerance for the change in the loss function. Default is None, in which case the reconstruction will terminate once the maximum number of iterations have been performed.
kwargs – Miscellaneous keyword arguments. See notes for details.

Notes

Many options can be specified through kwargs. These include:

Projector
The projector class to use.

absorbances
If use_absorbances is set to True, these absorbances will be used instead of ones calculated from the diode.

preconditioner_cutoff
The cutoff to use when computing the SIRT preconditioner. Default value is 0.1, which will lead to a roughly ellipsoidal mask.

weights_cutoff
The cutoff to use when computing the SIRT weights. Default value is 0.1, which will clip some projection edges.

BasisSet
The basis set class to use. If not provided TrivialBasis will be used for absorbances and GaussianKernels for other data.

basis_set_kwargs
Keyword arguments for BasisSet.

ResidualCalculator
The residual calculator class to use. If not provided, then GradientResidualCalculator will be used.

residual_calculator_kwargs
Keyword arguments for ResidualCalculator.

LossFunction
The loss function class to use. If not provided SquaredLoss will be used.

loss_function_kwargs
Keyword arguments for LossFunction.

Regularizers
A list of dictionaries with three entries, a name (str), a regularizer object, and a regularization weight (float); used by loss_function.add_regularizer().

Optimizer
The optimizer class to use. If not provided GradientDescent will be used. By default, the keyword argument nestorov_weight is set to 0.95, and enforce_non_negativity is True

optimizer_kwargs
Keyword arguments for Optimizer.

mumott.pipelines.reconstruction.run_sirt(data_container, use_absorbances=True, use_gpu=False, maxiter=20, **kwargs)[source]¶

A reconstruction pipeline for the SIRT algorithm, which uses a gradient preconditioner and a set of weights for the projections to achieve fast convergence. Generally, one varies the number of iterations until a good reconstruction is obtained.

Advanced users may wish to also modify the preconditioner_cutoff and weights_cutoff keyword arguments.

Parameters

data_container (DataContainer) – The DataContainer from loading the data set of interest.
use_absorbances (bool) – If True, the reconstruction will use the absorbances calculated from the diode, or absorbances provided via a keyword argument. If False, the data in data_container.data will be used.
use_gpu (bool) – Whether to use GPU resources in computing the projections. Default is False. If set to True, the method will use SAXSProjectorCUDA.
maxiter (int) – Maximum number of iterations for the gradient descent solution.
kwargs – Miscellaneous keyword arguments. See notes for details.

Notes

Many options can be specified through kwargs. These include:

Projector
The projector class to use.

preconditioner_cutoff
The cutoff to use when computing the SIRT preconditioner. Default value is 0.1, which will lead to a roughly ellipsoidal mask.

weights_cutoff
The cutoff to use when computing the SIRT weights. Default value is 0.1, which will clip some projection edges.

absorbances
If use_absorbances is set to True, these absorbances will be used instead of ones calculated from the diode.

BasisSet
The basis set class to use. If not provided TrivialBasis will be used for absorbances and GaussianKernels for other data.

basis_set_kwargs
Keyword arguments for BasisSet.

no_tqdm
Used to avoid a tqdm progress bar in the optimizer.

mumott.pipelines.reconstruction.run_sigtt(data_container, use_gpu=False, maxiter=20, ftol=0.01, regularization_weight=0.0001, **kwargs)[source]¶

A reconstruction pipeline for the SIGTT algorithm, which uses a gradient and a regularizer to accomplish reconstruction.

Parameters

data_container (DataContainer) – The DataContainer from loading the data set of interest.
use_gpu (bool) – Whether to use GPU resources in computing the projections. Default is False. If set to True, the method will use SAXSProjectorCUDA.
maxiter (int) – Maximum number of iterations for the gradient descent solution.
ftol (float) – Tolerance for the change in the loss function. Default is None, in which case the reconstruction will terminate once the maximum number of iterations have been performed.
regularization_weight (float) – Regularization weight for the default Laplacian regularizer. Ignored if a loss function is provided.
kwargs – Miscellaneous keyword arguments. See notes for details.

Notes

Many options can be specified through kwargs. Miscellaneous ones are passed to the optimizer. Specific keywords include:

Projector
The projector class to use.

BasisSet
The basis set class to use. If not provided SphericalHarmonics will be used.

basis_set_kwargs
Keyword arguments for BasisSet.

ResidualCalculator
The residual_calculator class to use. If not provided GradientResidualCalculator will be used.

residual_calculator_kwargs
Keyword arguments for ResidualCalculator.

LossFunction
The loss function class to use. If not provided SquaredLoss will be used.

loss_function_kwargs
Keyword arguments for LossFunction.

Regularizers
A list of dictionaries with three entries, a name (str), a regularizer object, and a regularization weight (float); used by loss_function.add_regularizer(). By default, a Laplacian with the weight regularization_weight will be used. If other regularizers are specified, this will be overridden.

Optimizer
The optimizer class to use. If not provided LBFGS will be used.

optimizer_kwargs
Keyword arguments for Optimizer.

Alignment¶

mumott.pipelines.alignment.run_cross_correlation_alignment(data_container, ignored_subset=None, projection_cropping=(slice(None, None, None), slice(None, None, None)), reconstruction_pipeline=<function run_mitra>, reconstruction_pipeline_kwargs=None, use_gpu=False, use_absorbances=True, maxiter=20, upsampling=1, shift_tolerance=None, shift_cutoff=None, relative_sample_size=1.0, relaxation_weight=0.0, center_of_mass_shift_weight=0.0, align_j=True, align_k=True)[source]¶

A pipeline for alignment using the phase cross-correlation method as implemented by scikit-image.

For details on the cross-correlation algorithm, see this article by Guizar-Sicairos et al., (2008). Briefly, the algorithm calculates the cross-correlation between a reference image (the data) and the corresponding projection of a reconstruction, and finds the shift that would result in maximal correlation between the two. It supports large upsampling factors with very little computational overhead.

This implementation applies this algorithm to a randomly sampled subset of the projections in each iteration, and adds to this two smoothing terms – a stochastic relaxation term, and a shift toward the center of mass of the reconstruction. These terms are added partly to reduce the determinism in the algorithm, and partly to improve the performance when no upsampling is used.

The relaxation term is given by

\[d(x_i) = \text{sgn}(x_i) \cdot \text{max} (0, \vert x_i \vert - \vert \mathcal{N}(\overline{\mu}(x), \sigma(x)) \vert)\]

where \(x_i\) is a given offset and \(\mathcal{N}(\mu, \sigma)\) is a random variable from a normal distribution with mean \(\mu\) and standard deviation \(\sigma\). \(x_i\) is then updated by

\[x_i \leftarrow x_i + \lambda \cdot \text{sign}(d(x_i)) \cdot \text{max}(1, \vert d(x_i) \vert)\]

where \(\lambda\) is the relaxation_weight.

The shift toward the center of mass is given by

\[t(x_i) = \mathbf{v_i} \cdot (\mathbf{v_i}(\text{CoM}(P_i) + x_i)_j - \text{CoM}(R))\]

where \(\mathbf{v_i}\) is the three-dimensional basis vector that maps out \(x_i\). This expression assumes that the basis vectors of the two shift directions are orthogonal, but the general expression is similar. The term \(t(x_i)\) is then used to update \(x_i\) similarly to \(d(x_i)\).

Parameters

data_container (DataContainer) – The data container from loading the data set of interest. Note that the offset factors in data_container.geometry will be modified during the alignment.
ignored_subset (Optional[Set[int]]) – A subset of projection numbers which will not have their alignment modified. The subset is still used in the reconstruction.
projection_cropping (Tuple[slice, slice]) – A tuple of two slices (slice), which specify the cropping of the projection and data. For example, to clip the first and last 5 pixels in each direction, set this parameter to (slice(5, -5), slice(5, -5)).
reconstruction_pipeline (Callable) – A callable, typically from the reconstruction pipelines, that performs the reconstruction at each alignment iteration. Must return a dictionary with a entry labelled 'result', which has an entry labelled 'x' containing the reconstruction. Additionally, it must expose a 'weights' entry containing the weights used during the reconstruction as well as a Projector object under the keyword 'projector'. If the pipeline supports the use_absorbances keyword argument, then an absorbances entry must also be exposed. If the pipeline supports using multi-channel data, absorbances, then a basis set object must be available under 'basis_set'.
reconstruction_pipeline_kwargs (Optional[Dict[str, any]]) – Keyword arguments to pass to reconstruction_pipeline. If 'data_container' or 'use_gpu' are set as keys, they will override data_container and use_gpu
use_gpu (bool) – Whether to use GPU resources in computing the reconstruction. Default is False. Will be overridden if set in reconstruction_pipeline_kwargs.
use_absorbances (bool) – Whether to use the absorbances to compute the reconstruction and align the projections. Default is True. Will be overridden if set in reconstruction_pipeline_kwargs.
maxiter (int) – Maximum number of iterations for the alignment.
upsampling (int) – Upsampling factor during alignment. If used, any masking in data_container.weights will be ignored, but projection_clipping will still be used. The suggested range of use is [1, 20].
shift_tolerance (Optional[float]) – Tolerance for the the maximal shift distance of each iteration of the alignment. The alignment will terminate when the maximal shift falls below this value. The maximal shift is the largest Euclidean distance any one projection is shifted by. Default value is 1 / upsampling.
shift_cutoff (Optional[float]) – Largest permissible shift due to cross-correlation in each iteration, as measured by the Euclidean distance. Larger shifts will be rescaled so as to not exceed this value. Default value is 5 / upsampling.
relative_sample_size (float) – Fraction of projections to align in each iteration. At each alignment iteration, ceil(number_of_projections * relative_sample_size) will be randomly selected for alignment. If set to 1, all projections will be aligned at each iteration.
relaxation_weight (float) – A relaxation parameter for stochastic relaxation; the larger this weight is, the more shifts will tend toward the mean shift in each direction. The relaxation step size in each direction at each iteration cannot be larger than this weight. This is \(\lambda\) in the expression given above.
center_of_mass_shift_weight (float) – A parameter that controls the tendency for the projection center of mass to be shifted toward the reconstruction center of mass. The relaxation step size in each direction at each iteration cannot be larger than this weight.
align_j (bool) – Whether to align in the j direction. Default is True.
align_k (bool) – Whether to align in the k direction. Default is True.

Return type

Dict[str, Any]

Returns

A dictionary with three entries for inspection:

reconstruction: The reconstruction used in the last alignment step.
projections: Projections of the reconstruction.
reference: The reference image derived from the data used to align the projections.

mumott.pipelines.alignment.shift_center_of_reconstruction(geometry, shift_vector=(0.0, 0.0, 0.0))[source]¶

This utility function will shift the offsets in the geometry based on a three-dimensional shift vector. Use this to reposition the reconstruction within the volume.

Parameters

geometry (Geometry) – A Geometry instance. Its j_offsets and k_offsets are modified in-place.
shift_vector (Tuple[float]) – A tuple that indicates the direction and magnitude of the desired shift in (x, y, z), in units of voxels.

Return type

None

Example

>>> import numpy as np
>>> from scipy.spatial.transform import Rotation
>>> from mumott.core.geometry import Geometry, GeometryTuple
>>> from mumott.pipelines.alignment import shift_center_of_reconstruction
>>> geo = Geometry()
>>> geo.append(GeometryTuple(rotation=np.eye(3), j_offset=0., k_offset=0.))
>>> geo.append(GeometryTuple(rotation=Rotation.from_euler('y', np.pi/4).as_matrix(),
                             j_offset=0.,
                             k_offset=0.))
>>> geo.append(GeometryTuple(rotation=Rotation.from_euler('z', np.pi/4).as_matrix(),
                             j_offset=0.,
                             k_offset=0.))
>>> print(geo.j_direction_0)
[1. 0. 0.]
>>> print(geo.k_direction_0)
[0. 0. 1.]
>>> shift_center_of_reconstruction(geo, shift_vector=(-5., 1.7, 3.))
>>> print(geo[0].j_offset, geo[0].k_offset)
-5.0 3.0
>>> print(geo[1].j_offset, geo[1].k_offset)
-1.4142... 5.6568...
>>> print(geo[2].j_offset, geo[2].k_offset)
-4.737... 3.0