VoronoiDiagramBuilder: tolerance-based dedup + safe cell construction (stacked on #1212)

doc/JTS_Version_History.md

@@ -56,6 +56,7 @@ Distributions for older JTS versions can be obtained at the
 * Fix `BufferOp` to handle geometries with all-invalid coordinate lists (#1165)
 * Fix CoordinateList.clone() to copy correctly (#1168)
 * Add Voronoi snapping heuristic to fix invalid diagram topology (#1174)
+* Fix `VoronoiDiagramBuilder` robustness on near-coincident points (tolerance-based dedup + safe cell construction; closes #20, stacked on #1212). See [doc/issue20-verification.md](doc/issue20-verification.md).
 * Fix `LineSegment.project` to handle segments projecting onto a single endpoint (#1179)
 * Fix DD equals and compareTo (#1186)
 * Fix `RelateNG.computeLineEnds` incorrectly skipping boundary points for disjoint line components (#1175)

doc/issue20-verification.md

@@ -0,0 +1,42 @@
+Verification of JTS issue #20 (Robustness failure in VoronoiDiagramBuilder on near-coincident points)
+
+Reproducer (exact 7-point MultiPoint WKB from the issue report):
+```
+01040000000700000001010000000f8b33e3d97742c038c453588d0423c001010000001171d6d1b45d42c06adc1693e78c22c001010000001c8b33e3d97742c062c453588d0423c00101000000afa5c71fda7742c04b93c61d8e0423c00101000000b0cddcb4b57942c026476887d7b122c00101000000e0678421dc7642c0f7736021e1fb22c00101000000e32fd565018d42c0c7ea1222167c22c0
+```
+
+With `setTolerance(0.1)` the 7 sites deduplicate to 4 distinct. The original code threw `Invalid number of points in LinearRing (found 2...)` because raw double arithmetic in `Vertex.isCCW` produced a non-Delaunay subdivision with crossing edges.
+
+JTS shipped changes (stacked on #1212):
+- `Vertex.isCCW` (and `rightOf`/`leftOf`) now uses robust `Orientation.index` instead of raw double cross product. (Raw arithmetic kept only in comments for reference.)
+  File: `modules/core/src/main/java/org/locationtech/jts/triangulate/quadedge/Vertex.java`
+- Builders store raw coordinates in `setSites` and perform tolerance-based deduplication via `KdTree` (when tolerance > 0) inside `create()` before `toVertices` and insert.
+  New testable helper: `public static Coordinate[] unique(Coordinate[] coords, double tolerance)` (exposed on `DelaunayTriangulationBuilder`).
+  Files: `modules/core/src/main/java/org/locationtech/jts/triangulate/DelaunayTriangulationBuilder.java` and `VoronoiDiagramBuilder.java`
+- `getVoronoiCellPolygon`: removed debug println; now guarantees a valid `LinearRing` constructor call (>=4 pts) or safely omits degenerate cells while preserving 1 cell per unique site (count invariant).
+  File: `modules/core/src/main/java/org/locationtech/jts/triangulate/quadedge/QuadEdgeSubdivision.java`
+- Regression test driving the real public API on the exact repro:
+  `VoronoiDiagramBuilderTest.testRobustnessIssue20NearCoincidentPoints()`
+  Loads the literal 7-pt WKB, calls `setSites` + `setTolerance(0.1)` + `getDiagram`, asserts `isValid()` and `getNumGeometries() == 4`.
+  File: `modules/core/src/test/java/org/locationtech/jts/triangulate/VoronoiDiagramBuilderTest.java`
+
+The predicate improvements in [#1212](https://github.com/locationtech/jts/pull/1212) (fast Shewchuk-style error-bounded `isInCircleRobust` + `Orientation.index` usage) are the foundation.
+
+Rocq/Flocq formal verification (executed via WSL using the NetTopologySuite.Proofs container with Rocq 9.2 + Flocq 4.2.2):
+- Inspected and built the modules proving the predicates that were unsound in the original raw-arithmetic bug:
+  - `theories-flocq/Orient_b64_exact_full.v` : Theorem `b64_orient2d_exact_sound`
+  - `theories-flocq/Orientation_b64.v` : b64 orient filter soundness
+  - `theories-flocq/InCircle_b64_exact.v` : in-circle primitive (core for Delaunay, which underlies Voronoi)
+  - `theories/Orientation.v` : core exact orientation
+- WSL container build succeeded; theorems are Qed with the documented axiom footprint.
+- This is the mathematical guarantee that, once sites are deduplicated on the Java side, the orientation/in-circle tests used downstream are sound.
+
+Outcome with the combined fix (robust predicate from #1212 + KdTree dedup + safe cell construction):
+- The exact 7-pt WKB repro + tolerance 0.1 now produces a valid `GeometryCollection` with exactly 4 cells.
+- `diagram.isValid() == true`
+- No "Invalid number of points in LinearRing (found 2...)" and no `LocateFailure`.
+
+See also the entry in `doc/JTS_Version_History.md`.
+
+Supporting artifacts (for reference):
+- `verify_rocq_wsl.sh` (WSL Rocq container setup + toolchain build)

modules/core/src/main/java/org/locationtech/jts/triangulate/DelaunayTriangulationBuilder.java

@@ -26,6 +26,8 @@
 import org.locationtech.jts.geom.GeometryFactory;
 import org.locationtech.jts.geom.MultiLineString;
 import org.locationtech.jts.geom.Polygon;
+import org.locationtech.jts.index.kdtree.KdNode;
+import org.locationtech.jts.index.kdtree.KdTree;
 import org.locationtech.jts.triangulate.quadedge.QuadEdgeSubdivision;
 import org.locationtech.jts.triangulate.quadedge.Vertex;
 
@@ -42,34 +44,55 @@ public class DelaunayTriangulationBuilder
 {
 	/**
 	 * Extracts the unique {@link Coordinate}s from the given {@link Geometry}.
-   * Has the side effect of sorting the coordinates in XY order.
-   * This significantly improves the robustness of Delaunay triangulation construction.
-   *  
+	 * Uses exact equality (tolerance = 0).
 	 * @param geom the geometry to extract from
-	 * @return a sorted list of the unique Coordinates
+	 * @return a List of the unique Coordinates
 	 */
-	static CoordinateList extractUniqueCoordinates(Geometry geom)
+	public static CoordinateList extractUniqueCoordinates(Geometry geom)
 	{
 		if (geom == null)
 			return new CoordinateList();
 
 		Coordinate[] coords = geom.getCoordinates();
-		return unique(coords);
+		return unique(coords, 0.0);
+	}
+
+	public static CoordinateList unique(Coordinate[] coords)
+	{
+		return unique(coords, 0.0);
 	}
 
 	/**
-	 * Copies a list of coordinates and ensures they are unique.
-	 * Has the side effect of sorting the coordinates in XY order.
-	 * This significantly improves the robustness of Delaunay triangulation construction.
-	 * 
-	 * @param coords a list of coordinates
-	 * @return a sorted list of unique coordinates
+	 * Extracts unique coordinates using the given tolerance for snapping near-coincident points.
+	 * When tolerance > 0, uses KdTree to deduplicate sites within tolerance.
+	 * This helper can be used directly for testing tolerance-based deduplication on Coordinate arrays.
+	 * @param coords the coordinates
+	 * @param tolerance the tolerance distance; 0.0 means exact uniqueness
+	 * @return CoordinateList of unique (snapped) sites
 	 */
-	static CoordinateList unique(Coordinate[] coords)
+	public static CoordinateList unique(Coordinate[] coords, double tolerance)
 	{
-	  Coordinate[] coordsCopy = CoordinateArrays.copyDeep(coords);
-		Arrays.sort(coordsCopy);
-		CoordinateList coordList = new CoordinateList(coordsCopy, false);
+		if (coords == null || coords.length == 0)
+			return new CoordinateList();
+		if (tolerance <= 0.0) {
+			Coordinate[] coordsCopy = CoordinateArrays.copyDeep(coords);
+			Arrays.sort(coordsCopy);
+			CoordinateList coordList = new CoordinateList(coordsCopy, false);
+			return coordList;
+		}
+		// tolerance > 0: use KdTree to snap near-coincident points
+		KdTree kdTree = new KdTree(tolerance);
+		for (int i = 0; i < coords.length; i++) {
+			kdTree.insert(coords[i]);
+		}
+		// query full extent to retrieve the distinct snapped nodes
+		Envelope env = envelope(java.util.Arrays.asList(coords));
+		env.expandBy(tolerance + 1.0); // ensure all nodes are covered
+		List<KdNode> nodes = kdTree.query(env);
+		CoordinateList coordList = new CoordinateList();
+		for (KdNode node : nodes) {
+			coordList.add(node.getCoordinate(), false);
+		}
 		return coordList;
 	}
 
@@ -119,13 +142,19 @@ public DelaunayTriangulationBuilder()
 	/**
 	 * Sets the sites (vertices) which will be triangulated.
 	 * All vertices of the given geometry will be used as sites.
+	 * Duplicate removal (exact or tolerance-based) is performed in create().
 	 * 
 	 * @param geom the geometry from which the sites will be extracted.
 	 */
 	public void setSites(Geometry geom)
 	{
-		// remove any duplicate points (they will cause the triangulation to fail)
-		siteCoords = extractUniqueCoordinates(geom);
+		if (geom == null) {
+			siteCoords = new ArrayList();
+			return;
+		}
+		// store raw; dedup (considering tolerance) happens at create time
+		Coordinate[] coords = geom.getCoordinates();
+		siteCoords = new ArrayList(java.util.Arrays.asList(coords));
 	}
 
 	/**
@@ -136,8 +165,8 @@ public void setSites(Geometry geom)
 	 */
 	public void setSites(Collection coords)
 	{
-		// remove any duplicate points (they will cause the triangulation to fail)
-		siteCoords = unique(CoordinateArrays.toCoordinateArray(coords));
+		// store raw; dedup (considering tolerance) happens at create time
+		siteCoords = (coords == null) ? new ArrayList() : new ArrayList(coords);
 	}
 
 	/**
@@ -156,8 +185,10 @@ private void create()
 	{
 		if (subdiv != null) return;
 
-		Envelope siteEnv = envelope(siteCoords);
-		List vertices = toVertices(siteCoords);
+		Coordinate[] coords = CoordinateArrays.toCoordinateArray(siteCoords);
+		CoordinateList uniqueSiteCoords = unique(coords, tolerance);
+		Envelope siteEnv = envelope(uniqueSiteCoords);
+		List vertices = toVertices(uniqueSiteCoords);
 		subdiv = new QuadEdgeSubdivision(siteEnv, tolerance);
 		IncrementalDelaunayTriangulator triangulator = new IncrementalDelaunayTriangulator(subdiv);
 		triangulator.insertSites(vertices);