diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 3a57c542..ee88ccd7 100644 --- a/CHANGELOG.rst +++ b/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 diff --git a/neurom/core/soma.py b/neurom/core/soma.py index 3451547a..9c6efaea 100755 --- a/neurom/core/soma.py +++ b/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) ' % (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) ' % (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. diff --git a/tests/core/test_soma.py b/tests/core/test_soma.py index 517aeba0..aeeae238 100644 --- a/tests/core/test_soma.py +++ b/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])