Make line intersection computation more robust (#937)

libgeos · Feb 26, 2024 · 6d0c831 · 6d0c831
1 parent abd1766
commit 6d0c831
Show file tree

Hide file tree

Showing 7 changed files with 82 additions and 71 deletions.
diff --git a/NEWS b/NEWS
@@ -2,6 +2,7 @@
 20xx-xx-xx
 
 - Fixes/Improvements:
+  - Intersection: change to using DoubleDouble computation to improve robustness (GH-937, Martin Davis)
   - Fix build on Illumus (GH-971)
   - PointOnSurface crashes with a collection containing a empty linestring (GH-1002, Paul Ramsey)
   - Fix IsSimpleOp for MultiPoint with empty element (GH-1005, Martin Davis)

diff --git a/src/algorithm/Intersection.cpp b/src/algorithm/Intersection.cpp
@@ -15,6 +15,7 @@
 #include <cmath>
 #include <vector>
 
+#include <geos/algorithm/CGAlgorithmsDD.h>
 #include <geos/algorithm/Intersection.h>
 
 namespace geos {
@@ -26,59 +27,7 @@ geom::Coordinate
 Intersection::intersection(const geom::Coordinate& p1, const geom::Coordinate& p2,
                            const geom::Coordinate& q1, const geom::Coordinate& q2)
 {
-    double minX0 = p1.x < p2.x ? p1.x : p2.x;
-    double minY0 = p1.y < p2.y ? p1.y : p2.y;
-    double maxX0 = p1.x > p2.x ? p1.x : p2.x;
-    double maxY0 = p1.y > p2.y ? p1.y : p2.y;
-
-    double minX1 = q1.x < q2.x ? q1.x : q2.x;
-    double minY1 = q1.y < q2.y ? q1.y : q2.y;
-    double maxX1 = q1.x > q2.x ? q1.x : q2.x;
-    double maxY1 = q1.y > q2.y ? q1.y : q2.y;
-
-    double intMinX = minX0 > minX1 ? minX0 : minX1;
-    double intMaxX = maxX0 < maxX1 ? maxX0 : maxX1;
-    double intMinY = minY0 > minY1 ? minY0 : minY1;
-    double intMaxY = maxY0 < maxY1 ? maxY0 : maxY1;
-
-    double midx = (intMinX + intMaxX) / 2.0;
-    double midy = (intMinY + intMaxY) / 2.0;
-
-    // condition ordinate values by subtracting midpoint
-    double p1x = p1.x - midx;
-    double p1y = p1.y - midy;
-    double p2x = p2.x - midx;
-    double p2y = p2.y - midy;
-    double q1x = q1.x - midx;
-    double q1y = q1.y - midy;
-    double q2x = q2.x - midx;
-    double q2y = q2.y - midy;
-
-    // unrolled computation using homogeneous coordinates eqn
-    double px = p1y - p2y;
-    double py = p2x - p1x;
-    double pw = p1x * p2y - p2x * p1y;
-
-    double qx = q1y - q2y;
-    double qy = q2x - q1x;
-    double qw = q1x * q2y - q2x * q1y;
-
-    double x = py * qw - qy * pw;
-    double y = qx * pw - px * qw;
-    double w = px * qy - qx * py;
-
-    double xInt = x/w;
-    double yInt = y/w;
-    geom::Coordinate rv;
-    // check for parallel lines
-    if (!std::isfinite(xInt) || !std::isfinite(yInt)) {
-        rv.setNull();
-        return rv;
-    }
-    // de-condition intersection point
-    rv.x = xInt + midx;
-    rv.y = yInt + midy;
-    return rv;
+    return CGAlgorithmsDD::intersection(p1, p2, q1, q2);
 }
 
 

diff --git a/tests/unit/algorithm/RobustLineIntersectionTest.cpp b/tests/unit/algorithm/RobustLineIntersectionTest.cpp
@@ -305,43 +305,37 @@ void object::test<2>
 // but apparently works now.
 // Possibly normalization has fixed this?
 //
-#if 0 // fails: finds 1 intersection rather then two
 template<>
 template<>
 void object::test<3>
 ()
 {
     std::vector<Coordinate> intPt;
-    intPt.push_back(Coordinate(2089426.5233462777, 1180182.3877339689));
-    intPt.push_back(Coordinate(2085646.6891757075, 1195618.7333999649));
+    intPt.push_back(Coordinate(2087600.4716727887, 1187639.7426241424));
 
     checkIntersection(
         "LINESTRING ( 2089426.5233462777 1180182.3877339689, 2085646.6891757075 1195618.7333999649 )",
         "LINESTRING ( 1889281.8148903656 1997547.0560044837, 2259977.3672235999 483675.17050843034 )",
-        2,
+        1,
         intPt, 0);
 }
-#endif // fails
 
-#if 0 // fails: the intersection point doesn't match
 // 4 - Outside envelope using HCoordinate method (testCmp5CaseWKT)
 template<>
 template<>
 void object::test<4>
 ()
 {
     std::vector<Coordinate> intPt;
-    intPt.push_back(Coordinate(4348437.0557510145, 5552597.375203926));
+    intPt.push_back(Coordinate(4348440.8493873989, 5552599.2720221197));
 
     checkIntersection(
         "LINESTRING (4348433.262114629 5552595.478385733, 4348440.849387404 5552599.272022122 )",
         "LINESTRING (4348433.26211463  5552595.47838573,  4348440.8493874   5552599.27202212  )",
         1,
         intPt, 0);
 }
-#endif // fails
 
-#if 0 // fails: the intersection point doesn't match
 // 5 - Result of this test should be the same as the WKT one!
 //     (testCmp5CaseRaw)
 template<>
@@ -356,11 +350,10 @@ void object::test<5>
     pt.push_back(Coordinate(4348440.8493874,   5552599.27202212));
 
     std::vector<Coordinate> intPt;
-    intPt.push_back(Coordinate(4348437.0557510145, 5552597.375203926));
+    intPt.push_back(Coordinate(4348440.8493873989, 5552599.2720221197));
 
     checkIntersection(pt, 1, intPt, 0);
 }
-#endif // fails
 
 /**
  * Test involving two non-almost-parallel lines.
@@ -453,7 +446,7 @@ void object::test<10>
         "LINESTRING (-215.22279674875 -158.65425425385, -218.1208801283 -160.68343590235)",
         1,
         "POINT ( -215.22279674875 -158.65425425385 )",
-        0);
+        1.0e-10);
 }
 
 /**

diff --git a/tests/unit/operation/relate/RelateOpTest.cpp b/tests/unit/operation/relate/RelateOpTest.cpp
@@ -84,4 +84,28 @@ void object::test<2> ()
     "FF10F0102"    );
 }
 
+// see https://github.com/locationtech/jts/issues/396
+// testContainsNoding
+template<>
+template<>
+void object::test<4> ()
+{
+    checkRelate(
+        "LINESTRING (1 0, 0 2, 0 0, 2 2)",
+        "LINESTRING (0 0, 2 2)",
+        "101F00FF2"    );
+}
+
+// see https://github.com/libgeos/geos/issues/933
+// testContainsNoding
+template<>
+template<>
+void object::test<5> ()
+{
+    checkRelate(
+        "MULTILINESTRING ((0 0, 1 1), (0.5 0.5, 1 0.1, -1 0.1))",
+        "LINESTRING (0 0, 1 1)",
+        "1F1000FF2"    );
+}
+
 } // namespace tut
diff --git a/tests/xmltester/tests/general/TestRelateLL.xml b/tests/xmltester/tests/general/TestRelateLL.xml
@@ -1,5 +1,4 @@
 <run>
-  <precisionModel scale="1.0" offsetx="0.0" offsety="0.0"/>
 
 <case>
   <desc>LL - disjoint, non-overlapping envelopes</desc>
@@ -308,7 +307,7 @@
 </case>
 
 <case>
-<desc>LA - closed multiline / empty line</desc>
+<desc>LL - closed multiline / empty line</desc>
   <a>
     MULTILINESTRING ((0 0, 0 1), (0 1, 1 1, 1 0, 0 0))
   </a>
@@ -320,4 +319,50 @@
   </test>
 </case>
 
+<case>
+<desc>LL - test intersection node computation (see https://github.com/locationtech/jts/issues/396)</desc>
+  <a>
+    LINESTRING (1 0, 0 2, 0 0, 2 2)
+  </a>
+  <b>
+    LINESTRING (0 0, 2 2)
+  </b>
+  <test>
+    <op name="relate" arg1="A" arg2="B" arg3="101F00FF2">true</op>
+  </test>
+  <test><op name="contains"   arg1="A" arg2="B"> true  </op></test>
+  <test><op name="coveredBy"  arg1="A" arg2="B"> false </op></test>
+  <test><op name="covers"     arg1="A" arg2="B"> true  </op></test>
+  <test><op name="crosses"    arg1="A" arg2="B"> false </op></test>
+  <test><op name="disjoint"   arg1="A" arg2="B"> false </op></test>
+  <test><op name="equalsTopo" arg1="A" arg2="B"> false </op></test>
+  <test><op name="intersects" arg1="A" arg2="B"> true  </op></test>
+  <test><op name="overlaps"   arg1="A" arg2="B"> false </op></test>
+  <test><op name="touches"    arg1="A" arg2="B"> false </op></test>
+  <test><op name="within"     arg1="A" arg2="B"> false </op></test>
+</case>
+
+<case>
+<desc>LL - test intersection node computation (see https://github.com/libgeos/geos/issues/933)</desc>
+  <a>
+    MULTILINESTRING ((0 0, 1 1), (0.5 0.5, 1 0.1, -1 0.1))
+  </a>
+  <b>
+    LINESTRING (0 0, 1 1)
+  </b>
+  <test>
+    <op name="relate" arg1="A" arg2="B" arg3="1F1000FF2">true</op>
+  </test>
+  <test><op name="contains"   arg1="A" arg2="B"> true  </op></test>
+  <test><op name="coveredBy"  arg1="A" arg2="B"> false </op></test>
+  <test><op name="covers"     arg1="A" arg2="B"> true  </op></test>
+  <test><op name="crosses"    arg1="A" arg2="B"> false </op></test>
+  <test><op name="disjoint"   arg1="A" arg2="B"> false </op></test>
+  <test><op name="equalsTopo" arg1="A" arg2="B"> false </op></test>
+  <test><op name="intersects" arg1="A" arg2="B"> true  </op></test>
+  <test><op name="overlaps"   arg1="A" arg2="B"> false </op></test>
+  <test><op name="touches"    arg1="A" arg2="B"> false </op></test>
+  <test><op name="within"     arg1="A" arg2="B"> false </op></test>
+</case>
+
 </run>
diff --git a/tests/xmltester/tests/issue/issue-geos-398.xml b/tests/xmltester/tests/issue/issue-geos-398.xml
@@ -32,9 +32,8 @@ POLYGON ((50.0044466166182886 0.0,
   </b>
 <test>
  <op name="union" arg1="A" arg2="B">
-MULTIPOLYGON (((0 25.950779066386122, 0 26.860827855779647, 0 29.83879230192533, 0 35.73921397193218, 60 13.183894681853793, 60 10.39852836763784, 60 7.802134559422377, 60 6.657099879646016, 60 6.510515132098641, 60 0, 55.31611974965083 0, 50.00444661661829 0, 43.31611974965083 0, 0 0, 0 14.034613342373582, 0 17.92266612923108, 0 21.587486526024364, 0 25.950779066386122),
-  (0 21.587486526024364, 34.025852439655786 6.898140262297274, 34.02585243965579 6.898140262297273, 0 21.587486526024364)),
-  ((13.44557253233479 36, 60 36, 60 16.794451829809802, 60 16.36440115550932, 60 14.043996030454757, 2.9187843276549756 36, 11.8945390820011 36, 13.44557253233479 36)))
+MULTIPOLYGON (((0 26.860827855779647, 0 29.83879230192533, 0 35.73921397193218, 60 13.183894681853793, 60 10.39852836763784, 60 7.802134559422377, 60 6.657099879646016, 60 6.510515132098641, 60 0, 55.31611974965083 0, 50.00444661661829 0, 43.31611974965083 0, 0 0, 0 14.034613342373582, 0 17.92266612923108, 0 21.587486526024364, 0 25.950779066386122, 0 26.860827855779647)), 
+  ((60 36, 60 16.794451829809802, 60 16.36440115550932, 60 14.043996030454757, 2.9187843276549756 36, 11.8945390820011 36, 13.44557253233479 36, 60 36)))
  </op>
 </test>
 </case>

diff --git a/tests/xmltester/tests/robust/overlay/TestOverlay-pg-list.xml b/tests/xmltester/tests/robust/overlay/TestOverlay-pg-list.xml
@@ -47,7 +47,7 @@ MULTIPOLYGON (((742591.4178018438 5087769.04526206, 742591.417801844 5087769.045
 
 <test>
 <op name="intersection" arg1="A" arg2="B">
-MULTIPOLYGON (((434608.9999999999 2012571.9999999907, 434608.99999999977 2012571.9999999907, 434505.9999999997 2013031.9999999907, 434889.4374999999 2013081.8073453507, 434506 2013031.99999999, 434608.9999999999 2012571.9999999907)), ((434990.9999999999 2013094.99999999, 435264.9999999997 2013496.9999999898, 435311.2499999999 2013535.9783236892, 435265 2013496.99999999, 434991 2013094.99999999, 434990.9999999999 2013094.99999999)), ((436813.9999999997 2014014.99999999, 436802.9999999997 2014278.9999999893, 436803 2014278.9999999893, 436814 2014014.99999999, 436813.9999999997 2014014.99999999)), ((437427.6587183307 2014442.90461996, 437488.99999999977 2014458.99999999, 437372.9999999997 2016223.9999999898, 437686.7499999999 2016676.926282041, 437373 2016223.99999999, 437489 2014458.99999999, 437427.6587183307 2014442.90461996)))
+MULTIPOLYGON (((437372.9999999997 2016223.9999999898, 437686.7499999997 2016676.9262820408, 437373 2016223.99999999, 437489 2014458.99999999, 437427.6587183306 2014442.90461996, 437488.99999999977 2014458.99999999, 437372.9999999997 2016223.9999999898)), ((434990.9999999999 2013094.99999999, 435264.9999999997 2013496.9999999898, 435311.2499999997 2013535.9783236892, 435265 2013496.99999999, 434991 2013094.99999999, 434990.9999999999 2013094.99999999)), ((434505.9999999997 2013031.9999999907, 434889.4374999997 2013081.8073453507, 434506 2013031.99999999, 434608.9999999999 2012571.9999999907, 434608.99999999977 2012571.9999999907, 434505.9999999997 2013031.9999999907)), ((436802.9999999997 2014278.9999999893, 436803 2014278.9999999893, 436814 2014014.99999999, 436813.9999999997 2014014.99999999, 436802.9999999997 2014278.9999999893)))
 </op>
 </test>