Methods

mumott distinguishes basis sets, projectors, and residual calculators, which can be combined to implement different methods.

• Basis sets
  • SphericalHarmonics
  • TrivialBasis
  • GaussianKernels
  • NearestNeighbor
• Projectors
  • SAXSProjector
  • SAXSProjectorCUDA
• Residual calculators
  • GradientResidualCalculator
  • ZHTTResidualCalculator

Utilities

mumott.methods.utilities.preconditioning.get_largest_eigenvalue(basis_set, projector, weights=None, preconditioner=None, niter=5, seed=None)

Calculate the largest eigenvalue of the matrix \(A^{T}*A\) using the power method. Here, \(A\) is the linear forward model defined be the input Projector and BasisSet. The largest eigenvalue can be used to set a safe step size for various optimizers.

Parameters

• basis_set (BasisSet) – basis set object

• projector (Projector) – projector object

• weights (Optional[ndarray[Any, dtype[float]]]) – Weights which will be applied to the residual. Default is None.

• preconditioner (Optional[ndarray[Any, dtype[float]]]) – A preconditioner which will be applied to the gradient. Default is None.

• niter (int) – number of iterations. Default is 5

• seed (Optional[int]) – Seed for random generation as starting state. Used for testing.

• Returns. –

• ------- – An estimate of the matrix norm (largest singular value)

Return type

float

mumott.methods.utilities.preconditioning.get_sirt_preconditioner(projector, cutoff=0.1)

Retrieves the SIRT preconditioner, which can be used together with the SIRT weights to condition the gradient of tomographic calculations for faster convergence and scaling of the step size for gradient descent.

Notes

The preconditioner normalizes the gradient according to the number of data points that map to each voxel in the computation of the projection adjoint. This preconditioner scales and conditions the gradient for better convergence. It is best combined with the SIRT weights, which normalize the residual for the number of voxels. When used together, they condition the reconstruction sufficiently well that a gradient descent optimizer with step size unity can arrive at a good solution. Other gradient-based solvers can also benefit from this preconditioning.

In addition, the calculation of these preconditioning terms makes it easy to identify regions of the volume or projections that are rarely probed, allowing them to be masked from the solution altogether.

If the projection operation is written in sparse matrix form as

\[P_{ij} X_{j} = Y_i\]

where \(P_{ij}\) is the projection matrix, \(X_j\) is a vector of voxels, and \(Y_i\) is the projection, then the preconditioner can be understood as

\[C_{jj} = \frac{I(n_j)}{\sum_i P_{ij}}\]

where \(I(n_j)\) is the identity matrix of the same size as \(X_j\). Similarly, the weights (of get_sirt_weights()) are computed as

\[W_{ii} = \frac{I(n_i)}{\sum_j P_{ij}}.\]

Here, any singularities in the system (e.g., where \(\sum_j P_{ij} = 0\)) can be masked out by setting the corresponding weight to zero. We thus end up with a weighted least-squares system

\[\text{argmin}_X(\Vert W_{ii}(P_{ij}X_{j} - D_{i})\Vert^2_2)\]

where \(D_{i}\) is some data, which we can solve iteratively by preconditioned gradient descent,

\[X_j^{k + 1} = X_j^k - C_{jj}P_{ji}^TW_{ii}(P_ij X_j^k - D_i)\]

As mentioned, we can add additional regularization terms, and because the preconditioning scales the problem appropriately, computing an optimal step size is not a requirement, although it can speed up the solution. This establishes a very flexible system, where quasi-Newton solvers such as LBFGS can be seamlessly combined with less restrictive gradient descent methods.

A good discussion of the algorithmic properties of SIRT can be found in this article by Gregor et al., while this article by van der Sluis et al. discusses SIRT as a least-squares solver in comparison to the conjugate gradient (CG) method.

Parameters

• projector (Projector) – A projector object which is used to calculate the weights. The computation of the weights is based on the geometry attached to the projector.

• cutoff (float) – The minimal number of rays that need to map to a voxel for it to be considered valid. Default is 0.1. Invalid voxels will be masked from the preconditioner.

Return type

ndarray[Any, dtype[float]]

