From 668c4bce68c0043f5f2bf6da6d8d481264e391f9 Mon Sep 17 00:00:00 2001 From: Martin Davis Date: Tue, 4 Jul 2023 08:47:10 -0700 Subject: [PATCH] Make intersection computation more robust (#989) --- .../jts/algorithm/Intersection.java | 17 +++++++ .../algorithm/RobustLineIntersectionTest.java | 8 +-- .../jts/operation/relate/ContainsTest.java | 41 +++++++++++++--- .../jts/operation/relate/RelateTest.java | 32 +++++++++--- .../testxml/general/TestRelateLL.xml | 49 ++++++++++++++++++- .../general/TestUnaryUnionFloating.xml | 3 +- 6 files changed, 128 insertions(+), 22 deletions(-) diff --git a/modules/core/src/main/java/org/locationtech/jts/algorithm/Intersection.java b/modules/core/src/main/java/org/locationtech/jts/algorithm/Intersection.java index cc86506199..34ac57675b 100644 --- a/modules/core/src/main/java/org/locationtech/jts/algorithm/Intersection.java +++ b/modules/core/src/main/java/org/locationtech/jts/algorithm/Intersection.java @@ -44,6 +44,23 @@ public class Intersection { * @see CGAlgorithmsDD#intersection(Coordinate, Coordinate, Coordinate, Coordinate) */ public static Coordinate intersection(Coordinate p1, Coordinate p2, Coordinate q1, Coordinate q2) { + return CGAlgorithmsDD.intersection(p1, p2, q1, q2); + } + + /** + * Compute intersection of two lines, using a floating-point algorithm. + * This is less accurate than {@link CGAlgorithmsDD#intersection(Coordinate, Coordinate, Coordinate, Coordinate)}. + * It has caused spatial predicate failures in some cases. + * This is kept for testing purposes. + * + * @param p1 an endpoint of line 1 + * @param p2 an endpoint of line 1 + * @param q1 an endpoint of line 2 + * @param q2 an endpoint of line 2 + * @return the intersection point between the lines, if there is one, + * or null if the lines are parallel or collinear + */ + private static Coordinate intersectionFP(Coordinate p1, Coordinate p2, Coordinate q1, Coordinate q2) { // compute midpoint of "kernel envelope" double minX0 = p1.x < p2.x ? p1.x : p2.x; double minY0 = p1.y < p2.y ? p1.y : p2.y; diff --git a/modules/core/src/test/java/org/locationtech/jts/algorithm/RobustLineIntersectionTest.java b/modules/core/src/test/java/org/locationtech/jts/algorithm/RobustLineIntersectionTest.java index 08f7485878..15419ca95c 100644 --- a/modules/core/src/test/java/org/locationtech/jts/algorithm/RobustLineIntersectionTest.java +++ b/modules/core/src/test/java/org/locationtech/jts/algorithm/RobustLineIntersectionTest.java @@ -74,7 +74,7 @@ public void testCentralEndpointHeuristicFailure2() "LINESTRING (-58.00593335955 -1.43739086465, -513.86101637525 -457.29247388035)", "LINESTRING (-215.22279674875 -158.65425425385, -218.1208801283 -160.68343590235)", 1, - "POINT ( -215.22279674875 -158.65425425385 )", + "POINT ( -215.22279674875003 -158.65425425385004 )", 0); } @@ -192,7 +192,7 @@ public void testDaveSkeaCase() "LINESTRING ( 1889281.8148903656 1997547.0560044837, 2259977.3672235999 483675.17050843034 )", 1, new Coordinate[] { - new Coordinate(2087536.6062609926, 1187900.560566967), + new Coordinate(2087600.4716727887, 1187639.7426241424), }, 0); } @@ -209,7 +209,7 @@ public void testCmp5CaseWKT() "LINESTRING (4348433.26211463 5552595.47838573, 4348440.8493874 5552599.27202212 )", 1, new Coordinate[] { - new Coordinate(4348440.8493874, 5552599.27202212), + new Coordinate(4348440.849387399, 5552599.27202212), }, 0); } @@ -230,7 +230,7 @@ public void testCmp5CaseRaw() new Coordinate(4348440.8493874, 5552599.27202212) }, 1, new Coordinate[] { - new Coordinate(4348440.8493874, 5552599.27202212), + new Coordinate(4348440.849387399, 5552599.27202212), }, 0); } diff --git a/modules/core/src/test/java/org/locationtech/jts/operation/relate/ContainsTest.java b/modules/core/src/test/java/org/locationtech/jts/operation/relate/ContainsTest.java index 2ca6069608..0cc2a459f1 100644 --- a/modules/core/src/test/java/org/locationtech/jts/operation/relate/ContainsTest.java +++ b/modules/core/src/test/java/org/locationtech/jts/operation/relate/ContainsTest.java @@ -16,8 +16,8 @@ import org.locationtech.jts.geom.GeometryFactory; import org.locationtech.jts.io.WKTReader; -import junit.framework.TestCase; import junit.textui.TestRunner; +import test.jts.GeometryTestCase; /** @@ -27,9 +27,9 @@ * @version 1.7 */ public class ContainsTest - extends TestCase + extends GeometryTestCase { - public static void main(String args[]) { + public static void main(String[] args) { TestRunner.run(ContainsTest.class); } @@ -45,6 +45,8 @@ public ContainsTest(String name) * From GEOS #572. * A case where B is contained in A, but * the JTS relate algorithm fails to compute this correctly. + * when using an FP intersection algorithm. + * This case works when using CGAlgorithmsDD#intersection(Coordinate, Coordinate, Coordinate, Coordinate). * * The cause is that the long segment in A nodes the single-segment line in B. * The node location cannot be computed precisely. @@ -63,8 +65,35 @@ public void testContainsIncorrect() { String a = "LINESTRING (1 0, 0 2, 0 0, 2 2)"; String b = "LINESTRING (0 0, 2 2)"; - - // for now assert this as false, although it should be true - assertTrue(! a.contains(b)); + checkContains(a, b); + } + + /** + * From GEOS #933. + * A case where B is contained in A, but + * the JTS relate algorithm fails to compute this correctly. + * when using an FP intersection algorithm. + * This case works when using CGAlgorithmsDD#intersection(Coordinate, Coordinate, Coordinate, Coordinate). + */ + public void testContainsGEOS933() + throws Exception + { + String a = "MULTILINESTRING ((0 0, 1 1), (0.5 0.5, 1 0.1, -1 0.1))"; + String b = "LINESTRING (0 0, 1 1)"; + checkContains(a, b); + } + + private void checkContains(String wktA, String wktB) { + Geometry geomA = read(wktA); + Geometry geomB = read(wktB); + boolean actual = geomA.contains(geomB); + assertTrue(actual); + } + + private void checkContainsError(String wktA, String wktB) { + Geometry geomA = read(wktA); + Geometry geomB = read(wktB); + boolean actual = geomA.contains(geomB); + assertFalse(actual); } } diff --git a/modules/core/src/test/java/org/locationtech/jts/operation/relate/RelateTest.java b/modules/core/src/test/java/org/locationtech/jts/operation/relate/RelateTest.java index d3f8fee011..8b7ee3e872 100644 --- a/modules/core/src/test/java/org/locationtech/jts/operation/relate/RelateTest.java +++ b/modules/core/src/test/java/org/locationtech/jts/operation/relate/RelateTest.java @@ -42,21 +42,37 @@ public RelateTest(String name) } /** - * From GEOS #572 + * From https://github.com/locationtech/jts/issues/396 * - * The cause is that the longer line nodes the single-segment line. - * The node then tests as not lying precisely on the original longer line. + * The original failure is caused by the intersection computed + * during noding not lying exactly on each original line segment. + * This is due to numerical error in the FP intersection algorithm. + * This is fixed by using DD intersection calculation. */ - public void testContainsIncorrectIMMatrix() + public void testContainsNoding() { String a = "LINESTRING (1 0, 0 2, 0 0, 2 2)"; String b = "LINESTRING (0 0, 2 2)"; - // actual matrix is 001F001F2 - // true matrix should be 101F00FF2 - runRelateTest(a, b, "001F001F2" ); + runRelateTest(a, b, "101F00FF2" ); } + /** + * From GEOS https://github.com/libgeos/geos/issues/933 + * + * The original failure is caused by the intersection computed + * during noding not lying exactly on each original line segment. + * This is due to numerical error in the FP intersection algorithm. + * This is fixed by using DD intersection calculation. + */ + public void testContainsNoding2() + { + String a = "MULTILINESTRING ((0 0, 1 1), (0.5 0.5, 1 0.1, -1 0.1))"; + String b = "LINESTRING (0 0, 1 1)"; + + runRelateTest(a, b, "1F1000FF2" ); + } + /** * Tests case where segments intersect properly, but computed intersection point * snaps to a boundary endpoint due to roundoff. @@ -93,6 +109,6 @@ void runRelateTest(String wkt1, String wkt2, String expectedIM) IntersectionMatrix im = RelateOp.relate(g1, g2); String imStr = im.toString(); //System.out.println(imStr); - assertTrue(im.matches(expectedIM)); + assertEquals(expectedIM, imStr); } } diff --git a/modules/tests/src/test/resources/testxml/general/TestRelateLL.xml b/modules/tests/src/test/resources/testxml/general/TestRelateLL.xml index 4c05b66642..4a4fb4289d 100644 --- a/modules/tests/src/test/resources/testxml/general/TestRelateLL.xml +++ b/modules/tests/src/test/resources/testxml/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/modules/tests/src/test/resources/testxml/general/TestUnaryUnionFloating.xml b/modules/tests/src/test/resources/testxml/general/TestUnaryUnionFloating.xml index a8bc222ad3..fd6af522e7 100644 --- a/modules/tests/src/test/resources/testxml/general/TestUnaryUnionFloating.xml +++ b/modules/tests/src/test/resources/testxml/general/TestUnaryUnionFloating.xml @@ -12,8 +12,7 @@ - POLYGON ((0.2942036115049298 2.226702215615205, -3 -2, -6 900, 700 900, 699.6114719806972 899.5000276853219, 700 899.5, 700 860, 670.2204017222961 861.6785046602191, 0.3 -0.4, 0.2942036115049298 2.226702215615205)) - + POLYGON ((-3 -2, -6 900, 700 900, 699.6114719806972 899.5000276853219, 700 899.5, 700 860, 670.2204017222961 861.6785046602191, 0.3 -0.4, 0.29420361150493 2.226702215615145, -3 -2))