Make intersection computation more robust (#989)

locationtech · Jul 4, 2023 · 668c4bc · 668c4bc
1 parent d6a280c
commit 668c4bc
Show file tree

Hide file tree

Showing 6 changed files with 128 additions and 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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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/modules/tests/src/test/resources/testxml/general/TestUnaryUnionFloating.xml b/modules/tests/src/test/resources/testxml/general/TestUnaryUnionFloating.xml
@@ -12,8 +12,7 @@
   </a>
 <test>
   <op name="union" arg1="A">
-    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))
-  </op>
+  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))  </op>
 </test>
 </case>