Fix: Convert apical_sections to new id (#60)

neurots/generate/grower.py

@@ -39,6 +39,7 @@
 from neurots.preprocess import preprocess_inputs
 from neurots.utils import NeuroTSError
 from neurots.utils import convert_from_legacy_neurite_type
+from neurots.utils import point_to_section_segment
 
 L = logging.getLogger(__name__)
 
@@ -169,8 +170,20 @@ def grow(self):
     def _post_grow(self):
         """Actions after the morphology has been grown and before its diametrization."""
         # ensures section.id are consistent with morphio loader
+        apical_points = [
+            self.neuron.sections[apical_section].points[-1] if apical_section is not None else None
+            for apical_section in self.apical_sections
+        ]
+
         self.neuron = Morphology(self.neuron)
 
+        self.apical_sections = [
+            point_to_section_segment(self.neuron, apical_point)[0]
+            if apical_point is not None
+            else None
+            for apical_point in apical_points
+        ]
+
     def _init_diametrizer(self, external_diametrizer=None):
         """Set a diametrizer function."""
         if (

neurots/utils.py

@@ -19,6 +19,7 @@
 from copy import deepcopy
 
 import numpy as np
+from neurom import COLS
 
 PIA_DIRECTION = [0.0, 1.0, 0.0]
 
@@ -80,3 +81,28 @@ def convert_from_legacy_neurite_type(data):
                     data[key][i] = "basal_dendrite"
 
     return data
+
+
+def point_to_section_segment(neuron, point, rtol=1e-05, atol=1e-08):
+    """Find section and segment that matches the point (also in morph_tool.spatial).
+
+    Only the first point found with the *exact* same coordinates as the point argument is considered
+
+    Args:
+        neuron (morphio.Morphology): neuron object
+        point (point): value of the point to find in the h5 file
+        rtol, atol (floats): precision of np.isclose
+
+    Returns:
+        Tuple: (NeuroM/MorphIO section ID, point ID) of the point the matches the input coordinates.
+        Since NeuroM v2, section ids of NeuroM and MorphIO are the same excluding soma.
+    """
+    for section in neuron.iter():
+        points = section.points
+        offset = np.where(
+            np.isclose(points[:, COLS.XYZ], point[COLS.XYZ], rtol=rtol, atol=atol).all(axis=1)
+        )
+        if offset[0].size:
+            return section.id, offset[0][0]
+
+    raise ValueError(f"Cannot find point in morphology that matches: {point}")

tests/test_neuron_functional.py

@@ -288,7 +288,7 @@ def test_breaker_of_tmd_algo():
     N = NeuronGrower(input_distributions=distributions, input_parameters=params)
     n = N.grow()
 
-    assert_array_equal(N.apical_sections, [33])
+    assert_array_equal(N.apical_sections, [11])
     assert_array_almost_equal(
         n.sections[118].points[-1],
         np.array([-220.93813, -21.49141, -55.93323]),
@@ -304,7 +304,7 @@ def test_breaker_of_tmd_algo():
     N = NeuronGrower(input_distributions=distributions, input_parameters=params, rng_or_seed=rng)
     n = N.grow()
 
-    assert_array_equal(N.apical_sections, [33])
+    assert_array_equal(N.apical_sections, [11])
     assert_array_almost_equal(
         n.sections[118].points[-1],
         np.array([-220.93813, -21.49141, -55.93323]),

tests/test_soma_grower.py

@@ -482,7 +482,7 @@ def test_apical_points():
     grower.grow()
 
     apicals = grower.apical_sections
-    expected = [60]
+    expected = [47]
     assert_array_equal(apicals, expected)
 
     # Apical point not found so keep the last bifurcation
@@ -497,7 +497,7 @@ def test_apical_points():
     grower.grow()
 
     apicals = grower.apical_sections
-    expected = [11]
+    expected = [23]
     assert_array_equal(apicals, expected)
 
 

tests/test_utils.py

@@ -23,6 +23,7 @@
 import dictdiffer  # pylint: disable=import-error
 import numpy as np
 import pytest
+from morphio import Morphology
 
 from neurots import utils
 
@@ -91,3 +92,24 @@ def test_convert_from_legacy_neurite_type():
 
     data_converted = utils.convert_from_legacy_neurite_type(data_legacy)
     assert data_converted == data
+
+
+def test_point_to_section_segment():
+    neuron = Morphology(DATA / "dummy_neuron.asc")
+
+    section, segment = utils.point_to_section_segment(neuron, [0.0, 15.279950, 0.0])
+    assert section == 63
+    assert segment == 0
+
+    with pytest.raises(ValueError):
+        utils.point_to_section_segment(neuron, [24, 0, 0])
+
+    section, segment = utils.point_to_section_segment(
+        neuron, [0.0, 15.2, 0.0], rtol=1e-1, atol=1e-1
+    )
+    assert section == 63
+    assert segment == 0
+
+    section, segment = utils.point_to_section_segment(neuron, [0.0, 15.28, 0.0])
+    assert section == 63
+    assert segment == 0