Repository: BlueBrain/NeuroM

Issue:
Bug in principal direction extent
There is a bug in the calculation of the principal direction extents:

https://github.com/BlueBrain/NeuroM/blob/fdf9f78be79508e34c0349936e233c184279c955/neurom/morphmath.py#L462-L496

At the point where the most negative value is subtracted from the entire extent:
```python
         if scalar_projs[0] < 0.: 
             extent -= scalar_projs[0] 
```
It shouldn't subtract from the entire extent, rather from the i-th component `extent[i]`.

Failing test example with points on a circle with radius `0.5`, and center at `0.0`:
```python
from neurom.morphmath import principal_direction_extent
from numpy import testing as npt

def test_principal_direction_extent():

    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],
        [-4.9e-01,  8.2e-02,  0.0e+00],
        [-4.9e-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(
        principal_direction_extent(circle_points),
        [1., 1., 0.],
    )
```
Which returns instead `[1.49, 1.0, 0.0]`.

Thanks, @lidakanari for finding that there is something wrong with the principal direction extent function.

@arnaudon @adrien-berchet 

Discussion:
eleftherioszisis: @lidakanari , could you please check that with these changes you get the expected results?
codecov-commenter: # [Codecov](https://codecov.io/gh/BlueBrain/NeuroM/pull/1008?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=BlueBrain) Report
> Merging [#1008](https://codecov.io/gh/BlueBrain/NeuroM/pull/1008?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=BlueBrain) (b13e977) into [master](https://codecov.io/gh/BlueBrain/NeuroM/commit/fdf9f78be79508e34c0349936e233c184279c955?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=BlueBrain) (fdf9f78) will **not change** coverage.
> The diff coverage is `100.00%`.

> :exclamation: Current head b13e977 differs from pull request most recent head def17d2. Consider uploading reports for the commit def17d2 to get more accurate results



```diff
@@            Coverage Diff            @@
##            master     #1008   +/-   ##
=========================================
  Coverage   100.00%   100.00%           
=========================================
  Files           36        36           
  Lines         2355      2349    -6     
=========================================
- Hits          2355      2349    -6     
```

eleftherioszisis: > The function seems to return correct values for the tested cases. However, the signature of what it returns changes with respect to the previous one, as it now returns a vector of the same dimension as the input data. i.e. 2d point cloud returned 3d vector before, now it returns 2d.

Indeed, the dimension of the extents depends on the dimension of the points. Therefore, for 2D points it would return a 2D vector of extents, and for 3D a 3D respectively. For NeuroM purposes, only 3D points are used, so this is not really an issue.
eleftherioszisis: > Ok, that should be good. Also something to keep in mind, currently the extents are **not** sorted (i.e. largest to smallest). If this is by choice, it would be good to mention in documentation because it might be confusing.

That's true indeed. There was no requirement that they are sorted from max to min, which imho renders the `direction` argument a bit useless in the respective feature. 

Should I order them?

Patch:
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 540b7e2..80313f9 100644
--- a/CHANGELOG.rst
+++ b/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
diff --git a/neurom/morphmath.py b/neurom/morphmath.py
index b01c516..41e090d 100644
--- a/neurom/morphmath.py
+++ b/neurom/morphmath.py
@@ -462,17 +462,17 @@ section_length = path_distance
 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)


Test patch:
diff --git a/tests/features/test_get_features.py b/tests/features/test_get_features.py
index 6661b89..a78efa7 100644
--- a/tests/features/test_get_features.py
+++ b/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)
diff --git a/tests/test_morphmath.py b/tests/test_morphmath.py
index 0188ccf..6119f7b 100644
--- a/tests/test_morphmath.py
+++ b/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,
+    )