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

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 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.

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.

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.

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=array([[[[1., 0., 0.]]]]))[source]¶

A small container class for probed coordinates for the :class:` BasisSet <mumott.methods.BasisSet>` being used together with the Method.

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).

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.

Numba version of John transform.

mumott.core.john_transform.john_transform(projection, volume, unit_vector_p, unit_vector_j, unit_vector_k, step_size, projection_offsets, sampling_kernel, kernel_offsets)[source]¶

Frontend for performing the John transform with parallel CPU computing. The number of threads is regulated with numba.set_num_threads(). All arrays must be row-major and contiguous, also called C-contiguous.

Parameters

projection (ndarray[Any, dtype[Float]]) – A 3-dimensional numpy array where the projection is stored. The last index contains the data channels for multi-channel data. For scalar data, use 1 as the size of the last dimension.
volume (ndarray[Any, dtype[Float]]) – The volume to be projected, with 4 dimensions. The last index is the same as for projection.
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 projection.
unit_vector_k (ndarray[Any, dtype[Float]]) – The other direction for the pixels of projection.
step_size (Float) – The size of each step for the quadrature of the integral.
projection_offsets (ndarray[Any, dtype[Float]]) – The j, k offsets that align the projections.
sampling_kernel (Optional[ndarray[Any, dtype[Float]]]) – A kernel for uppsampling or convolution of the projection. Should be normalized so that it sums to unity.
kernel_offsets (Optional[ndarray[Any, dtype[Float]]]) – Offsets for each entry of the kernel.

Return type

None

Notes

The computation performed by this function may be written as

\[p(J, K)_i = \sum_{s=0}^{N} \sum_{j=0}^{M} d \cdot w_j \cdot V_i(\lfloor \mathbf{r}_j + d s \cdot \mathbf{v} \rfloor)\]

where \(p(J, K)_i\) is projection[J, K, i], \(d\) is step_size, \(N\) is the total number of steps, \(M\) is sampling_kernel.size, \(w_j\) is sampling_kernel[j], \(V_i\) is volume[..., i], and \(\mathbf{v}\) is unit_vector_p. mathbf{r}_j is the starting position for the ray weighted by sampling_kernel[j], 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 the sum of projection_offsets[0] and kernel_offsets[0, j], \(o_K\) is the sum of projection_offsets[1] and kernel_offsets[1, j], \(\mathbf{u}\) is unit_vector_j, \(mathbf{w}\) is unit_vector_k, \(J_\text{max}\) and :math`K_text{max}` are projection.shape[0] and projection.shape[1] 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(projection, volume, unit_vector_p, unit_vector_j, unit_vector_k, step_size, projection_offsets, sampling_kernel, kernel_offsets)[source]¶

Frontend for performing the John transform adjoint (back-projection) with parallel CPU computation. The number of threads is regulated by numba.set_num_threads.

Parameters

projection (ndarray[Any, dtype[Float]]) – A 3-dimensional numpy array for which the adjoint is calculated. The last index contains the data channels for multi-channel data. For scalar data, use 1 as the size of the last dimension.
volume (ndarray[Any, dtype[Float]]) – Output variable modified in-place. The volume where the adjoint is stored, with 5 dimensions. The last index is the same as the last index in projection. The first index is for each thread of multi-threaded computation, so volume.shape[0] must equal numba.get_num_threads. After computing, it is necessary to reduce over the first index to obtain the total result.
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 projection.
unit_vector_k (ndarray[Any, dtype[Float]]) – The other direction for the pixels of projection.
step_size (ndarray[Any, dtype[Float]]) – The size of each step for the quadrature of the integral.
projection_offsets (ndarray[Any, dtype[Float]]) – The j, k offsets that align the projections.
sampling_kernel (Optional[ndarray[Any, dtype[Float]]]) – A kernel for upsampling or convolution of the projection. Should be normalized so that it sums to unity.
kernel_offsets (Optional[ndarray[Any, dtype[Float]]]) – Offsets for each entry of the kernel.

Return type

None

Notes

The computation performed by this function may be written as

\[V_i(\mathbf{x}) = \sum_{s=0}^{N} \sum_{j=0}^{M} d \cdot w_j \cdot p(J, K)_i \delta_{\lfloor\mathbf{r}_j + d s \cdot \mathbf{v} \rfloor}^{\mathbf{x}}\]

where \(d_i\) is distributor[i], \(N\) is the total number of steps, \(d\) is step_size. \(M\) is sampling_kernel.size, \(w_j\) is sampling_kernel[j], \(V_i\) is volume[:, i], and \(\mathbf{v}\) is projection_vector. \(\mathbf{x}\) is understood to run over every triplet of non-negative integers in the domain of \(V_i\). mathbf{r}_j is the starting position for the ray weighted by sampling_kernel[j], 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 the sum of projection_offsets[0] and kernel_offsets[0, j], \(o_K\) is the sum of projection_offsets[1] and kernel_offsets[1, j], \(\mathbf{u}\) is unit_vector_j, \(mathbf{w}\) is unit_vector_k, \(J_\text{max}\) and :math`K_text{max}` are projection.shape[0] and projection.shape[1] 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_bilinear.john_transform_adjoint_bilinear(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

For mathematical details, see the :meth:<~.john_transform.john_transform_adjoint.

mumott.core.john_transform_bilinear.john_transform_bilinear(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

For mathematical details, see the john_transform().

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[Any, dtype[TypeVar(_ScalarType_co, bound= generic, covariant=True)]]) – 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[Any, dtype[TypeVar(_ScalarType_co, bound= generic, covariant=True)]]) – 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[Any, dtype[TypeVar(_ScalarType_co, bound= generic, covariant=True)]]) – The direction of projection in Cartesian coordinates.
unit_vector_j (ndarray[Any, dtype[TypeVar(_ScalarType_co, bound= generic, covariant=True)]]) – One of the directions for the pixels of projection.
unit_vector_k (ndarray[Any, dtype[TypeVar(_ScalarType_co, bound= generic, covariant=True)]]) – The other direction for the pixels of projection.
offsets_j (ndarray[Any, dtype[TypeVar(_ScalarType_co, bound= generic, covariant=True)]]) – Offsets which align projections in the direction of j
offsets_k (ndarray[Any, dtype[TypeVar(_ScalarType_co, bound= generic, covariant=True)]]) – Offsets which align projections in the direction of k.

Notes

For mathematical details, see the :meth:<~.john_transform.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[Any, dtype[TypeVar(_ScalarType_co, bound= generic, covariant=True)]]) – 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[Any, dtype[TypeVar(_ScalarType_co, bound= generic, covariant=True)]]) – 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[Any, dtype[TypeVar(_ScalarType_co, bound= generic, covariant=True)]]) – The direction of projection in Cartesian coordinates.
unit_vector_j (ndarray[Any, dtype[TypeVar(_ScalarType_co, bound= generic, covariant=True)]]) – One of the directions for the pixels of projection.
unit_vector_k (ndarray[Any, dtype[TypeVar(_ScalarType_co, bound= generic, covariant=True)]]) – The other direction for the pixels of projection.
offsets_j (ndarray[Any, dtype[TypeVar(_ScalarType_co, bound= generic, covariant=True)]]) – Offsets which align projections in the direction of j
offsets_k (ndarray[Any, dtype[TypeVar(_ScalarType_co, bound= generic, covariant=True)]]) – Offsets which align projections in the direction of k.

Notes

For mathematical details, see the :meth:<~.john_transform.john_transform.

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 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)