Add method to Soma class to check wether it overlaps a point or not (#…

…956)

CHANGELOG.rst

@@ -3,6 +3,7 @@ Changelog
 
 Version 3.0.1
 -------------
+- Add method to Soma class to check wether it overlaps a point or not.
 - Ensure ``neurom.morphmath.angle_between_vectors`` always return 0 when the vectors are equal.
 
 Version 3.0.0

neurom/core/soma.py

@@ -84,6 +84,15 @@ def volume(self):
         warnings.warn('Approximating soma volume by a sphere. {}'.format(self))
         return 4. / 3 * math.pi * self.radius ** 3
 
+    def overlaps(self, points, exclude_boundary=False):
+        """Check that the given points are located inside the soma."""
+        points = np.atleast_2d(np.asarray(points, dtype=np.float64))
+        if exclude_boundary:
+            mask = np.linalg.norm(points - self.center, axis=1) < self.radius
+        else:
+            mask = np.linalg.norm(points - self.center, axis=1) <= self.radius
+        return mask
+
 
 class SomaSinglePoint(Soma):
     """Type A: 1point soma.
@@ -149,6 +158,29 @@ def __str__(self):
         return ('SomaCylinders(%s) <center: %s, virtual radius: %s>' %
                 (repr(self.points), self.center, self.radius))
 
+    def overlaps(self, points, exclude_boundary=False):
+        """Check that the given points are located inside the soma."""
+        points = np.atleast_2d(np.asarray(points, dtype=np.float64))
+        mask = np.ones(len(points)).astype(bool)
+        for p1, p2 in zip(self.points[:-1], self.points[1:]):
+            vec = p2[COLS.XYZ] - p1[COLS.XYZ]
+            vec_norm = np.linalg.norm(vec)
+            dot = (points[mask] - p1[COLS.XYZ]).dot(vec) / vec_norm
+
+            cross = np.linalg.norm(np.cross(vec, points[mask]), axis=1) / vec_norm
+            dot_clipped = np.clip(dot / vec_norm, a_min=0, a_max=1)
+            radii = p1[COLS.R] * (1 - dot_clipped) + p2[COLS.R] * dot_clipped
+
+            if exclude_boundary:
+                in_cylinder = (dot > 0) & (dot < vec_norm) & (cross < radii)
+            else:
+                in_cylinder = (dot >= 0) & (dot <= vec_norm) & (cross <= radii)
+            mask[np.where(mask)] = ~in_cylinder
+            if not mask.any():
+                break
+
+        return ~mask
+
 
 class SomaNeuromorphoThreePointCylinders(SomaCylinders):
     """NeuroMorpho compatible soma.
@@ -227,6 +259,56 @@ def __str__(self):
         return ('SomaSimpleContour(%s) <center: %s, radius: %s>' %
                 (repr(self.points), self.center, self.radius))
 
