Morphology level radial distance features use the soma as reference p…

…oint (#1030)

Use soma center as a reference point for morphology-level radial distance features.

Radial distance features on the morphology level were using
as reference point the root of each neurite. While this is ok
when per-neurite features are calculated, it is incorrect when
the radial distance features are calculated for the entire cell.

This change renders soma as the reference point for calculating
the radial distances, not the root of each neurite. When the same
features are used per-neurite the old behavior of using the neurite
root is still the same.

Co-authored-by: Adrien Berchet <[email protected]>

neurom/features/morphology.py

@@ -59,7 +59,7 @@
 from neurom.core.types import NeuriteType
 from neurom.exceptions import NeuroMError
 from neurom.features import feature, NameSpace, neurite as nf, section as sf
-from neurom.utils import str_to_plane
+from neurom.utils import str_to_plane, flatten
 from neurom.morphmath import convex_hull
 
 
@@ -115,15 +115,77 @@ def soma_radius(morph):
 @feature(shape=())
 def max_radial_distance(morph, origin=None, neurite_type=NeuriteType.all, use_subtrees=False):
     """Get the maximum radial distances of the termination sections."""
+    origin = morph.soma.center if origin is None else origin
+
     term_radial_distances = _map_neurites(
         partial(nf.max_radial_distance, origin=origin),
-        morph,
-        neurite_type=neurite_type,
-        use_subtrees=use_subtrees,
+        morph, neurite_type, use_subtrees
     )
     return max(term_radial_distances) if term_radial_distances else 0.
 
 
+@feature(shape=(...,))
+def section_radial_distances(morph, origin=None, neurite_type=NeuriteType.all, use_subtrees=False):
+    """Section radial distances.
+
+    The iterator_type can be used to select only terminal sections (ileaf)
+    or only bifurcations (ibifurcation_point).
+    """
+    origin = morph.soma.center if origin is None else origin
+
+    return list(flatten(_map_neurites(
+        partial(nf.section_radial_distances, origin=origin),
+        morph=morph,
+        neurite_type=neurite_type,
+        use_subtrees=use_subtrees,
+    )))
+
+
+@feature(shape=(...,))
+def section_term_radial_distances(
+    morph, origin=None, neurite_type=NeuriteType.all, use_subtrees=False
+):
+    """Get the radial distances of the termination sections."""
+    origin = morph.soma.center if origin is None else origin
+
+    return list(flatten(_map_neurites(
+        partial(nf.section_term_radial_distances, origin=origin),
+        morph=morph,
+        neurite_type=neurite_type,
+        use_subtrees=use_subtrees
+    )))
+
+
+@feature(shape=(...,))
+def section_bif_radial_distances(
+    morph, origin=None, neurite_type=NeuriteType.all, use_subtrees=False
+):
+    """Get the radial distances of the bifurcation sections."""
+    origin = morph.soma.center if origin is None else origin
+
+    return list(flatten(_map_neurites(
+        partial(nf.section_bif_radial_distances, origin=origin),
+        morph=morph,
+        neurite_type=neurite_type,
+        use_subtrees=use_subtrees,
+    )))
+
+
+@feature(shape=(...,))
+def segment_radial_distances(
+    morph, origin=None, neurite_type=NeuriteType.all, use_subtrees=False
+):
+    """Ger the radial distances of the segments."""
+    origin = morph.soma.center if origin is None else origin
+
+    return list(flatten(_map_neurites(
+        partial(nf.segment_radial_distances, origin=origin),
+        morph=morph,
+        neurite_type=neurite_type,
+        use_subtrees=use_subtrees,
+    )))
+
+
 @feature(shape=(...,))
 def number_of_sections_per_neurite(morph, neurite_type=NeuriteType.all, use_subtrees=False):
     """List of numbers of sections per neurite."""

neurom/features/neurite.py

@@ -81,15 +81,6 @@ def homogeneous_filter(section):
     return list(map(fun, filter(filt, iterator_type(neurite.root_node))))
 
 
-@feature(shape=())
-def max_radial_distance(neurite, origin=None, section_type=NeuriteType.all):
-    """Get the maximum radial distances of the termination sections."""
-    term_radial_distances = section_term_radial_distances(
-        neurite, origin=origin, section_type=section_type
-    )
-    return max(term_radial_distances) if term_radial_distances else 0.
-
-
 @feature(shape=())
 def number_of_segments(neurite, section_type=NeuriteType.all):
     """Number of segments."""
@@ -423,34 +414,48 @@ def diameter_power_relations(neurite, method='first', section_type=NeuriteType.a
     )
 
 
+def _radial_distances(neurite, origin, iterator_type, section_type):
+
+    if origin is None:
+        origin = neurite.root_node.points[0]
+
+    return _map_sections(
+        partial(sf.section_radial_distance, origin=origin),
+        neurite=neurite,
+        iterator_type=iterator_type,
+        section_type=section_type
+    )
+
+
 @feature(shape=(...,))
-def section_radial_distances(
-    neurite, origin=None, iterator_type=Section.ipreorder, section_type=NeuriteType.all
-):
+def section_radial_distances(neurite, origin=None, section_type=NeuriteType.all):
     """Section radial distances.
 
     The iterator_type can be used to select only terminal sections (ileaf)
     or only bifurcations (ibifurcation_point).
     """
-    pos = neurite.root_node.points[0] if origin is None else origin
-    return _map_sections(partial(sf.section_radial_distance, origin=pos),
-                         neurite,
-                         iterator_type,
-                         section_type=section_type)
+    return _radial_distances(neurite, origin, Section.ipreorder, section_type)
 
 
 @feature(shape=(...,))
 def section_term_radial_distances(neurite, origin=None, section_type=NeuriteType.all):
     """Get the radial distances of the termination sections."""
-    return section_radial_distances(neurite, origin, Section.ileaf, section_type=section_type)
+    return _radial_distances(neurite, origin, Section.ileaf, section_type)
+
+
+@feature(shape=())
+def max_radial_distance(neurite, origin=None, section_type=NeuriteType.all):
+    """Get the maximum radial distances of the termination sections."""
+    term_radial_distances = section_term_radial_distances(
+        neurite, origin=origin, section_type=section_type
+    )
+    return max(term_radial_distances) if term_radial_distances else 0.
 
 
 @feature(shape=(...,))
 def section_bif_radial_distances(neurite, origin=None, section_type=NeuriteType.all):
     """Get the radial distances of the bf sections."""
-    return section_radial_distances(
-        neurite, origin, Section.ibifurcation_point, section_type=section_type
-    )
+    return _radial_distances(neurite, origin, Section.ibifurcation_point, section_type)
 
 
 @feature(shape=(...,))

tests/features/test_get_features.py

@@ -112,23 +112,23 @@ def test_number_of_sections():
 def test_max_radial_distance():
     assert_allclose(
         features.get('max_radial_distance', POP),
-        [99.58945832, 94.43342439, 1053.77939245])
+        [99.62086,   94.43019, 1072.9137])
     assert_allclose(
         features.get('max_radial_distance', POP, neurite_type=NeuriteType.all),
-        [99.58945832, 94.43342439, 1053.77939245])
+        [99.62086,   94.43019, 1072.9137])
     assert_allclose(
         features.get('max_radial_distance', POP, neurite_type=NeuriteType.axon),
-        [82.442545, 82.442545, 1053.779392])
+        [82.52528,   82.44438, 1072.9137])
     assert_allclose(
         features.get('max_radial_distance', POP, neurite_type=NeuriteType.basal_dendrite),
-        [94.43342563, 94.43342439, 207.56977859])
+        [94.36033,  94.43019, 209.92587])
 
     assert_allclose(
-        features.get('max_radial_distance', NRN), 99.58945832)
+        features.get('max_radial_distance', NRN), 99.62086)
     assert_allclose(
-        features.get('max_radial_distance', NRN, neurite_type=NeuriteType.all), 99.58945832)
+        features.get('max_radial_distance', NRN, neurite_type=NeuriteType.all), 99.62086)
     assert_allclose(features.get(
-        'max_radial_distance', NRN, neurite_type=NeuriteType.apical_dendrite), 99.589458)
+        'max_radial_distance', NRN, neurite_type=NeuriteType.apical_dendrite), 99.62086)
 
     assert_allclose(
         features.get('max_radial_distance', NRN.neurites),
@@ -771,7 +771,19 @@ def test_section_strahler_orders():
 
 
 def test_section_bif_radial_distances():
+
+    # the feature applied on morph calculates radial distance from soma
     trm_rads = features.get('section_bif_radial_distances', NRN, neurite_type=nm.AXON)
+
+    assert_allclose(
+        trm_rads,
+        [8.92228 , 16.825268, 23.152378, 30.262894, 36.71048 ,
+        44.049297, 52.00228 , 59.510105, 66.33529 , 74.134636]
+    )
+
+    # the feature applied per neurite calculates radial distance from root
+    trm_rads = features.get('section_bif_radial_distances', NRN.neurites[3])
+
     assert_allclose(trm_rads,
                     [8.842008561870646,
                      16.7440421479104,
@@ -786,7 +798,25 @@ def test_section_bif_radial_distances():
 
 
 def test_section_term_radial_distances():
+
     trm_rads = features.get('section_term_radial_distances', NRN, neurite_type=nm.APICAL_DENDRITE)
+
+    print(trm_rads)
+    assert_allclose(trm_rads,
+                [16.258472,
+                 26.040075,
+                 33.35425 ,
+                 42.755745,
+                 52.41365 ,
+                 59.476284,
+                 67.11225 ,
+                 80.00984 ,
+                 87.13672 ,
+                 97.284706,
+                 99.62086 ])
+
+    apical = NRN.neurites[0]
+    trm_rads = features.get('section_term_radial_distances', apical, section_type=nm.APICAL_DENDRITE)
     assert_allclose(trm_rads,
                     [16.22099879395879,
                      25.992977561564082,