Source code for monalisa.geometry._volume_element_center_out_radial3

"""Port of src/geom3/bmVolumeElement_center_out_radial3.m."""

import numpy as np

from monalisa._mat import msize
from monalisa.geometry._volume_element_voronoi_center_out_radial2 import _radial_dr
from monalisa.trajectory._traj_format import traj_format
from monalisa.trajectory._traj_line_reshape import traj_line_reshape
from monalisa.trajectory._traj_squared_norm import traj_squared_norm




def volume_element_center_out_radial3(t):
    """Volume elements of a 3D center-out radial trajectory (uniform solid angle)."""
    if msize(t, 1) != 3:
        raise ValueError("The trajectory must be 3Dim")
    t, _, formatedIndex, formatedWeight = traj_format(t, nargout=4)
    centerFlag = np.linalg.norm(t[:, 0]) < 100 * np.finfo(np.float64).eps
    if centerFlag:
        t = t[:, 1:]
    t = traj_line_reshape(t)
    N, nLine = msize(t, 2), msize(t, 3)
    t3 = np.reshape(t, (3, N, nLine), order="F")
    dr, _, r_1 = _radial_dr(t3, centerFlag)
    ds = np.reshape(4 * np.pi * traj_squared_norm(t) / nLine, (N, nLine), order="F")
    v = np.ravel(dr, order="F") * np.ravel(ds, order="F")
    if centerFlag:
        dr_0 = np.mean(r_1 / 2)
        v = np.concatenate([[(4 / 3) * np.pi * dr_0**3], v])
    idx = np.ravel(formatedIndex).astype(np.int64) - 1
    return (v[idx] * np.ravel(formatedWeight))[None, :]

↑