Port JTS Concave Hull

NEWS.md

@@ -3,6 +3,7 @@
 
 - New things:
   - OffsetCurve (GH-530, Paul Ramsey/Martin Davis)
+  - ConcaveHull (GH-549, Paul Ramsey/Martin Davis)
 
 - Fixes/Improvements:
   - Fix unaryUnion to avoid segfault with empty polygon (GH-501, Mike Taves)

capi/geos_c.cpp

@@ -468,6 +468,15 @@ extern "C" {
         return GEOSConvexHull_r(handle, g);
     }
 
+    Geometry*
+    GEOSConcaveHull(const Geometry* g,
+                    double ratio,
+                    unsigned int allowHoles)
+
+    {
+        return GEOSConcaveHull_r(handle, g, ratio, allowHoles);
+    }
+
     Geometry*
     GEOSMinimumRotatedRectangle(const Geometry* g)
     {

capi/geos_c.h.in

@@ -741,6 +741,13 @@ extern GEOSGeometry GEOS_DLL *GEOSConvexHull_r(
     GEOSContextHandle_t handle,
     const GEOSGeometry* g);
 
+/** \see GEOSConcaveHull */
+extern GEOSGeometry GEOS_DLL *GEOSConcaveHull_r(
+    GEOSContextHandle_t handle,
+    const GEOSGeometry* g,
+    double ratio,
+    unsigned int allowHoles);
+
 /** \see GEOSMinimumRotatedRectangle */
 extern GEOSGeometry GEOS_DLL *GEOSMinimumRotatedRectangle_r(
     GEOSContextHandle_t handle,
@@ -2575,6 +2582,31 @@ extern GEOSGeometry GEOS_DLL *GEOSBoundary(const GEOSGeometry* g);
 */
 extern GEOSGeometry GEOS_DLL *GEOSConvexHull(const GEOSGeometry* g);
 
+/**
+* Returns "concave hull" of a geometry. The concave hull is fully
+* contained within the convex hull and also contains all the
+* points of the input, but in a smaller area.
+* The area ratio is the ratio
+* of the area of the convex hull and the concave hull. Frequently
+* used to convert a multi-point into a polygonal area.
+* that contains all the points in the input Geometry
+* \param g The input geometry
+* \param ratio The ratio value, between 0 and 1.
+* \param allowHoles When non-zero, the polygonal output may contain holes.
+* \return A newly allocated geometry of the concave hull. NULL on exception.
+* The area ratio is the ratio of the concave hull area to the convex hull area.
+* 1 produces the convex hull; 0 produces maximum concaveness.
+* The Length Ratio is a fraction determining the length of the longest
+* edge in the computed hull.  1 produces the convex hull;
+* 0 produces a hull with maximum concaveness
+* Caller is responsible for freeing with GEOSGeom_destroy().
+* \see geos::algorithm::hull::ConcaveHull
+*/
+extern GEOSGeometry GEOS_DLL *GEOSConcaveHull(
+    const GEOSGeometry* g,
+    double ratio,
+    unsigned int allowHoles);
+
 /**
 * Returns the minimum rotated rectangular POLYGON which encloses
 * the input geometry. The rectangle has width equal to the

capi/geos_ts_c.cpp

@@ -55,6 +55,7 @@
 #include <geos/algorithm/construct/LargestEmptyCircle.h>
 #include <geos/algorithm/distance/DiscreteHausdorffDistance.h>
 #include <geos/algorithm/distance/DiscreteFrechetDistance.h>
+#include <geos/algorithm/hull/ConcaveHull.h>
 #include <geos/simplify/DouglasPeuckerSimplifier.h>
 #include <geos/simplify/TopologyPreservingSimplifier.h>
 #include <geos/noding/GeometryNoder.h>
@@ -168,6 +169,7 @@ using geos::io::GeoJSONWriter;
 
 using geos::algorithm::distance::DiscreteFrechetDistance;
 using geos::algorithm::distance::DiscreteHausdorffDistance;
+using geos::algorithm::hull::ConcaveHull;
 
 using geos::operation::buffer::BufferBuilder;
 using geos::operation::buffer::BufferParameters;
@@ -1211,6 +1213,21 @@ extern "C" {
         });
     }
 
+    Geometry*
+    GEOSConcaveHull_r(GEOSContextHandle_t extHandle,
+        const Geometry* g1,
+        double ratio,
+        unsigned int allowHoles)
+    {
+        return execute(extHandle, [&]() {
+            ConcaveHull hull(g1);
+            hull.setMaximumEdgeLengthRatio(ratio);
+            hull.setHolesAllowed(allowHoles);
+            std::unique_ptr<Geometry> g3 = hull.getHull();
+            g3->setSRID(g1->getSRID());
+            return g3.release();
+        });
+    }
 
     Geometry*
     GEOSMinimumRotatedRectangle_r(GEOSContextHandle_t extHandle, const Geometry* g)

include/geos/algorithm/hull/ConcaveHull.h

@@ -0,0 +1,252 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2021 Paul Ramsey <[email protected]>
+ *
+ * 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
+
+#include <geos/geom/Triangle.h>
+#include <geos/triangulate/tri/Tri.h>
+#include <geos/triangulate/tri/TriList.h>
+#include <geos/triangulate/quadedge/TriangleVisitor.h>
+#include <geos/algorithm/hull/HullTri.h>
+
+#include <queue>
+#include <deque>
+
+namespace geos {
+namespace geom {
+class Coordinate;
+class Geometry;
+class GeometryFactory;
+}
+namespace triangulate {
+namespace quadedge {
+class Quadedge;
+class QuadEdgeSubdivision;
+}
+}
+}
+
+using geos::geom::Coordinate;
+using geos::geom::Geometry;
+using geos::geom::GeometryFactory;
+using geos::geom::Triangle;
+using geos::triangulate::quadedge::QuadEdge;
+using geos::triangulate::quadedge::QuadEdgeSubdivision;
+using geos::triangulate::quadedge::TriangleVisitor;
+using geos::triangulate::tri::Tri;
+using geos::triangulate::tri::TriList;
+
+namespace geos {
+namespace algorithm { // geos::algorithm
+namespace hull {      // geos::algorithm::hull
+
+
+typedef std::priority_queue<HullTri*, std::vector<HullTri*>, HullTri::HullTriCompare> HullTriQueue;
+
+
+/**
+* Constructs a concave hull of a set of points.
+* The hull is constructed by eroding the Delaunay Triangulation of the points
+* until specified target criteria are reached.
+* The target criteria are:
+*
+*  * Maximum Edge Length Ratio - determine the Maximum Edge Length
+*    as a fraction of the difference between the longest and shortest edge lengths
+*    in the Delaunay Triangulation. This normalizes the Maximum Edge Length to be scale-independent.
+*  * Maximum Area Ratio - the ratio of the concave hull area to the convex
+*    hull area will be no larger than this value
+*
+* The preferred criterium is the Maximum Edge Length Ratio, since it is
+* scale-free and local (so that no assumption needs to be made about the
+* total amount of concavity present.
+*
+* Other length criteria can be used by setting the Maximum Edge Length.
+* For example, use a length relative  to the longest edge length
+* in the Minimum Spanning Tree of the point set.
+* Or, use a length derived from the uniformGridEdgeLength() value.
+*
+* The computed hull is always a single connected geom::Polygon
+* (unless it is degenerate, in which case it will be a geom::Point or a geom::LineString).
+* This constraint may cause the concave hull to fail to meet the target criteria.
+*
+* Optionally the concave hull can be allowed to contain holes.
+* Note that this may be substantially slower than not permitting holes,
+* and it can produce results of lower quality.
+*
+* @author mdavis
+*/
+class GEOS_DLL ConcaveHull {
+
+public:
+
+    ConcaveHull(const Geometry* geom)
+        : inputGeometry(geom)
+        , maxEdgeLength(0.0)
+        , maxEdgeLengthRatio(-1.0)
+        , isHolesAllowed(false)
+        , geomFactory(geom->getFactory())
+        {};
+
+    /**
+    * Computes the approximate edge length of
+    * a uniform square grid having the same number of
+    * points as a geometry and the same area as its convex hull.
+    * This value can be used to determine a suitable length threshold value
+    * for computing a concave hull.
+    * A value from 2 to 4 times the uniform grid length
+    * seems to produce reasonable results.
+    *
+    * @param geom a geometry
+    * @return the approximate uniform grid length
+    */
+    static double uniformEdgeLength(const Geometry* geom);
+
+    /**
+    * Computes the concave hull of the vertices in a geometry
+    * using the target criteria of maximum edge length.
+    *
+    * @param geom the input geometry
+    * @param maxLength the target maximum edge length
+    * @return the concave hull
+    */
+    static std::unique_ptr<Geometry> concaveHullByLength(
+        const Geometry* geom, double maxLength);
+
+    static std::unique_ptr<Geometry> concaveHullByLength(
+        const Geometry* geom, double maxLength, bool isHolesAllowed);
+
+    /**
+    * Computes the concave hull of the vertices in a geometry
+    * using the target criteria of maximum edge length ratio.
+    * The edge length ratio is a fraction of the length difference
+    * between the longest and shortest edges
+    * in the Delaunay Triangulation of the input points.
+    *
+    * @param geom the input geometry
+    * @param lengthRatio the target edge length factor
+    * @return the concave hull
+    */
+    static std::unique_ptr<Geometry> concaveHullByLengthRatio(
+        const Geometry* geom, double lengthRatio);
+
+    /**
+    * Computes the concave hull of the vertices in a geometry
+    * using the target criterion of maximum edge length factor,
+    * and optionally allowing holes.
+    * The edge length factor is a fraction of the length difference
+    * between the longest and shortest edges
+    * in the Delaunay Triangulation of the input points.
+    *
+    * @param geom the input geometry
+    * @param lengthRatio the target maximum edge length
+    * @param isHolesAllowed whether holes are allowed in the result
+    * @return the concave hull
+    */
+    static std::unique_ptr<Geometry> concaveHullByLengthRatio(
+        const Geometry* geom, double lengthRatio, bool isHolesAllowed);
+
+    /**
+    * Sets the target maximum edge length for the concave hull.
+    * The length value must be zero or greater.
+    *
+    *  * The value 0.0 produces the concave hull of smallest area
+    *    that is still connected.
+    *  * Larger values produce less concave results.
+    *    A value equal or greater than the longest Delaunay Triangulation edge length
+    *    produces the convex hull.
+    *
+    * The uniformGridEdgeLength value may be used as
+    * the basis for estimating an appropriate target maximum edge length.
+    *
+    * @param edgeLength a non-negative length
+    */
+    void setMaximumEdgeLength(double edgeLength);
+
+    /**
+    * Sets the target maximum edge length ratio for the concave hull.
+    * The edge length ratio is a fraction of the length delta
+    * between the longest and shortest edges
+    * in the Delaunay Triangulation of the input points.
+    * A value of 1.0 produces the convex hull.
+    * A value of 0.0 produces a concave hull of minimum area
+    * that is still connected.
+    *
+    * @param edgeLengthRatio a length ratio value between 0 and 1
+    */
+    void setMaximumEdgeLengthRatio(double edgeLengthRatio);
+
+    /**
+    * Sets whether holes are allowed in the concave hull polygon.
+    *
+    * @param holesAllowed true if holes are allowed in the result
+    */
+    void setHolesAllowed(bool holesAllowed);
+
+    /**
+    * Gets the computed concave hull.
+    *
+    * @return the concave hull
+    */
+    std::unique_ptr<Geometry> getHull();
+
+
+private:
+
+    // Members
+    const Geometry* inputGeometry;
+    double maxEdgeLength;
+    double maxEdgeLengthRatio;
+    bool isHolesAllowed;
+    const GeometryFactory* geomFactory;
+
+    static double computeTargetEdgeLength(
+        TriList<HullTri>& triList,
+        double edgeLengthFactor);
+
+    void computeHull(TriList<HullTri>& triList);
+    void computeHullBorder(TriList<HullTri>& triList);
+    void createBorderQueue(HullTriQueue& queue, TriList<HullTri>& triList);
+
+    /**
+    * Adds a Tri to the queue.
+    * Only add tris with a single border edge.
+    * The ordering size is the length of the border edge.
+    *
+    * @param tri the Tri to add
+    * @param queue the priority queue
+    */
+    void addBorderTri(HullTri* tri, HullTriQueue& queue);
+    bool isBelowLengthThreshold(const HullTri* tri) const;
+    void computeHullHoles(TriList<HullTri>& triList);
+
+    static std::vector<HullTri*> findCandidateHoles(
+        TriList<HullTri>& triList, double minEdgeLen);
+
+    void removeHole(TriList<HullTri>& triList, HullTri* triHole);
+
+    bool isRemovableBorder(const HullTri* tri) const;
+    bool isRemovableHole(const HullTri* tri) const;
+
+    std::unique_ptr<Geometry> toGeometry(
+        TriList<HullTri>& triList,
+        const GeometryFactory* factory);
+
+
+};
+
+
+} // geos::algorithm::hull
+} // geos::algorithm
+} // geos
+