+    def overlaps(self, points, exclude_boundary=False):
+        """Check that the given points are located inside the soma.
+
+        The contour is supposed to be in the plane XY, the Z component is ignored.
+        """
+        # pylint: disable=too-many-locals
+        points = np.atleast_2d(np.asarray(points, dtype=np.float64))
+
+        # Convert points to angles from the center
+        relative_pts = points - self.center
+        pt_angles = np.arctan2(relative_pts[:, COLS.Y], relative_pts[:, COLS.X])
+
+        # Convert soma points to angles from the center
+        relative_soma_pts = self.points[:, COLS.XYZ] - self.center
+        soma_angles = np.arctan2(relative_soma_pts[:, COLS.Y], relative_soma_pts[:, COLS.X])
+
+        # Order the soma points by ascending angles
+        soma_angle_order = np.argsort(soma_angles)
+        ordered_soma_angles = soma_angles[soma_angle_order]
+        ordered_relative_soma_pts = relative_soma_pts[soma_angle_order]
+
+        # Find the two soma points which form the segment crossed by the one from the center
+        # to the point
+        angles = np.atleast_2d(pt_angles).T - ordered_soma_angles
+        closest_indices = np.argmin(np.abs(angles), axis=1)
+        neighbors = np.ones_like(closest_indices)
+        neighbors[angles[np.arange(len(closest_indices)), closest_indices] < 0] = -1
+        signs = (neighbors == 1) * 2. - 1.
+        neighbors[
+            (closest_indices >= len(relative_soma_pts) - 1)
+            & (neighbors == 1)
+        ] = -len(relative_soma_pts) + 1
+
+        # Compute the cross product and multiply by neighbors to get the same result as if all
+        # vectors were clockwise
+        cross_z = np.cross(
+            (
+                ordered_relative_soma_pts[closest_indices + neighbors]
+                - ordered_relative_soma_pts[closest_indices]
+            ),
+            relative_pts - ordered_relative_soma_pts[closest_indices],
+        )[:, COLS.Z] * signs
+
+        if exclude_boundary:
+            interior_side = (cross_z > 0)
+        else:
+            interior_side = (cross_z >= 0)
+
+        return interior_side
+
 
 def make_soma(morphio_soma):
     """Make a soma object from a MorphIO soma.

tests/core/test_soma.py

@@ -206,3 +206,65 @@ def test_Soma_Cylinders():
 
     assert list(s.center) == [0., 0., 0.]
     assert_almost_equal(s.area, 444.288293851) # cone area, not including bottom
+
+
+def test_soma_overlaps():
+    # Test with spherical soma
+    sm = load_morphology(StringIO(u"""1 1 11 22 33 44 -1"""), reader='swc').soma
+    assert isinstance(sm, soma.SomaSinglePoint)
+    points = [
+        [11, 22, 33],  # at soma center
+        [11, 22, 33 + 44],  # on the edge of the soma
+        [100, 100, 100],  # outside the soma
+    ]
+    np.testing.assert_array_equal(sm.overlaps(points), [True, True, False])
+    np.testing.assert_array_equal(sm.overlaps(points, exclude_boundary=True), [True, False, False])
+
+    # Test with cynlindrical soma
+    sm = load_morphology(StringIO(u"""
+                1 1 0 0 -10 40 -1
+                2 1 0 0   0 40  1
+                3 1 0 0  10 40  2"""), reader='swc').soma
+    assert isinstance(sm, soma.SomaCylinders)
+    points = [
+        [0, 0, -20],  # on the axis of the cylinder but outside it
+        [0, 0, -10],  # on the axis of the cylinder and on it's edge
+        [0, 0, -5],  # on the axis of the cylinder and inside it
+        [10, 10, 5],  # inside a cylinder
+        [100, 0, 0],  # outside all cylinders
+    ]
+    np.testing.assert_array_equal(sm.overlaps(points), [False, True, True, True, False])
+    np.testing.assert_array_equal(sm.overlaps(points, exclude_boundary=True), [False, False, True, True, False])
+
+    # Test with all points in soma for coverage
+    sm = load_morphology(StringIO(u"""
+                1 1 0 0 -10 40 -1
+                2 1 0 0   0 40  1
+                3 1 0 0  10 40  2"""), reader='swc').soma
+    assert isinstance(sm, soma.SomaCylinders)
+    points = [
+        [0, 0, -10],  # on the axis of the cylinder and on it's edge
+        [0, 0, -5],  # on the axis of the cylinder and inside it
+        [10, 10, 5],  # inside a cylinder
+    ]
+    np.testing.assert_array_equal(sm.overlaps(points), [True, True, True])
+    np.testing.assert_array_equal(sm.overlaps(points, exclude_boundary=True), [False, True, True])
+
+    # Test with contour soma
+    sm = load_morphology(StringIO(u"""
+                ((CellBody)
+                    (1 0 0 1)
+                    (1 1 0 1)
+                    (-1 1 0 1)
+                    (-1 0 0 1)) """), reader='asc').soma
+    assert isinstance(sm, soma.SomaSimpleContour)
+    points = [
+        [0, 0.5, 0],  # on the center of the soma
+        [1, 0.5, 0],  # on an edge of the soma
+        [1, 0, 0],  # on a corner of the soma
+        [0.25, 0.75, 0],  # inside the soma
+        [-0.9, 0.55, 0],  # inside the soma with closest index equal to the last soma point
+        [2, 3, 0],  # outside the soma
+    ]
+    np.testing.assert_array_equal(sm.overlaps(points), [True, True, True, True, True, False])
+    np.testing.assert_array_equal(sm.overlaps(points, exclude_boundary=True), [True, False, False, True, True, False])