Add polygon triangulation algorithms and Tri data structure (#775)

* Add polygon triangulation algorithms and Tri data structure

Signed-off-by: Martin Davis <[email protected]>

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

@@ -35,6 +35,7 @@
 import org.locationtech.jts.index.strtree.Boundable;
 import org.locationtech.jts.index.strtree.GeometryItemDistance;
 import org.locationtech.jts.index.strtree.STRtree;
+//import org.locationtech.jts.triangulatepoly.VertexSequencePackedRtree;
 
 
 public class SpatialIndexFunctions
@@ -249,4 +250,20 @@ public static Geometry monotoneChains(Geometry geom) {
     }
     return geom.getFactory().buildGeometry(lines);
   }
+
+  /*
+  public static Geometry sprTreeBounds(Geometry geom)
+  {
+    Coordinate[] pts = geom.getCoordinates();
+    VertexSequencePackedRtree index = new VertexSequencePackedRtree(pts);
+    Envelope[] bounds = index.getBounds();
+    Geometry[] polys = new Geometry[bounds.length];
+    int i = 0;
+    for (Envelope env : bounds) {
+      polys[i++] = geom.getFactory().toGeometry(env);
+    }
+    return geom.getFactory().createGeometryCollection(polys);
+  }
+  */
+
 }

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

@@ -0,0 +1,32 @@
+/*
+ * Copyright (c) 2021 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.jtstest.function;
+
+import org.locationtech.jts.geom.Geometry;
+import org.locationtech.jts.triangulate.polygon.ConstrainedDelaunayTriangulator;
+import org.locationtech.jts.triangulate.polygon.PolygonTriangulator;
+
+
+public class TriangulatePolyFunctions 
+{
+  public static Geometry triangulate(Geometry geom)
+  {
+    return PolygonTriangulator.triangulate(geom);
+  }
+
+  public static Geometry constrainedDelaunay(Geometry geom)
+  {
+    return ConstrainedDelaunayTriangulator.triangulate(geom);
+  }
+
+
+}

modules/app/src/main/java/org/locationtech/jtstest/geomfunction/GeometryFunctionRegistry.java

@@ -65,6 +65,7 @@
 import org.locationtech.jtstest.function.SpatialIndexFunctions;
 import org.locationtech.jtstest.function.SpatialPredicateFunctions;
 import org.locationtech.jtstest.function.TriangleFunctions;
+import org.locationtech.jtstest.function.TriangulatePolyFunctions;
 import org.locationtech.jtstest.function.TriangulationFunctions;
 import org.locationtech.jtstest.function.ValidationFunctions;
 import org.locationtech.jtstest.function.WriterFunctions;
@@ -131,6 +132,7 @@ public static GeometryFunctionRegistry createTestBuilderRegistry()
     funcRegistry.add(SnappingFunctions.class);
     funcRegistry.add(SortingFunctions.class);
     funcRegistry.add(TriangulationFunctions.class);
+    funcRegistry.add(TriangulatePolyFunctions.class);
     funcRegistry.add(TriangleFunctions.class);
     funcRegistry.add(ValidationFunctions.class);
     funcRegistry.add(WriterFunctions.class);

modules/core/src/main/java/org/locationtech/jts/geom/Triangle.java

@@ -33,12 +33,9 @@ public class Triangle
    * Note: this implementation is not robust for angles very close to 90
    * degrees.
    * 
-   * @param a
-   *          a vertex of the triangle
-   * @param b
-   *          a vertex of the triangle
-   * @param c
-   *          a vertex of the triangle
+   * @param a a vertex of the triangle
+   * @param b a vertex of the triangle
+   * @param c a vertex of the triangle
    * @return true if the triangle is acute
    */
   public static boolean isAcute(Coordinate a, Coordinate b, Coordinate c)
@@ -51,6 +48,38 @@ public static boolean isAcute(Coordinate a, Coordinate b, Coordinate c)
       return false;
     return true;
   }