mumott.methods.utilities.preconditioning.get_sirt_weights(projector, cutoff=0.1)

Retrieves the SIRT weights, which can be used together with the SIRT preconditioner to weight the residual of tomographic calculations for faster convergence and scaling of the step size for gradient descent.

Notes

See get_sirt_preconditioner() for theoretical details.

Parameters

• projector (Projector) – A projector object which is used to calculate the weights. The computation of the weights is based on the geometry attached to the projector.

• cutoff (float) – The minimal number of voxels that need to map to a point for it to be considered valid. Default is 0.1. Invalid pixels will be masked.

Return type

ndarray[Any, dtype[float]]

mumott.methods.utilities.preconditioning.get_tensor_sirt_preconditioner(projector, basis_set, cutoff=0.1)

Retrieves the SIRT preconditioner for tensor representations, which can be used together with the tensor SIRT weights to condition the gradient of tensor tomographic calculations for faster convergence and scaling of the step size for gradient descent.

Notes

The preconditioner is computed similarly as in the scalar case, except the underlying system of equations is

\[P_{ij} X_{jk} U_{kl} = Y_{il}\]

where \(X_{jk}\) is a vector of tensor-valued voxels, and \(Y_{il}\) is the projection into measurement space. The preconditioner then corresponds to

\[C_{jjkk} = \frac{I(n_j) \otimes I(n_k)}{\sum_{i,l} P_{ij} U_{kl}},\]

where \(I(n_j)\) is the identity matrix of the same size as \(X_j\). Similarly, the weights (of get_tensor_sirt_weights()) are computed as

\[W_{iill} = \frac{I(n_j) \otimes I(n_k)}{\sum_{j, k} P_{ij} U_{kl}}.\]

Here, any singularities in the system (e.g., where \(\sum_j P_{ij} = 0\)) can be masked out by setting the corresponding weight to zero. We thus end up with a weighted least-squares system

\[\text{argmin}_X(\Vert W_{iill}(P_{ij} X_{jk} U_{kl} - D_{il})\Vert^2_2),\]

where \(D_{il}\) is some data. This problem can be solved iteratively by preconditioned gradient descent,

\[X_j^{k + 1} = X_j^k - C_{jjkk}P_{ji}^TW_{iill}(P_ij X_{jk}^k U_{kl} - D_{il}).\]

Parameters

• projector (Projector) – A projector object which is used to calculate the weights. The computation of the weights is based on the geometry attached to the projector.

• basis_set (BasisSet) – A basis set object which is used to calculate the weights. The tensor-space-to-detector-space transform uses the basis set. Should use a local representation.

• cutoff (float) – The minimal number of rays that need to map to a voxel for it to be considered valid. Default is 0.1. Invalid voxels will be masked from the preconditioner.

Return type

ndarray[Any, dtype[float]]

mumott.methods.utilities.preconditioning.get_tensor_sirt_weights(projector, basis_set, cutoff=0.1)

Retrieves the tensor SIRT weights, which can be used together with the SIRT preconditioner to weight the residual of tensor tomographic calculations for faster convergence and scaling of the step size for gradient descent.

Notes

See get_tensor_sirt_preconditioner() for theoretical details.

Parameters

• projector (Projector) – A projector object which is used to calculate the weights. The computation of the weights is based on the geometry attached to the projector.

• basis_set (BasisSet) – A basis set object which is used to calculate the weights. The tensor-space-to-detector-space transform uses the basis set. Should use a local representation.

• cutoff (float) – The minimal number of voxels that need to map to a point for it to be considered valid. Default is 0.1. Invalid pixels will be masked.

Return type

ndarray[Any, dtype[float]]

mumott.methods.utilities.fiber_fit.degree_of_symmetry_map(coefficients, ell_max, resolution=10, filter=None)

Make a longitude-latitude map of the degree of symmetry. Can be used to make illustrating plots and to decide which filter is more appropriate.

Parameters

• coefficients (ndarray[Any, dtype[float]]) – Spherical harmonic coefficients of a single voxel.

• ell_max (int) – Largest order \(\ell\) used in the fitting. If coefficients contains higher orders, it will be truncated.

• resolution (int) – Number or angular steps along a half-circle used in the search for the optimal axis.

