Add Voronoi snapping heuristic to fix invalid diagram topology (#1174)

Author: dr-jts

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

@@ -1,5 +1,5 @@
 /*
- * Copyright (c) 2016 Vivid Solutions.
+ * Copyright (c) 2026 Martin Davis.
  *
  * All rights reserved. This program and the accompanying materials
  * are made available under the terms of the Eclipse Public License 2.0
@@ -17,11 +17,14 @@
 
 import org.locationtech.jts.geom.Coordinate;
 import org.locationtech.jts.geom.CoordinateArrays;
+import org.locationtech.jts.geom.CoordinateList;
 import org.locationtech.jts.geom.Envelope;
 import org.locationtech.jts.geom.Geometry;
 import org.locationtech.jts.geom.GeometryCollection;
 import org.locationtech.jts.geom.GeometryFactory;
+import org.locationtech.jts.geom.LinearRing;
 import org.locationtech.jts.geom.Polygon;
+import org.locationtech.jts.noding.snap.SnappingPointIndex;
 import org.locationtech.jts.triangulate.quadedge.QuadEdgeSubdivision;
 
 
@@ -44,7 +47,13 @@
  */
 public class VoronoiDiagramBuilder 
 {
-	private Collection siteCoords;
+  /**
+   * A very small factor which detects short Voronoi cell segments 
+   * which might be caused by nearly-cocircular site circumcentres.
+   */
+	private static final double SHORT_SEG_TOLERANCE_FACTOR = 1.0e-10;
+	
+  private Collection siteCoords;
 	private double tolerance = 0.0;
 	private QuadEdgeSubdivision subdiv = null;
 	private Envelope clipEnv = null;
@@ -160,39 +169,141 @@ public Geometry getDiagram(GeometryFactory geomFact)
 		create();
 		Geometry polys = subdiv.getVoronoiDiagram(geomFact);
 
+    //System.out.println(polys);
+		//System.out.println( subdiv.getTriangles(true, geomFact) );
+		
 		/*
-    System.out.println(polys);
-		Geometry tris = subdiv.getTriangles(true, geomFact);
-		System.out.println(tris);
 		if (! subdiv.isFrameDelaunay()) {
-		  throw new IllegalStateException("Triangulation frame is not Delaunay");
-		}
+      throw new IllegalStateException("Triangulation frame is not Delaunay");
+    }
 		//*/
-		
-		//-- clip polys to diagramEnv
-		return clipGeometryCollection(polys, diagramEnv);
-	}
 
+    Geometry polysClean = clean(polys);
+    //Geometry polysClean = polys;   //  TESTING ONLY
+    
+		//-- clip cell polygons to diagram boundary
+		return clipGeometryCollection(polysClean, diagramEnv);
+	}
+	
   private static Geometry clipGeometryCollection(Geometry geom, Envelope clipEnv)
-	{
-		Geometry clipPoly = geom.getFactory().toGeometry(clipEnv);
-		List clipped = new ArrayList();
-		for (int i = 0; i < geom.getNumGeometries(); i++) {
-			Geometry g = geom.getGeometryN(i);
-			Geometry result = null;
-			// don't clip unless necessary
-			if (clipEnv.contains(g.getEnvelopeInternal()))
-					result = g;
-			else if (clipEnv.intersects(g.getEnvelopeInternal())) {
-				result = clipPoly.intersection(g);
-				// keep vertex key info
-				result.setUserData(g.getUserData());
-			}
+  {
+    Geometry clipPoly = geom.getFactory().toGeometry(clipEnv);
+    List<Geometry> clipped = new ArrayList<Geometry>();
+    for (int i = 0; i < geom.getNumGeometries(); i++) {
+      Geometry g = geom.getGeometryN(i);
+      Geometry result = null;
+      // don't clip unless necessary
+      if (clipEnv.contains(g.getEnvelopeInternal()))
+          result = g;
+      else if (clipEnv.intersects(g.getEnvelopeInternal())) {
+        result = clipPoly.intersection(g);
+        // keep vertex key info
+        result.setUserData(g.getUserData());
+      }
 
-			if (result != null && ! result.isEmpty()) {
-				clipped.add(result);
-			}
+      if (result != null && ! result.isEmpty()) {
+        clipped.add(result);
+      }
+    }
+    return geom.getFactory().createGeometryCollection(GeometryFactory.toGeometryArray(clipped));
+  }
+  
+	/**
+	 * Cleans diagram polygons to fix invalid topology caused by robustness errors,
+	 * 
+	 * @param polys a GeometryCollection containing the raw polygons for the diagram
+	 * @return the clean polygons
+	 */
+  private Geometry clean(Geometry polys) {
+    /**
+		 * Check for a diagram polygon with a very short edge which is invalid.
+		 * This can indicate invalid diagram topology caused by nearly cocircular input points.
+     * This is an efficient test which should not trigger on most typical datasets.
+     * 
+     * If found, snap the polygons to fix the topology.
+		 * This is a heuristic fix, but should generally restore correct topology
+		 * with very little effect on the diagram geometry.
+		 * 
+		 * See https://github.com/locationtech/jts/issues/1171
+		 */
+		double segmentLenTolerance = SHORT_SEG_TOLERANCE_FACTOR * diagramEnv.getDiameter();
+		if (hasInvalidPolygonWithShortEdge(polys, segmentLenTolerance)) {
+		  //System.out.println("SNAPPING!");
+		  Geometry polysSnap = snap(polys, segmentLenTolerance);
+		  return polysSnap;
 		}
-		return geom.getFactory().createGeometryCollection(GeometryFactory.toGeometryArray(clipped));
-	}
+    return polys;
+  }
+
+  /**
+   * Tests for a polygon with a very short edge which is invalid.
+   * This check is efficient for valid input,
+   * since that is unlikely to contain very short edges.
+   * 
+   * @param polys
+   * @param segmentLenTolerance
+   * @return true if a short edge in an invalid polygon is found
+   */
+  private static boolean hasInvalidPolygonWithShortEdge(Geometry polys, double segmentLenTolerance) {
+    for (int i = 0; i < polys.getNumGeometries(); i++) {
+      Polygon poly = (Polygon) polys.getGeometryN(i);
+      if (hasShortSegment(poly, segmentLenTolerance)) {
+        if (! poly.isValid())
+          return true;
+      }
+    }
+    return false;
+  }
+  
+  /**
+   * Tests if a polygon shell contains a short edge.
+   * 
+   * @param poly a polygon
+   * @param segmentLenTolerance the minimum segment length
+   * @return true if the polygon has a short edge
+   */
+  private static boolean hasShortSegment(Polygon poly, double segmentLenTolerance) {
+    LinearRing ring = poly.getExteriorRing();
+    Coordinate prev = ring.getCoordinateN(0);
+    for (int i = 1; i < ring.getNumPoints(); i++) {
+      Coordinate p = ring.getCoordinateN(i);
+      if (p.distance(prev) < segmentLenTolerance)
+        return true;
+      prev = p;
+    }
+    return false;
+  }
+
+  /**
+   * Snaps the vertices of a collection of polygons to eliminate short edges.
+   * Using a snapping map for all vertices ensures that adjacent polygons 
+   * match after snapping.
+   * 
+   * @param polys a GeometryCollection of single-ring polygons
+   * @param snapTolerance the snapping tolerance
+   * @return a GeometryCollection of snapped polygons
+   */
+  private static Geometry snap(Geometry polys, double snapTolerance) {
+    SnappingPointIndex snapMap = new SnappingPointIndex(snapTolerance);
+    List<Polygon> polysSnap = new ArrayList<Polygon>();
+    for (int i = 0; i < polys.getNumGeometries(); i++) {
+      Polygon polySnap = snapPolygon((Polygon) polys.getGeometryN(i), snapMap);
+      polysSnap.add(polySnap);
+    }
+    GeometryFactory geomFact = polys.getFactory();
+    return geomFact.createGeometryCollection(GeometryFactory.toGeometryArray(polysSnap));  
+  }
+
+  private static Polygon snapPolygon(Polygon poly, SnappingPointIndex snapMap) {
+    CoordinateList ptsSnap = new CoordinateList();
+    //-- voronoi polygons do not contain holes
+    Coordinate[] pts = poly.getExteriorRing().getCoordinates();
+    for (Coordinate pt : pts) {
+      Coordinate snapPt = snapMap.snap(pt);
+      ptsSnap.add(snapPt.copy(), false);
+    }
+    Polygon polySnap = poly.getFactory().createPolygon(ptsSnap.toCoordinateArray());
+    return polySnap;
+  }
+
 }

