Improved Code for Segments Intersect (#11553)

TarasKuzyo · timhoffm · commit 4d1380f7915a · 2019-03-02T23:54:06.000+01:00
* added function to check if a point lies on a segment

* check if collinear segments have common points in segments_intersect

* test cases added for Path.intersects_path

* PEP8 check

* added a test case for a range of slopes of intersecting paths

* changed the way segments are compared provided their slopes and intercepts match

* typo fix

* simplified logic of segments_intersect when handling the infinite slope case

* implemented isclose for path intersect check; added more test cases

* removed debug include

* flake8 fixes

* typo fix

Co-Authored-By: TarasKuzyo &lt;kuzyo.taras@gmail.com&gt;

* reduced the number of checks to 30 angles and 4 eps values
diff --git a/lib/matplotlib/tests/test_path.py b/lib/matplotlib/tests/test_path.py
@@ -273,6 +273,83 @@ def test_path_deepcopy():
     copy.deepcopy(path2)
 
 
+def test_path_intersect_path():
+    # test for the range of intersection angles
+    base_angles = np.array([0, 15, 30, 45, 60, 75, 90, 105, 120, 135])
+    angles = np.concatenate([base_angles, base_angles + 1, base_angles - 1])
+    eps_array = [1e-5, 1e-8, 1e-10, 1e-12]
+
+    for phi in angles:
+
+        transform = transforms.Affine2D().rotate(np.deg2rad(phi))
+
+        # a and b intersect at angle phi
+        a = Path([(-2, 0), (2, 0)])
+        b = transform.transform_path(a)
+        assert a.intersects_path(b) and b.intersects_path(a)
+
+        # a and b touch at angle phi at (0, 0)
+        a = Path([(0, 0), (2, 0)])
+        b = transform.transform_path(a)
+        assert a.intersects_path(b) and b.intersects_path(a)
+
+        # a and b are orthogonal and intersect at (0, 3)
+        a = transform.transform_path(Path([(0, 1), (0, 3)]))
+        b = transform.transform_path(Path([(1, 3), (0, 3)]))
+        assert a.intersects_path(b) and b.intersects_path(a)
+
+        # a and b are collinear and intersect at (0, 3)
+        a = transform.transform_path(Path([(0, 1), (0, 3)]))
+        b = transform.transform_path(Path([(0, 5), (0, 3)]))
+        assert a.intersects_path(b) and b.intersects_path(a)
+
+        # self-intersect
+        assert a.intersects_path(a)
+
+        # a contains b
+        a = transform.transform_path(Path([(0, 0), (5, 5)]))
+        b = transform.transform_path(Path([(1, 1), (3, 3)]))
+        assert a.intersects_path(b) and b.intersects_path(a)
+
+        # a and b are collinear but do not intersect
+        a = transform.transform_path(Path([(0, 1), (0, 5)]))
+        b = transform.transform_path(Path([(3, 0), (3, 3)]))
+        assert not a.intersects_path(b) and not b.intersects_path(a)
+
+        # a and b are on the same line but do not intersect
+        a = transform.transform_path(Path([(0, 1), (0, 5)]))
+        b = transform.transform_path(Path([(0, 6), (0, 7)]))
+        assert not a.intersects_path(b) and not b.intersects_path(a)
+
+        # Note: 1e-13 is the absolute tolerance error used for
+        # `isclose` function from src/_path.h
+
+        # a and b are parallel but do not touch
+        for eps in eps_array:
+            a = transform.transform_path(Path([(0, 1), (0, 5)]))
+            b = transform.transform_path(Path([(0 + eps, 1), (0 + eps, 5)]))
+            assert not a.intersects_path(b) and not b.intersects_path(a)
+
+        # a and b are on the same line but do not intersect (really close)
+        for eps in eps_array:
+            a = transform.transform_path(Path([(0, 1), (0, 5)]))
+            b = transform.transform_path(Path([(0, 5 + eps), (0, 7)]))
+            assert not a.intersects_path(b) and not b.intersects_path(a)
+
+        # a and b are on the same line and intersect (really close)
+        for eps in eps_array:
+            a = transform.transform_path(Path([(0, 1), (0, 5)]))
+            b = transform.transform_path(Path([(0, 5 - eps), (0, 7)]))
+            assert a.intersects_path(b) and b.intersects_path(a)
+
+        # b is the same as a but with an extra point
+        a = transform.transform_path(Path([(0, 1), (0, 5)]))
+        b = transform.transform_path(Path([(0, 1), (0, 2), (0, 5)]))
+        assert a.intersects_path(b) and b.intersects_path(a)
+
+    return
+
+
 @pytest.mark.parametrize('offset', range(-720, 361, 45))
 def test_full_arc(offset):
     low = offset
diff --git a/src/_path.h b/src/_path.h
@@ -813,6 +813,14 @@ int count_bboxes_overlapping_bbox(agg::rect_d &a, BBoxArray &bboxes)
     return count;
 }
 
+
+inline bool isclose(double a, double b, double rtol, double atol)
+{
+    // as per python's math.isclose
+    return fabs(a-b) <= fmax(rtol * fmax(fabs(a), fabs(b)), atol);
+}
+
+
 inline bool segments_intersect(const double &x1,
                                const double &y1,
                                const double &x2,
@@ -822,8 +830,27 @@ inline bool segments_intersect(const double &x1,
                                const double &x4,
                                const double &y4)
 {
+    // relative and absolute tolerance values are chosen empirically
+    // it looks the atol value matters here bacause of round-off errors
+    const double rtol = 1e-10;
+    const double atol = 1e-13;
+    // determinant
     double den = ((y4 - y3) * (x2 - x1)) - ((x4 - x3) * (y2 - y1));
-    if (den == 0.0) {
+
+    if (isclose(den, 0.0, rtol, atol)) {  // collinear segments
+        if (x1 == x2 && x2 == x3) { // segments have infinite slope (vertical lines)
+                                    // and lie on the same line
+            return (fmin(y1, y2) <= fmin(y3, y4) && fmin(y3, y4) <= fmax(y1, y2)) ||
+                   (fmin(y3, y4) <= fmin(y1, y2) && fmin(y1, y2) <= fmax(y3, y4));
+        }
+        else {
+            double intercept = (y1*x2 - y2*x1)*(x4 - x3) - (y3*x4 - y4*x3)*(x1 - x2);
+            if (isclose(intercept, 0.0, rtol, atol)) { // segments lie on the same line
+                return (fmin(x1, x2) <= fmin(x3, x4) && fmin(x3, x4) <= fmax(x1, x2)) ||
+                       (fmin(x3, x4) <= fmin(x1, x2) && fmin(x1, x2) <= fmax(x3, x4));
+            }
+        }
+       
         return false;
     }
 
@@ -833,7 +860,10 @@ inline bool segments_intersect(const double &x1,
     double u1 = n1 / den;
     double u2 = n2 / den;
 
-    return (u1 >= 0.0 && u1 <= 1.0 && u2 >= 0.0 && u2 <= 1.0);
+    return ((u1 > 0.0 || isclose(u1, 0.0, rtol, atol)) && 
+            (u1 < 1.0 || isclose(u1, 1.0, rtol, atol)) && 
+            (u2 > 0.0 || isclose(u2, 0.0, rtol, atol)) && 
+            (u2 < 1.0 || isclose(u2, 1.0, rtol, atol)));
 }
 
 template <class PathIterator1, class PathIterator2>
