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
) – TheDataContainer
holding the data set of interest.use_gpu (
bool
) – Whether to use GPU resources in computing the projections. Default isFalse
. If set toTrue
, the method will useSAXSProjectorCUDA
.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
- Returns
A dictionary with the entry
'result'
, with an entry'x'
, containing a filtered back projection reconstruction of the absorptivity calculated from thediode
.
- 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
) – TheDataContainer
from loading the data set of interest.use_absorbances (
bool
) – IfTrue
, the reconstruction will use the absorbances calculated from the diode, or absorbances provided via a keyword argument. IfFalse
, the data indata_container.data
will be used.use_sirt_weights (
bool
) – IfTrue
(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 isFalse
, which meansSAXSProjector
. If set toTrue
, the method will useSAXSProjectorCUDA
.maxiter (
int
) – Maximum number of iterations for the gradient descent solution.ftol (
Optional
[float
]) – Tolerance for the change in the loss function. Default isNone
, 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 toTrue
, 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 andGaussianKernels
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 byloss_function.add_regularizer()
.- Optimizer
The optimizer class to use. If not provided
GradientDescent
will be used. By default, the keyword argumentnestorov_weight
is set to0.95
, andenforce_non_negativity
isTrue
- 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
andweights_cutoff
keyword arguments.- Parameters
data_container (
DataContainer
) – TheDataContainer
from loading the data set of interest.use_absorbances (
bool
) – IfTrue
, the reconstruction will use the absorbances calculated from the diode, or absorbances provided via a keyword argument. IfFalse
, the data indata_container.data
will be used.use_gpu (
bool
) – Whether to use GPU resources in computing the projections. Default isFalse
. If set toTrue
, the method will useSAXSProjectorCUDA
.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 isFalse
.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 toTrue
, 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 andGaussianKernels
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
) – TheDataContainer
from loading the data set of interest.use_gpu (
bool
) – Whether to use GPU resources in computing the projections. Default isFalse
. If set toTrue
, the method will useSAXSProjectorCUDA
.maxiter (
int
) – Maximum number of iterations for the gradient descent solution.ftol (
float
) – Tolerance for the change in the loss function. Default isNone
, in which case the reconstruction will terminate once the maximum number of iterations have been performed.regularization_weight (
float
) – Regularization weight for the defaultLaplacian
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 byloss_function.add_regularizer()
. By default, aLaplacian
with the weightregularization_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
) – TheDataContainer
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 isFalse
. If set toTrue
, the method will useSAXSProjectorCUDA
.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
) – TheDataContainer
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 isFalse
. If set toTrue
, the method will useSAXSProjectorCUDA
.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, aGaussianKernels
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 gradientl1_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, aGaussianKernels
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 asGaussianKernels
. Not used if set toNone
, which is the default settingz.step_size (
float
) – Step size for each iteration in the reconstruction. Default value is1.
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 withmaxiter // update_frequency
rows where the entries in each column areiteration, 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 gradientstep_size (
float
) – Step size for L1 optimization. Step size for L2 part of optimization isstep_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, aGaussianKernels
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 asGaussianKernels
. Not used if set toNone
, 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 defaultBasisSet
. Not used if a customBasisSet
is provided.basis_set (
Optional
[BasisSet
]) – A customBasisSet
instance. By default, aGaussianKernels
withenforce_sparsity=True
andsparsity_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 gradientl1_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, aGaussianKernels
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 asGaussianKernels
. Not used if set toNone
, which is the default settingz.step_size (
float
) – Step size for each iteration in the reconstruction. Default value is1.
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 withmaxiter // update_frequency
rows where the entries in each column areiteration, 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 isstep_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, aGaussianKernels
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 asGaussianKernels
. Not used if set toNone
, 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 indata_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 theprojection
anddata
. 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
) – Acallable
, 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 theuse_absorbances
keyword argument, then anabsorbances
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 toreconstruction_pipeline
. If'data_container'
or'use_gpu'
are set as keys, they will overridedata_container
anduse_gpu
use_gpu (
bool
) – Whether to use GPU resources in computing the reconstruction. Default isFalse
. Will be overridden if set inreconstruction_pipeline_kwargs
.use_absorbances (
bool
) – Whether to use the absorbances to compute the reconstruction and align the projections. Default isTrue
. Will be overridden if set inreconstruction_pipeline_kwargs
.maxiter (
int
) – Maximum number of iterations for the alignment.upsampling (
int
) – Upsampling factor during alignment. If used, any masking indata_container.weights
will be ignored, butprojection_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 is1 / 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 is5 / 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 to1
, 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 thej
direction. Default isTrue
.align_k (
bool
) – Whether to align in thek
direction. Default isTrue
.
- Return type
- 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 indata_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
- 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 functionshifts_to_vector()
transforms the shifts in projection space into shifts in 3D space based on the detector geometry. The functionshifts_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 theDataContainer
geometry and data, whilesinogram_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
) – TheProjectionStack
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. Seeget_transmittivities
for details.transmittivity_cutoff (
tuple
[float
,float
]) – The cutoffs to use for the transmittivity calculation. Seeget_transmittivities
for details.
- Return type
- 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
- 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 thegeometry
based on a three-dimensional shift vector. Use this to reposition the reconstruction within the volume.- Parameters
- Return type
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