Add LargestEmptyCircle handling of polygon obstacles (#939)

capi/geos_c.h.in

@@ -4129,8 +4129,7 @@ extern GEOSGeometry GEOS_DLL *GEOSMaximumInscribedCircle(
 * Constructs the "largest empty circle" (LEC) for a set of obstacle geometries
 * and within a polygonal boundary,
 * with accuracy to to a specified distance tolerance.
-* The obstacles are point and line geometries.
-* Polygonal obstacles are treated as linear features.
+* The obstacles may be any collection of points, lines and polygons.
 * The LEC is the largest circle whose interior does not intersect with any obstacle.
 * and which has its **center** inside the given boundary.
 * If no boundary is provided, the
@@ -4139,15 +4138,21 @@ extern GEOSGeometry GEOS_DLL *GEOSMaximumInscribedCircle(
 * the obstacles (up to the given distance tolerance).
 * The LEC is determined by the center point and a point indicating the circle radius
 * (which will lie on an obstacle).
+*
+* To compute an LEC which lies **wholly** within a polygonal boundary, 
+* include the boundary of the polygon(s) as a linear obstacle.
+*
 * The implementation uses a successive-approximation technique over a grid of square cells covering the obstacles and boundary.
 * The grid is refined using a branch-and-bound algorithm.  Point containment and distance are computed in a performant
 * way by using spatial indexes.
+*
 * Returns the LEC radius as a two-point linestring, with the start point at the center of the inscribed circle and the end
 * on the boundary of the circle.
+*
 * \param obstacles The geometries that the LEC must not cross
 * \param boundary The area to contain the LEC center (may be null or empty)
 * \param tolerance Stop the algorithm when the search area is smaller than this tolerance
-* \return A newly allocated geometry of the LEC radius. NULL on exception.
+* \return A newly allocated geometry of the LEC radius line. NULL on exception.
 * Caller is responsible for freeing with GEOSGeom_destroy().
 * \see geos::algorithm::construct::LargestEmptyCircle
 *

include/geos/algorithm/construct/IndexedDistanceToPoint.h

@@ -0,0 +1,88 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2023 Martin Davis <[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.
+ *
+ **********************************************************************
+ *
+ * Last port: algorithm/construct/IndexedDistanceToPoint.java
+ * https://github.com/locationtech/jts/commit/d92f783163d9440fcc10c729143787bf7b9fe8f9
+ *
+ **********************************************************************/
+
+#pragma once
+
+#include <geos/geom/CoordinateSequence.h>
+#include <geos/geom/Geometry.h>
+#include <geos/geom/Point.h>
+#include <geos/algorithm/construct/IndexedPointInPolygonsLocator.h>
+#include <geos/operation/distance/IndexedFacetDistance.h>
+
+using geos::geom::Geometry;
+using geos::geom::Point;
+using geos::operation::distance::IndexedFacetDistance;
+
+namespace geos {
+namespace algorithm { // geos::algorithm
+namespace construct { // geos::algorithm::construct
+
+/**
+ * \brief Computes the distance between a point and a geometry
+ * (which may be a collection containing any type of geometry).
+ * 
+ * Also computes the pair of points containing the input
+ * point and the nearest point on the geometry.
+ * 
+ * \author Martin Davis
+ */
+class GEOS_DLL IndexedDistanceToPoint {
+
+public:
+    /** 
+     * \brief Creates an instance to find the distance from points to a geometry.
+     * 
+     * \param geom the geometry to compute distances to
+     */
+    IndexedDistanceToPoint(const Geometry& geom);
+
+    /**
+     * \brief Computes the distance from the base geometry to the given point.
+     *
+     * \param pt the point to compute the distance to
+     *
+     * \return the computed distance
+     */
+    double distance(const Point& pt);
+
+    /**
+     * \brief Computes the nearest point on the geometry to the point.
+     * 
+     * The first location lies on the geometry, 
+     * and the second location is the provided point.
+     *
+     * \param pt the point to compute the nearest point for
+     *
+     * \return the points that are nearest
+     */
+    std::unique_ptr<geom::CoordinateSequence> nearestPoints(const Point& pt);
+
+private:
+    void init();
+
+    bool isInArea(const Point& pt);
+
+    //-- members
+    const Geometry& targetGeometry;
+    std::unique_ptr<operation::distance::IndexedFacetDistance> facetDistance;
+    std::unique_ptr<IndexedPointInPolygonsLocator> ptLocator;
+
+};
+
+}}}

include/geos/algorithm/construct/IndexedPointInPolygonsLocator.h

@@ -0,0 +1,79 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2023 Martin Davis <[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.
+ *
+ **********************************************************************
+ *
+ * Last port: algorithm/construct/IndexedDistanceToPoint.java
+ * https://github.com/locationtech/jts/commit/d92f783163d9440fcc10c729143787bf7b9fe8f9
+ *
+ **********************************************************************/
+
+#pragma once
+
+#include <geos/geom/Coordinate.h>
+#include <geos/geom/Geometry.h>
+#include <geos/geom/Location.h>
+#include <geos/index/strtree/TemplateSTRtree.h>
+#include <geos/algorithm/locate/IndexedPointInAreaLocator.h>
+
+using geos::geom::Geometry;
+using geos::geom::CoordinateXY;
+using geos::geom::Location;
+using geos::index::strtree::TemplateSTRtree;
+using geos::algorithm::locate::IndexedPointInAreaLocator;
+
+namespace geos {
+namespace algorithm { // geos::algorithm
+namespace construct { // geos::algorithm::construct
+
+/**
+ * \brief Determines the location of a point in the polygonal elements of a geometry.
+ * 
+ * Uses spatial indexing to provide efficient performance.
+ * 
+ * \author Martin Davis
+ */
+class GEOS_DLL IndexedPointInPolygonsLocator {
+
+public:
+    /** 
+     * \brief Creates an instance to locate a point in polygonal elements.
+     * 
+     * \param geom the geometry to locate in
+     */
+    IndexedPointInPolygonsLocator(const Geometry& geom);
+
+    /** \brief
+     * Determines the [Location](@ref geom::Location) of a point in 
+     * the polygonal elements of a
+     * [Geometry](@ref geom::Geometry).
+     *
+     * @param p the point to test
+     * @return the location of the point in the geometry
+     */
+    Location locate(const CoordinateXY* /*const*/ p);
+
+private:
+    void init();
+
+    // Declare type as noncopyable
+    IndexedPointInPolygonsLocator(const IndexedPointInPolygonsLocator& other) = delete;
+    IndexedPointInPolygonsLocator& operator=(const IndexedPointInPolygonsLocator& rhs) = delete;
+
+    //-- members
+    const Geometry& geom;
+    bool isInitialized;
+    TemplateSTRtree<IndexedPointInAreaLocator*> index;
+    std::vector<std::unique_ptr<IndexedPointInAreaLocator>> locators;
+};
+
+}}}

include/geos/algorithm/construct/LargestEmptyCircle.h

@@ -22,14 +22,13 @@
 #include <geos/geom/Coordinate.h>
 #include <geos/geom/Point.h>
 #include <geos/geom/Envelope.h>
+#include <geos/algorithm/construct/IndexedDistanceToPoint.h>
 #include <geos/algorithm/locate/IndexedPointInAreaLocator.h>
 #include <geos/operation/distance/IndexedFacetDistance.h>
 
 #include <memory>
 #include <queue>
 
-
-
 namespace geos {
 namespace geom {
 class Coordinate;
@@ -39,29 +38,33 @@ class GeometryFactory;
 class LineString;
 class Point;
 }
-namespace operation {
-namespace distance {
-class IndexedFacetDistance;
-}
-}
 }
 
+using geos::operation::distance::IndexedFacetDistance;
 
 namespace geos {
 namespace algorithm { // geos::algorithm
 namespace construct { // geos::algorithm::construct
 
 /**
 * Constructs the Largest Empty Circle for a set of obstacle geometries,
-* up to a specified tolerance. The obstacles are point and line geometries.
+* up to a specified tolerance. 
+* The obstacles may be any combination of point, linear and polygonal geometries.
 *
-* The Largest Empty Circle is the largest circle which has its center
-* in the convex hull of the obstacles (the boundary), and whose
-* interior does not intersect with any obstacle. The circle center
-* is the point in the interior of the boundary which has the
-* farthest distance from the obstacles (up to tolerance).
-* The circle is determined by the center point and a point lying
-* on an obstacle indicating the circle radius.
+* The Largest Empty Circle (LEC) is the largest circle 
+* whose interior does not intersect with any obstacle
+* and whose center lies within a polygonal boundary.
+* The circle center is the point in the interior of the boundary 
+* which has the farthest distance from the obstacles 
+* (up to the accuracy of the distance tolerance).
+* The circle itself is determined by the center point
+* and a point lying on an obstacle determining the circle radius.
+* 
+* The polygonal boundary may be supplied explicitly.
+* If it is not specified the convex hull of the obstacles is used as the boundary.
+* 
+* To compute an LEC which lies wholly within
+* a polygonal boundary, include the boundary of the polygon(s) as an obstacle.
 *
 * The implementation uses a successive-approximation technique
 * over a grid of square cells covering the obstacles and boundary.
@@ -77,19 +80,36 @@ class GEOS_DLL LargestEmptyCircle {
 
     /**
     * Creates a new instance of a Largest Empty Circle construction.
+    * The obstacles may be any collection of points, lines and polygons.
+    * The constructed circle center lies within the convex hull of the obstacles.
     *
-    * @param p_obstacles a geometry representing the obstacles (points and lines)
+    * @param p_obstacles a geometry representing the obstacles
     * @param p_tolerance the distance tolerance for computing the circle center point
     */
     LargestEmptyCircle(const geom::Geometry* p_obstacles, double p_tolerance);
+
+    /**
+    * Creates a new instance of a Largest Empty Circle construction,
+    * interior-disjoint to a set of obstacle geometries 
+    * and having its center within a polygonal boundary.
+    * The obstacles may be any collection of points, lines and polygons.
+    * If the boundary is null or empty the convex hull
+    * of the obstacles is used as the boundary.
+    *
+    * @param p_obstacles a geometry representing the obstacles
+    * @param p_boundary a polygonal geometry to contain the LEC center
+    * @param p_tolerance the distance tolerance for computing the circle center point
+    */
     LargestEmptyCircle(const geom::Geometry* p_obstacles, const geom::Geometry* p_boundary, double p_tolerance);
+
     ~LargestEmptyCircle() = default;
 
     /**
     * Computes the center point of the Largest Empty Circle
     * within a set of obstacles, up to a given tolerance distance.
+    * The obstacles may be any collection of points, lines and polygons.
     *
-    * @param p_obstacles a geometry representing the obstacles (points and lines)
+    * @param p_obstacles a geometry representing the obstacles
     * @param p_tolerance the distance tolerance for computing the center point
     * @return the center point of the Largest Empty Circle
     */
@@ -98,8 +118,9 @@ class GEOS_DLL LargestEmptyCircle {
     /**
     * Computes a radius line of the Largest Empty Circle
     * within a set of obstacles, up to a given distance tolerance.
+    * The obstacles may be any collection of points, lines and polygons.
     *
-    * @param p_obstacles a geometry representing the obstacles (points and lines)
+    * @param p_obstacles a geometry representing the obstacles
     * @param p_tolerance the distance tolerance for computing the center point
     * @return a line from the center of the circle to a point on the edge
     */
@@ -118,10 +139,10 @@ class GEOS_DLL LargestEmptyCircle {
     std::unique_ptr<geom::Geometry> boundary;
     const geom::GeometryFactory* factory;
     geom::Envelope gridEnv;
-    operation::distance::IndexedFacetDistance obstacleDistance;
     bool done;
-    std::unique_ptr<algorithm::locate::IndexedPointInAreaLocator> ptLocator;
-    std::unique_ptr<operation::distance::IndexedFacetDistance> boundaryDistance;
+    std::unique_ptr<algorithm::locate::IndexedPointInAreaLocator> boundaryPtLocater;
+    IndexedDistanceToPoint obstacleDistance;
+    std::unique_ptr<IndexedFacetDistance> boundaryDistance;
     geom::CoordinateXY centerPt;
     geom::CoordinateXY radiusPt;
 

include/geos/geom/util/PolygonalExtracter.h

@@ -0,0 +1,54 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2001-2002 Vivid Solutions Inc.
+ * Copyright (C) 2006 Refractions Research Inc.
+ *
+ * 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/export.h>
+#include <geos/geom/Geometry.h>
+#include <vector>
+
+namespace geos {
+namespace geom { // geos.geom
+namespace util { // geos.geom.util
+
+/**
+ * \brief Extracts the polygonal (Polygon and MultiPolygon) elements from a Geometry.
+ */
+class GEOS_DLL PolygonalExtracter {
+
+public:
+
+    /**
+     * Pushes the polygonal (Polygon and MultiPolygon) elements from a geometry into
+     * the provided vector.
+     * 
+     * @param geom the geometry to extract from
+     * @param polys the vector to add the polygonal elements to
+     */
+    static void getPolygonals(const Geometry& geom, std::vector<const Geometry*>& polys);
+
+private:
+
+    static void getPolygonals(const Geometry* geom, std::vector<const Geometry*>& polys);
+
+    // Declare type as noncopyable
+    PolygonalExtracter(const PolygonalExtracter& other) = delete;
+    PolygonalExtracter& operator=(const PolygonalExtracter& rhs) = delete;
+};
+
+} // namespace geos.geom.util
+} // namespace geos.geom
+} // namespace geos
+