modules/core/src/test/java/org/locationtech/jts/triangulate/VoronoiChecker.java

@@ -0,0 +1,71 @@
+/*
+ * Copyright (c) 2026 Martin Davis.
+ *
+ * All rights reserved. This program and the accompanying materials
+ * are made available under the terms of the Eclipse Public License 2.0
+ * and Eclipse Distribution License v. 1.0 which accompanies this distribution.
+ * The Eclipse Public License is available at http://www.eclipse.org/legal/epl-v20.html
+ * and the Eclipse Distribution License is available at
+ *
+ * http://www.eclipse.org/org/documents/edl-v10.php.
+ */
+package org.locationtech.jts.triangulate;
+
+import org.locationtech.jts.algorithm.Orientation;
+import org.locationtech.jts.geom.Coordinate;
+import org.locationtech.jts.geom.Geometry;
+import org.locationtech.jts.geom.Polygon;
+
+public class VoronoiChecker {
+  
+  public static boolean isValid(Geometry voronoiDiagram) {
+    boolean isValid = voronoiDiagram.isValid();
+    if (! isValid) return false;
+    
+    boolean isConvex = isConvex(voronoiDiagram);
+    if (! isConvex) return false;
+    
+    boolean isNonoverlapping = isNonOverlapping(voronoiDiagram);
+    if (! isNonoverlapping) return false;
+    
+    return true;
+  }
+  
+  private static boolean isConvex(Geometry voronoiDiagram) {
+    Geometry union = voronoiDiagram.union();
+    if (! (union instanceof Polygon)) {
+      return false;
+    }
+    return isConvex((Polygon) union);
+  }
+  
+  private static boolean isConvex(Polygon poly) {
+    Coordinate[] pts = poly.getCoordinates();
+    for (int i = 0; i < pts.length - 1; i++) {
+      int iprev = i - 1;
+      if (iprev < 0) iprev = pts.length - 2;
+      int inext = i + 1;
+      //-- orientation must be CLOCKWISE or COLLINEAR
+      boolean isConvex = Orientation.COUNTERCLOCKWISE != Orientation.index(pts[iprev], pts[i], pts[inext]);
+      if (! isConvex)
+        return false;
+    }
+    return true;
+  }
+  
+  private static final String INTERIOR_INTERSECTS = "T********";
+  
+  private static boolean isNonOverlapping(Geometry result) {
+    int n = result.getNumGeometries();
+    for (int i1 = 0; i1 < n; i1++) {
+      Polygon poly1 = (Polygon) result.getGeometryN(i1);
+      for (int i2 = i1 + 1; i2 < n; i2++) {
+        Polygon poly2 = (Polygon) result.getGeometryN(i2);
+        boolean isOverlapping = poly1.relate(poly2, INTERIOR_INTERSECTS);
+        if (isOverlapping)
+          return false;
+      }
+    }
+    return true;
+  }
+}