Pipelines¶

A number of pipelines are included in mumott, which pre-define common workflows.

The standard reconstruction pipelines include:

Filtered back-projection (FBP): Standard implementation of the inverse Radon formula for reconstruction of absorptivity from diode data.
Modular Iterative Tomographic Reconstruction Algorithm (MITRA): Reliable tensor-tomographic reconstruction algorithm which uses Gaussian kernels, normalizing weights and preconditioners, and momentum-accelerated gradient descent. Highly configurable.
Simultaneous Iterative Reconstruction Tomography (SIRT): Standard tomographic algorithm which reconstructs through preconditioned gradient descent. Configured mainly by modifying the maximum number of iterations.
Spherical Integral Geometric Tensor Tomography (SIGTT): Tensor tomography using a spherical harmonic representation, a quasi-Newton solver, and Laplacian regularization. Highly configurable.
Discrete Directions (DD). Bins reciprocal space directions on the reciprocal space sphere and solves for each of them separately using SIRT. Easy to get started with.

The asynchronous pipelines include:

Tensor SIRT: Asynchronous implementation of SIRT for tensor tomography, using weights and a preconditioner which account for the representation used.
Momentum Total variation Reconstruction (MOTR): Similar to Tensor SIRT, but uses Nestorov momentum, and optionally an L1 and total variation regularizer.
Robust And Denoised Tensor Tomography (RADTT): Uses a Huber norm for robust regression, which makes it less sensitive to outliers and noise than other pipelines. Includes total variation regularization.

The asynchronous pipelines have sparse counterparts, which function similarly, but are more efficient with respect to memory at the expense of only optimizing a few basis functions for each projection.

The alignment pipelines are:

Phase matching alignment: An alignment pipeline that uses cross-correlation to align data.
Optical flow alignment: Uses the optical flow algorithm described in [Odstrcil2019] to align data in multiple steps.

Reconstruction¶

mumott.pipelines.filtered_back_projection.run_fbp(data_container, use_gpu=False, fbp_axis='inner', filter_type='Ram-Lak', **kwargs)[source]¶

This pipeline is used to compute the filtered back projection of the absorbances calculated from the diode. This allows for a quick, one-step solution to the problem of scalar tomography.

Parameters

data_container (DataContainer) – The DataContainer holding 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.
fbp_axis (str) – Default is 'inner', the value depends on how the sample is mounted to the holder. Typically, the inner axis is the rotation axis while the 'outer' axis refers to the tilt axis.
filter_type (str) – Default is 'Ram-Lak', a high-pass filter. Other options are 'Hamming', 'Hann', 'Shepp-Logan' and 'cosine'.
kwargs – Miscellaneous keyword arguments. See notes for details.

Notes

Three possible kwargs can be provided:

Projector
The projector class to use.

normalization_percentile
The normalization percentile to use for the transmittivity calculation. See get_transmittivities for details.

transmittivity_cutoff
The cutoffs to use for the transmittivity calculation. See get_transmittivities for details.

Return type: dict
Returns: A dictionary with the entry 'result', with an entry 'x', containing a filtered back projection reconstruction of the absorptivity calculated from the diode.

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 SAXSProjector. 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, enforce_non_negativity=False, **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.
enforce_non_negativity (bool) – Enforces strict positivity on all the coefficients. Should only be used with local or scalar representations. Default value is False.
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.

mumott.pipelines.reconstruction.run_discrete_directions(data_container, directions, use_gpu=False, maxiter=20, no_tqdm=False)[source]¶

A reconstruction pipeline for the discrete directions algorithm, which is similar the the algorithms first descibed in [Schaff2015].

Parameters

data_container (DataContainer) – The DataContainer from loading the data set of interest.
directions (ndarray[Any, dtype[float]]) – A N by 3 Numpy array of unit-vectors descibing a grid covering the half unit sphere.
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.
no_tqdm (bool) – Flag whether ot not to print a progress bar for the reconstruction.

mumott.pipelines.reconstruction.run_group_lasso(data_container, regularization_parameter, step_size_parameter=None, x0=None, basis_set=None, ell_max=8, use_gpu=False, maxiter=100, enforce_non_negativity=False, no_tqdm=False)[source]¶

A reconstruction pipeline to do least squares reconstructions regularized with the group-lasso regularizer and solved with the Iterative Soft-Thresholding Algorithm (ISTA), a proximal gradient decent method. This reconstruction automatically masks out voxels with zero scattering but needs the regularization weight as input.

Parameters

