Core¶

class mumott.core.geometry.Geometry(filename=None)[source]¶

Stores information about the system geometry. Instances of this class are used by DataContainer and ProjectionStack to maintain geometry information. They can be stored as a file using the write() method. This allows one to (re)create a Geometry instance from an earlier and overwrite the geometry information read by a DataContainer instance. This is useful, for example, in the context of alignment.

Parameters: filename (Optional[str]) – Name of file from which to read geometry information. Defaults to None, in which case the instance is created with default parameters.

append(value)[source]¶

Appends projection-wise geometry data provided as a GeometryTuple.

Return type: None

delete_projections()[source]¶

Delete all projections.

Return type: None

property detector_angles: NDArray(float)¶: Azimuthal angles of detector segments in radians. One-dimensional sequence of center positions

property detector_direction_origin: ndarray[float]¶: The direction at which the angle on the detector is zero, when no experimental rotation is applied.

property detector_direction_positive_90: ndarray[float]¶: Rotation matrices corresponding to each projection of data.

property full_circle_covered: bool¶: Whether the azimuthal bins cover a half-circle of the detector (False) or the full circle (True).

property hash_inner_angles: str¶: A blake2b hash of inner_angle.

property hash_inner_axes: str¶: A blake2b hash of inner_axes.

property hash_j_offsets: str¶: A blake2b hash of j_offsets_as_array.

property hash_k_offsets: str¶: A blake2b hash of k_offsets_as_array.

property hash_outer_angles: str¶: A blake2b hash of outer_anglesy.

property hash_outer_axes: str¶: A blake2b hash of outer_axes.

property hash_rotations: str¶: A blake2b hash of rotations_as_array.

property inner_angles: list[float]¶: Rotation angles for inner rotation, in radians.

property inner_angles_as_array: ndarray[float]¶: Rotation angles for inner rotations, in radians, as an array.

property inner_axes: list[numpy.ndarray[float]]¶: Inner rotation axes. All axes can be set at once using a single array with three entries.

property inner_axes_as_array: ndarray[float]¶: Inner rotation axes as an array.

insert(key, value)[source]¶

Inserts projection-wise data handed via a GeometryTuple.

Return type: None

property j_direction_0: ndarray[float]¶: The direction corresponding to the first index in each projection when no experimental rotation is applied.

property j_offsets: list[float]¶: Offsets to align projection in the direction j.

property j_offsets_as_array: ndarray[float]¶: Offsets to align projection in the direction j as an array.

property k_direction_0: ndarray[float]¶: The direction corresponding to the second index in each projection when no experimental rotation is applied.

property k_offsets: list[float]¶: Offsets to align projection in the direction k.

property k_offsets_as_array: ndarray[float]¶: Offsets to align projection in the direction k as an array.

property outer_angles: list[float]¶: Rotation angles for inner rotation, in radians.

property outer_angles_as_array: ndarray[float]¶: Rotation angles for inner rotations, in radians, as an array.

property outer_axes: list[numpy.ndarray[float]]¶: Inner rotation axes. All axes can be set at once using a single array with three entries.

property outer_axes_as_array: ndarray[float]¶: Outer rotation axes as an array.

property p_direction_0: ndarray[float]¶: The projection direction when no experimental rotation is applied.

property probed_coordinates: ProbedCoordinates¶

An array of 3-vectors with the (x, y, z)-coordinates on the reciprocal space map probed by the method. Structured as (N, K, 3, 3), where N is the number of projections, K is the number of detector segments, the second-to-last axis contains start-, mid-, and endpoints, and the last axis contains the (x, y, z)-coordinates.

Notes

The number of detector segments is len(geometry.detecor_angles)*len(geometry.two_theta) i.e. the product of the number of two_theta bins times the number of azimuthal bins. As a default, only on two theta bin is used. When several two_theta bins are used, the second index corresponds to a raveled array, where the azimuthal is the fast index and two theta is the slow index.

property projection_shape: ndarray[Any, dtype[int]]¶: 2D shape of the raster-scan. 1st element is the number of steps in the j-direction and the second is the number of steps in the k-direction.

read(filename)[source]¶

Method for reading the current state of a Geometry instance from file.

Parameters: filename (str) – Name of input file.
Return type: None

property reconstruction_rotations: list[numpy.ndarray[float]]¶: list of rotation matrices sequentially applied to the reconstruction geometry of the system.

rotate_reconstruction(A=None, axis=None, angle=None)[source]¶

Rotates the reconstruction geometry. The given rotation matrix will modify the rotation matrix of each projection by multiplication from the right, such that

