From 6d0c831dfbbb861cf3445257a322bbb21600471e Mon Sep 17 00:00:00 2001 From: Mike Taves Date: Mon, 26 Feb 2024 22:02:51 +1300 Subject: [PATCH] Make line intersection computation more robust (#937) --- NEWS | 1 + src/algorithm/Intersection.cpp | 55 +------------------ .../algorithm/RobustLineIntersectionTest.cpp | 17 ++---- tests/unit/operation/relate/RelateOpTest.cpp | 24 ++++++++ .../xmltester/tests/general/TestRelateLL.xml | 49 ++++++++++++++++- .../xmltester/tests/issue/issue-geos-398.xml | 5 +- .../robust/overlay/TestOverlay-pg-list.xml | 2 +- 7 files changed, 82 insertions(+), 71 deletions(-) diff --git a/NEWS b/NEWS index a51f014697..9c3e38b9ec 100644 --- 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 index 455bd6048e..f62025f16a 100644 --- a/src/algorithm/Intersection.cpp +++ b/src/algorithm/Intersection.cpp @@ -15,6 +15,7 @@ #include #include +#include #include 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 index 8ee79e2945..595dbd79f3 100644 --- a/tests/unit/algorithm/RobustLineIntersectionTest.cpp +++ b/tests/unit/algorithm/RobustLineIntersectionTest.cpp @@ -305,25 +305,21 @@ 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 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<> @@ -331,7 +327,7 @@ void object::test<4> () { std::vector 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 )", @@ -339,9 +335,7 @@ void object::test<4> 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 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 index 711aa6850a..787a41ad7a 100644 --- 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 index 4c05b66642..4a4fb4289d 100644 --- a/tests/xmltester/tests/general/TestRelateLL.xml +++ b/tests/xmltester/tests/general/TestRelateLL.xml @@ -1,5 +1,4 @@ - LL - disjoint, non-overlapping envelopes @@ -308,7 +307,7 @@ -LA - closed multiline / empty line +LL - closed multiline / empty line MULTILINESTRING ((0 0, 0 1), (0 1, 1 1, 1 0, 0 0)) @@ -320,4 +319,50 @@ + +LL - test intersection node computation (see https://github.com/locationtech/jts/issues/396) + + LINESTRING (1 0, 0 2, 0 0, 2 2) + + + LINESTRING (0 0, 2 2) + + + true + + true + false + true + false + false + false + true + false + false + false + + + +LL - test intersection node computation (see https://github.com/libgeos/geos/issues/933) + + MULTILINESTRING ((0 0, 1 1), (0.5 0.5, 1 0.1, -1 0.1)) + + + LINESTRING (0 0, 1 1) + + + true + + true + false + true + false + false + false + true + false + false + false + + diff --git a/tests/xmltester/tests/issue/issue-geos-398.xml b/tests/xmltester/tests/issue/issue-geos-398.xml index 3143f09e04..edb3984ef0 100644 --- 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, -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))) diff --git a/tests/xmltester/tests/robust/overlay/TestOverlay-pg-list.xml b/tests/xmltester/tests/robust/overlay/TestOverlay-pg-list.xml index 66f488591c..f5f9e69ee9 100644 --- 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 -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)))