Fix IncrementalDelaunayTriangulator to ensure triangulation boundary …

…is convex (#1004)

modules/app/src/main/java/org/locationtech/jtstest/function/TriangleFunctions.java

@@ -52,6 +52,15 @@ public static double circumradius(Geometry g)
       return Triangle.circumradius(pts[0], pts[1], pts[2]);
   }
 
+  public static Geometry circumcircle(Geometry g, int quadSegs)
+  {
+    Coordinate[] pts = trianglePts(g);
+    Coordinate cc = Triangle.circumcentreDD(pts[0], pts[1], pts[2]);
+    Geometry ccPt = g.getFactory().createPoint(cc);
+    double cr = Triangle.circumradius(pts[0], pts[1], pts[2]);
+    return ccPt.buffer(cr, quadSegs);
+  }
+
   public static Geometry circumcentreDD(Geometry g)
   {
     return GeometryMapper.map(g, 

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

@@ -15,6 +15,8 @@
 import java.util.Collection;
 import java.util.Iterator;
 
+import org.locationtech.jts.algorithm.Orientation;
+import org.locationtech.jts.geom.Coordinate;
 import org.locationtech.jts.triangulate.quadedge.LocateFailureException;
 import org.locationtech.jts.triangulate.quadedge.QuadEdge;
 import org.locationtech.jts.triangulate.quadedge.QuadEdgeSubdivision;
@@ -32,6 +34,7 @@ public class IncrementalDelaunayTriangulator
 {
 	private QuadEdgeSubdivision subdiv;
 	private boolean isUsingTolerance = false;
+  private boolean isForceConvex = true;
 
 	/**
 	 * Creates a new triangulator using the given {@link QuadEdgeSubdivision}.
@@ -46,6 +49,22 @@ public IncrementalDelaunayTriangulator(QuadEdgeSubdivision subdiv) {
 
 	}
 
+  /**
+   * Sets whether the triangulation is forced to have a convex boundary. Because
+   * of the use of a finite-size frame, this condition requires special logic to
+   * enforce. The default is true, since this is a requirement for some uses of
+   * Delaunay Triangulations (such as Concave Hull generation). However, forcing
+   * the triangulation boundary to be convex may cause the overall frame
+   * triangulation to be non-Delaunay. This can cause a problem for Voronoi
+   * generation, so the logic can be disabled via this method.
+   * 
+   * @param isForceConvex true if the triangulation boundary is forced to be
+   *                      convex
+   */
+  public void forceConvex(boolean isForceConvex) {
+    this.isForceConvex = isForceConvex;
+  }
+
 	/**
 	 * Inserts all sites in a collection. The inserted vertices <b>MUST</b> be
 	 * unique up to the provided tolerance value. (i.e. no two vertices should be
@@ -105,19 +124,89 @@ else if (subdiv.isOnEdge(e, v.getCoordinate())) {
 			e = base.oPrev();
 		} while (e.lNext() != startEdge);
 
-		// Examine suspect edges to ensure that the Delaunay condition
-		// is satisfied.
+		/**
+		 * Examine suspect edges to ensure that the Delaunay condition is satisfied.
+		 * If it is not, flip the edge and continue scanning.
+		 * 
+		 * Since the frame is not infinitely far away,
+		 * edges which touch the frame or are adjacent to it require special logic
+		 * to ensure the inner triangulation maintains a convex boundary.
+		 */
 		do {
-			QuadEdge t = e.oPrev();
-			if (t.dest().rightOf(e) && v.isInCircle(e.orig(), t.dest(), e.dest())) {
-				QuadEdge.swap(e);
-				e = e.oPrev();
-			} else if (e.oNext() == startEdge) {
-				return base; // no more suspect edges.
-			} else {
-				e = e.oNext().lPrev();
-			}
-		} while (true);
+       //-- general case - flip if vertex is in circumcircle
+      QuadEdge t = e.oPrev();
+      boolean doFlip = t.dest().rightOf(e) && v.isInCircle(e.orig(), t.dest(), e.dest());
+
+      if (isForceConvex) {
+        //-- special cases to ensure triangulation boundary is convex
+        if (isConcaveBoundary(e)) {
+          //-- flip if the triangulation boundary is concave
+          doFlip = true;
+        }
+        else if (isBetweenFrameAndInserted(e, v)) {
+          //-- don't flip if edge lies between the inserted vertex and a frame vertex
+          doFlip = false;
+        }
+      }
+
+      if (doFlip) {
+        //-- flip the edge within its quadrilateral
+        QuadEdge.swap(e);
+        e = e.oPrev();
+        continue;
+      }
+
+      if (e.oNext() == startEdge) {
+        return base; // no more suspect edges.
+      } else {
+        e = e.oNext().lPrev();
+      }
+    } while (true);
 	}
 
+	/**
+	 * Tests if a edge touching a frame vertex 
+	 * creates a concavity in the triangulation boundary.
+	 * 
+	 * @param e the edge to test
+	 * @return true if the triangulation boundary is concave at the edge
+	 */
+  private boolean isConcaveBoundary(QuadEdge e) {
+    if (subdiv.isFrameVertex(e.dest())) {
+      return isConcaveAtOrigin(e);
+    }
+    if (subdiv.isFrameVertex(e.orig())) {
+      return isConcaveAtOrigin(e.sym());
+    }
+    return false;
+  }
+
+  /**
+   * Tests if the quadrilateral surrounding an edge is concave at the edge origin.
+   * Used to determine if the triangulation boundary has a concavity.
+   * @param e
+   * @return
+   */
+  private static boolean isConcaveAtOrigin(QuadEdge e) {
+    Coordinate p = e.orig().getCoordinate();
+    Coordinate pp = e.oPrev().dest().getCoordinate();
+    Coordinate pn = e.oNext().dest().getCoordinate();
+    boolean isConcave = Orientation.COUNTERCLOCKWISE == Orientation.index(pp, pn, p);
+    return isConcave;
+  }
+
+  /**
+   * Edges whose adjacent triangles contain
+   * a frame vertex and the inserted vertex must not be flipped.
+   *
+   * @param e the edge to test
+   * @param vInsert the inserted vertex
+   * @return true if the edge is between the frame and inserted vertex
+   */
+  private boolean isBetweenFrameAndInserted(QuadEdge e, Vertex vInsert) {
+    Vertex v1 = e.oNext().dest();
+    Vertex v2 = e.oPrev().dest();
+    return (v1 == vInsert && subdiv.isFrameVertex(v2))
+        || (v2 == vInsert && subdiv.isFrameVertex(v1));
+  }
 }

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

