From 73369fab557499419f5d645cc37374dd12551191 Mon Sep 17 00:00:00 2001 From: iback Date: Thu, 27 Aug 2026 12:45:10 +0000 Subject: [PATCH] feat(core): add point-search, boundary-normal and vertebra-angle helpers Upstreams three helpers that sibling repos had each re-implemented locally, and removes the duplicate vector helpers that already existed inside TPTBox. np_utils: - np_index / np_find_closest_point_index: row lookup and nearest-point search over an (N, D) point array. The nearest-point search tries an exact match first for integer (voxel) arrays, then falls back to scipy's cdist. - np_compute_boundary_normals: per-voxel surface normals where one label touches another, via a 6-connectivity contact test and the Sobel gradient of a Gaussian-smoothed mask. Complements np_compute_surface, which has no normals, and calculate_pca_normal_np, which returns a single axis for a whole segmentation. - np_unit_vector / np_angle_between: canonical implementations. angles.py and ray_casting.py both carried their own copy; both now alias these, and angles.py keeps its historic public names. poi_fun.vertebra_direction: - get_vert_direction_angles: angle between each of a vertebra's P/I/R direction vectors and the corresponding global axis. Useful as a QC signal for direction landmarks that are not orthonormal. Four existing call sites now use the shared point search instead of hand-rolling cdist/argmin: pixel_based_point_finder (get_nearest_neighbor and the point-set projection loop) and endplates (two norm(axis=1) distance computations). All four were verified equivalent to their previous behaviour over randomized inputs. Note on the source material: the upstream fast_cdist was documented as "much faster than scipy.cdist" but benchmarks 2-6x slower, since scipy's is compiled C while the einsum form materialises the difference array. Only its exact-match fast path was adopted, on top of cdist. Co-Authored-By: Claude Opus 5 --- TPTBox/core/np_utils.py | 128 ++++++++++++++++++ .../core/poi_fun/pixel_based_point_finder.py | 15 +- TPTBox/core/poi_fun/ray_casting.py | 5 +- TPTBox/core/poi_fun/vertebra_direction.py | 27 ++++ TPTBox/spine/spinestats/angles.py | 37 +---- TPTBox/spine/spinestats/poi_fun/endplates.py | 7 +- unit_tests/test_nputils.py | 87 ++++++++++++ unit_tests/test_poi_ops.py | 50 +++++++ 8 files changed, 306 insertions(+), 50 deletions(-) diff --git a/TPTBox/core/np_utils.py b/TPTBox/core/np_utils.py index c7200f48..ecadbb72 100755 --- a/TPTBox/core/np_utils.py +++ b/TPTBox/core/np_utils.py @@ -20,10 +20,13 @@ from scipy.ndimage import ( binary_erosion, center_of_mass, + convolve, distance_transform_edt, gaussian_filter, generate_binary_structure, + sobel, ) +from scipy.spatial.distance import cdist from skimage.measure import euler_number as _euler_number from skimage.measure import label as _label @@ -893,6 +896,131 @@ def np_point_coordinates( return surface_points +def np_unit_vector(vector: np.ndarray) -> np.ndarray: + """Returns the unit vector of the input vector. + + Args: + vector (np.ndarray): Any non-zero numeric array. + + Returns: + np.ndarray: Array with the same direction as ``vector`` but unit length. + """ + return vector / np.linalg.norm(vector) + + +def np_angle_between(v1, v2, degrees: bool = False) -> float: + """Calculates the angle between two vectors. + + Args: + v1: The first vector. + v2: The second vector. + degrees (bool, optional): Return the angle in degrees instead of radians. Defaults to False. + + Returns: + float: The angle between ``v1`` and ``v2``. + + Examples: + >>> np_angle_between((1, 0, 0), (0, 1, 0)) + 1.5707963267948966 + >>> np_angle_between((1, 0, 0), (0, 1, 0), degrees=True) + 90.0 + >>> np_angle_between((1, 0, 0), (-1, 0, 0)) + 3.141592653589793 + """ + rad = np.arccos(np.clip(np.dot(np_unit_vector(v1), np_unit_vector(v2)), -1.0, 1.0)) + return float(np.degrees(rad)) if degrees else float(rad) + + +def np_index(arr: np.ndarray, entry) -> np.ndarray: + """Finds the rows of a point array that exactly equal a given point. + + Args: + arr (np.ndarray): Array of shape ``(N, D)`` holding ``N`` points. + entry: The point to look for, broadcastable to shape ``(D,)``. + + Returns: + np.ndarray: 1-D integer array with the indices of every matching row. + Empty if the point is not present. + """ + arr = np.asarray(arr) + assert arr.ndim == 2, f"expected a (N, D) point array, got shape {arr.shape}" + return np.flatnonzero((arr == np.asarray(entry)).all(axis=1)) + + +def np_find_closest_point_index(point_arr: np.ndarray, point) -> int: + """Finds the index of the point in ``point_arr`` closest to ``point``. + + For integer point arrays (i.e. voxel coordinates) an exact match is tried first, which is + the common case when the query already sits on the grid; only if that fails is the full + distance search run. ``scipy``'s compiled ``cdist`` is used for the search -- a hand-rolled + ``einsum`` over ``point_arr - point`` is roughly 2-6x slower because it materialises the + difference array. + + Args: + point_arr (np.ndarray): Array of shape ``(N, D)`` holding the candidate points. + point: The query point, broadcastable to shape ``(D,)``. + + Returns: + int: Index into ``point_arr`` of the nearest point. When several points are equally + close, the lowest index is returned. + """ + arr = np.asarray(point_arr) + assert arr.ndim == 2, f"expected a (N, D) point array, got shape {arr.shape}" + assert len(arr) != 0, "cannot search an empty point array" + p = np.asarray(point) + if np.issubdtype(arr.dtype, np.integer): + # Only meaningful for grids; rounding a float query would silently snap it to integers. + matches = np_index(arr, np.round(p).astype(arr.dtype)) + if len(matches) != 0: + return int(matches[0]) + return int(np.argmin(cdist([p], arr)[0])) + + +def np_compute_boundary_normals( + arr: UINTARRAY, + label: int | Sequence[int], + other_label: int | Sequence[int], + sigma: float = 1.0, +) -> tuple[np.ndarray, np.ndarray]: + """Computes per-voxel surface normals where one label touches another. + + Voxels of ``label`` that share a face with ``other_label`` form the interface. The normal at + each of them is the gradient of a Gaussian-smoothed ``label`` mask, so it points *into* + ``label`` (the direction of increasing mask density); negate it to point outwards. + + Unlike :func:`TPTBox.core.poi_fun.ray_casting.calculate_pca_normal_np`, which returns a single + principal axis for a whole segmentation, this returns a normal per interface voxel. + + Args: + arr (UINTARRAY): 3-dimensional label array. + label (int | Sequence[int]): Label(s) whose interface voxels are returned. + other_label (int | Sequence[int]): Label(s) that ``label`` must touch to count as interface. + sigma (float, optional): Standard deviation of the Gaussian applied before the gradient. + Larger values give smoother, less noisy normals. Defaults to 1.0. + + Returns: + tuple[np.ndarray, np.ndarray]: ``(coords, normals)``, both of shape ``(N, 3)``. ``coords`` + are the integer voxel coordinates of the interface, ``normals`` the matching unit vectors. + Both are empty when the two labels do not touch. + """ + assert arr.ndim == 3, arr.ndim + mask = np_isin(arr, label) + other = np_isin(arr, other_label) + # 6-connectivity kernel: the six face neighbours, so a purely diagonal contact does not count. + kernel = np.zeros((3, 3, 3), dtype=np.uint8) + kernel[1, 1, 0] = kernel[1, 1, 2] = 1 + kernel[1, 0, 1] = kernel[1, 2, 1] = 1 + kernel[0, 1, 1] = kernel[2, 1, 1] = 1 + touching = mask & (convolve(other.astype(np.uint8), kernel, mode="constant") > 0) + coords = np.argwhere(touching) + if len(coords) == 0: + return coords, np.zeros((0, 3), dtype=float) + smoothed = gaussian_filter(mask.astype(float), sigma=sigma) + gradient = np.stack([sobel(smoothed, axis=a) for a in range(3)], axis=-1) + gradient /= np.linalg.norm(gradient, axis=-1, keepdims=True) + 1e-8 + return coords, gradient[touching] + + def np_connected_components( arr: UINTARRAY, label_ref: LABEL_REFERENCE | None = None, diff --git a/TPTBox/core/poi_fun/pixel_based_point_finder.py b/TPTBox/core/poi_fun/pixel_based_point_finder.py index b5e6ba63..a8a6b911 100644 --- a/TPTBox/core/poi_fun/pixel_based_point_finder.py +++ b/TPTBox/core/poi_fun/pixel_based_point_finder.py @@ -7,6 +7,7 @@ from scipy.spatial.distance import cdist from TPTBox import NII, POI, Logger_Interface, Print_Logger +from TPTBox.core.np_utils import np_find_closest_point_index from TPTBox.core.poi_fun._help import to_local_np from TPTBox.core.poi_fun.vertebra_direction import _get_sub_array_by_direction, get_direction, get_vert_direction_matrix from TPTBox.core.vert_constants import COORDINATE, DIRECTIONS, Location @@ -32,13 +33,9 @@ def get_nearest_neighbor( 1-D integer array ``(x, y, z)`` of the voxel in ``sr_msk`` with label ``region_label`` that minimises the Euclidean distance to ``p``. """ - if len(p.shape) == 1: - p = np.expand_dims(p, 1) - locs = np.where(sr_msk == region_label) - locs_array = np.array(list(locs)).T - distances = cdist(p.T, locs_array) - - return locs_array[distances.argmin()] + p = np.asarray(p).reshape(-1) + locs_array = np.array(list(np.where(sr_msk == region_label))).T + return locs_array[np_find_closest_point_index(locs_array, p)] def max_distance_ray_cast_pixel_level( @@ -248,9 +245,7 @@ def project_pois_onto_set_of_points(poi: POI, point_set: list[COORDINATE]) -> PO point_arr = np.asarray(point_set) for r, s, c in poi.items(): - distance_to_point = cdist_to_point(c, point_arr) - new_coord = point_arr[np.argmin(distance_to_point)] - poi_n[r, s] = new_coord + poi_n[r, s] = point_arr[np_find_closest_point_index(point_arr, c)] return poi_n diff --git a/TPTBox/core/poi_fun/ray_casting.py b/TPTBox/core/poi_fun/ray_casting.py index e645ac62..2807c389 100644 --- a/TPTBox/core/poi_fun/ray_casting.py +++ b/TPTBox/core/poi_fun/ray_casting.py @@ -8,6 +8,7 @@ from sklearn.decomposition import PCA from TPTBox import NII, POI, Print_Logger, Vertebra_Instance +from TPTBox.core.np_utils import np_unit_vector from TPTBox.core.poi_fun._help import sacrum_w_o_arcus, to_local_np from TPTBox.core.poi_fun.pixel_based_point_finder import get_direction from TPTBox.core.vert_constants import COORDINATE, DIRECTIONS, Location @@ -16,9 +17,7 @@ _log = Print_Logger() -def unit_vector(vector: np.ndarray) -> np.ndarray: - """Returns the unit vector of the vector.""" - return vector / np.linalg.norm(vector) +unit_vector = np_unit_vector # @njit(fastmath=True) diff --git a/TPTBox/core/poi_fun/vertebra_direction.py b/TPTBox/core/poi_fun/vertebra_direction.py index babec962..7a510819 100644 --- a/TPTBox/core/poi_fun/vertebra_direction.py +++ b/TPTBox/core/poi_fun/vertebra_direction.py @@ -8,6 +8,7 @@ from numpy.linalg import norm from TPTBox import NII, POI, Print_Logger, calc_poi_from_subreg_vert +from TPTBox.core.np_utils import np_angle_between from TPTBox.core.poi_fun._help import make_spine_plot, sacrum_w_o_direction from TPTBox.core.vert_constants import DIRECTIONS, Location, Vertebra_Instance, _plane_dict, never_called @@ -409,6 +410,32 @@ def get_vert_direction_matrix(poi: POI, vert_id: int, to_pir: bool = False) -> t return to_vert_orient, from_vert_orient +def get_vert_direction_angles(poi: POI, vert_id: int, to_pir: bool = False, degrees: bool = True) -> tuple[float, float, float]: + """Return how far a vertebra's local frame is tilted from the global PIR axes. + + Each of the vertebra's Posterior/Inferior/Right direction vectors is compared with the + corresponding global axis, giving one angle per axis. Useful as a quality check: a healthy + frame stays within a few degrees of orthogonal to its neighbours, whereas an implausibly + large angle (e.g. a "posterior" direction more than 90 degrees from global posterior) + indicates that the direction landmarks are wrong. + + Args: + poi: ``POI`` object with pre-computed vertebra direction landmarks. + vert_id: Vertebra identifier (integer label). + to_pir: Whether to convert the POI to isotropic PIR space before computing. + Defaults to ``False``. + degrees: Return the angles in degrees rather than radians. Defaults to ``True``. + + Returns: + Tuple of three angles ``(posterior, inferior, right)`` between the vertebra's direction + vectors and the global PIR axes. + """ + directions = get_vert_direction_PIR(poi, vert_id=vert_id, to_pir=to_pir) + global_pir = (np.array([1, 0, 0]), np.array([0, 1, 0]), np.array([0, 0, 1])) + a, b, c = (np_angle_between(v, g, degrees=degrees) for v, g in zip(directions, global_pir)) + return a, b, c + + def calc_center_spinal_cord( poi: POI, subreg: NII, diff --git a/TPTBox/spine/spinestats/angles.py b/TPTBox/spine/spinestats/angles.py index ce2fda4d..3a6af2d1 100644 --- a/TPTBox/spine/spinestats/angles.py +++ b/TPTBox/spine/spinestats/angles.py @@ -10,6 +10,7 @@ from TPTBox import POI, Image_Reference from TPTBox.core.compat import zip_strict from TPTBox.core.nii_wrapper import to_nii +from TPTBox.core.np_utils import np_angle_between, np_unit_vector from TPTBox.core.vert_constants import DIRECTIONS, Location, Vertebra_Instance from TPTBox.spine.snapshot2D.snapshot_modular import Snapshot_Frame, create_snapshot @@ -158,39 +159,9 @@ def get_stop_vert(self, poi) -> Vertebra_Instance: } -def unit_vector(vector: np.ndarray) -> np.ndarray: - """Return the unit vector of the input vector. - - Args: - vector: Any non-zero numeric array. - - Returns: - Array with the same direction as ``vector`` but unit length. - """ - return vector / np.linalg.norm(vector) - - -def angle_between(v1, v2) -> float: - """Calculates the angle in radians between two vectors. - - Args: - v1 (tuple): The first vector. - v2 (tuple): The second vector. - - Returns: - float: The angle in radians between vectors 'v1' and 'v2'. - - Examples: - >>> angle_between((1, 0, 0), (0, 1, 0)) - 1.5707963267948966 - >>> angle_between((1, 0, 0), (1, 0, 0)) - 0.0 - >>> angle_between((1, 0, 0), (-1, 0, 0)) - 3.141592653589793 - """ - v1_u = unit_vector(v1) - v2_u = unit_vector(v2) - return np.arccos(np.clip(np.dot(v1_u, v2_u), -1.0, 1.0)) +# Canonical implementations live in np_utils; re-exported here under their historic names. +unit_vector = np_unit_vector +angle_between = np_angle_between def get_to_space(a, b, c) -> tuple[np.ndarray, np.ndarray]: diff --git a/TPTBox/spine/spinestats/poi_fun/endplates.py b/TPTBox/spine/spinestats/poi_fun/endplates.py index 19294948..48c87302 100644 --- a/TPTBox/spine/spinestats/poi_fun/endplates.py +++ b/TPTBox/spine/spinestats/poi_fun/endplates.py @@ -8,6 +8,7 @@ from stl.mesh import Mesh from TPTBox import NII, POI, Location, Logger_Interface, Print_Logger +from TPTBox.core.poi_fun.pixel_based_point_finder import cdist_to_point from TPTBox.core.vert_constants import Vertebra_Instance _log = Print_Logger() @@ -67,9 +68,7 @@ def _ray_cast_to_mesh(mesh: Mesh | trimesh.Trimesh, origin: np.ndarray, directio if len(locations) == 0: return None # closest hit - # d = np.linalg.norm(locations - origin, axis=1) - # return locations[np.argmin(d)] - d = np.linalg.norm(locations - origin, axis=1) + d = cdist_to_point(origin, locations) return origin + d.mean() * direction # numpy-stl fallback ts = [] @@ -133,7 +132,7 @@ def _local_curvature_grid( """ world_pts = poi.local_to_global_arr(voxel_pts_full) casted_world = np.asarray(poi.local_to_global(tuple(casted_full)), dtype=float) - dists = np.linalg.norm(world_pts - casted_world, axis=1) + dists = cdist_to_point(casted_world, world_pts) neighborhood = world_pts[dists <= radius] if len(neighborhood) < 3: return 0.0 diff --git a/unit_tests/test_nputils.py b/unit_tests/test_nputils.py index 4477540f..8a3edca6 100755 --- a/unit_tests/test_nputils.py +++ b/unit_tests/test_nputils.py @@ -463,6 +463,93 @@ def test_np_binary_fill_holes_and_set_inter_labels_based_on_majority(self): self.assertTrue(np.array_equal(filled, expected)) +class Test_point_helpers(unittest.TestCase): + def test_np_index_single_and_multiple_hits(self): + arr = np.array([[1, 2, 3], [4, 5, 6], [1, 2, 3]]) + self.assertTrue(np.array_equal(np_utils.np_index(arr, [4, 5, 6]), [1])) + self.assertTrue(np.array_equal(np_utils.np_index(arr, [1, 2, 3]), [0, 2])) + + def test_np_index_no_hit_is_empty(self): + arr = np.array([[1, 2, 3], [4, 5, 6]]) + idx = np_utils.np_index(arr, [7, 8, 9]) + self.assertEqual(len(idx), 0) + + def test_np_find_closest_point_index_exact_match_returns_first(self): + arr = np.array([[0, 0, 0], [5, 5, 5], [0, 0, 0]]) + # both row 0 and row 2 match exactly; the lowest index wins + self.assertEqual(np_utils.np_find_closest_point_index(arr, [0, 0, 0]), 0) + + def test_np_find_closest_point_index_matches_bruteforce(self): + rng = np.random.default_rng(42) + for dtype in (np.int32, np.float64): + for _ in range(repeats): + arr = (rng.random((200, 3)) * 50).astype(dtype) + point = (rng.random(3) * 50).astype(dtype) + expected = int(np.argmin(np.linalg.norm(arr - point, axis=1))) + got = np_utils.np_find_closest_point_index(arr, point) + # ties may resolve to a different index, so compare distances not indices + self.assertAlmostEqual( + float(np.linalg.norm(arr[got] - point)), + float(np.linalg.norm(arr[expected] - point)), + places=6, + ) + + def test_np_find_closest_point_index_rejects_empty(self): + with self.assertRaises(AssertionError): + np_utils.np_find_closest_point_index(np.zeros((0, 3), dtype=int), [0, 0, 0]) + + +class Test_boundary_normals(unittest.TestCase): + def make_two_slabs(self): + """Label 1 fills x < 6, label 2 fills x >= 6, so the interface is the plane x == 5/6.""" + arr = np.zeros((12, 12, 12), dtype=np.uint8) + arr[:6] = 1 + arr[6:] = 2 + return arr + + def test_interface_voxels_are_on_the_boundary_of_the_first_label(self): + arr = self.make_two_slabs() + coords, normals = np_utils.np_compute_boundary_normals(arr, 1, 2) + self.assertEqual(len(coords), 12 * 12) + self.assertTrue((coords[:, 0] == 5).all()) + self.assertTrue((arr[tuple(coords.T)] == 1).all()) + self.assertEqual(len(normals), len(coords)) + + def test_normals_are_unit_length_and_point_into_the_label(self): + arr = self.make_two_slabs() + _, normals = np_utils.np_compute_boundary_normals(arr, 1, 2) + self.assertTrue(np.allclose(np.linalg.norm(normals, axis=1), 1.0, atol=1e-3)) + # label 1 lies towards -x, and the normal follows increasing mask density + self.assertTrue(np.allclose(normals[:, 0], -1.0, atol=1e-3)) + self.assertTrue(np.allclose(normals[:, 1:], 0.0, atol=1e-3)) + + def test_labels_that_do_not_touch_give_empty_result(self): + arr = np.zeros((10, 10, 10), dtype=np.uint8) + arr[1:3] = 1 + arr[7:9] = 2 + coords, normals = np_utils.np_compute_boundary_normals(arr, 1, 2) + self.assertEqual(coords.shape, (0, 3)) + self.assertEqual(normals.shape, (0, 3)) + + def test_diagonal_only_contact_does_not_count(self): + arr = np.zeros((6, 6, 6), dtype=np.uint8) + arr[2, 2, 2] = 1 + arr[3, 3, 3] = 2 # shares only a corner, not a face + coords, _ = np_utils.np_compute_boundary_normals(arr, 1, 2) + self.assertEqual(len(coords), 0) + + +class Test_vector_helpers(unittest.TestCase): + def test_np_unit_vector(self): + self.assertTrue(np.allclose(np_utils.np_unit_vector(np.array([3.0, 4.0, 0.0])), [0.6, 0.8, 0.0])) + + def test_np_angle_between(self): + self.assertAlmostEqual(np_utils.np_angle_between((1, 0, 0), (0, 1, 0)), np.pi / 2) + self.assertAlmostEqual(np_utils.np_angle_between((1, 0, 0), (1, 0, 0)), 0.0) + self.assertAlmostEqual(np_utils.np_angle_between((1, 0, 0), (-1, 0, 0)), np.pi) + self.assertAlmostEqual(np_utils.np_angle_between((1, 0, 0), (0, 1, 0), degrees=True), 90.0) + + if __name__ == "__main__": unittest.main() diff --git a/unit_tests/test_poi_ops.py b/unit_tests/test_poi_ops.py index 6db0e2fd..17cbb445 100644 --- a/unit_tests/test_poi_ops.py +++ b/unit_tests/test_poi_ops.py @@ -352,5 +352,55 @@ def test_subtract_disjoint(self): self.assertEqual(len(result), len(p1)) +class Test_vert_direction_angles(unittest.TestCase): + def make_poi(self, frame: np.ndarray, corpus=(30.0, 30.0, 30.0)): + from TPTBox import POI, Location + + corpus = np.asarray(corpus, dtype=float) + poi = POI({}, orientation=("P", "I", "R"), zoom=(1, 1, 1), shape=(80, 80, 80)) + poi[23, Location.Vertebra_Corpus] = tuple(corpus) + poi[23, Location.Vertebra_Direction_Posterior] = tuple(corpus + 10 * frame[:, 0]) + poi[23, Location.Vertebra_Direction_Inferior] = tuple(corpus + 10 * frame[:, 1]) + poi[23, Location.Vertebra_Direction_Right] = tuple(corpus + 10 * frame[:, 2]) + return poi + + def test_frame_aligned_to_global_axes_has_zero_angles(self): + from TPTBox.core.poi_fun.vertebra_direction import get_vert_direction_angles + + poi = self.make_poi(np.eye(3)) + angles = get_vert_direction_angles(poi, vert_id=23, to_pir=False, degrees=True) + self.assertEqual(len(angles), 3) + for a in angles: + self.assertAlmostEqual(a, 0.0, places=5) + + def test_radians_and_degrees_agree(self): + from TPTBox.core.poi_fun.vertebra_direction import get_vert_direction_angles + + frame = np.array([[0.0, 0.0, 1.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]]) + poi = self.make_poi(frame) + deg = get_vert_direction_angles(poi, vert_id=23, to_pir=False, degrees=True) + rad = get_vert_direction_angles(poi, vert_id=23, to_pir=False, degrees=False) + for d, r in zip(deg, rad): + self.assertAlmostEqual(d, float(np.degrees(r)), places=6) + + def test_angles_are_consistent_with_the_direction_matrix(self): + """The angle triple must match the columns of ``from_vert_orient`` it is derived from.""" + from TPTBox.core.np_utils import np_angle_between + from TPTBox.core.poi_fun.vertebra_direction import get_vert_direction_angles, get_vert_direction_matrix + + rng = np.random.default_rng(7) + for _ in range(repeats): + frame, _ = np.linalg.qr(rng.random((3, 3))) + if np.linalg.det(frame) < 0: + frame[:, 0] *= -1 + poi = self.make_poi(frame) + _, from_vert_orient = get_vert_direction_matrix(poi, vert_id=23, to_pir=False) + angles = get_vert_direction_angles(poi, vert_id=23, to_pir=False, degrees=True) + for axis, angle in enumerate(angles): + global_axis = np.eye(3)[axis] + expected = np_angle_between(from_vert_orient[:, axis], global_axis, degrees=True) + self.assertAlmostEqual(angle, expected, places=6) + + if __name__ == "__main__": unittest.main()