Repository: BlueBrain/NeuroM

Issue:
Raise more informative error when a morph has no soma
Some morphology files have no soma, which NeuroM is able to load. But when one compute features on such morphology, one may get obscure errors. For example:
```python
nm.get("trunk_angles", neuron, neurite_type=nm.core.types.NeuriteType.basal_dendrite)
```
will raise the following error on a morphology without soma:
```python
File neurom/features/__init__.py:150, in get(feature_name, obj, **kwargs)
    135 def get(feature_name, obj, **kwargs):
    136     """Obtain a feature from a set of morphology objects.
    137 
    138     Features can be either Neurite, Morphology or Population features. For Neurite features see
   (...)
    148         List|Number: feature value as a list or a single number.
    149     """
--> 150     return _get_feature_value_and_func(feature_name, obj, **kwargs)[0]

File neurom/features/__init__.py:110, in _get_feature_value_and_func(feature_name, obj, **kwargs)
    108 if feature_name in _MORPHOLOGY_FEATURES:
    109     feature_ = _MORPHOLOGY_FEATURES[feature_name]
--> 110     res = feature_(obj, **kwargs)
    111 elif feature_name in _NEURITE_FEATURES:
    112     feature_ = _NEURITE_FEATURES[feature_name]

File neurom/features/morphology.py:202, in trunk_angles(morph, neurite_type, coords_only, sort_along, consecutive_only)
    172 @feature(shape=(...,))
    173 def trunk_angles(
    174     morph,
   (...)
    178     consecutive_only=True,
    179 ):
    180     """Calculate the angles between all the trunks of the morph.
    181 
    182     By default, the angles are defined on the x-y plane and the trees are sorted from the y axis
   (...)
    200             only the angle with the next trunk is returned for each trunk.
    201     """
--> 202     vectors = np.array(trunk_vectors(morph, neurite_type=neurite_type))
    203     # In order to avoid the failure of the process in case the neurite_type does not exist
    204     if len(vectors) == 0:

File neurom/features/morphology.py:169, in trunk_vectors(morph, neurite_type)
    166 def vector_to_root_node(neurite):
    167     return morphmath.vector(neurite.root_node.points[0], morph.soma.center)
--> 169 return _map_neurites(vector_to_root_node, morph, neurite_type)

File neurom/features/morphology.py:67, in _map_neurites(function, morph, neurite_type)
     66 def _map_neurites(function, morph, neurite_type):
---> 67     return list(
     68         iter_neurites(morph, mapfun=function, filt=is_type(neurite_type))
     69     )

File neurom/features/morphology.py:167, in trunk_vectors.<locals>.vector_to_root_node(neurite)
    166 def vector_to_root_node(neurite):
--> 167     return morphmath.vector(neurite.root_node.points[0], morph.soma.center)

File neurom/morphmath.py:60, in vector(p1, p2)
     50 def vector(p1, p2):
     51     """Compute vector between two 3D points.
     52 
     53     Args:
   (...)
     58         3-vector from p1 - p2
     59     """
---> 60     return np.subtract(p1[COLS.XYZ], p2[COLS.XYZ])

TypeError: 'NoneType' object is not subscriptable
```
which is not really helpful for users.
I suggest we add some checkers to raise more understandable errors.

Discussion:

Patch:
diff --git a/neurom/features/morphology.py b/neurom/features/morphology.py
index 69098cb..4b2c4de 100644
--- a/neurom/features/morphology.py
+++ b/neurom/features/morphology.py
@@ -63,6 +63,13 @@ from neurom.morphmath import convex_hull
 feature = partial(feature, namespace=NameSpace.NEURON)
 
 