\[R_i' = R_i A\]

where \(R_i\) is the rotation matrix of projection \(i\) and \(A\) is the rotation matrix. For each projection, the system vectors are then rotated by

\[v_i = (R_i')^T v = A^T R_i^T v\]

where \(v\) corresponds to e.g., p_direction_0.

Notes

It is not possible to directly modify rotations after adding a reconstruction rotation.

Parameters

A (Optional[ndarray[float]]) – A 3-by-3 rotation matrix. If not given, then axis and angle must be provided.
axis (Optional[ndarray[float]]) – An axis, given as a unit length 3-vector, about which the rotation is defined. Not used if A is provided.
angle (Optional[ndarray[float]]) – The angle in radians of the rotation about axis. Not used if A is provided.

rotate_system_vectors(A=None, axis=None, angle=None)[source]¶

Rotates the system vectors. The given rotation matrix will modify the system vectors by

\[v' = A v\]

where \(v\) is a system vector, e.g., p_direction_0, and \(A\) is the rotation matrix. For each projection, the system vectors are then rotated by

\[v_i = R_i^T A v\]

where \(R_i\) corresponds to rotations for projection \(i\).

Notes

It is not possible to directly modify the system vectors after adding a system rotation.

Parameters

A (Optional[ndarray[float]]) – A 3-by-3 rotation matrix. If not given, then axis and angle must be provided.
axis (Optional[ndarray[float]]) – An axis, given as a 3-vector, about which a rotation can be defined. Not used if A is provided.
angle (Optional[ndarray[float]]) – The angle in radians of the rotation about axis. Not used if A is provided.

property rotations: list[numpy.ndarray[float]]¶: Rotation matrices for the experimental rotation corresponding to each projection of data.

property rotations_as_array: ndarray[float]¶: Rotation matrices corresponding to each projection of data as an array.

property system_rotations: list[numpy.ndarray[float]]¶: list of rotation matrices sequentially applied to the basis vectors of the system.

property two_theta: NDArray(float)¶: Scattering angle in radians. Can be list of angles, if multiple radial bins are used.

property volume_shape: ndarray[Any, dtype[int]]¶: 3D shape of the reconstruction voxel array. 1st element is the number of points along the x-direction. 2nd is y and 3rd is z.

write(filename)[source]¶

Method for writing the current state of a Geometry instance to file.

Notes

Any rotations in reconstruction_rotations and system_rotations will be applied to the rotations and system vectors respectively prior to writing.

Parameters: filename (str) – Name of output file.
Return type: None

class mumott.core.geometry.GeometryTuple(rotation: np.ndarray[float] = array([[1., 0., 0.], [0., 1., 0.], [0., 0., 1.]]), j_offset: float = 0.0, k_offset: float = 0.0, inner_angle: float = None, outer_angle: float = None, inner_axis: np.ndarray[float] = None, outer_axis: np.ndarray[float] = None)[source]¶

Tuple for passing and returning projection-wise geometry information. This is a helper class used by Geometry.

rotation¶: Rotation matrix. If the angle and axis arguments are given, this matrix should correspond to the matrix given by R_outer @ R_inner, where R_inner is defined by a rotation by inner_angle about inner_axis, and similarly for R_outer.

j_offset¶: Offset to align projection in the direction j.

k_offset¶: Offset to align projection in the direction k.

inner_angle¶: Angle of rotation about inner_axis in radians.

outer_angle¶: Angle of rotation about outer_axis in radians.

inner_axis¶: Inner rotation axis.

outer_axis¶: Outer rotation axis.

class mumott.core.projection_stack.ProjectionStack[source]¶

Instances of this class contain data, geometry and other pertinent information for a series of measurements. The individual measurements are stored as Projection objects. The latter are accessible via list-like operations, which enables, for example, iteration over measurements but also retrieval of individual measurements by index, in-place modification or deletion.

The geometry information (i.e., rotations and offsets for each projection) is accessible via the geometry attribute. Data, diode readouts, and weights can be retrieved as contiguous arrays via the properties data, diode, and weights, respectively.

Example

The following code snippet illustrates how individual measurements can be accessed via list operations. For demonstration, here, we use default (“empty”) projections. In practice the individual measurements are read from a data file via the DataContainer class, which makes them available via the DataContainer.projections attribute.

First we create an empty projection stack.

>>> from mumott.core.projection_stack import Projection, ProjectionStack
>>> projection_stack = ProjectionStack()

Next we create a projection and attach it to the projection stack. In order to be able to distinguish this projection during this example, we assign it a Projection.j_offset of 0.5.

>>> projection = Projection(j_offset=0.5)
>>> projection_stack.append(projection)

The geometry information can now be accessed via the projection stack in several different but equivalent ways, including via the original projection object,

>>> print(projection.j_offset)

via indexing projection_stack

>>> print(projection_stack[0].geometry.j_offset)

or by indexing the respective geometry property of the projection stack itself.

>>> print(projection_stack.geometry.j_offsets[0])

We can modify the geometry parameters via any of these properties with identical outcome. For example,

>>> projection_stack[0].j_offset = -0.2
>>> print(projection.j_offset,
          projection_stack[0].geometry.j_offset,
          projection_stack.geometry.j_offsets[0])
-0.2 -0.2 -0.2

Next consider a situation where several projections are included in the projection stack.

>>> projection_stack.append(Projection(j_offset=0.1))
>>> projection_stack.append(Projection(j_offset=-0.34))
>>> projection_stack.append(Projection(j_offset=0.23))
>>> projection_stack.append(Projection(j_offset=0.78))
>>> print(projection_stack.geometry.j_offsets)
[-0.2, 0.1, -0.34, 0.23, 0.78]

The summary of the projection stack includes hashes for the data, the diode readout, and the weights. This allows one to get a quick indication for whether the content of these fields has changed.

>>> print(projection_stack)
--------------------------------------------------------------------------
                                  ProjectionStack
--------------------------------------------------------------------------
hash_data          : ...

We could, for example, decide to remove an individual projection as we might have realized that the data from that measurement was corrupted.

>>> del projection_stack[1]
>>> print(projection_stack)
--------------------------------------------------------------------------
                                  ProjectionStack
--------------------------------------------------------------------------
hash_data          : ...

From the output it is readily apparent that the content of the data field has changed as a result of this operation.

Finally, note that we can also loop over the projection stack, for example, to print the projections.

>>> for projection in projection_stack:
>>>     print(projection)
...

append(projection)[source]¶

Appends a measurement in the form of a Projection object. Once a projection is attached to a projection_stack, the geometry information of the projection will be synchronized with the geometry information of the projection_stack (see geometry).

Parameters: projection (Projection) – Projection object to be appended.
Return type: None

property data: ndarray[Any, dtype[_ScalarType_co]]¶: Scattering data, structured (n, j, k, w), where n is the projection number, j is the pixel in the j-direction, k is the pixel in the k-direction, and w is the detector segment. Before the reconstruction, this should be normalized by the diode. This may already have been done prior to loading the data.

property diode: ndarray[Any, dtype[_ScalarType_co]]¶: The diode readout, used to normalize the data. Can be blank if data is already normalized. The diode value should not be normalized per projection, i.e., it is distinct from the transmission value used in standard tomography.

property geometry: Geometry¶: Contains geometry information for each projection as well as information about the geometry of the whole system.

property hash_data: str¶: A hash of data.

property hash_diode: str¶: A sha1 hash of diode.

property hash_weights: str¶: A sha1 hash of weights.

index_by_key(key)[source]¶: Returns an index from a key.

insert(k, projection)[source]¶

Inserts a projection at a particular index, increasing the indices of all subsequent projections by 1.

Return type: None

property weights: ndarray[Any, dtype[_ScalarType_co]]¶: Weights applied multiplicatively during optimization. A value of 0 means mask, a value of 1 means no weighting, and other values means weighting each data point either less (weights < 1) or more (weights > 1) than a weight of 1.

class mumott.core.projection_stack.Projection(data=None, diode=None, weights=None, rotation=array([[1., 0., 0.], [0., 1., 0.], [0., 0., 1.]]), j_offset=0.0, k_offset=0.0, inner_angle=None, outer_angle=None, inner_axis=None, outer_axis=None)[source]¶

Instances of this class contain data and metadata from a single measurement. Typically they are appended to a ProjectionStack object.

Parameters

data (Optional[ndarray[Any, dtype[float]]]) – Data from measurement, structured into 3 dimensions representing the two scanning directions and the detector angle.
diode (Optional[ndarray[Any, dtype[float]]]) – Diode or transmission data from measurement, structured into 2 dimensions representing the two scanning directions.
weights (Optional[ndarray[Any, dtype[float]]]) – Weights or masking information, represented as a number between 0 and 1. 0 means mask, 1 means do not mask. Structured the same way as data.
rotation (ndarray[Any, dtype[float]]) – 3-by-3 rotation matrix, representing the rotation of the sample in the laboratory coordinate system.
j_offset (float) – The offset needed to align the projection in the j-direction.
k_offset (float) – The offset needed to align the projection in the k-direction.
inner_angle (Optional[float]) – Angle of rotation about inner_axis in radians.
outer_angle (Optional[float]) – Angle of rotation about outer_axis in radians.
inner_axis (Optional[ndarray[Any, dtype[float]]]) – Inner rotation axis.
outer_axis (Optional[ndarray[Any, dtype[float]]]) – Outer rotation axis.

attach_to_stack(projection_stack, index)[source]¶: Used to attach the projection to a projection_stack. This method should not be called by users.

property attached¶: Returns true if projection is attached to a ProjectionStack object.

property data: ndarray[Any, dtype[_ScalarType_co]]¶: Scattering data, structured (j, k, w), where j is the pixel in the j-direction, k is the pixel in the k-direction, and w is the detector segment. Before the reconstruction, the data should be normalized by the diode. This may already have been done prior to loading the data.

detach_from_stack()[source]¶: Used to detach the projection from a projection stack. This method should not be called by users.

property diode: ndarray[Any, dtype[float64]]¶: The diode readout, used to normalize the data. Can be blank if data is already normalized.

property geometry: GeometryTuple¶: Returns geometry information as a named tuple.

property hash_data: str¶: A hash of data.

property hash_diode: str¶: A sha1 hash of diode.

property hash_weights: str¶: A sha1 hash of weights.

property inner_angle: float¶: Rotation angle about inner axis.

property inner_axis: ndarray[Any, dtype[float]]¶: Rotation angle about inner axis.

property j_offset: float64¶: The offset needed to align the projection in the j-direction.

property k_offset: float64¶: The offset needed to align the projection in the k-direction.

property outer_angle: float¶: Rotation angle about inner axis.

property outer_axis: ndarray[Any, dtype[float]]¶: Rotation angle about inner axis.

property rotation: ndarray[Any, dtype[float64]]¶: 3-by-3 rotation matrix, representing the rotation of the sample in the laboratory coordinate system.

property weights: ndarray[Any, dtype[_ScalarType_co]]¶: Weights to be applied multiplicatively during optimization. A value of 0 means mask, a value of 1 means no weighting, and other values means weighting each data point either less (weights < 1) or more (weights > 1) than a weight of 1.

class mumott.core.spherical_harmonic_mapper.SphericalHarmonicMapper(ell_max=2, polar_resolution=16, azimuthal_resolution=32, polar_zero=0, azimuthal_zero=0, enforce_friedel_symmetry=True)[source]¶

Helper class for visualizing and analyzing spherical harmonics. Using this class, one can obtain the amplitudes for a given set of even-ordered harmonics and apply the Funk-Radon transform. These can then be plotted, analyzed, and so forth. In addition, the class allows one to represent functions in terms of spherical harmonics, if they are given as functions of the azimuthal and polar angles of the class instance.

Parameters

ell_max (int, optional) – Maximum order of the spherical harmonics. Default is 2.
polar_resolution (int, optional) – Number of samples in the polar direction. Default is 16.
azimuthal_resolution (int, optional) – Number of samples in the azimuthal direction. Default is 32.
polar_zero (float, optional) – The polar angle of the spherical harmonic coordinate system’s pole, relative to the reference coordinate system. Default is 0.
azimuthal_zero (float, optional) – The azimuthal angle of the spherical harmonic coordinate system’s pole, relative to a reference coordinate system. Default is 0.
enforce_friedel_symmetry (bool) – If set to True, Friedel symmetry will be enforced, using the assumption that points on opposite sides of the sphere are equivalent.

Example

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> S = SphericalHarmonicMapper(ell_max=4, polar_resolution=16, azimuthal_resolution=8)
>>> S.ell_indices
array([0, 2, 2, 2, 2, 2, 4, 4, ...])
>>> S.emm_indices
array([ 0, -2, -1,  0,  1,  2, -4, -3, ...])
>>> a = S.get_amplitudes(np.random.rand(S.ell_indices.size) - 0.5)
>>> plt.pcolormesh(S.phi * np.sqrt(1 / 2), # basic cylindrical equal-area projection
                   np.cos(S.theta) / np.sqrt(1 / 2),
                   a[:-1, :-1])
...

property azimuthal_zero: ndarray[Any, dtype[float]]¶: The azimuthal angle of the spherical harmonic pole, relative to a fixed reference system.

property coordinates: Tuple[ndarray[Any, dtype[float]], ndarray[Any, dtype[float]], ndarray[Any, dtype[float]]]¶: The X, Y, Z coordinates that the amplitudes are mapped to.

property ell_indices: ndarray[Any, dtype[int]]¶: The orders of the harmonics calculated by the class instance.

property ell_max: int¶: Maximum order of the spherical harmonics.

property emm_indices: ndarray[Any, dtype[int]]¶: The degrees of the harmonics calculated by the class instance.

get_amplitudes(coefficients, apply_funk_transform=False)[source]¶

Returns the amplitudes of a set of spherical harmonics. For sorting of the coefficients, see the orders and degrees attributes.

Parameters

coefficients (ndarray[Any, dtype[float]]) – The coefficients for which the amplitudes are to be calculated. Size must be equal to (ell_max + 1) * (ell_max / 2 + 1).
apply_funk_fransform – Whether to apply the Funk-Radon transform to the coefficients. This is useful for orientation analysis in some cases. Default is False.

Return type

ndarray[Any, dtype[float]]

get_harmonic_coefficients(amplitude)[source]¶

Returns the spherical harmonic coefficients for the given amplitudes. Can be used with amplitudes from another instance of SphericalHarmonicParameters with a different orientation in order to solve for the rotated coefficients. The accuracy of the representation depends on the maximum order and the polar and azimuthal resolution.

Parameters: amplitude (ndarray[Any, dtype[float]]) – The amplitude of the spherical function to be represented, as a function of theta and phi.
Return type: ndarray[Any, dtype[float]]

property map: ndarray[Any, dtype[float]]¶: Map between amplitude and harmonics.

property phi: ndarray[Any, dtype[float]]¶: The azimuthal angle to which the amplitude is mapped.

property polar_zero: ndarray[Any, dtype[float]]¶: The polar angle of the spherical harmonic pole, relative to a fixed reference system.

property theta: ndarray[Any, dtype[float]]¶: The polar angle to which the amplitude is mapped.

property unit_vectors¶: Probed coordinates object used by BasisSets to calculate function maps.

update_zeros(polar_zero, azimuthal_zero)[source]¶

Changes the orientation of the coordinate system.

Parameters

polar_zero (float) – The new polar angle at which the pole should be, relative to a fixed reference system.
azimuthal_zero (float) – The new azimuthal angle at which the pole should be, relative to a fixed reference system.

Return type

None

class mumott.core.probed_coordinates.ProbedCoordinates(vector=<factory>, great_circle_offset=<factory>)[source]¶

A small container class for probed coordinates for the :class:` BasisSet <mumott.methods.BasisSet>`.

Parameters

vector (NDArray[float]) – The coordinates on the sphere probed at each detector segment by the experimental method. Should be structured (N, M, I, 3) where N is the number of projections, M is the number of detector segments, I is the number of points on the detector to be integrated over, and the last index gives the (x, y, z) components of the coordinates. I can be as small as 1, but the array still needs to be structured in this way explicitly. By default, the value will be np.array([1., 0., 0.]).reshape(1, 1, 1, 3).
great_circle_offset (np.ndarray[float]) – The vector which offsets the probed coordinate vectors from lying on a great circle. Must have the same number of dimensions as vector, and be broadcastable (dimensions where the value would be repeated may be 1). This vector is used when interpolating coordinates that lie on a small circle.

property to_spherical: tuple¶: Returns spherical coordinates of vector, in the order (radius, polar angle, azimuthal angle).

mumott.core.hashing.list_to_hash(list_to_hash, hashing_algorithm='blake2b')[source]¶

Function which takes a list containing a set of objects and automatically generates a deterministic hash for them.

Parameters

list_to_hash (list) – List of a set of objects of various types, see notes for a complete list.
hashing_algorithm (str) – The hashing algorithm to use. Can be any algorithm name in hashlib.algorithms_available. Default is 'blake2b'.

Return type

str

Example

The following code snippets illustrate hashing lists that will work, and ones that will not work.

Works: A list of an integer, an array, a dictionary with valid types, and a None.

>>> from mumott.core.hashing import list_to_hash
>>> print(list_to_hash([1, np.array((1, 3, 5)), dict(val=1, string='abc'), None]))
2a949c...

Does not work: an array containing a None, due to the dtype being object.

>>> print(list_to_hash([np.array([None])]))
Traceback (most recent call last):
...
TypeError: Hash of dtype `object` is not deterministic, cannot hash [None]

Does not work: a generator expression, which is an unknown object.

>>> print(list_to_hash([(a for a in [1, 2, 3])]))
Traceback (most recent call last):
...
TypeError: Cannot hash unknown object: <generator object...

Notes

float-type objects are rounded to five significant digits in the mantissa before hashing. This is necessary to obtain semi-deterministic hashes that obey a subset of fuzzy equality for float comparison. There are edge cases where equality can fail due to rounding errors, but these should be extremely rare.

Supported entry types in list_to_hash:

int: Cast to string. Works along with similar numpy types.
float: Mantissa rounded to five significant digits and concatenated with exponent. Works along with similar numpy types.
complex: Real and imaginary parts concatenated and treated like float. Works along with similar numpy types.
str: Automatically given 'utf-8' encoding.
bytes: Cast to string.
None: Ignored.
np.ndarray: Provided dtype is not object, hence arrays of None are not allowed.
list, tuple: Provided they can be cast to allowed, i.e. non-ragged np.ndarray
dict: Assuming entries are allowed types. Keys and entries are concatenated. If an entry is None, the key is added to the hash while the entry is ignored.

mumott.core.john_transform.john_transform(field, projections, unit_vector_p, unit_vector_j, unit_vector_k, offsets_j, offsets_k, float_type='float64')[source]¶

Frontend for performing the John transform with parallel CPU computing, using an algorithm akin to mumott.core.john_transform_cuda().

Parameters

field (ndarray[Any, dtype[float]]) – The field to be projected, with 4 dimensions. The last index should have the same size as the last index of projections.
projections (ndarray[Any, dtype[float]]) – A 4-dimensional numpy array where the projections are stored. The first index runs over the different projection directions.
unit_vector_p (ndarray[Any, dtype[float]]) – The direction of projection in Cartesian coordinates.
unit_vector_j (ndarray[Any, dtype[float]]) – One of the directions for the pixels of projections`.
unit_vector_k (ndarray[Any, dtype[float]]) – The other direction for the pixels of projections.
offsets_j (ndarray[Any, dtype[float]]) – Offsets which align projections in the direction of j
offsets_k (ndarray[Any, dtype[float]]) – Offsets which align projections in the direction of k.
float_type (str) – Whether to use ‘float64’ (default) or ‘float32’. The argument should be supplied as a string. The types of field and projections must match this type.

Return type

callable

Notes

The computation performed by this function may be written as

\[p(I, J, K)_i = \sum_{s=0}^{N-1} d \cdot \sum_{t=0}^3 t_w V_i(\lfloor \mathbf{r}_j + d s \cdot \mathbf{v} \rfloor + \mathbf{t})\]

where \(p(I, J, K)_i\) is projection[I, J, K, i], \(V_i\) is volume[..., i], and \(\mathbf{v}\) is unit_vector_p[I]. \(N\) is the number of voxels in the maximal direction of unit_vector_p[I]. \(d\) is the step length, and is given by the norm of \(\Vert \mathbf{v}_I \Vert / \vert \max(\mathbf{v})) \vert\). \(t_w\) and \(\mathbf{t}\) are weights and index shifts, respectively. The latter are necessary to perform bilinear interpolation between the four closest voxels in the two directions orthogonal to the maximal direction of \(\mathbf{v}\). \(\mathbf{r}_j\) is the starting position for the ray, and is given by

\[\begin{split}\mathbf{r}_j(J, K) = (J + 0.5 + o_J - 0.5 \cdot J_\text{max}) \cdot \mathbf{u} + \\ (K + 0.5 + o_K - 0.5 \cdot K_\text{max}) \cdot \mathbf{w} + \\ (\Delta_p) \cdot \mathbf{v} + 0.5 \cdot \mathbf{r}_\text{max} \\\end{split}\]

where \(o_J\) is offsets_j[I], \(o_K\) is offsets_k[I], \(\mathbf{u}\) is unit_vector_j[I], \(\mathbf{w}\) is unit_vector_k[I], \(J_\text{max}\) and \(K_\text{max}\) are projections.shape[1] and projections.shape[2] respectively, and \(\mathbf{r}_\text{max}\) is volume.shape[:3]. \(\Delta_p\) is an additional offset that places the starting position at the edge of the volume.

mumott.core.john_transform.john_transform_adjoint(field, projections, unit_vector_p, unit_vector_j, unit_vector_k, offsets_j, offsets_k, float_type='float64')[source]¶

Frontend for performing the adjoint of the John transform with parallel CPU computing, using an algorithm akin to mumott.core.john_transform_cuda().

Parameters

field (ndarray[Any, dtype[float]]) – The field into which the adjoint is projected, with 4 dimensions. The last index should have the same size as the last index of projections.
projections (ndarray[Any, dtype[float]]) – The projections from which the adjoint is calculated. The first index runs over the different projection directions.
unit_vector_p (ndarray[Any, dtype[float]]) – The direction of projection in Cartesian coordinates.
unit_vector_j (ndarray[Any, dtype[float]]) – One of the directions for the pixels of projections.
unit_vector_k (ndarray[Any, dtype[float]]) – The other direction for the pixels of projections.
offsets_j (ndarray[Any, dtype[float]]) – Offsets which align projections in the direction of j
offsets_k (ndarray[Any, dtype[float]]) – Offsets which align projections in the direction of k.
float_type (str) – Whether to use ‘float64’ (default) or ‘float32’. The argument should be supplied as a string. The types of field and projections must match this type.

Return type

callable

Notes

The computation performed by this function may be written as

\[V_i(\mathbf{x}) = \sum_{s=0}^{N-1} \cdot \sum_{t = 0}^{4} t_w p_i(s, \mathbf{x} \cdot \mathbf{u} + \Delta_J + t_J, \mathbf{x} \cdot \mathbf{w} + \Delta_K + t_K)\]

\(N\) is the total number of projections. \(V_i\) is volume[:, i], and \(\mathbf{v}\) is projection_vector[s]. \(\mathbf{u}\) is unit_vector_j[s], \(\mathbf{w}\) is unit_vector_k[s]. \(\Delta_J\) and \(\Delta_K\) are additional offsets based on the unit vectors, shapes, and offsets, which align the centers of the projection and volume, so that the intersection of each ray is correctly computed. \(t_W\) are weights for bilinear interpolation between the four pixels nearest to each ray; \(t_J\) and \(t_K\) are the index offsets necessary to perform this interpolation.

mumott.core.john_transform_cuda.john_transform_adjoint_cuda(field, projections, unit_vector_p, unit_vector_j, unit_vector_k, offsets_j, offsets_k)[source]¶

Frontend for performing the adjoint of the John transform with parallel GPU computing.

Parameters

field (ndarray[float]) – The field into which the adjoint is projected, with 4 dimensions. The last index should have the same size as the last index of projections. Can be either a numpy.ndarray, or a device array that implements the CUDA array interface. If a device array is given, no copying to device is needed.
projections (ndarray[float]) – The projections from which the adjoint is calculated. The first index runs over the different projection directions. Can be either a numpy.ndarray, or a device array that implements the CUDA array interface. If a device array is given, no copying or synchronization is needed.
unit_vector_p (ndarray[float]) – The direction of projection in Cartesian coordinates.
unit_vector_j (ndarray[float]) – One of the directions for the pixels of projection.
unit_vector_k (ndarray[float]) – The other direction for the pixels of projection.
offsets_j (ndarray[float]) – Offsets which align projections in the direction of j
offsets_k (ndarray[float]) – Offsets which align projections in the direction of k.

Notes

For mathematical details, see john_transform_adjoint().

mumott.core.john_transform_cuda.john_transform_cuda(field, projections, unit_vector_p, unit_vector_j, unit_vector_k, offsets_j, offsets_k)[source]¶

Frontend for performing the John transform with parallel GPU computing.

Parameters

field (ndarray[float]) – The field to be projected, with 4 dimensions. The last index should have the same size as the last index of projections. Can be either a numpy.ndarray, or a device array that implements the CUDA array interface. If a device array is given, no copying to device is needed.
projections (ndarray[float]) – A 4-dimensional numpy array where the projections are stored. The first index runs over the different projection directions. Can be either a numpy.ndarray, or a device array that implements the CUDA array interface. If a device array is given, no copying or synchronization is needed.
unit_vector_p (ndarray[float]) – The direction of projection in Cartesian coordinates.
unit_vector_j (ndarray[float]) – One of the directions for the pixels of projection.
unit_vector_k (ndarray[float]) – The other direction for the pixels of projection.
offsets_j (ndarray[float]) – Offsets which align projections in the direction of j
offsets_k (ndarray[float]) – Offsets which align projections in the direction of k.

Notes

For mathematical details, see john_transform().

mumott.core.john_transform_sparse_cuda.john_transform_adjoint_sparse_cuda(field, projections, sparse_matrix, unit_vector_p, unit_vector_j, unit_vector_k, offsets_j, offsets_k)[source]¶

Frontend for performing the adjoint of the John transform with parallel GPU computing, with a sparse index for projection-to-field mapping.

Parameters

field (ndarray[float]) – The field into which the adjoint is projected, with 4 dimensions. The last index should have the same size as the number of tensor channels. Can be either a numpy.ndarray, or a device array that implements the CUDA array interface. If a device array is given, no copying to device is needed.
projections (ndarray[float]) – The projections from which the adjoint is calculated. Should have the same shape as the input data matrix. The first index runs over the different projection directions. Can be either a numpy.ndarray, or a device array that implements the CUDA array interface. If a device array is given, no copying or synchronization is needed.
sparse_matrix (tuple) – Contains three np.ndarrays, indices, indptr, and data. indices[indptr[i, j]] to indices[indptr[i, j+1]] where i is the projection channel index and j is the projection channel index contains the corresponding volume channel indices. data[indptr[i, j]] contains the weights for each index.
unit_vector_p (ndarray[float]) – The direction of projection in Cartesian coordinates.
unit_vector_j (ndarray[float]) – One of the directions for the pixels of projections.
unit_vector_k (ndarray[float]) – The other direction for the pixels of projections.
offsets_j (ndarray[float]) – Offsets which align projections in the direction of j
offsets_k (ndarray[float]) – Offsets which align projections in the direction of k.

Notes

For mathematical details, see john_transform_adjoint().

mumott.core.john_transform_sparse_cuda.john_transform_sparse_cuda(field, projections, sparse_matrix, unit_vector_p, unit_vector_j, unit_vector_k, offsets_j, offsets_k)[source]¶

Frontend for performing the John transform with parallel GPU computing, with a sparse matrix for projection-to-field mapping.

Parameters

field (ndarray[float]) – The field to be projected, with 4 dimensions. The last index should have the same size as the number of tensor channels. Can be either a numpy.ndarray, or a device array that implements the CUDA array interface. If a device array is given, no copying to device is needed.
projections (ndarray[float]) – A 4-dimensional numpy array where the projections are stored. Should have the same shape as the input data matrix. The first index runs over the different projection directions. Can be either a numpy.ndarray, or a device array that implements the CUDA array interface. If a device array is given, no copying or synchronization is needed.
sparse_matrix (tuple) – Contains three np.ndarrays, indices, indptr, and data. indices[indptr[i, j]] to indices[indptr[i, j+1]] where i is the projection channel index and j is the projection channel index contains the corresponding volume channel indices. data[indptr[i, j]] contains the weights for each index.
unit_vector_p (ndarray[float]) – The direction of projection in Cartesian coordinates.
unit_vector_j (ndarray[float]) – One of the directions for the pixels of projections.
unit_vector_k (ndarray[float]) – The other direction for the pixels of projections.
offsets_j (ndarray[float]) – Offsets which align projections in the direction of j
offsets_k (ndarray[float]) – Offsets which align projections in the direction of k

Notes

For mathematical details, see john_transform().

This transform also encodes the detector-to-tensor transform using a sparse matrix.

mumott.core.cuda_kernels.cuda_difference(shape)[source]¶

Computes a CUDA kernel for the difference of 2 4D arrays, e.g. a gradient and a value

Parameters: shape (tuple[int]) – A 4-tuple of the shape of the two gradients.
Returns: A CUDA callable which takes an old_gradient input/output, and a new_gradient input. The sum is stored in-place in old_gradient.

mumott.core.cuda_kernels.cuda_framewise_contraction(shape, rows, columns)[source]¶

Computes a CUDA kernel for the framewise contraction of a tensor field and a matrix stack: out[i, j, k, g] = sum_h(field[i, j, k, h], matrix[i, g, h]). In numpy.einsum notation, this would be 'ijkh, igh -> ijkg'.

Parameters

shape (tuple[int]) – A 3-tuple giving the shapes of the first three dimensions of the field ((i, j, k) dimensions).
rows (int) – The number of rows in the matrix/output vector length (g dimension).
columns (int) – The number of columns in the matrix/input vector length (h dimension).

Returns

A CUDA callable which takes field and matrix inputs, and out output.

mumott.core.cuda_kernels.cuda_framewise_contraction_adjoint(shape, rows, columns)[source]¶

Computes a CUDA kernel for the adjoint of the framewise contraction of a tensor field and a matrix stack: out[i, j, k, h] = sum_h(field[i, j, k, g], matrix[i, g, h]). In numpy.einsum notation, this would be 'ijkg, igh -> ijkh'.

Parameters

shape (tuple[int]) – A 3-tuple giving the shapes of the first three dimensions of the field ((i, j, k) dimensions).
rows (int) – The number of rows in the matrix/input vector length (g dimension).
columns (int) – The number of columns in the matrix/output vector length (h dimension).

Returns

A CUDA callable which takes field and matrix inputs, and out output.

mumott.core.cuda_kernels.cuda_l1_gradient(shape, weight=0.0001)[source]¶

Compiles a CUDA kernel for the gradient of an L1 regularizer.

Parameters

shape (tuple[int]) – The shape of the coefficients to which the gradient will ultimately be applied.
weight (float) – The weight of the L1 gradient.

Returns

A compiled CUDA callable which takes two arrays, an input coefficients array and an output gradient array, both of shape shape and dtype float32.

mumott.core.cuda_kernels.cuda_lower_bound(shape, lower_bound=0.0)[source]¶

Compiles a CUDA kernel for the enforcement of a lower bound to a 4-dimensional field. The computation is field[i, j, k, h] = max(field[i, j, k, h], lower_bound).

Parameters: shape (tuple[int]) – The shape of the coefficients to threshold with the lower bound.
Returns: A compiled CUDA callable which takes one array, an input/output field.

mumott.core.cuda_kernels.cuda_rescale(shape, momentum=0.9)[source]¶

Compiles a CUDA kernel for the rescaling of a gradient with a momentum term, or similar rescaling of a 4-dimensional array with a scalar.

Parameters

shape (tuple[int]) – The shape of the coefficients to which the gradient will ultimately be applied.
momentum (float) – The momentum weight, from 0 to 1. Default is 0.9.

Returns

A compiled CUDA callable which takes one array, an input/output gradient.

mumott.core.cuda_kernels.cuda_rescale_array(shape)[source]¶

Compiles a CUDA kernel for the rescaling of a gradient with a momentum term, or similar rescaling of a 4-dimensional array with another 4-dimensional array.

Parameters: shape (tuple[int]) – The shape of the coefficients to which the gradient will ultimately be applied.
Returns: A compiled CUDA callable which takes one array, an input/output gradient.

mumott.core.cuda_kernels.cuda_scaled_difference(shape)[source]¶

Compiles a CUDA kernel for a ‘scaled difference’, i.e., a -= b * c for 3 4-dimensional arrays, e.g. a data, gradient, and preconditioner array.

Parameters: shape (tuple[int]) – The shape of a and b as a 4-tuple.
Returns: A CUDA callable which takes 3 inputs, a gradient, value, and scaling. The output is stored in value. All inputs must be 4D arrays with shape shape.

mumott.core.cuda_kernels.cuda_sum(shape)[source]¶

Computes a CUDA kernel for the summation of 2 4D arrays, e.g. 2 gradients.

Parameters: shape (tuple[int]) – A 4-tuple of the shape of the two gradients.
Returns: A CUDA callable which takes an old_gradient input/output, and a new_gradient input. The sum is stored in-place in old_gradient.

mumott.core.cuda_kernels.cuda_tv_gradient(shape, weight=0.0001)[source]¶

Compiles a CUDA kernel for the gradient of a Total Variation regularizer. Gradient values at edges are sets to 0.

Parameters

shape (tuple[int]) – The shape of the coefficients to which the gradient will ultimately be applied.
weight (float) – The weight of the TV gradient.

Returns

A compiled CUDA callable which takes two arrays, an input coefficients array and an output gradient array, both of shape shape and dtype float32. The gradient array will have the value of the TV gradient added to it.

mumott.core.cuda_kernels.cuda_weighted_difference(shape)[source]¶

Compiles a CUDA kernel for a ‘weighted difference’, i.e. a = (a - b) * c. For example, a could be an approximation of b, and c could be the weight to assign to the residual of a and b.

Parameters: shape (tuple[int]) – The shape of a, b, and c.
Returns: A CUDA callable that takes 3 inputs, data, value, and weights, and stores the output in value. The difference is computed as ((value * weights) - data * weights).

mumott.core.cuda_kernels.cuda_weighted_sign(shape, delta=0.0)[source]¶

Compiles a CUDA kernel for a ‘weighted sign’, i.e. a = sgn(a - b) * c. For example, a could be an approximation of b, and c could be the weight to assign to the residual of a and b.

If delta is set to be greater than 0, then this function will return (a - b) * c / (2 * delta) if abs(a - b) < delta.

Parameters

shape (tuple[int]) – The shape of a, b, and c.
delta (float) – Threshold at which to switch from sign to actual difference.

Returns

A CUDA callable that takes 3 inputs, data, value, and weights, and stores the output in value.

mumott.core.wigner_d_utilities.calculate_sph_coefficients_rotated_around_z(input_array, angle, ell_list, output_array=None)[source]¶

Calculate the spherical harmonic coefficients after a rotation around +z.

Parameters

input_array (ndarray[Any, dtype[float]]) – spherical harmonic coefficients
angle (ndarray[Any, dtype[float]]) – either an array with the same shape as input_array.shape[:-1] or a scalar. If a scalar is given all coefficient lists are rotated by the same angle.
ell_list (List[int]) – list of \(\ell\) values to use. If None (default) we assume that all even orders are inluded until some \(\ell_\mathrm{max}\), which will be calculated from the shape of the coefficients array.
d_matrices – list of Wigner (small) d-matrices corresponding to a 90 degree rotation around x. If None (default) pre-calculated matrix elements are used.
output_array (Optional[ndarray[Any, dtype[float]]]) – numpy array with same shape as input_array. If None (default) a new array is initialized. If output_array is given, the calculations are carried out in place.

Return type

ndarray[Any, dtype[float]]

Returns

numpy array with same shape as input_array.

mumott.core.wigner_d_utilities.calculate_sph_coefficients_rotated_around_z_derived_wrt_the_angle(input_array, angle, ell_list, output_array=None)[source]¶

Calculate the angular derivatives of spherical harmonic coefficients with respect at the specified angle, where the angle refers to the rotation around +z.

Parameters

input_array (ndarray[Any, dtype[float]]) – spherical harmonic coefficients
angle (List[float]) – either an array with the same shape as input_array.shape[:-1] or a scalar. If a scalar is given all coefficient lists are rotated by the same angle.
ell_list (List[int]) – list of \(\ell\) values to use. If None (default) we assume that all even orders are inluded until some \(\ell_\mathrm{max}\), which will be calculated from the shape of the coefficients array.
d_matrices – list of Wigner (small) d-matrices corresponding to a 90 degree rotation about x. If None (default) pre-calculated matrix elements are used.
output_array (Optional[ndarray[Any, dtype[float]]]) – numpy array with same shape as input_array. If None (default) a new array is initialized. If output_array is given, the calculations are carried out in place.

Return type

ndarray[Any, dtype[float]]

Returns

numpy array with same shape as input_array.

mumott.core.wigner_d_utilities.calculate_sph_coefficients_rotated_by_90_degrees_around_negative_x(input_array, ell_list, d_matrices, output_array=None)[source]¶

Calculate the spherical harmonic coefficients after a rotation by 90 degrees around -x.

Parameters

input_array (ndarray[Any, dtype[float]]) – spherical harmonic coefficients
ell_list (List[int]) – list of \(\ell\) values to use. If None (default) we assume that all even orders are inluded until some \(\ell_\mathrm{max}\), which will be calculated from the shape of the coefficients array.
d_matrices (List[ndarray[Any, dtype[float]]]) – list of Wigner (small) d-matrices corresponding to a 90 degree rotation around x. If None (default) pre-calculated matrix elements are used.
output_array (Optional[ndarray[Any, dtype[float]]]) – numpy array with same shape as input_array. If None (default) a new array is initialized. If output_array is given, the calculations are carried out in place.

Return type

ndarray[Any, dtype[float]]

Returns

numpy array with same shape as input_array.

mumott.core.wigner_d_utilities.calculate_sph_coefficients_rotated_by_90_degrees_around_positive_x(input_array, ell_list, d_matrices, output_array=None)[source]¶

Calculate the spherical harmonic coefficients after a rotation by 90 degrees around +x.

Parameters

input_array (ndarray[Any, dtype[float]]) – spherical harmonic coefficients
ell_list (List[int]) – list of \(\ell\) values to use. If None (default) we assume that all even orders are inluded until some \(\ell_\mathrm{max}\), which will be calculated from the shape of the coefficients array.
d_matrices (List[ndarray[Any, dtype[float]]]) – list of Wigner (small) d-matrices corresponding to a 90 degree rotation around x. If None (default) pre-calculated matrix elements are used.
output_array (Optional[ndarray[Any, dtype[float]]]) – numpy array with same shape as input_array. If None (default) a new array is initialized. If output_array is given, the calculations are carried out in place.

Return type

ndarray[Any, dtype[float]]

Returns

numpy array with the same shape as input_array.

mumott.core.wigner_d_utilities.calculate_sph_coefficients_rotated_by_euler_angles(input_array, Psi, Theta, Phi, ell_list=None, d_matrices=None, output_array=None)[source]¶

Calculate the spherical harmonics coefficients after a rotation specified by a set of Euler angles. The Euler angles need to be given in ‘zyz’ format. A rotation with (0, Theta, Phi) will move the z-axis into the coordinates (Theta, Phi).

Parameters

input_array (ndarray[Any, dtype[float]]) – array, the last dimension of which runs over the spherical harmonics coefficients.
Psi (ndarray[Any, dtype[float]]) – First Euler angle. Initial rotation about the z-axis. Can be either a scalar or a numpy array with shape input_array.shape[:-1]. If None the rotation is skipped.
Theta (ndarray[Any, dtype[float]]) – Second Euler angle. Rotation about the y-axis. Can be either a scalar or a numpy array with shape input_array.shape[:-1]. If None the rotation is skipped.
Phi (ndarray[Any, dtype[float]]) – Third Euler angle. Final rotation about the z-axis. Can be either a scalar or a numpy array with shape input_array.shape[:-1]. If None the rotation is skipped.
ell_list (Optional[List[int]]) – list of \(\ell\) values to use. If None (default) we assume that all even orders are inluded until some \(\ell_\mathrm{max}\), which will be calculated from the shape of the coefficients array.
d_matrices (Optional[List[ndarray[Any, dtype[TypeVar(_ScalarType_co, bound= generic, covariant=True)]]]]) – list of Wigner (small) d-matrices corresponding to a 90 degree rotation around x. If None (default) pre-calculated matrix elements are used.
output_array (Optional[ndarray[Any, dtype[float]]]) – numpy array with same shape as input_array. If None (default) a new array is initialized. If output_array is given, the calculations are carried out in place.

Return type

ndarray[Any, dtype[float]]

Returns

numpy array with same shape as input_array.

mumott.core.wigner_d_utilities.load_d_matrices(ell_max)[source]¶

Load Wigner (small) d-matrices corresponding to a rotation by 90 degrees around the x-axis up to order ell_max. The sign and normalization conventions are those adopted throughout mumott.

Parameters: ell_max (int) – maximum order of the Wigner (small) d-matrices to return
Return type: List[ndarray[Any, dtype[float]]]
Returns: list of Wigner (small) d-matrices as numpy arrays with shape (l*2+1)x(l*2+1)