data_container (DataContainer) – The DataContainer containing the data set of interest.
regularization_parameter (float) – Scalar weight of the regularization term. Should be optimized by performing reconstructions for a range of possible values.
step_size_parameter (Optional[float]) – Step size parameter of the reconstruction. If no value is given, a largest-safe value is estimated.
x0 (Optional[ndarray[Any, dtype[float]]]) – Starting guess for the solution. By default (None) the coefficients are initialized with zeros.
basis_set (Optional[BasisSet]) – Optionally a basis set can be specified. By default (None) SphericalHarmonics is used.
ell_max (int) – If no basis set is given, this is the maximum spherical harmonics order used in the generated basis set.
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) – Number of iterations for the ISTA optimization. No stopping rules are implemented.
enforce_non_negativity (bool) – Whether or not to enforce non-negativitu of the solution coefficients.
no_tqdm (bool) – Flag whether ot not to print a progress bar for the reconstruction.

mumott.pipelines.async_pipelines.run_tensor_sirt(data_container, maxiter=10, basis_set=None, update_frequency=5)[source]¶

An asynchronous implementation of the tensor SIRT algorithm. This approach uses only asynchronous, in-place operations on the GPU, and is therefore much faster than standard pipelines, as well as more memory-efficient.

Parameters

data_container (DataContainer) – The data.
maxiter (int) – Maximum number of iterations.
basis_set (Optional[BasisSet]) – User-provided basis set to be used, if desired. By default, a GaussianKernels basis set is used.
update_frequency (int) – Synchronization and norm reduction progress printing frequency. If set too small, the optimization will be slower due to frequent host-device synchronization. This effect can be seen by noting the iterations per second on the progress bar. The printed norm does not account for the total variation regularization.

Returns

Dictionary with reconstruction, projector, basis_set entries.

mumott.pipelines.async_pipelines.run_motr(data_container, maxiter=10, momentum=0.9, l1_weight=0.001, tv_weight=0.001, basis_set=None, lower_bound=None, step_size=1.0, update_frequency=5)[source]¶

MOTR (MOmentum Total variation Reconstruction) pipeline that uses asynchronous GPU (device) computations only during the reconstruction. This leads to a large speedup compared to standard pipelines which synchronize with the CPU several times per iteration. However, this implementation is less modular than standard pipelines.

This pipeline uses Nestorov Momentum in combination with total variation and L1 regularization, as well as the Tensor SIRT preconditioner-weight pair, which normalize the gradient based on the projection geometry and basis set.

This is also a relatively efficient implementation with respect to device memory, as all arithmetic operations are done in-place.

Parameters

data_container (DataContainer) – The container for the data.
maxiter (int) – Maximum number of iterations.
momentum (float) – Momentum for the Nestorov gradient
l1_weight (float) – Weight for L1 regularization.
tv_weight (float) – Weight for total variation regularization.
basis_set (Optional[BasisSet]) – User-provided basis set to be used, if desired. By default, a GaussianKernels basis set is used.
lower_bound (Optional[float]) – Lower bound to threshold coefficients at in each iteration. Can be used to e.g. enforce non-negativity for local basis sets such as GaussianKernels. Not used if set to None, which is the default settingz.
step_size (float) – Step size for each iteration in the reconstruction. Default value is 1. which should be a suitable value in the vast majority of cases, due to the normalizing effect of the weight-preconditioner pair used.
update_frequency (int) – Synchronization and norm reduction progress printing frequency. If set too small, the optimization will be slower due to frequent host-device synchronization. This effect can be seen by noting the iterations per second on the progress bar. The printed norm does not account for the total variation regularization.

Returns

Dictionary with reconstruction, projector, basis_set, loss_curve entries. loss_curve is an array with maxiter // update_frequency rows where the entries in each column are iteration, loss, residual_norm, tv_norm.

mumott.pipelines.async_pipelines.run_radtt(data_container, maxiter=10, momentum=0.9, step_size=1.0, delta=1.0, tv_weight=0.001, basis_set=None, lower_bound=0.0, update_frequency=5)[source]¶

RADTT (Robust And Denoised Tensor Tomography) pipeline that uses asynchronous GPU (device) computations only during the reconstruction. This leads to a large speedup compared to standard pipelines which synchronize with the CPU several times per iteration. However, this implementation is less modular than standard pipelines.

This pipeline uses Nestorov accelerated gradient descent, with a Huber loss function and total variation regularization. This reconstruction approach requires a large number of iterations, but is robust to outliers in the data.

This is also a relatively efficient implementation with respect to device memory, as all arithmetic operations are done in-place.