• filter (Optional[str]) – Weighing of different orders used to calculate the degree of symmetry. Possible values are “ramp” and “square”. By default no filter is applied (None).

Return type

Tuple[ndarray[Any, dtype[float]]]

Returns

• dos – Map of the calculated degree of symmetry.

• theta – Polar angle coordinates of the map.

• phi – Azimuthal angle coordinates of the map.

mumott.methods.utilities.fiber_fit.find_approximate_symmetry_axis(coefficients, ell_max, resolution=10, filter=None)

Find the axis of highest apparent symmetry voxel-by-voxel for a voxel map of sperical harmonics coefficients. As a default, the measure of degree of symmetry is a power of the function rotationally averaged around the given axis.

Parameters

• coefficients (ndarray[Any, dtype[float]]) – Voxel map of spherical harmonic coefficients.

• ell_max (int) – Largest order \(\ell\) used in the fitting. If coefficients contains higher orders, it will be truncated.

• resolution (int) – Number or angular steps along a half-circle, used in the search for the optimal axis.

• filter (Optional[str]) – Weighing of different orders used to calculate the degree of symmetry. Possible values are “ramp” and “square”. By default no filter is applied (None).

Return type

Tuple[ndarray[Any, dtype[float]]]

Returns

• optimal_zonal_coeffs – Zonal harmonics coefficients in the frame-of-reference corresponding to the axis of highest degree of symmetry.

• optimal_theta – Voxel-by-voxel polar angles of the axis with highest degree of symmetry.

• optimal_phi – Voxel-by-voxel azimuthal angles of the axis with highest degree of symmetry.

mumott.methods.utilities.fiber_fit.symmetric_part_along_given_direction(coefficients, theta, phi, ell_max)

Find the zonal harmonic coefficients along the given input directions. This can be used if eigenvector analysis is used to generate the symmetry directions used for a further zonal-harmonics refinement step.

Parameters

• coefficients (ndarray[Any, dtype[float]]) – Voxel map of spherical harmonic coefficients.

• theta (ndarray[Any, dtype[float]]) – Voxel map of polar angles.

• phi (ndarray[Any, dtype[float]]) – Voxel map of azimuthal angles.

• ell_max (int) – Largest order \(\ell\) used in the fitting. If coefficients contains higher orders, it will be truncated.

Return type

ndarray[Any, dtype[float]]

Returns

Zonal harmonics coefficients in the frame of reference corresponding to the axis of highest degree of symmetry.

Utilities for efficient tensor algebra.

mumott.methods.utilities.tensor_operations.framewise_contraction(matrices, data, out)

This function is a more efficient implementation of the series of matrix multiplications which is carried out the expression np.einsum('ijk, inmj -> inmk', matrices, data, out=out), using numba.jit, for a four-dimensional data array. Used to linearly transform data from one representation to another.

Parameters

• matrices (ndarray[float]) – Stack of matrices with shape (i, j, k) where i is the stacking dimension.

• data (ndarray[float]) – Stack of data with shape (i, n, m, ..., j) where i is the stacking dimension and j is the representation dimension.

• out (ndarray[float]) – Output array with shape (i, n, m, ..., k) where i is the stacking dimension and k is the representation dimension. Will be modified in-place and must be contiguous.

Return type

None

Notes

All three arrays will be cast to out.dtype.

mumott.methods.utilities.tensor_operations.framewise_contraction_transpose(matrices, data, out)

This function is a more efficient implementation of the series of matrix multiplications which is carried out the expression np.einsum('ijk, inmk -> inmj', matrices, data, out=out), using numba.jit, for a four-dimensional data array. Used to linearly transform data from one representation to another.

Parameters

• matrices (ndarray[float]) – Stack of matrices with shape (i, j, k) where i is the stacking dimension.

• data (ndarray[float]) – Stack of data with shape (i, n, m, ..., k) where i is the stacking dimension and k is the representation dimension.

• out (ndarray[float]) – Output array with shape (i, n, m, ..., j) where i is the stacking dimension and j is the representation dimension. Will be modified in-place and must be contiguous.

Return type

ndarray[float]

Notes

All three arrays will be cast to out.dtype.