Fix principal direction extent

CHANGELOG.rst

@@ -4,6 +4,7 @@ Changelog
 Version 3.2.0
 -------------
 
+- Fix ``neurom.morphmath.principal_direction_extent`` to calculate correctly the pca extent.
 - Fix ``neurom.features.neurite.segment_taper_rates`` to return signed taper rates.
 - Fix warning system so that it doesn't change the pre-existing warnings configuration
 - Fix ``neurom.features.bifurcation.partition_asymmetry`` Uylings variant to not throw

neurom/morphmath.py

@@ -462,17 +462,17 @@ def sphere_area(r):
 def principal_direction_extent(points):
     """Calculate the extent of a set of 3D points.
 
-    The extent is defined as the maximum distance between
-    the projections on the principal directions of the covariance matrix
-    of the points.
+    The extent is defined as the maximum distance between the projections on the principal
+    directions of the covariance matrix of the points.
 
-    Parameter:
-        points : a 2D numpy array of points
+    Args:
+        points : a 2D numpy array of points with 2 or 3 columns for (x, y, z)
 
     Returns:
         extents : the extents for each of the eigenvectors of the cov matrix
-        eigs : eigenvalues of the covariance matrix
-        eigv : respective eigenvectors of the covariance matrix
+
+    Note:
+        Direction extents are not ordered from largest to smallest.
     """
     # pca can be biased by duplicate points
     points = np.unique(points, axis=0)
@@ -483,14 +483,8 @@ def principal_direction_extent(points):
     # principal components
     _, eigv = pca(points)
 
-    extent = np.zeros(3)
-
-    for i in range(eigv.shape[1]):
-        # orthogonal projection onto the direction of the v component
-        scalar_projs = np.sort(np.array([np.dot(p, eigv[:, i]) for p in points]))
-        extent[i] = scalar_projs[-1]
-
-        if scalar_projs[0] < 0.:
-            extent -= scalar_projs[0]
+    # for each eigenvector calculate the scalar projection of the points on it (n_points, n_eigv)
+    scalar_projections = points.dot(eigv)
 
-    return extent
+    # and return the range of the projections (abs(max - min)) along each column (eigenvector)
+    return np.ptp(scalar_projections, axis=0)

tests/features/test_get_features.py

@@ -793,18 +793,18 @@ def test_section_term_radial_distances():
 def test_principal_direction_extents():
     m = nm.load_morphology(SWC_PATH / 'simple.swc')
     principal_dir = features.get('principal_direction_extents', m)
-    assert_allclose(principal_dir, [14.736052694538641, 12.105102672688004])
+    assert_allclose(principal_dir, [10.99514 , 10.997688])
 
     # test with a realistic morphology
     m = nm.load_morphology(DATA_PATH / 'h5/v1' / 'bio_neuron-000.h5')
     p_ref = [
-        1672.969491,
-        142.437047,
-        224.607978,
-        415.50613,
-        429.830081,
-        165.954097,
-        346.832825,
+        1210.569727,
+        38.493958,
+        147.098687,
+        288.226628,
+        330.166506,
+        152.396521,
+        293.913857
     ]
     p = features.get('principal_direction_extents', m)
     assert_allclose(p, p_ref, rtol=1e-6)

tests/test_morphmath.py

@@ -29,6 +29,7 @@
 from math import fabs, pi, sqrt
 
 import numpy as np
+from numpy import testing as npt
 from neurom import morphmath as mm
 from neurom.core.dataformat import Point
 from numpy.random import uniform
@@ -532,3 +533,59 @@ def test_spherical_coordinates():
 
         new_elevation, new_azimuth = mm.spherical_from_vector(vect)
         assert np.allclose([elevation, azimuth], [new_elevation, new_azimuth])
+
+
+def test_principal_direction_extent():
+
+    # test with points on a circle with radius 0.5, and center at 0.0
+    circle_points = np.array([
+        [ 5.0e-01,  0.0e+00,  0.0e+00],
+        [ 4.7e-01,  1.6e-01,  0.0e+00],
+        [ 3.9e-01,  3.1e-01,  0.0e+00],
+        [ 2.7e-01,  4.2e-01,  0.0e+00],
+        [ 1.2e-01,  4.8e-01,  0.0e+00],
+        [-4.1e-02,  5.0e-01,  0.0e+00],
+        [-2.0e-01,  4.6e-01,  0.0e+00],
+        [-3.4e-01,  3.7e-01,  0.0e+00],
+        [-4.4e-01,  2.4e-01,  0.0e+00],
+        [-5.0e-01,  8.2e-02,  0.0e+00],
+        [-5.0e-01, -8.2e-02,  0.0e+00],
+        [-4.4e-01, -2.4e-01,  0.0e+00],
+        [-3.4e-01, -3.7e-01,  0.0e+00],
+        [-2.0e-01, -4.6e-01,  0.0e+00],
+        [-4.1e-02, -5.0e-01,  0.0e+00],
+        [ 1.2e-01, -4.8e-01,  0.0e+00],
+        [ 2.7e-01, -4.2e-01,  0.0e+00],
+        [ 3.9e-01, -3.1e-01,  0.0e+00],
+        [ 4.7e-01, -1.6e-01,  0.0e+00],
+        [ 5.0e-01, -1.2e-16,  0.0e+00]
+    ])
+
+    npt.assert_allclose(
+        mm.principal_direction_extent(circle_points),
+        [1., 1., 0.], atol=1e-6,
+    )
+
+    # extent should be invariant to translations
+    npt.assert_allclose(
+        mm.principal_direction_extent(circle_points + 100.),
+        [1., 1., 0.], atol=1e-6,
+    )
+    npt.assert_allclose(
+        mm.principal_direction_extent(circle_points - 100.),
+        [1., 1., 0.], atol=1e-6,
+    )
+
+    cross_3D_points = np.array([
+        [-5.2, 0.0, 0.0],
+        [ 4.8, 0.0, 0.0],
+        [ 0.0,-1.3, 0.0],
+        [ 0.0, 4.7, 0.0],
+        [ 0.0, 0.0,-11.2],
+        [ 0.0, 0.0, 0.8],
+    ])
+
+    npt.assert_allclose(
+        sorted(mm.principal_direction_extent(cross_3D_points)),
+        [6.0, 10.0, 12.0], atol=0.1,
+    )