Add path length option to sholl analysis

neurom/apps/cli.py

@@ -76,6 +76,8 @@ def view(input_file, is_3d, plane, backend, realistic_diameters):
 
     plot(load_morphology(input_file))
     if is_matplotlib:
+        if not is_3d:
+            plt.axis('equal')
         plt.show()
 
 

neurom/check/morphtree.py

@@ -53,7 +53,7 @@ def is_monotonic(neurite, tol):
             if sec[point_id + 1][COLS.R] > sec[point_id][COLS.R] + tol:
                 return False
         # Check that section boundary points satisfy monotonicity
-        if(node.parent is not None and
+        if (node.parent is not None and
            sec[0][COLS.R] > node.parent.points[-1][COLS.R] + tol):
             return False
 

neurom/features/morphology.py

@@ -50,10 +50,11 @@
 import numpy as np
 
 from neurom import morphmath
-from neurom.core.morphology import iter_neurites, iter_segments, Morphology
+from neurom.core.morphology import iter_neurites, iter_segments, Morphology, iter_sections
 from neurom.core.types import tree_type_checker as is_type
 from neurom.core.dataformat import COLS
 from neurom.core.types import NeuriteType
+from neurom.core.morphology import Section
 from neurom.exceptions import NeuroMError
 from neurom.features import feature, NameSpace, neurite as nf, section as sf
 from neurom.utils import str_to_plane
@@ -457,7 +458,9 @@ def neurite_volume_density(morph, neurite_type=NeuriteType.all):
 
 
 @feature(shape=(...,))
-def sholl_crossings(morph, neurite_type=NeuriteType.all, center=None, radii=None):
+def sholl_crossings(
+    morph, neurite_type=NeuriteType.all, center=None, radii=None, distance_type="radial"
+):
     """Calculate crossings of neurites.
 
     Args:
@@ -466,6 +469,7 @@ def sholl_crossings(morph, neurite_type=NeuriteType.all, center=None, radii=None
         center(Point): center point, if None then soma center is taken
         radii(iterable of floats): radii for which crossings will be counted,
             if None then soma radius is taken
+        distance_type(str): either `radial` or `path` as distance measure to use
 
     Returns:
         Array of same length as radii, with a count of the number of crossings
@@ -480,14 +484,25 @@ def sholl_crossings(morph, neurite_type=NeuriteType.all, center=None, radii=None
     """
     def _count_crossings(neurite, radius):
         """Used to count_crossings of segments in neurite with radius."""
-        r2 = radius ** 2
         count = 0
-        for start, end in iter_segments(neurite):
-            start_dist2, end_dist2 = (morphmath.point_dist2(center, start),
-                                      morphmath.point_dist2(center, end))
 
-            count += int(start_dist2 <= r2 <= end_dist2 or
-                         end_dist2 <= r2 <= start_dist2)
+        if distance_type == 'radial':
+            r2 = radius ** 2
+            for start, end in iter_segments(neurite):
+                start_dist2, end_dist2 = (morphmath.point_dist2(center, start),
+                                          morphmath.point_dist2(center, end))
+
+                count += int(start_dist2 <= r2 <= end_dist2 or
+                             end_dist2 <= r2 <= start_dist2)
+
+        if distance_type == 'path':
+            for section in iter_sections(neurite):
+                base_length = sf.section_path_length(section) - section.length
+                path_dists = base_length + sf.segment_lengths(section, prepend_zero=True)
+
+                forward = np.logical_and(path_dists[:-1] <= radius, path_dists[1:] >= radius)
+                backward = np.logical_and(path_dists[:-1] >= radius, path_dists[1:] <= radius)
+                count += sum(np.logical_or(forward, backward))
 
         return count
 
@@ -505,7 +520,9 @@ def _count_crossings(neurite, radius):
 
 
 @feature(shape=(...,))
-def sholl_frequency(morph, neurite_type=NeuriteType.all, step_size=10, bins=None):
+def sholl_frequency(
+    morph, neurite_type=NeuriteType.all, step_size=10, bins=None, distance_type='euclidean'
+):
     """Perform Sholl frequency calculations on a morph.
 
     Args:
@@ -514,6 +531,7 @@ def sholl_frequency(morph, neurite_type=NeuriteType.all, step_size=10, bins=None
         step_size(float): step size between Sholl radii
         bins(iterable of floats): custom binning to use for the Sholl radii. If None, it uses
         intervals of step_size between min and max radii of ``morphologies``.
+        distance_type(str): either `euclidean` or `path` as distance measure to use
 
     Note:
         Given a morphology, the soma center is used for the concentric circles,
@@ -528,19 +546,34 @@ def sholl_frequency(morph, neurite_type=NeuriteType.all, step_size=10, bins=None
     neurite_filter = is_type(neurite_type)
 
     if bins is None:
-        min_soma_edge = morph.soma.radius
 
-        max_radius_per_neurite = [
-            np.max(np.linalg.norm(n.points[:, COLS.XYZ] - morph.soma.center, axis=1))
-            for n in morph.neurites if neurite_filter(n)
-        ]
-
-        if not max_radius_per_neurite:
-            return []
-
-        bins = np.arange(min_soma_edge, min_soma_edge + max(max_radius_per_neurite), step_size)
-
-    return sholl_crossings(morph, neurite_type, morph.soma.center, bins)
+        if distance_type == "euclidean":
+            max_radius_per_neurite = [
+                np.max(np.linalg.norm(n.points[:, COLS.XYZ] - morph.soma.center, axis=1))
+                for n in morph.neurites
+                if neurite_filter(n)
+            ]
+            if not max_radius_per_neurite:
+                return []
+            min_distance = morph.soma.radius
+            max_distance = min_distance + max(max_radius_per_neurite)
+
+        if distance_type == "path":
+            max_path_distance_per_neurite = [
+                max(sf.section_path_length(section) for section in iter_sections(n, Section.ileaf))
+                for n in morph.neurites
+                if neurite_filter(n)
+            ]
+            if not max_path_distance_per_neurite:
+                return []
+            min_distance = 0
+            max_distance = max(max_path_distance_per_neurite)
+
+        bins = np.arange(min_distance, max_distance, step_size)
+
+    return sholl_crossings(
+        morph, neurite_type, morph.soma.center, bins, distance_type=distance_type
+    )
 
 
 def _extent_along_axis(morph, axis, neurite_type):

setup.py

@@ -55,8 +55,6 @@
     extras_require={
         'plotly': ['plotly>=3.6.0', 'psutil>=5.5.1'],  # for plotly image saving
         'docs': [
-            'Jinja2<3.1',  # New release 3.1.0 breaks sphinx-bluebrain-theme
-            'sphinx',
             'sphinx-bluebrain-theme',
             'sphinx-autorun',
         ],

tests/check/test_morphology_checks.py

@@ -387,7 +387,7 @@ def test__bool__():
 
 def test_has_multifurcation():
     m = load_morphology(StringIO(u"""
-	((CellBody) (0 0 0 2))
+	((CellBody) (-1 0 0 2) (1 0 0 2))
 ( (Color Blue)
   (Axon)
   (0 5 0 2)
@@ -418,11 +418,7 @@ def test_has_multifurcation():
 
 def test_has_unifurcation():
     m = load_morphology(StringIO(u"""
-("CellBody"
- (Color Red)
- (CellBody)
- (0 0 0 2)
- )
+((CellBody) (-1 0 0 2) (1 0 0 2))
 
  ((Dendrite)
   (0 0 0 2)

tests/check/test_morphtree.py

@@ -41,8 +41,7 @@ def _generate_back_track_tree(n, dev):
     points = np.array(dev) + np.array([1, 3 if n == 0 else -3, 0])
 
     m = load_morphology(StringIO(u"""
-    ((CellBody)
-     (0 0 0 0.4))
+    ((CellBody) (-1 0 0 2) (1 0 0 2))
 
     ((Dendrite)
     (0 0 0 0.4)

tests/core/test_iter.py

@@ -200,8 +200,7 @@ def test_iter_segments_pop():
 
 def test_iter_segments_section():
     sec = load_morphology(StringIO(u"""
-	                      ((CellBody)
-	                       (0 0 0 2))
+	                      ((CellBody) (-1 0 0 2) (1 0 0 2))
 
 	                      ((Dendrite)
                           (1 2 3 8)

tests/features/test_get_features.py

@@ -402,7 +402,7 @@ def test_segment_meander_angles():
 
 
 def test_segment_meander_angles_single_section():
-    m = nm.load_morphology(StringIO(u"""((CellBody) (0 0 0 0))
+    m = nm.load_morphology(StringIO(u"""((CellBody) (-1 0 0 2) (1 0 0 2))
                                       ((Dendrite)
                                        (0 0 0 2)
                                        (1 0 0 2)
@@ -736,6 +736,7 @@ def test_sholl_frequency():
 
     # check that if there is no neurite of a specific type, an empty list is returned
     assert features.get('sholl_frequency', m, neurite_type=NeuriteType.axon) == []
+    assert features.get('sholl_frequency', m, neurite_type=NeuriteType.axon, distance_type='path') == []
 
 
 def test_bifurcation_partitions():

tests/features/test_morphology.py

@@ -506,6 +506,16 @@ def test_sholl_analysis_custom():
     assert (list(morphology.sholl_crossings(morph_A, center=center, radii=radii)) ==
             [2, 2, 2, 2, 2, 2, 2, 2])
 
+    assert list(morphology.sholl_frequency(morph_A)) == [2, 2, 2, 2, 2, 2, 2, 2]
+
+    assert (
+        list(morphology.sholl_frequency(morph_A, distance_type='path')) == [2, 2, 2, 2, 2, 2, 2, 2]
+    )
+
+    assert list(
+        morphology.sholl_crossings(morph_A, center=center, radii=radii, distance_type='path')
+    ) == [2, 2, 2, 2, 2, 2, 2, 2]
+
     morph_B = load_swc("""\
  1 1   0   0  0 1. -1
  2 3   0   0  0 1.  1
@@ -522,8 +532,16 @@ def test_sholl_analysis_custom():
 13 4 -51  -5  0 1.  9
 14 4 -51 -10  0 1.  9
                        """)
-    assert (list(morphology.sholl_crossings(morph_B, center=center, radii=radii)) ==
-            [2, 2, 2, 10, 10, 0, 0, 0])
+    assert (list(
+        morphology.sholl_crossings(morph_B, center=center, radii=radii)
+    ) == [2, 2, 2, 10, 10, 0, 0, 0])
+
+    assert (list(
+        morphology.sholl_crossings(morph_B, center=center, radii=radii, distance_type='path')
+    ) == [2, 2, 2, 10, 10, 0, 0, 0])
+
+    assert list(morphology.sholl_frequency(morph_B)) == [2, 2, 2, 2, 10, 10]
+    assert list(morphology.sholl_frequency(morph_B, distance_type='path')) == [2, 2, 2, 2, 10, 10]
 
     morph_C = load_swc("""\
  1 1   0   0  0 1. -1
@@ -544,6 +562,15 @@ def test_sholl_analysis_custom():
     assert (list(morphology.sholl_crossings(morph_C, center=center, radii=radii)) ==
             [2, 2, 2, 2, 2, 2, 10, 10])
 
+    assert (list(
+        morphology.sholl_crossings(morph_C, center=center, radii=radii, distance_type='path')
+    ) == [2, 2, 2, 2, 2, 2, 10, 10])
+
+    assert list(morphology.sholl_frequency(morph_C)) == [2, 2, 2, 2, 2, 2, 2, 10, 10]
+    assert (list(
+        morphology.sholl_frequency(morph_C, distance_type='path')
+    ) == [2, 2, 2, 2, 2, 2, 2, 10, 10])
+
 
 def test_extent_along_axis():
     morph = load_swc("""