+
+  /**
+   * Tests whether a triangle is oriented counter-clockwise.
+   * 
+   * @param a a vertex of the triangle
+   * @param b a vertex of the triangle
+   * @param c a vertex of the triangle
+   * @return true if the triangle orientation is counter-clockwise
+   */
+  public static boolean isCCW(Coordinate a, Coordinate b, Coordinate c)
+  {
+    return Orientation.COUNTERCLOCKWISE == Orientation.index(a, b, c);
+  }
+
+  /**
+   * Tests whether a triangle intersects a point.
+   * 
+   * @param a a vertex of the triangle
+   * @param b a vertex of the triangle
+   * @param c a vertex of the triangle
+   * @param p the point to test
+   * @return true if the triangle intersects the point
+   */
+  public static boolean intersects(Coordinate a, Coordinate b, Coordinate c, Coordinate p)
+  {
+    int exteriorIndex = isCCW(a, b, c) ? 
+        Orientation.CLOCKWISE : Orientation.COUNTERCLOCKWISE;
+    if (exteriorIndex == Orientation.index(a, b, p)) return false;
+    if (exteriorIndex == Orientation.index(b, c, p)) return false;
+    if (exteriorIndex == Orientation.index(c, a, p)) return false;
+    return true;
+  }
 
   /**
    * Computes the line which is the perpendicular bisector of the line segment
@@ -437,8 +466,8 @@ public static double interpolateZ(Coordinate p, Coordinate v0, Coordinate v1,
     double u = (-c * dx + a * dy) / det;
     double z = v0.getZ() + t * (v1.getZ() - v0.getZ()) + u * (v2.getZ() - v0.getZ());
     return z;
-  }
-
+  }  
+  
   /**
    * The coordinates of the vertices of the triangle
    */
@@ -487,9 +516,18 @@ public Coordinate inCentre()
    */
   public boolean isAcute()
   {
-    return isAcute(this.p0, this.p1, this.p2);
+    return isAcute(p0, p1, p2);
   }
 