+def _assert_soma_center(morph):
+    if morph.soma.center is None:
+        raise NeuroMError(
+            f"The morphology named '{morph.name}' has no soma so the feature can not be computed."
+        )
+
+
 def _map_neurites(function, morph, neurite_type):
     return list(
         iter_neurites(morph, mapfun=function, filt=is_type(neurite_type))
@@ -132,6 +139,8 @@ def trunk_origin_azimuths(morph, neurite_type=NeuriteType.all):
 
     The range of the azimuth angle [-pi, pi] radians
     """
+    _assert_soma_center(morph)
+
     def azimuth(neurite):
         """Azimuth of a neurite trunk."""
         return morphmath.azimuth_from_vector(
@@ -151,6 +160,8 @@ def trunk_origin_elevations(morph, neurite_type=NeuriteType.all):
 
     The range of the elevation angle [-pi/2, pi/2] radians
     """
+    _assert_soma_center(morph)
+
     def elevation(neurite):
         """Elevation of a section."""
         return morphmath.elevation_from_vector(
@@ -163,6 +174,8 @@ def trunk_origin_elevations(morph, neurite_type=NeuriteType.all):
 @feature(shape=(...,))
 def trunk_vectors(morph, neurite_type=NeuriteType.all):
     """Calculate the vectors between all the trunks of the morphology and the soma center."""
+    _assert_soma_center(morph)
+
     def vector_to_root_node(neurite):
         return morphmath.vector(neurite.root_node.points[0], morph.soma.center)
 
@@ -496,6 +509,7 @@ def sholl_crossings(morph, neurite_type=NeuriteType.all, center=None, radii=None
             '`sholl_crossings` input error. If `center` or `radii` is not set then `morph` is ' \
             'expected to be an instance of Morphology and have a soma.'
         if center is None:
+            _assert_soma_center(morph)
             center = morph.soma.center
         if radii is None:
             radii = [morph.soma.radius]
@@ -525,6 +539,7 @@ def sholl_frequency(morph, neurite_type=NeuriteType.all, step_size=10, bins=None
         If a `neurite_type` is specified and there are no trees corresponding to it, an empty
         list will be returned.
     """
+    _assert_soma_center(morph)
     neurite_filter = is_type(neurite_type)
 
     if bins is None:
@@ -701,6 +716,7 @@ def length_fraction_above_soma(morph, neurite_type=NeuriteType.all, up="Y"):
     Returns:
         The fraction of neurite length that lies on the right of the soma along the given axis.
     """
+    _assert_soma_center(morph)
     axis = up.upper()
 
     if axis not in {"X", "Y", "Z"}:
diff --git a/neurom/features/population.py b/neurom/features/population.py
index 98cbd4d..dfd0a7e 100644
--- a/neurom/features/population.py
+++ b/neurom/features/population.py
@@ -47,6 +47,7 @@ from neurom.core.dataformat import COLS
 from neurom.core.types import NeuriteType
 from neurom.core.types import tree_type_checker as is_type
 from neurom.features import feature, NameSpace
+from neurom.features.morphology import _assert_soma_center
 from neurom.features.morphology import sholl_crossings
 
 feature = partial(feature, namespace=NameSpace.POPULATION)
@@ -78,7 +79,8 @@ def sholl_frequency(morphs, neurite_type=NeuriteType.all, step_size=10, bins=Non
                         for n in m.neurites if neurite_filter(n))
         bins = np.arange(min_soma_edge, min_soma_edge + max_radii, step_size)
 
-    return np.array([
-        sholl_crossings(m, neurite_type, m.soma.center, bins)
-        for m in morphs
-    ]).sum(axis=0)
+    def _sholl_crossings(morph):
+        _assert_soma_center(morph)
+        return sholl_crossings(morph, neurite_type, morph.soma.center, bins)
+
+    return np.array([_sholl_crossings(m) for m in morphs]).sum(axis=0)


Test patch:
diff --git a/tests/features/test_morphology.py b/tests/features/test_morphology.py
index 2337c83..3bd49d5 100644
--- a/tests/features/test_morphology.py
+++ b/tests/features/test_morphology.py
@@ -33,6 +33,7 @@ import warnings
 from io import StringIO
 from pathlib import Path
 
+import morphio
 import numpy as np
 import pytest
 from morphio import PointLevel, SectionType
@@ -44,9 +45,9 @@ from numpy.testing import assert_array_equal
 
 from neurom import morphmath
 from neurom import NeuriteType, load_morphology, AXON, BASAL_DENDRITE
-from neurom.core import Morphology
+from neurom.core import Morphology, Population
 from neurom.exceptions import NeuroMError
-from neurom.features import morphology, section
+from neurom.features import morphology, population, section
 
 
 DATA_PATH = Path(__file__).parent.parent / 'data'
@@ -725,3 +726,26 @@ def test_unique_projected_points():
         morphology._unique_projected_points(morph, "airplane", NeuriteType.all)
 
     assert len(morphology._unique_projected_points(morph, "yz", NeuriteType.apical_dendrite)) == 0
+
+
+def test_missing_soma():
+    NRN_missing_soma = load_morphology(SWC_PATH / 'Single_apical_no_soma.swc')
+
+    with pytest.raises(NeuroMError):
+        morphology.trunk_origin_elevations(NRN_missing_soma)
+    with pytest.raises(NeuroMError):
+        morphology.trunk_origin_azimuths(NRN_missing_soma)
+    with pytest.raises(NeuroMError):
+        morphology.trunk_origin_elevations(NRN_missing_soma)
+    with pytest.raises(NeuroMError):
+        morphology.trunk_vectors(NRN_missing_soma)
+    with pytest.raises(NeuroMError):
+        morphology.sholl_crossings(NRN_missing_soma)
+    with pytest.raises(NeuroMError):
+        morphology.sholl_frequency(NRN_missing_soma)
+    with pytest.raises(NeuroMError):
+        morphology.length_fraction_above_soma(NRN_missing_soma)
+
+    POP_missing_soma = Population([NRN_missing_soma])
+    with pytest.raises(NeuroMError):
+        population.sholl_frequency(POP_missing_soma)