Port locationtech/jts#775, a constrained

delaunay triangulation algorithm.

Adds PolygonTriangulator::triangulate() for a naive ear-clipped
triangulation, and ConstrainedDelaunayTriangulator::triangulate()
for a more aesthetically pleasing polygon triangulation.

Access for the CAPI is via GEOSConstrainedDelaunayTriangulation()

.github/workflows/ci.yml

@@ -60,15 +60,6 @@ jobs:
           cmake: 3.17.*,
           os: ubuntu-20.04,
         }
-        - {
-          compiler: clang++,
-          build_type: Release,
-          cxxstd: 11,
-          arch: 64,
-          packages: 'clang',
-          cmake: 3.8.*,
-          os: ubuntu-16.04
-        }
         - {
           compiler: clang++,
           build_type: Release,

NEWS

@@ -15,6 +15,9 @@ Changes in 3.10.0
           (Paul Ramsey, Martin Davis)
   - CAPI: GEOSWKBWriter_getFlavor, GEOSWKBWriter_setFlavor support
           outputting ISO or Extended WKB flavors (#466, Paul Ramsey)
+  - CAPI: GEOSConstrainedDelaunayTriangulation, builds a constrained triangulation
+          of an input Polygon or MultiPolygon, returning a GeometryCollection(Polygon)
+          of the triangles.
 
 - Fixes/Improvements:
   - Preserve ordering of lines in overlay results (Martin Davis)

capi/geos_c.cpp

@@ -1624,6 +1624,12 @@ extern "C" {
         return GEOSDelaunayTriangulation_r(handle, g, tolerance, onlyEdges);
     }
 
+    Geometry*
+    GEOSConstrainedDelaunayTriangulation(const Geometry* g)
+    {
+        return GEOSConstrainedDelaunayTriangulation_r(handle, g);
+    }
+
     Geometry*
     GEOSVoronoiDiagram(const Geometry* g, const Geometry* env, double tolerance, int onlyEdges)
     {

capi/geos_c.h.in

@@ -941,6 +941,11 @@ extern GEOSGeometry GEOS_DLL * GEOSDelaunayTriangulation_r(
     double tolerance,
     int onlyEdges);
 
+/** \see GEOSConstrainedDelaunayTriangulation */
+extern GEOSGeometry GEOS_DLL * GEOSConstrainedDelaunayTriangulation_r(
+    GEOSContextHandle_t handle,
+    const GEOSGeometry *g);
+
 /** \see GEOSVoronoiDiagram */
 extern GEOSGeometry GEOS_DLL * GEOSVoronoiDiagram_r(
     GEOSContextHandle_t extHandle,
@@ -2948,6 +2953,18 @@ extern GEOSGeometry GEOS_DLL * GEOSDelaunayTriangulation(
     double tolerance,
     int onlyEdges);
 
+/**
+* Return a constrained Delaunay triangulation of the vertices of the
+* given polygon(s).
+* For non-polygonal inputs, returns an empty geometry collection.
+*
+* \param g the input geometry whose rings will be used as input
+* \return A newly allocated geometry. NULL on exception.
+* Caller is responsible for freeing with GEOSGeom_destroy().
+*/
+extern GEOSGeometry GEOS_DLL * GEOSConstrainedDelaunayTriangulation(
+    const GEOSGeometry *g);
+
 /**
 * Returns the Voronoi polygons of the vertices of the given geometry.
 *

capi/geos_ts_c.cpp

@@ -85,6 +85,7 @@
 #include <geos/linearref/LengthIndexedLine.h>
 #include <geos/triangulate/DelaunayTriangulationBuilder.h>
 #include <geos/triangulate/VoronoiDiagramBuilder.h>
+#include <geos/triangulate/polygon/ConstrainedDelaunayTriangulator.h>
 #include <geos/util.h>
 #include <geos/util/IllegalArgumentException.h>
 #include <geos/util/Interrupt.h>
@@ -3706,6 +3707,16 @@ extern "C" {
         });
     }
 
+    Geometry*
+    GEOSConstrainedDelaunayTriangulation_r(GEOSContextHandle_t extHandle, const Geometry* g1)
+    {
+        using geos::triangulate::polygon::ConstrainedDelaunayTriangulator;
+
+        return execute(extHandle, [&]() -> Geometry* {
+            return ConstrainedDelaunayTriangulator::triangulate(g1).release();
+        });
+    }
+
     Geometry*
     GEOSVoronoiDiagram_r(GEOSContextHandle_t extHandle, const Geometry* g1, const Geometry* env, double tolerance,
                          int onlyEdges)

include/geos/geom/Coordinate.h

@@ -99,6 +99,8 @@ class GEOS_DLL Coordinate {
 
     bool equals2D(const Coordinate& other) const;
 
+    bool equals2D(const Coordinate& other, double tolerance) const;
+
     /// 2D only
     bool equals(const Coordinate& other) const;
 

include/geos/geom/Coordinate.inl

@@ -68,6 +68,18 @@ Coordinate::equals2D(const Coordinate& other) const
     return true;
 }
 
+INLINE bool
+Coordinate::equals2D(const Coordinate& other, double tolerance) const
+{
+    if (std::abs(x - other.x) > tolerance) {
+        return false;
+    }
+    if (std::abs(y - other.y) > tolerance) {
+        return false;
+    }
+    return true;
+}
+
 INLINE bool
 Coordinate::equals(const Coordinate& other) const
 {

include/geos/geom/CoordinateArraySequence.h

@@ -123,6 +123,8 @@ class GEOS_DLL CoordinateArraySequence : public CoordinateSequence {
 
     void expandEnvelope(Envelope& env) const override;
 
+    void closeRing();
+
     std::size_t getDimension() const override;
 
     void apply_rw(const CoordinateFilter* filter) override;

include/geos/geom/Envelope.h

@@ -522,6 +522,11 @@ GEOS_DLL bool operator==(const Envelope& a, const Envelope& b);
 
 GEOS_DLL bool operator!=(const Envelope& a, const Envelope& b);
 
+/// Strict weak ordering operator for Envelope
+/// This is the C++ equivalent of JTS's compareTo
+GEOS_DLL bool operator< (const Envelope& a, const Envelope& b);
+
+
 } // namespace geos::geom
 } // namespace geos
 

include/geos/geom/GeometryFactory.h

@@ -280,6 +280,9 @@ class GEOS_DLL GeometryFactory {
     std::unique_ptr<Polygon> createPolygon(std::unique_ptr<LinearRing> && shell,
                                            std::vector<std::unique_ptr<LinearRing>> && holes) const;
 
+    /// Construct a Polygon from a Coordinate vector, taking ownership of the vector
+    std::unique_ptr<Polygon> createPolygon(std::vector<Coordinate> && coords) const;
+
     /// Construct a Polygon with a deep-copy of given arguments
     Polygon* createPolygon(const LinearRing& shell,
                            const std::vector<LinearRing*>& holes) const;

include/geos/geom/Triangle.h

@@ -33,11 +33,10 @@ class GEOS_DLL Triangle {
     Coordinate p0, p1, p2;
 
     Triangle(const Coordinate& nP0, const Coordinate& nP1, const Coordinate& nP2)
-        :
-        p0(nP0),
-        p1(nP1),
-        p2(nP2)
-    {}
+        : p0(nP0)
+        , p1(nP1)
+        , p2(nP2) {}
+
 
     /** \brief
      * The inCentre of a triangle is the point which is equidistant
@@ -69,10 +68,71 @@ class GEOS_DLL Triangle {
     void circumcentre(Coordinate& resultPoint);
     void circumcentreDD(Coordinate& resultPoint);
 
+    /** Computes the circumcentre of a triangle. */
+    static const Coordinate circumcentre(
+        const Coordinate& p0, const Coordinate& p1, const Coordinate& p2);
+
     bool isIsoceles();
 
-    /// Computes the circumcentre of a triangle.
-    static const Coordinate circumcentre(const Coordinate& p0, const Coordinate& p1, const Coordinate& p2);
+    /**
+    * Tests whether a triangle is acute. A triangle is acute if all interior
+    * angles are acute. This is a strict test - right triangles will return
+    * <tt>false</tt>. A triangle which is not acute is either right or obtuse.
+    * <p>
+    * 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
+    * @return true if the triangle is acute
+    */
+    static bool isAcute(const Coordinate& a, const Coordinate& b, const Coordinate& c);
+
+    /**
+    * 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
+    */
+    static bool isCCW(const Coordinate& a, const Coordinate& b, const Coordinate& 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
+    */
+    static bool intersects(const Coordinate& a, const Coordinate& b, const Coordinate& c,
+        const Coordinate& p);
+
+
+    /**
+    * Tests whether a triangle intersects a point.
+    * @param p the point to test
+    * @return true if the triangle intersects the point
+    */
+    bool intersects(const Coordinate& p) { return intersects(p0, p1, p2, p); };
+
+    /**
+    * Tests whether this triangle is oriented counter-clockwise.
+    * @return true if the triangle orientation is counter-clockwise
+    */
+    bool isCCW() { return isCCW(p0, p1, p2); };
+
+    /**
+    * Tests whether this triangle is acute.
+    * @return true if this triangle is acute
+    */
+    bool isAcute() { return isAcute(p0, p1, p2); };
+
+
 
 private:
 

include/geos/triangulate/polygon/ConstrainedDelaunayTriangulator.h

@@ -0,0 +1,102 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2020 Paul Ramsey <pramsey@cleverelephant.ca>
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU Lesser General Public Licence as published
+ * by the Free Software Foundation.
+ * See the COPYING file for more information.
+ *
+ **********************************************************************/
+
+#pragma once
+
+
+// Forward declarations
+namespace geos {
+namespace geom {
+class Geometry;
+class GeometryFactory;
+class Polygon;
+}
+namespace triangulate {
+namespace tri {
+class TriList;
+}
+}
+}
+
+using geos::geom::Geometry;
+using geos::geom::GeometryFactory;
+using geos::geom::Polygon;
+using geos::triangulate::tri::TriList;
+
+
+namespace geos {
+namespace triangulate {
+namespace polygon {
+
+
+/**
+ * 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.
+ */
+class GEOS_DLL ConstrainedDelaunayTriangulator {
+
+private:
+
+    // Members
+    const Geometry* inputGeom;
+    const GeometryFactory* 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
+    */
+    void triangulatePolygon(const Polygon* poly, TriList& triList);
+
+    std::unique_ptr<Geometry> compute();
+
+
+public:
+
+    /**
+    * Constructs a new triangulator.
+    *
+    * @param p_inputGeom the input geometry
+    */
+    ConstrainedDelaunayTriangulator(const Geometry* p_inputGeom)
+        : inputGeom(p_inputGeom)
+        , geomFact(p_inputGeom->getFactory())
+    {}
+
+    /**
+    * 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
+    */
+    static std::unique_ptr<Geometry> triangulate(const Geometry* geom);
+
+
+
+
+};
+
+
+
+} // namespace geos.triangulate.polygon
+} // namespace geos.triangulate
+} // namespace geos
+