+  /**
+   * Tests whether this triangle is oriented counter-clockwise.
+   * 
+   * @return true if the triangle orientation is counter-clockwise
+   */
+  public boolean isCCW() {
+    return isCCW(p0, p1, p2);
+  }
+
   /**
    * Computes the circumcentre of this triangle. The circumcentre is the centre
    * of the circumcircle, the smallest circle which encloses the triangle. It is
@@ -507,7 +545,7 @@ public boolean isAcute()
    */
   public Coordinate circumcentre()
   {
-    return circumcentre(this.p0, this.p1, this.p2);
+    return circumcentre(p0, p1, p2);
   }
 
   /**
@@ -522,7 +560,7 @@ public Coordinate circumcentre()
    */
   public Coordinate centroid()
   {
-    return centroid(this.p0, this.p1, this.p2);
+    return centroid(p0, p1, p2);
   }
 
   /**
@@ -532,7 +570,7 @@ public Coordinate centroid()
    */
   public double longestSideLength()
   {
-    return longestSideLength(this.p0, this.p1, this.p2);
+    return longestSideLength(p0, p1, p2);
   }
 
   /**
@@ -545,7 +583,7 @@ public double longestSideLength()
    */
   public double area()
   {
-    return area(this.p0, this.p1, this.p2);
+    return area(p0, p1, p2);
   }
 
   /**
@@ -563,7 +601,7 @@ public double area()
    */
   public double signedArea()
   {
-    return signedArea(this.p0, this.p1, this.p2);
+    return signedArea(p0, p1, p2);
   }
 
   /**
@@ -574,7 +612,7 @@ public double signedArea()
    */
   public double area3D()
   {
-    return area3D(this.p0, this.p1, this.p2);
+    return area3D(p0, p1, p2);
   }
 
   /**

modules/core/src/main/java/org/locationtech/jts/math/MathUtil.java

@@ -47,6 +47,31 @@ public static int clamp(int x, int min, int max)
     return x;
   }
 
+  /**
+   * Clamps an integer to a given maximum limit.
+   * 
+   * @param x the value to clamp
+   * @param max the maximum value
+   * @return the clamped value
+   */
+  public static int clampMax(int x, int max)
+  {
+    if (x > max) return max;
+    return x;
+  }
+
+  /**
+   * Computes the ceiling function of the dividend of two integers.
+   * 
+   * @param num the numerator
+   * @param denom the denominator
+   * @return the ceiling of num / denom
+   */
+  public static int ceil(int num, int denom) {
+    int div = num / denom;
+    return div * denom >= num ? div : div + 1;
+  }
+
   private static final double LOG_10 = Math.log(10);
 
   /**

modules/core/src/main/java/org/locationtech/jts/triangulate/polygon/ConstrainedDelaunayTriangulator.java

@@ -0,0 +1,95 @@
+/*
+ * Copyright (c) 2021 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.polygon;
+
+import java.util.ArrayList;
+import java.util.List;
+
+import org.locationtech.jts.geom.Coordinate;
+import org.locationtech.jts.geom.Geometry;
+import org.locationtech.jts.geom.GeometryFactory;
+import org.locationtech.jts.geom.Polygon;
+import org.locationtech.jts.geom.util.PolygonExtracter;
+import org.locationtech.jts.triangulate.tri.Tri;
+import org.locationtech.jts.triangulate.tri.TriangulationBuilder;
+
+/**
+ * Computes the Constrained Delaunay Triangulation of polygons.
+ * The Constrained Delaunay Triangulation of a polygon is a set of triangles
+ * covering the polygon, with the maximum total interior angle over all 
+ * possible triangulations.  It provides the "best quality" triangulation
+ * of the polygon.
+ * <p>
+ * Holes are supported.
+ */
+public class ConstrainedDelaunayTriangulator {
+
+  /**
+   * Computes the Constrained Delaunay Triangulation of each polygon element in a geometry.
+   * 
+   * @param geom the input geometry
+   * @return a GeometryCollection of the computed triangle polygons
+   */
+  public static Geometry triangulate(Geometry geom) {
+    ConstrainedDelaunayTriangulator cdt = new ConstrainedDelaunayTriangulator(geom);
+    return cdt.compute();
+  }
+
+  private final GeometryFactory geomFact;
+  private final Geometry inputGeom;
+
+  /**
+   * Constructs a new Constrained Delaunay triangulator.
+   * 
+   * @param inputGeom the input geometry
+   */
+  public ConstrainedDelaunayTriangulator(Geometry inputGeom) {
+    geomFact = new GeometryFactory();
+    this.inputGeom = inputGeom;
+  }
+
+  private Geometry compute() {
+    List<Polygon> polys = PolygonExtracter.getPolygons(inputGeom);
+    List<Tri> triList = new ArrayList<Tri>();
+    for (Polygon poly : polys) {
+      List<Tri> polyTriList = triangulatePolygon(poly);
+      triList.addAll(polyTriList);
+    }
+    return Tri.toGeometry(triList, geomFact);
+  }
+
+  /**
+   * Computes the triangulation of a single polygon
+   * and returns it as a list of {@link Tri}s.
+   * 
+   * @param poly the input polygon
+   * @return list of Tris forming the triangulation
+   */
+  List<Tri> triangulatePolygon(Polygon poly) {
+    /**
+     * Normalize to ensure that shell and holes have canonical orientation.
+     * 
+     * TODO: perhaps better to just correct orientation of rings?
+     */
+    Polygon polyNorm = (Polygon) poly.norm();
+    Coordinate[] polyShell = PolygonHoleJoiner.join(polyNorm);
+    List<Tri> triList = PolygonEarClipper.triangulate(polyShell);
+
+    //long start = System.currentTimeMillis();
+    TriangulationBuilder.build(triList);
+    TriDelaunayImprover.improve(triList);
+    //System.out.println("swap used: " + (System.currentTimeMillis() - start) + " milliseconds");
+
+    return triList;
+  }
+
+}