@@ -125,9 +125,14 @@ private void create()
 		List vertices = DelaunayTriangulationBuilder.toVertices(siteCoords);
 		subdiv = new QuadEdgeSubdivision(diagramEnv, tolerance);
 		IncrementalDelaunayTriangulator triangulator = new IncrementalDelaunayTriangulator(subdiv);
+		/**
+		 * Avoid creating very narrow triangles along triangulation boundary.
+		 * These otherwise can cause malformed Voronoi cells.
+		 */
+		triangulator.forceConvex(false);
 		triangulator.insertSites(vertices);
 	}
-	
+  
 	/**
 	 * Gets the {@link QuadEdgeSubdivision} which models the computed diagram.
 	 * 
@@ -155,11 +160,20 @@ public Geometry getDiagram(GeometryFactory geomFact)
 		create();
 		Geometry polys = subdiv.getVoronoiDiagram(geomFact);
 
-		// clip polys to diagramEnv
+		/*
+    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");
+		}
+		//*/
+
+		//-- clip polys to diagramEnv
 		return clipGeometryCollection(polys, diagramEnv);
 	}
-	
-	private static Geometry clipGeometryCollection(Geometry geom, Envelope clipEnv)
+
+  private static Geometry clipGeometryCollection(Geometry geom, Envelope clipEnv)
 	{
 		Geometry clipPoly = geom.getFactory().toGeometry(clipEnv);
 		List clipped = new ArrayList();

modules/core/src/main/java/org/locationtech/jts/triangulate/quadedge/QuadEdgeSubdivision.java

@@ -80,7 +80,7 @@ public static void getTriangleEdges(QuadEdge startQE, QuadEdge[] triEdge) {
 
 	private final static double EDGE_COINCIDENCE_TOL_FACTOR = 1000;
 
-  private static final double FRAME_SIZE_FACTOR = 100.0;
+  private static final double FRAME_SIZE_FACTOR = 10.0;
 
 	// debugging only - preserve current subdiv statically
 	// private static QuadEdgeSubdivision currentSubdiv;
@@ -617,6 +617,24 @@ public List getPrimaryEdges(boolean includeFrame) {
 		return edges;
 	}
 
+  /**
+   * Gets the edges which touch frame vertices. The returned edges are oriented so
+   * that their origin is a frame vertex.
+   * 
+   * @return the edges which touch the frame
+   */
+  public List<QuadEdge> getFrameEdges() {
+    List<QuadEdge> edges = getPrimaryEdges(true);
+    List<QuadEdge> frameEdges = new ArrayList<QuadEdge>();
+    for (QuadEdge e : edges) {
+      if (isFrameEdge(e)) {
+        QuadEdge fe = isFrameVertex(e.orig()) ? e : e.sym();
+        frameEdges.add(fe);
+      }
+    }
+    return frameEdges;
+  }
+
   /**
    * A TriangleVisitor which computes and sets the 
    * circumcentre as the origin of the dual 
@@ -867,6 +885,25 @@ public Geometry getTriangles(GeometryFactory geomFact) {
 		return geomFact.createGeometryCollection(tris);
 	}
 
+	/**
+   * Gets the geometry for the triangles in a triangulated subdivision as a {@link GeometryCollection}
+   * of triangular {@link Polygon}s, optionally including the frame triangles.
+   * 
+ 	 * @param includeFrame true if the frame triangles should be included
+ 	 * @param geomFact the GeometryFactory to use
+   * @return a GeometryCollection of triangular Polygons
+	 */
+  public Geometry getTriangles(boolean includeFrame, GeometryFactory geomFact) {
+    List triPtsList = getTriangleCoordinates(includeFrame);
+    Polygon[] tris = new Polygon[triPtsList.size()];
+    int i = 0;
+    for (Iterator it = triPtsList.iterator(); it.hasNext();) {
+      Coordinate[] triPt = (Coordinate[]) it.next();
+      tris[i++] = geomFact.createPolygon(geomFact.createLinearRing(triPt));
+    }
+    return geomFact.createGeometryCollection(tris);
+  }
+
 	/**
 	 * Gets the cells in the Voronoi diagram for this triangulation.
 	 * The cells are returned as a {@link GeometryCollection} of {@link Polygon}s
@@ -957,4 +994,61 @@ public Polygon getVoronoiCellPolygon(QuadEdge qe, GeometryFactory geomFact)
     return cellPoly;
   }
 
+  /**
+   * Tests whether a subdivision is a valid Delaunay Triangulation.
+   * This is the case iff every edge is locally Delaunay, meaning that
+   * the apex of one adjacent triangle is not inside the circumcircle 
+   * of the other adjacent triangle.
+   * 
+   * @return true if the subdivision is Delaunay
+   */
+  public boolean isDelaunay() {
+    List<QuadEdge> edges = getPrimaryEdges(true);
+    for (QuadEdge e : edges) {
+      Vertex a0 = e.oPrev().dest();
+      Vertex a1 = e.oNext().dest();
+      boolean isDelaunay = ! a1.isInCircle(e.orig(), a0, e.dest());
+      if (! isDelaunay) {
+        /*
+        System.out.println(WKTWriter.toLineString(new Coordinate[] {
+            e.orig().getCoordinate(), a0.getCoordinate(), e.dest().getCoordinate()
+        }));
+        */
+        return false;
+      }
+    }
+    return true;
+  }
+
+  /**
+   * Tests whether the frame edges are Delaunay
+   * @return true if the frame edges are Delaunay
+   */
+  /*
+  public boolean isFrameDelaunay() {
+    List<QuadEdge> edges = getFrameEdges();
+    for (QuadEdge e : edges) {
+      Vertex a0 = e.oPrev().dest();
+      Vertex a1 = e.oNext().dest();
+      boolean isDelaunay = ! a1.isInCircle(e.orig(), a0, e.dest());
+      if (! isDelaunay) {
+
+        return false;
+      }
+    }
+    return true;
+  }
+  
+  public void makeFrameDelaunay() {
+    List<QuadEdge> edges = getFrameEdges();
+    for (QuadEdge e : edges) {
+      Vertex a0 = e.oPrev().dest();
+      Vertex a1 = e.oNext().dest();
+      boolean isDelaunay = ! a1.isInCircle(e.orig(), a0, e.dest());
+      if (! isDelaunay) {
+        QuadEdge.swap(e);
+      }
+    }
+  }
+  */
 }