Source code for mumott.output_handling.reconstruction_derived_quantities

from dataclasses import dataclass
import numpy as np


def get_sorted_eigenvectors(tensors_array: np.array):
    """ Caluculate eigenvectors and eigenvalues of an array of
    3 by 3 matrices and sort according to increasing eigenvalues.

    Parameters
    ----------
    tensors_array
        numpy array containing the 3 by 3 tensors with the tensor
        indicies as the two last.

    Returns
    -------
    eigenvalues
        numpy array containing the eigenvalues of the tensors where
        the last dimension indexes the eigenvalues. Smallest eigenvalues
        come first.
    eigenvectors
        numpy array containing the eigenvalues of the tensors where
        the last dimension indexes the eigenvalues and the second
        to last dimension is the vector index.
    """

    volume_shape = tensors_array.shape[:-2]

    # Compute and sort eigenvectors
    w, v = np.linalg.eigh(tensors_array.reshape(-1, 3, 3))
    sorting = np.argsort(w, axis=1).reshape(-1, 3, 1)
    v = v.transpose(0, 2, 1)
    v = np.take_along_axis(v, sorting, axis=1)
    v = v.transpose(0, 2, 1)
    v = v / np.sqrt(np.sum(v ** 2, axis=1).reshape(-1, 1, 3))
    eigenvectors = v.reshape(volume_shape + (3, 3,))
    eigenvalues = np.sort(w, axis=-1).reshape(volume_shape + (3,))

    # Flip eigenvectors to have a positive z-component
    for ii in range(3):
        whereflip = eigenvectors[..., 2, ii] < 0.0
        eigenvectors[whereflip, :, ii] = -eigenvectors[whereflip, :, ii]

    return eigenvalues, eigenvectors


[docs]@dataclass class ReconstructionDerivedQuantities: """ A number of useful quantities that have been computed from the coefficients of a reconstruction. Attributes ---------- volume_shape : tuple The shape if the reconstructed volume. mean_intensity : np.array The mean intensity of the reconstructed functions over the unit sphere for each voxel. fractional_anisostropy : np.array A measure of the relative amount of anisotropic scattering in each voxel. eigenvector_1 eigenvector_2 eigenvector_3 : np.array The eingenvectors of the second order tensor component of the function. Sorted by ascending eigenvalue. eigenvalue_1 eigenvalue_2 eigenvalue_3 : np.array The eingenvalues of the second order tensor component of the function. Sorted by ascending eigenvalue. second_moment_tensor : np.array The second-moment tensor, which includes both the zero-th and second-order parts of the reconstructed functions. """ volume_shape: tuple mean_intensity: np.array fractional_anisotropy: np.array eigenvector_1: np.array eigenvector_2: np.array eigenvector_3: np.array eigenvalue_1: np.array eigenvalue_2: np.array eigenvalue_3: np.array second_moment_tensor: np.array def save_to_disk(filename : str): raise NotImplementedError