Repository: matplotlib/matplotlib
Pull request: Improved Code for Segments Intersect
## PR Summary

Supposedly fixes #11527 (Inconsistent path intersection).

`Path.intersects_path` didn't work properly on
collinear paths intersection because underlying 
[`segments_intersect`](https://github.com/matplotlib/matplotlib/blob/master/src/_path.h#L815)
handled parallel segments as non-intersecting (ignoring possible segments' overlap).

This PR handles this case of collinear paths and makes `Path.intersects_path` return `True`
if paths have at least one common point.

Added an appropriate test case to https://github.com/matplotlib/matplotlib/blob/master/lib/matplotlib/tests/test_path.py

## PR Checklist

- [ ] Has Pytest style unit tests
- [ ] Code is PEP 8 compliant
- [ ] New features are documented, with examples if plot related
- [ ] Documentation is sphinx and numpydoc compliant
- [ ] Added an entry to doc/users/next_whats_new/ if major new feature (follow instructions in README.rst there)
- [ ] Documented in doc/api/api_changes.rst if API changed in a backward-incompatible way

<!--
Thank you so much for your PR!  To help us review your contribution, please
consider the following points:

- A development guide is available at https://matplotlib.org/devdocs/devel/index.html.

- Help with git and github is available at
  https://matplotlib.org/devel/gitwash/development_workflow.html.

- Do not create the PR out of master, but out of a separate branch.

- The PR title should summarize the changes, for example "Raise ValueError on
  non-numeric input to set_xlim".  Avoid non-descriptive titles such as
  "Addresses issue #8576".

- The summary should provide at least 1-2 sentences describing the pull request
  in detail (Why is this change required?  What problem does it solve?) and
  link to any relevant issues.

- If you are contributing fixes to docstrings, please pay attention to
  http://matplotlib.org/devel/documenting_mpl.html#formatting.  In particular,
  note the difference between using single backquotes, double backquotes, and
  asterisks in the markup.

We understand that PRs can sometimes be overwhelming, especially as the
reviews start coming in.  Please let us know if the reviews are unclear or
the recommended next step seems overly demanding, if you would like help in
addressing a reviewer's comments, or if you have been waiting too long to hear
back on your PR.
-->


Problem statement:
Inconsistent path intersection
### Bug report

**Bug summary**
If I two paths has a common segment, the result of `a.intersects_path(b)` is inconsistent. I think it should return `True` every time in this case.

**Code for reproduction**
```python
from matplotlib.path import Path

a = Path([(0,1),(0,3)]) 
b = Path([(0,0),(0,4)])   # larger than a
c = Path([(0,0),(0,2)])   # ends in a
d = Path([(0,1),(0,2),(0,3)])   # same as a, but has an extra point in a
e = Path([(0,0),(0,2),(0,4)])   # same as b, but has an extra point in a
f = Path([(0,0),(0,1),(0,3),(0,4)])  # contains a as a segment
g = Path([(0,1),(0,3),(1,3)])   # a and a non-collinear point

paths = [a,b,c,d,e,f, g]
for i in paths:
    print([1 if p.intersects_path(i) else 0 for p in paths])
```

**Actual outcome**
```
[0, 0, 0, 0, 1, 1, 1]
[0, 0, 0, 0, 0, 0, 1]
[0, 0, 0, 0, 0, 0, 0]
[0, 0, 0, 0, 1, 1, 1]
[1, 0, 0, 1, 0, 0, 1]
[1, 0, 0, 1, 0, 0, 1]
[1, 1, 0, 1, 1, 1, 1]
```

**Expected outcome**
True for all (since they all have the `[(0,1),(0,2)]` segment in common).

**Matplotlib version**
  * Operating system: Ubuntu 16.04
  * Matplotlib version: 2.2.2
  * Matplotlib backend (`print(matplotlib.get_backend())`): module://ipykernel.pylab.backend_inline
  * Python version: 3.6.5
  * Jupyter version (if applicable): 5.4.1

Hints from issue discussion:
That is interesting.  

I am not sure what the convention for 'intersection' with co-linear lines is.  If by "intersect" you mean "some points fall on one side and some points fall on the other" then there is a case that everything but 'g' should not be considered intersecting.
What does shapely do in these cases?  We should probably do the same thing as I trust that there are experts in computational geometry working on that project.
I expect that if two paths have at least one common point, intersects_path() should return True. It looks like shapely thinks the same:
```python
from shapely.geometry import LineString

a = LineString([(0,1),(0,3)]) 
b = LineString([(0,0),(0,4)])   # larger than a
c = LineString([(0,0),(0,2)])   # ends in a
d = LineString([(0,1),(0,2),(0,3)])   # same as a, but has an extra point in a
e = LineString([(0,0),(0,2),(0,4)])   # same as b, but has an extra point in a
f = LineString([(0,0),(0,1),(0,3),(0,4)])  # contains a as a segment
g = LineString([(0,1),(0,3),(1,3)])   # a and a not co-linear point

paths = [a,b,c,d,e,f, g]
for i in paths:
    print([1 if p.intersects(i) else 0 for p in paths])
```
returns
```
[1, 1, 1, 1, 1, 1, 1]
[1, 1, 1, 1, 1, 1, 1]
[1, 1, 1, 1, 1, 1, 1]
[1, 1, 1, 1, 1, 1, 1]
[1, 1, 1, 1, 1, 1, 1]
[1, 1, 1, 1, 1, 1, 1]
[1, 1, 1, 1, 1, 1, 1]
```

I can accept that if by "intersect" we mean two lines crossing each other in exactly one point, but this is seriously wrong:
```python
a = Path([(0,1),(0,2),(0,3)])   # co-linear points
b = Path([(0,1),(0,2),(1,3)])  # not co-linear points
>>> a.intersects_path(a)
False
>>> b.intersects_path(b)
True
``` 
Thanks for doing the legwork on getting those answers!

I suspect the bug lies someplace in https://github.com/matplotlib/matplotlib/blob/master/src/_path_wrapper.cpp or https://github.com/matplotlib/matplotlib/blob/master/src/_path.h

attn @mdboom @QuLogic 
It looks the issue is caused by the early exit condition in `segments_intersect`
https://github.com/matplotlib/matplotlib/blob/master/src/_path.h#L825

```
double den = ((y4 - y3) * (x2 - x1)) - ((x4 - x3) * (y2 - y1));
    if (den == 0.0) {
        return false;
    }
```

basically, it returns `false` if two segments have the same slopes, yet it doesn't take into account the case when segments overlap.

Commit message: Merge branch 'segments-intersect' of https://github.com/TarasKuzyo/matplotlib into segments-intersect

Patch to review:
diff --git a/lib/matplotlib/tests/test_path.py b/lib/matplotlib/tests/test_path.py
index c99701097d95..f9d16cf5dfd4 100644
--- a/lib/matplotlib/tests/test_path.py
+++ b/lib/matplotlib/tests/test_path.py
@@ -273,6 +273,82 @@ def test_path_deepcopy():
     copy.deepcopy(path2)
 
 
+def test_path_intersect_path():
+    # test for the range of intersection angles
+    for phi in np.linspace(0, 2*np.pi, 180):
+
+        transform = transforms.Affine2D().rotate(phi)
+
+        # a nd 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), np.rad2deg(phi)
+
+        # 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 power in range(5, 13):
+            eps = 10**(-power)
+            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 power in range(5, 13):
+            eps = 10**(-power)
+            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 power in range(5, 13):
+            eps = 10**(-power)
+            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
index c92321db2052..4a38aed5a621 100644
--- 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>


Code review comments:
diff_hunk: @@ -273,6 +273,82 @@ def test_path_deepcopy():
     copy.deepcopy(path2)
 
 
+def test_path_intersect_path():
+    # test for the range of intersection angles
+    for phi in np.linspace(0, 2*np.pi, 180):
path: lib/matplotlib/tests/test_path.py
text: @user1:
Going through 180 angles seems a bit much. Have you checked, that the test still runs fast (<< 1s)? Otherwise it should be sufficient to go through some selected angles like
~~~
angles = np.array([0, 30, 45, 90, ...])
angles = np.concatenate([angles, angles+1, angles-1])
~~~

@author:
@user1, it does take about 1s to complete. I will use the approach you've suggested.

@author:
So, I've reduced the total range of angles to 30 so it runs faster.