Parameters

data_container (DataContainer) – The container for the data.
maxiter (int) – Maximum number of iterations.
momentum (float) – Momentum for the Nestorov gradient
step_size (float) – Step size for L1 optimization. Step size for L2 part of optimization is step_size / (2 * delta). A good choice for the step size is typically around the same order of magnitude as each coefficient of the reconstruction.
delta (float) – Threshold for transition to L2 optimization. Should be small relative to data. Does not affect total variation.
tv_weight (float) – Weight for total variation regularization.
basis_set (Optional[BasisSet]) – User-provided basis set to be used, if desired. By default, a GaussianKernels basis set is used.
lower_bound (float) – Lower bound to threshold coefficients at in each iteration. Can be used to e.g. enforce non-negativity for local basis sets such as GaussianKernels. Not used if set to None, which is the default settings.
update_frequency (int) – Synchronization and norm reduction progress printing frequency. If set too small, the optimization will be slower due to frequent host-device synchronization. This effect can be seen by noting the iterations per second on the progress bar. The printed norm does not account for the total variation regularization.

Returns

Dictionary with reconstruction, projector, basis_set, loss_curve entries.

mumott.pipelines.sparse_pipelines.run_stsirt(data_container, maxiter=10, sparsity_count=3, basis_set=None, update_frequency=5)[source]¶

A sparse basis implementation of the tensor SIRT algorithm (Sparse Tensor SIRT). This approach uses only asynchronous, in-place operations on the GPU, and is therefore much faster than standard pipelines, as well as more memory-efficient. In exchange, it is less modular and requires a relatively sparse basis set.

Parameters

data_container (DataContainer) – The data.
maxiter (int) – Maximum number of iterations.
sparsity_count (int) – The maximum number of non-zero matrix elements per detector segment, if using the default BasisSet. Not used if a custom BasisSet is provided.
basis_set (Optional[BasisSet]) – A custom BasisSet instance. By default, a GaussianKernels with enforce_sparsity=True and sparsity_count=sparsity_count is used.
update_frequency (int) – Synchronization and norm reduction progress printing frequency. If set too small, the optimization will be slower due to frequent host-device synchronization. This effect can be seen by noting the iterations per second on the progress bar. The printed norm does not account for the total variation regularization.

Returns

Dictionary with reconstruction, projector, basis_set entries.

mumott.pipelines.sparse_pipelines.run_smotr(data_container, maxiter=10, sparsity_count=3, momentum=0.9, l1_weight=0.001, tv_weight=0.001, basis_set=None, update_frequency=5, lower_bound=None, step_size=1.0)[source]¶

SMOTR (Sparse MOmentum Total variation Reconstruction) pipeline that uses asynchronous GPU (device) computations only during the reconstruction. This leads to a large speedup compared to standard pipelines which synchronize with the CPU several times per iteration. However, this implementation is less modular than standard pipelines.

This pipeline uses Nestorov Momentum in combination with total variation and L1 regularization, as well as the Tensor SIRT preconditioner-weight pair.

This is also a highly efficient implementation with respect to device memory, as all arithmetic operations are done in-place.

Parameters

data_container (DataContainer) – The container for the data.
maxiter (int) – Maximum number of iterations.
sparsity_count (int) – The maximum number of non-zero matrix elements per detector segment.
momentum (float) – Momentum for the Nestorov gradient
l1_weight (float) – Weight for L1 regularization.
tv_weight (float) – Weight for total variation regularization.
basis_set (Optional[BasisSet]) – User-provided basis set to be used, if desired. By default, a GaussianKernels basis set is used.
update_frequency (int) – Synchronization and norm reduction progress printing frequency. If set too small, the optimization will be slower due to frequent host-device synchronization. This effect can be seen by noting the iterations per second on the progress bar. The printed norm does not account for the total variation regularization.
lower_bound (Optional[float]) – Lower bound to threshold coefficients at in each iteration. Can be used to e.g. enforce non-negativity for local basis sets such as GaussianKernels. Not used if set to None, which is the default settingz.
step_size (float) – Step size for each iteration in the reconstruction. Default value is 1. which should be a suitable value in the vast majority of cases, due to the normalizing effect of the weight-preconditioner pair used.

Returns

Dictionary with reconstruction, projector, basis_set, loss_curve entries. loss_curve is an array with maxiter // update_frequency rows where the entries in each column are iteration, loss, residual_norm, tv_norm.

