Merge pull request #154 from AllenInstitute/depth-field-crossing-fix

pia/wm boundary intersection check and fix

neuron_morphology/transforms/geometry.py

@@ -2,7 +2,7 @@
 """
 
 from typing import List, Tuple, Callable
-
+import warnings
 from shapely import geometry as geo
 
 
@@ -13,6 +13,52 @@ def get_ccw_vertices_from_two_lines(line1: List[Tuple], line2: List[Tuple]):
     """
     return get_ccw_vertices(get_vertices_from_two_lines(line1, line2))
 
+def prune_two_lines(line1: List[Tuple], line2: List[Tuple]):
+    """
+        check the boundary to avoid intersections with side lines
+
+        Parameters
+            ----------
+            line1, line2: List of coordinates describing two lines
+
+            Returns
+            -------
+            line1, line2: boundary pruned if needed
+    """
+    # validate the pia/wm does not cross the side lines
+    prune = True
+    while prune:
+        prune = False
+
+        # update lines
+        line1_str = geo.LineString(line1)
+        line2_str = geo.LineString(line2)
+
+        side1 = geo.LineString([line1[0], line2[-1]])
+        side2 = geo.LineString([line1[-1], line2[0]])
+
+        # prune the edge points if there are intersections of side and pia/wm
+        if side1.crosses(line1_str):
+            line1.pop(0)
+            prune = True
+
+        if side1.crosses(line2_str):
+            line2.pop(-1)
+            prune = True
+
+        if side2.crosses(line1_str):
+            line1.pop(-1)
+            prune = True
+
+        if side2.crosses(line2_str):
+            line2.pop(0)
+            prune = True
+
+        if prune:
+            warnings.warn(f"lines are modified \nline1: {line1}\nline2: {line2}", UserWarning)
+
+    return line1, line2
+
 def get_vertices_from_two_lines(line1: List[Tuple], line2: List[Tuple]):
     """
         Generates circular vertices from two lines
@@ -36,6 +82,8 @@ def get_vertices_from_two_lines(line1: List[Tuple], line2: List[Tuple]):
     if side1.crosses(side2):
         line2.reverse()
 
+    line1, line2 = prune_two_lines(line1, line2)
+
     vertices = line1 + line2 + [line1[0]]
     return vertices
 

tests/transforms/test_geometry.py

@@ -45,4 +45,13 @@ def test_ccw_both_cw(self):
         line1 = [(0, 1), (1, 1)]
         line2 = [(1, 0), (0, 0)]
         result = geo.get_ccw_vertices_from_two_lines(line1, line2)
-        self.assertEqual(result, [(0, 1), (0, 0), (1, 0), (1, 1), (0, 1)])
+        self.assertEqual(result, [(0, 1), (0, 0), (1, 0), (1, 1), (0, 1)])
+
+    def test_prune_two_lines(self):
+        """
+            line1 (1, 1), (0, 1), (0.5, 1.5) -> (1, 1), (0, 1)
+        """
+        line1 = [(1, 1), (0, 1), (0.5, 1.5)]
+        line2 = [(0, 0), (1, 0)]
+        line1_pruned, line2_pruned = geo.prune_two_lines(line1, line2)
+        self.assertEqual(line1_pruned, [(1, 1), (0, 1)])