Repository: BlueBrain/NeuroM

Issue:
`sholl_frequency` does not use soma center to calculate max radius
The soma center is not taken into account when calculating the distance to the points to find the max radius (L481):
https://github.com/BlueBrain/NeuroM/blob/e31eb63b522a3558233d821fcfbb178e6754aa56/neurom/features/morphology.py#L479-L485

Of course, for morphologies that are centered at (0, 0, 0) this does not make a difference, but it will mess up the rest.

Discussion:
mgeplf: Nice catch
mgeplf: Would it be worth adding a test where the soma isn't at the origin?
eleftherioszisis: > Would it be worth adding a test where the soma isn't at the origin?

Yes, good point. I added one.

Patch:
diff --git a/neurom/features/morphology.py b/neurom/features/morphology.py
index 2315e12..67583e9 100644
--- a/neurom/features/morphology.py
+++ b/neurom/features/morphology.py
@@ -478,8 +478,10 @@ def sholl_frequency(morph, neurite_type=NeuriteType.all, step_size=10, bins=None
 
     if bins is None:
         min_soma_edge = morph.soma.radius
-        max_radii = max(np.max(np.linalg.norm(n.points[:, COLS.XYZ], axis=1))
-                        for n in morph.neurites if neurite_filter(n))
+        max_radii = max(
+            np.max(np.linalg.norm(n.points[:, COLS.XYZ] - morph.soma.center, axis=1))
+            for n in morph.neurites if neurite_filter(n)
+        )
         bins = np.arange(min_soma_edge, min_soma_edge + max_radii, step_size)
 
     return sholl_crossings(morph, neurite_type, morph.soma.center, bins)


Test patch:
diff --git a/tests/features/test_get_features.py b/tests/features/test_get_features.py
index 6f20106..312edaf 100644
--- a/tests/features/test_get_features.py
+++ b/tests/features/test_get_features.py
@@ -706,6 +706,17 @@ def test_sholl_frequency():
     assert len(features.get('sholl_frequency', POP)) == 108
 
 
+    # check that the soma is taken into account for calculating max radius and num bins
+    m = nm.load_morphology(
+        """
+        1  1  -10  0  0    5.0 -1
+        2  3    0  0  0    0.1  1
+        3  3   10  0  0    0.1  2
+        """, reader="swc",
+    )
+
+    assert features.get('sholl_frequency', m, step_size=5.0) == [0, 1, 1, 1]
+
 def test_bifurcation_partitions():
     assert_allclose(features.get('bifurcation_partitions', POP)[:10],
                     [19., 17., 15., 13., 11., 9., 7., 5., 3., 1.])