Core¶
- class mumott.core.geometry.Geometry(filename=None)[source]¶
Stores information about the system geometry. Instances of this class are used by
DataContainer
andProjectionStack
to maintain geometry information. They can be stored as a file using thewrite()
method. This allows one to (re)create aGeometry
instance from an earlier and overwrite the geometry information read by aDataContainer
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 toNone
, in which case the instance is created with default parameters.
- append(value)[source]¶
Appends projection-wise geometry data provided as a
GeometryTuple
.- Return type
- 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_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_axes: str¶
A blake2b hash of
outer_axes
.
- property hash_rotations: str¶
A blake2b hash of
rotations_as_array
.
- 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.
- insert(key, value)[source]¶
Inserts projection-wise data handed via a
GeometryTuple
.- Return type
- 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_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_as_array: ndarray[float]¶
Offsets to align projection in the direction k as an array.
- 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 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)
, whereN
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.
- 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, thenaxis
andangle
must be provided.axis (
Optional
[ndarray
[float
]]) – An axis, given as a unit length 3-vector, about which the rotation is defined. Not used ifA
is provided.angle (
Optional
[ndarray
[float
]]) – The angle in radians of the rotation aboutaxis
. Not used ifA
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, thenaxis
andangle
must be provided.axis (
Optional
[ndarray
[float
]]) – An axis, given as a 3-vector, about which a rotation can be defined. Not used ifA
is provided.angle (
Optional
[ndarray
[float
]]) – The angle in radians of the rotation aboutaxis
. Not used ifA
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
andsystem_rotations
will be applied to therotations
and system vectors respectively prior to writing.
- 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
andaxis
arguments are given, this matrix should correspond to the matrix given by R_outer @ R_inner, where R_inner is defined by a rotation byinner_angle
aboutinner_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 propertiesdata
,diode
, andweights
, 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 theDataContainer.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
of0.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 (seegeometry
).- Parameters
projection (
Projection
) –Projection
object to be appended.- Return type
- property data: ndarray[Any, dtype[_ScalarType_co]]¶
Scattering data, structured
(n, j, k, w)
, wheren
is the projection number,j
is the pixel in the j-direction,k
is the pixel in the k-direction, andw
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.
- 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 between0
and1
.0
means mask,1
means do not mask. Structured the same way asdata
.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 aboutinner_axis
in radians.outer_angle (
Optional
[float
]) – Angle of rotation aboutouter_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)
, wherej
is the pixel in the j-direction,k
is the pixel in the k-direction, andw
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 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.
- 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 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
anddegrees
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
- 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.
- property polar_zero: ndarray[Any, dtype[float]]¶
The polar angle of the spherical harmonic pole, relative to a fixed reference system.
- property unit_vectors¶
Probed coordinates object used by BasisSets to calculate function maps.
- 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)
whereN
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 as1
, but the array still needs to be structured in this way explicitly. By default, the value will benp.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 be1
). This vector is used when interpolating coordinates that lie on a small circle.
- class mumott.core.simulator.Simulator(volume_mask, basis_set, seed=None, distance_radius=1.4283556979968262)[source]¶
Simulator for tensor tomography samples based on a geometry and a few sources with associated influence functions. Designed primarily for local representations. Polynomial representations, such as spherical harmonics, may require other residual functions which take the different frequency bands into account.
- Parameters
volume_mask (np.ndarray[int]) – A three-dimensional mask. The shape of the mask defines the shape of the entire simulated volume. The non-zero entries of the mask determine where the simulated sample is located. The non-zero entries should be mostly contiguous for good results.
basis_set (BasisSet) – The basis set used for the simulation. Ideally local representations such as
GaussianKernels
should be used. Do not modify the basis set after creating the simulator.seed (int) – Seed for the random number generator. Useful for generating consistent simulations. By default no seed is used.
distance_radius (float) – Radius for the balls used in determining interior distances. Usually, this value should not be changed, but it can be increased to take larger strides in the interior of the sample. Default value is
np.sqrt(2) * 1.01
.
- add_source(location=None, coefficients=None, influence_exponent=2, influence_scale_parameter=10)[source]¶
Add a source to your simulation.
- Parameters
location (
Optional
[tuple
[int
,int
,int
]]) – Location, in terms of indices, of the source. If not given, will be randomized weighted by inverse of distance to other source points.coefficients (
Optional
[ndarray
[float
]]) – Coefficients defining the source. If not given, will be randomized in the interval[0, 1]
.influence_exponent (
float
) – Exponent of the influence of the influence of the source. Default value is2
, giving a Gaussian.influence_scale_parameter (
float
) – Scale parameter of the influence of the influence of the source. Default value is10
.
- Return type
Notes
The equation for the source influence is
np.exp(-(d / (p * s)) ** p)
, whered
is the interior distance,s
is the scale parameter, andp
is the exponent.If a location is not given, a location will be searched for for no more than 1e6 iterations. For large and sparse samples, consider specifying locations manually.
- property basis_set: BasisSet¶
Basis set defining the representation used in the sample. Read-only property; do not modify.
- optimize(step_size=0.01, iterations=10, weighting_iterations=5, momentum=0.95, tv_weight=0.1, tv_exponent=1.0, tv_delta=0.01, residual_exponent=1.0, residual_delta=0.05, power_weight=0.01, power_exponent=1, power_delta=0.01, lower_bound=0.0)[source]¶
Optimizer for carrying out the simulation. Can be called repeatedly. Uses iteratively reweighted least squares, with weights calculated based on the Euclidean norm over each voxel.
- Parameters
step_size (
float
) – Step size for gradient descent.iterations (
int
) – Number of iterations for each gradient descent solution.weighting_iterations (
int
) – Number of reweighting iterations.momentum (
float
) – Nesterov momentum term.tv_weight (
float
) – Weight for the total variation term.tv_exponent (
float
) – Exponent for the total variation reweighting. Default is 1, which will approach a Euclidean norm considered.tv_delta (
float
) – Huber-style cutoff for the total variation factor normalization.residual_exponent (
float
) – Exponent for the residual norm reweighting.residual_delta (
float
) – Huber-style cutoff for the residual normalization.power_weight (
float
) – Weight for the power term.power_exponent (
float
) – Exponent for the power term.power_delta (
float
) – Huber-style cutoff for the power normalization.lower_bound (
float
) – Lower bound for coefficients. Coefficients will be clipped at these bounds at every reweighting.
- property sources: dict¶
Dictionary of source properties. They are given as arrays where the first index specifies the source, so that
len(array)
is the number of source points.Notes
The items are, in order:
- coefficients
The coefficients of each source point.
- distances
The interior distance from each support point to each point in the volume.
- scale_parameters
The scale parameter of each source point.
- locations
The location of each source point.
- exponents
The exponent of each source point.
- influences
The influence of each source point.
- 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
- Return type
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 thedtype
beingobject
.>>> 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 similarnumpy
types.str
Automatically given
'utf-8'
encoding.bytes
Cast to string.
None
Ignored.
np.ndarray
Provided
dtype
is notobject
, hence arrays ofNone
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 ofprojections
.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 ofprojections`
.unit_vector_k (
ndarray
[Any
,dtype
[float
]]) – The other direction for the pixels ofprojections
.offsets_j (
ndarray
[Any
,dtype
[float
]]) – Offsets which align projections in the direction of joffsets_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 offield
andprojections
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\) isvolume[..., i]
, and \(\mathbf{v}\) isunit_vector_p[I]
. \(N\) is the number of voxels in the maximal direction ofunit_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\) isoffsets_k[I]
, \(\mathbf{u}\) isunit_vector_j[I]
, \(\mathbf{w}\) isunit_vector_k[I]
, \(J_\text{max}\) and \(K_\text{max}\) areprojections.shape[1]
andprojections.shape[2]
respectively, and \(\mathbf{r}_\text{max}\) isvolume.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 ofprojections
.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 ofprojections
.unit_vector_k (
ndarray
[Any
,dtype
[float
]]) – The other direction for the pixels ofprojections
.offsets_j (
ndarray
[Any
,dtype
[float
]]) – Offsets which align projections in the direction of joffsets_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 offield
andprojections
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}\) isprojection_vector[s]
. \(\mathbf{u}\) isunit_vector_j[s]
, \(\mathbf{w}\) isunit_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 ofprojections
. 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 ofprojection
.unit_vector_k (
ndarray
[float
]) – The other direction for the pixels ofprojection
.offsets_j (
ndarray
[float
]) – Offsets which align projections in the direction of joffsets_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 ofprojections
. 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 ofprojection
.unit_vector_k (
ndarray
[float
]) – The other direction for the pixels ofprojection
.offsets_j (
ndarray
[float
]) – Offsets which align projections in the direction of joffsets_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 anumpy.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 anumpy.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 threenp.ndarrays
,indices
,indptr
, anddata
.indices[indptr[i, j]]
toindices[indptr[i, j+1]]
wherei
is the projection channel index andj
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 ofprojections
.unit_vector_k (
ndarray
[float
]) – The other direction for the pixels ofprojections
.offsets_j (
ndarray
[float
]) – Offsets which align projections in the direction of joffsets_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 anumpy.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 anumpy.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 threenp.ndarrays
,indices
,indptr
, anddata
.indices[indptr[i, j]]
toindices[indptr[i, j+1]]
wherei
is the projection channel index andj
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 ofprojections
.unit_vector_k (
ndarray
[float
]) – The other direction for the pixels ofprojections
.offsets_j (
ndarray
[float
]) – Offsets which align projections in the direction of joffsets_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
- 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])
. Innumpy.einsum
notation, this would be'ijkh, igh -> ijkg'
.- Parameters
- Returns
A CUDA callable which takes
field
andmatrix
inputs, andout
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])
. Innumpy.einsum
notation, this would be'ijkg, igh -> ijkh'
.- Parameters
- Returns
A CUDA callable which takes
field
andmatrix
inputs, andout
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
- Returns
A compiled CUDA callable which takes two arrays, an input
coefficients
array and an outputgradient
array, both of shapeshape
and dtypefloat32
.
- 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)
.
- 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.
- 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.
- 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.
- mumott.core.cuda_kernels.cuda_sum(shape)[source]¶
Computes a CUDA kernel for the summation of 2 4D arrays, e.g. 2 gradients.
- 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
- Returns
A compiled CUDA callable which takes two arrays, an input
coefficients
array and an outputgradient
array, both of shapeshape
and dtypefloat32
. Thegradient
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 ofb
, andc
could be the weight to assign to the residual ofa
andb
.
- 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 ofb
, andc
could be the weight to assign to the residual ofa
andb
.If
delta
is set to be greater than 0, then this function will return(a - b) * c / (2 * delta)
ifabs(a - b) < delta
.
- 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 coefficientsangle (
ndarray
[Any
,dtype
[float
]]) – either an array with the same shape asinput_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. IfNone
(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 asinput_array
. IfNone
(default) a new array is initialized. Ifoutput_array
is given, the calculations are carried out in place.
- Return type
- 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 coefficientsangle (
List
[float
]) – either an array with the same shape asinput_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. IfNone
(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 asinput_array
. IfNone
(default) a new array is initialized. Ifoutput_array
is given, the calculations are carried out in place.
- Return type
- 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 coefficientsell_list (
List
[int
]) – list of \(\ell\) values to use. IfNone
(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. IfNone
(default) pre-calculated matrix elements are used.output_array (
Optional
[ndarray
[Any
,dtype
[float
]]]) – numpy array with same shape asinput_array
. IfNone
(default) a new array is initialized. Ifoutput_array
is given, the calculations are carried out in place.
- Return type
- 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 coefficientsell_list (
List
[int
]) – list of \(\ell\) values to use. IfNone
(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. IfNone
(default) pre-calculated matrix elements are used.output_array (
Optional
[ndarray
[Any
,dtype
[float
]]]) – numpy array with same shape asinput_array
. IfNone
(default) a new array is initialized. Ifoutput_array
is given, the calculations are carried out in place.
- Return type
- 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 shapeinput_array.shape[:-1]
. IfNone
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 shapeinput_array.shape[:-1]
. IfNone
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 shapeinput_array.shape[:-1]
. IfNone
the rotation is skipped.ell_list (
Optional
[List
[int
]]) – list of \(\ell\) values to use. IfNone
(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. IfNone
(default) pre-calculated matrix elements are used.output_array (
Optional
[ndarray
[Any
,dtype
[float
]]]) – numpy array with same shape asinput_array
. IfNone
(default) a new array is initialized. Ifoutput_array
is given, the calculations are carried out in place.
- Return type
- Returns
numpy array with same shape as
input_array
.