mumott.pipelines.sparse_pipelines.run_spradtt(data_container, maxiter=10, sparsity_count=3, momentum=0.9, step_size=1.0, delta=1.0, tv_weight=0.001, basis_set=None, lower_bound=0.0, update_frequency=5)[source]¶

SPRADTT (SParse Robust And Denoised Tensor Tomography) pipeline that uses asynchronous GPU (device) computations only during the reconstruction. This leads to a large speedup compared to standard pipelines which synchronize with the CPU several times per iteration. However, this is only suitable for sparse representations and therefore this implementation is less modular than standard pipelines.

This pipeline uses Nestorov accelerated gradient descent, with a Huber loss function and total variation regularization. This reconstruction approach requires a large number of iterations, but is robust to outliers in the data.

This is also a highly efficient implementation with respect to device memory, as all arithmetic operations are done in-place.

Parameters

data_container (DataContainer) – The container for the data.
maxiter (int) – Maximum number of iterations.
sparsity_count (int) – The maximum number of non-zero matrix elements per detector segment.
momentum (float) – Momentum for the Nestorov gradient. Should be in the range [0., 1.].
step_size (float) – Step size for L1 optimization. Step size for L2 part of optimization is step_size / (2 * delta). A good choice for the step size is typically around the same order of magnitude as each coefficient of the reconstruction.
delta (float) – Threshold for transition to L2 optimization. Should be small relative to data. Does not affect total variation.
tv_weight (float) – Weight for total variation regularization.
basis_set (Optional[BasisSet]) – User-provided basis set to be used, if desired. By default, a GaussianKernels basis set is used.
lower_bound (float) – Lower bound to threshold coefficients at in each iteration. Can be used to e.g. enforce non-negativity for local basis sets such as GaussianKernels. Not used if set to None, which is the default settingz.
update_frequency (int) – Synchronization and norm reduction progress printing frequency. If set too small, the optimization will be slower due to frequent host-device synchronization. This effect can be seen by noting the iterations per second on the progress bar. The printed norm does not account for the total variation regularization.

Returns

Dictionary with reconstruction, projector, basis_set, loss_curve entries.

Alignment¶

This diagram shows an overview of how the alignment pipelines work. Dashed-outline steps are exclusive to the optical flow alignment.

mumott.pipelines.phase_matching_alignment.run_phase_matching_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.optical_flow_alignment.run_optical_flow_alignment(data_container, **kwargs)[source]¶

This pipeline implements the alignment algorithm described in [Odstrcil2019] The latter allows one to compute the shifts that are need to align the data according to the tomographical problem defined by the given geometry and projector. The alignment also relies on the related definition of the reconstruction volume and main axis of rotation (‘y’ axis). The procedure can be customized via the keyword arguments described under below (see Notes).

Parameters

data_container (DataContainer) – Input data.
volume (tuple[int, int, int]) – Geometry parameter: The size of the volume of the wanted tomogram. If not specified, deduced from information in data_container.
main_rot_axis (int) – Geometry parameter: The index of the main rotation axis (y axis on the real geometry) on the array. If not specified, deduced from information in data_container.
smooth_data (bool) – Data filtering: If True apply Gaussian filter to smoothen the data; default: False.
sigma_smooth (float) – Data filtering: The smoothing kernel related to the overall data smoothing; default: 0.
high_pass_filter (float) – Data filtering: The kernel for high pass filter in the gradient computation to avoid phase artifact; default: 0.01.
optimal_shift (np.ndarray[float]) – Shifts: The original shift; default: np.zeros((nb_projections, 2)).
rec_iteration (int) – Projector parameters: The number of iteration used to solve the tomogram from the projections; default: 20.
use_gpu (bool) – Projector parameters: Use GPU (True) or CPU (False) for the tomographic computation; default: False.
optimizer_kwargs (dict[str, Any]) – Projector parameters: Keyword arguments to pass on to the optimizer; default: dict(nestorov_weight = 0.6).
stop_max_iteration (int) – Alignment parameters: The maximum iteration allowed to find the shifts for the alignment; default: 20.
stop_min_correction (float) – Alignment parameters: The optimization is terminated if the correction (in pixel) drops below this value; default 0.01.
align_horizontal (bool) – Alignment parameters: Apply the horizontal alignment procedure; default: True.
align_vertical (bool) – Alignment parameters: Apply the vertical alignment procedure; default: True.
center_reconstruction (bool) – Alignment parameters: Shift the reconstructed tomogram to the center of the volume to avoid drift in the alignment procedure; default: True.

Return type

tuple[ndarray[float], ndarray[float], ndarray[float]]

Returns

A tuple comprising the shifts used for aligning the, the resulting (and aligned) projections obtained by projecting the reconstructed tomogram based on the aligned data, and the resulting tomogram reconstructed using the aligned data.

Example

The alignment procedure is simple to use.

>>> import numpy as np
>>> from mumott.data_handling import DataContainer
>>> from mumott.pipelines import optical_flow_alignment as ofa
>>> data_container = DataContainer('tests/test_fbp_data.h5')

We introduce some spurious offsets to this already-aligned data.

>>> length = len(data_container.geometry)
>>> data_container.geometry.j_offsets = np.arange(0., length) - length * 0.5
>>> data_container.geometry.k_offsets =         np.cos(np.arange(0., length) * np.pi / length)

We then perform the alignment with default parameters.

>>> shifts, sinogram_corr, tomogram_corr = ofa.run_optical_flow_alignment(data_container)
...

To use the alignment shifts, we have to translate them from the reconstruction (x, y, z)-coordinates into the projection (p, j, k) coordinates. The function shifts_to_vector() transforms the shifts in projection space into shifts in 3D space based on the detector geometry. The function shifts_to_geometry() transforms the shifts vector from the detector reference frame to the sample reference frame.

>>> j_offsets, k_offsets = ofa.shifts_to_geometry(data_container, shifts)
>>> data_container.geometry.j_offsets = j_offsets
>>> data_container.geometry.k_offsets = k_offsets

We can configure a variety of parameters and pass those to the alignment pipeline. For example, we can choose whether to align in the horizontal or vertical directions, whether to center the reconstruction, initial guesses for the shifts, the number of iterations for the reconstruction, the number of iterations for the alignment procedure, and so on.

A quick way to align very misaligned data is to use the line_vertical_alignment() function to obtain an initial guess. It can as well be used between alignment steps. Note that their is two version for the vertical It can also be used between alignment steps. There are two version for the vertical alignment. line_vertical_alignment() uses the DataContainer geometry and data, while sinogram_line_vertical_alignment() is based on a sinogram and the geometry. Note that the more strongly the object deviates from axial symmetry, the worse the vertical line alignment will be.

>>> initial_shift = ofa.line_vertical_alignment(data_container)

>>> alignment_param = dict(
...    optimal_shift=initial_shift,
...    rec_iteration=10,
...    stop_max_iteration=3,
...    align_horizontal=True,
...    align_vertical=True,
...    center_reconstruction=True,
...    optimizer_kwargs=dict(nestorov_weight=0.6))
>>> shifts, sinogram_corr, tomogram_corr = ofa.run_optical_flow_alignment(data_container,
...                                                                       **alignment_param)
>>> vertical_shift, index = ofa.sinogram_line_vertical_alignment(
...    sinogram_corr, ofa.get_alignment_geometry(data_container)[0])
>>> shifts[:,index] += vertical_shift # to update the shifts with vertical line alignment.

Utilities¶

mumott.pipelines.fbp_utilities.get_filtered_projections(projections, axis_string='inner', filter_type='Ram-Lak', normalization_percentile=99.9, transmittivity_cutoff=(None, None))[source]¶

Applies a high-pass filter to a selected subset of the absorbances for filtered back projection.

Parameters

projections (ProjectionStack) – The ProjectionStack to calculate the filtered projections from.
axis_string (str) – Default is 'inner', the value depends on how the sample is mounted to the holder. Typically, the inner axis is the rotation axis while the 'outer' axis refers to the tilt axis.
filter_string – Default is 'ram-lak', a high-pass filter. Other options are 'Hamming' and 'Hanning'.
normalization_percentile (float) – The normalization percentile to use for the transmittivity calculation. See get_transmittivities for details.
transmittivity_cutoff (tuple[float, float]) – The cutoffs to use for the transmittivity calculation. See get_transmittivities for details.

Return type

tuple[ndarray, ndarray]

Returns

A tuple containing the filtered subset of the absorbances and the index of the axis orthogonal to the inner or outer axis.

mumott.pipelines.utilities.alignment_geometry.get_alignment_geometry(data_container)[source]¶

Define the index of the main rotation axis relative to the data, and the related tomogram volume.

Parameters: data_container (DataContainer) – A :class:’DataContainer <mumott.data_handling.DataContainer>’ instance.
Return type: tuple[int, ndarray[int]]
Returns: a tuple comprising the index of the main rotation axis (0 or 1) and the volume of the tomogram related to the data and the main rotation axis

mumott.pipelines.utilities.alignment_geometry.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.utilities.alignment_geometry 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