/src/geos/include/geos/algorithm/construct/MaximumInscribedCircle.h

Source
/**********************************************************************
 *
 * GEOS - Geometry Engine Open Source
 * http://geos.osgeo.org
 *
 * Copyright (c) 2020 Martin Davis
 * 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.
 *
 **********************************************************************
 *
 * Last port: algorithm/construct/MaximumInscribedCircle.java
 * https://github.com/locationtech/jts/commit/f0b9a808bdf8a973de435f737e37b7a221e231cb
 *
 **********************************************************************/

#pragma once

#include <geos/geom/Coordinate.h>
#include <geos/geom/Point.h>
#include <geos/geom/Envelope.h>
#include <geos/algorithm/locate/IndexedPointInAreaLocator.h>
#include <geos/operation/distance/IndexedFacetDistance.h>

#include <memory>
#include <queue>



namespace geos {
namespace geom {
class Coordinate;
class Envelope;
class Geometry;
class GeometryFactory;
class LineString;
class Point;
}
}

namespace geos {
namespace algorithm { // geos::algorithm
namespace construct { // geos::algorithm::construct

/**
 * Constructs the Maximum Inscribed Circle for a
 * polygonal Geometry, up to a specified tolerance.
 * The Maximum Inscribed Circle is determined by a point in the interior of the area
 * which has the farthest distance from the area boundary,
 * along with a boundary point at that distance.
 *
 * In the context of geography the center of the Maximum Inscribed Circle
 * is known as the **Pole of Inaccessibility**.
 * A cartographic use case is to determine a suitable point
 * to place a map label within a polygon.
 *
 * The radius length of the Maximum Inscribed Circle is a
 * measure of how "narrow" a polygon is. It is the
 * distance at which the negative buffer becomes empty.
 *
 * The class supports testing whether a polygon is "narrower"
 * than a specified distance via
 * isRadiusWithin(Geometry, double) or
 * isRadiusWithin(double).
 * Testing for the maximum radius is generally much faster
 * than computing the actual radius value, since short-circuiting
 * is used to limit the approximation iterations.
 *
 * The class supports polygons with holes and multipolygons.
 *
 * The implementation uses a successive-approximation technique
 * over a grid of square cells covering the area geometry.
 * The grid is refined using a branch-and-bound algorithm.
 * Point containment and distance are computed in a performant
 * way by using spatial indexes.
 *
 * Future Enhancements
 *
 *   * Support a polygonal constraint on placement of center point,
 *     for example to produce circle-packing constructions,
 *     or support multiple labels.
 *
 * @author Martin Davis
 *
 */
class GEOS_DLL MaximumInscribedCircle {

    using IndexedPointInAreaLocator = geos::algorithm::locate::IndexedPointInAreaLocator;
    using IndexedFacetDistance = geos::operation::distance::IndexedFacetDistance;

public:

    MaximumInscribedCircle(const geom::Geometry* polygonal, double tolerance);
    ~MaximumInscribedCircle() = default;

    /**
    * Gets the center point of the maximum inscribed circle
    * (up to the tolerance distance).
    *
    * @return the center point of the maximum inscribed circle
    */
    std::unique_ptr<geom::Point> getCenter();

    /**
    * Gets a point defining the radius of the Maximum Inscribed Circle.
    * This is a point on the boundary which is
    * nearest to the computed center of the Maximum Inscribed Circle.
    * The line segment from the center to this point
    * is a radius of the constructed circle, and this point
    * lies on the boundary of the circle.
    *
    * @return a point defining the radius of the Maximum Inscribed Circle
    */
    std::unique_ptr<geom::Point> getRadiusPoint();

    /**
    * Gets a line representing a radius of the Largest Empty Circle.
    *
    * @return a 2-point line from the center of the circle to a point on the edge
    */
    std::unique_ptr<geom::LineString> getRadiusLine();

    /**
    * Tests if the radius of the maximum inscribed circle
    * is no longer than the specified distance.
    * This method determines the distance tolerance automatically
    * as a fraction of the maxRadius value.
    * After this method is called the center and radius
    * points provide locations demonstrating where
    * the radius exceeds the specified maximum.
    *
    * @param maxRadius the (non-negative) radius value to test
    * @return true if the max in-circle radius is no longer than the max radius
    */
    bool isRadiusWithin(double maxRadius);

    /**
    * Computes the center point of the Maximum Inscribed Circle
    * of a polygonal geometry, up to a given tolerance distance.
    *
    * @param polygonal a polygonal geometry
    * @param tolerance the distance tolerance for computing the center point
    * @return the center point of the maximum inscribed circle
    */
    static std::unique_ptr<geom::Point> getCenter(const geom::Geometry* polygonal, double tolerance);

    /**
    * Computes a radius line of the Maximum Inscribed Circle
    * of a polygonal geometry, up to a given tolerance distance.
    *
    * @param polygonal a polygonal geometry
    * @param tolerance the distance tolerance for computing the center point
    * @return a line from the center to a point on the circle
    */
    static std::unique_ptr<geom::LineString> getRadiusLine(const geom::Geometry* polygonal, double tolerance);

    /**
    * Tests if the radius of the maximum inscribed circle
    * is no longer than the specified distance.
    * This method determines the distance tolerance automatically
    * as a fraction of the maxRadius value.
    *
    * @param polygonal a polygonal geometry
    * @param maxRadius the radius value to test
    * @return true if the max in-circle radius is no longer than the max radius
    */
    static bool isRadiusWithin(const geom::Geometry* polygonal, double maxRadius);

    /**
     * Computes the maximum number of iterations allowed.
     * Uses a heuristic based on the area of the input geometry
     * and the tolerance distance.
     * The number of tolerance-sized cells that cover the input geometry area
     * is computed, times a safety factor.
     * This prevents massive numbers of iterations and created cells
     * for casees where the input geometry has extremely small area
     * (e.g. is very thin).
     *
     * @param geom the input geometry
     * @param toleranceDist the tolerance distance
     * @return the maximum number of iterations allowed
     */
    static std::size_t computeMaximumIterations(const geom::Geometry* geom, double toleranceDist);

private:

    /* private members */
    const geom::Geometry* inputGeom;
    std::unique_ptr<geom::Geometry> inputGeomBoundary;
    double tolerance;
    IndexedFacetDistance indexedDistance;
    IndexedPointInAreaLocator ptLocator;
    const geom::GeometryFactory* factory;
    bool done;
    geom::CoordinateXY centerPt;
    geom::CoordinateXY radiusPt;
    double maximumRadius = -1;

    /* private constant */
    static constexpr double MAX_RADIUS_FRACTION = 0.0001;

    /* private methods */
    double distanceToBoundary(const geom::Point& pt);
    double distanceToBoundary(double x, double y);
    void compute();
    void computeApproximation();
    void createResult(const geom::CoordinateXY& c, const geom::CoordinateXY& r);

    /* private class */
    class Cell {
    private:
        static constexpr double SQRT2 = 1.4142135623730951;
        double x;
        double y;
        double hSize;
        double distance;
        double maxDist;

    public:
        Cell(double p_x, double p_y, double p_hSize, double p_distanceToBoundary)
            : x(p_x)
            , y(p_y)
            , hSize(p_hSize)
            , distance(p_distanceToBoundary)
            , maxDist(p_distanceToBoundary+(p_hSize*SQRT2))
        {};

        geom::Envelope getEnvelope() const
        {
            geom::Envelope env(x-hSize, x+hSize, y-hSize, y+hSize);
            return env;
        }

        double getMaxDistance() const
        {
            return maxDist;
        }
        double getDistance() const
        {
            return distance;
        }
        double getHSize() const
        {
            return hSize;
        }
        double getX() const
        {
            return x;
        }
        double getY() const
        {
            return y;
        }

        bool operator< (const Cell& rhs) const
        {
            return maxDist <  rhs.maxDist;
        }

        bool operator> (const Cell& rhs) const
        {
            return maxDist >  rhs.maxDist;
        }

        bool operator==(const Cell& rhs) const
        {
            return maxDist == rhs.maxDist;
        }

        /**
         * The Cell priority queue is sorted by the natural order of maxDistance.
         * std::priority_queue sorts with largest first,
         * which is what is needed for this algorithm.
         */
        using CellQueue = std::priority_queue<Cell>;
    };

    void createInitialGrid(const geom::Envelope* env, Cell::CellQueue& cellQueue);
    Cell createInteriorPointCell(const geom::Geometry* geom);

};


} // geos::algorithm::construct
} // geos::algorithm
} // geos

Coverage Report

Created: 2025-12-18 06:34

Line	Count	Source
1		/**********************************************************************
2		*
3		* GEOS - Geometry Engine Open Source
4		* http://geos.osgeo.org
5		*
6		* Copyright (c) 2020 Martin Davis
7		* Copyright (C) 2020 Paul Ramsey <pramsey@cleverelephant.ca>
8		*
9		* This is free software; you can redistribute and/or modify it under
10		* the terms of the GNU Lesser General Public Licence as published
11		* by the Free Software Foundation.
12		* See the COPYING file for more information.
13		*
14		**********************************************************************
15		*
16		* Last port: algorithm/construct/MaximumInscribedCircle.java
17		* https://github.com/locationtech/jts/commit/f0b9a808bdf8a973de435f737e37b7a221e231cb
18		*
19		**********************************************************************/
20
21		#pragma once
22
23		#include <geos/geom/Coordinate.h>
24		#include <geos/geom/Point.h>
25		#include <geos/geom/Envelope.h>
26		#include <geos/algorithm/locate/IndexedPointInAreaLocator.h>
27		#include <geos/operation/distance/IndexedFacetDistance.h>
28
29		#include <memory>
30		#include <queue>
31
32
33
34		namespace geos {
35		namespace geom {
36		class Coordinate;
37		class Envelope;
38		class Geometry;
39		class GeometryFactory;
40		class LineString;
41		class Point;
42		}
43		}
44
45		namespace geos {
46		namespace algorithm { // geos::algorithm
47		namespace construct { // geos::algorithm::construct
48
49		/**
50		* Constructs the Maximum Inscribed Circle for a
51		* polygonal Geometry, up to a specified tolerance.
52		* The Maximum Inscribed Circle is determined by a point in the interior of the area
53		* which has the farthest distance from the area boundary,
54		* along with a boundary point at that distance.
55		*
56		* In the context of geography the center of the Maximum Inscribed Circle
57		* is known as the Pole of Inaccessibility.
58		* A cartographic use case is to determine a suitable point
59		* to place a map label within a polygon.
60		*
61		* The radius length of the Maximum Inscribed Circle is a
62		* measure of how "narrow" a polygon is. It is the
63		* distance at which the negative buffer becomes empty.
64		*
65		* The class supports testing whether a polygon is "narrower"
66		* than a specified distance via
67		* isRadiusWithin(Geometry, double) or
68		* isRadiusWithin(double).
69		* Testing for the maximum radius is generally much faster
70		* than computing the actual radius value, since short-circuiting
71		* is used to limit the approximation iterations.
72		*
73		* The class supports polygons with holes and multipolygons.
74		*
75		* The implementation uses a successive-approximation technique
76		* over a grid of square cells covering the area geometry.
77		* The grid is refined using a branch-and-bound algorithm.
78		* Point containment and distance are computed in a performant
79		* way by using spatial indexes.
80		*
81		* Future Enhancements
82		*
83		* * Support a polygonal constraint on placement of center point,
84		* for example to produce circle-packing constructions,
85		* or support multiple labels.
86		*
87		* @author Martin Davis
88		*
89		*/
90		class GEOS_DLL MaximumInscribedCircle {
91
92		using IndexedPointInAreaLocator = geos::algorithm::locate::IndexedPointInAreaLocator;
93		using IndexedFacetDistance = geos::operation::distance::IndexedFacetDistance;
94
95		public:
96
97		MaximumInscribedCircle(const geom::Geometry* polygonal, double tolerance);
98	0	~MaximumInscribedCircle() = default;
99
100		/**
101		* Gets the center point of the maximum inscribed circle
102		* (up to the tolerance distance).
103		*
104		* @return the center point of the maximum inscribed circle
105		*/
106		std::unique_ptr<geom::Point> getCenter();
107
108		/**
109		* Gets a point defining the radius of the Maximum Inscribed Circle.
110		* This is a point on the boundary which is
111		* nearest to the computed center of the Maximum Inscribed Circle.
112		* The line segment from the center to this point
113		* is a radius of the constructed circle, and this point
114		* lies on the boundary of the circle.
115		*
116		* @return a point defining the radius of the Maximum Inscribed Circle
117		*/
118		std::unique_ptr<geom::Point> getRadiusPoint();
119
120		/**
121		* Gets a line representing a radius of the Largest Empty Circle.
122		*
123		* @return a 2-point line from the center of the circle to a point on the edge
124		*/
125		std::unique_ptr<geom::LineString> getRadiusLine();
126
127		/**
128		* Tests if the radius of the maximum inscribed circle
129		* is no longer than the specified distance.
130		* This method determines the distance tolerance automatically
131		* as a fraction of the maxRadius value.
132		* After this method is called the center and radius
133		* points provide locations demonstrating where
134		* the radius exceeds the specified maximum.
135		*
136		* @param maxRadius the (non-negative) radius value to test
137		* @return true if the max in-circle radius is no longer than the max radius
138		*/
139		bool isRadiusWithin(double maxRadius);
140
141		/**
142		* Computes the center point of the Maximum Inscribed Circle
143		* of a polygonal geometry, up to a given tolerance distance.
144		*
145		* @param polygonal a polygonal geometry
146		* @param tolerance the distance tolerance for computing the center point
147		* @return the center point of the maximum inscribed circle
148		*/
149		static std::unique_ptr<geom::Point> getCenter(const geom::Geometry* polygonal, double tolerance);
150
151		/**
152		* Computes a radius line of the Maximum Inscribed Circle
153		* of a polygonal geometry, up to a given tolerance distance.
154		*
155		* @param polygonal a polygonal geometry
156		* @param tolerance the distance tolerance for computing the center point
157		* @return a line from the center to a point on the circle
158		*/
159		static std::unique_ptr<geom::LineString> getRadiusLine(const geom::Geometry* polygonal, double tolerance);
160
161		/**
162		* Tests if the radius of the maximum inscribed circle
163		* is no longer than the specified distance.
164		* This method determines the distance tolerance automatically
165		* as a fraction of the maxRadius value.
166		*
167		* @param polygonal a polygonal geometry
168		* @param maxRadius the radius value to test
169		* @return true if the max in-circle radius is no longer than the max radius
170		*/
171		static bool isRadiusWithin(const geom::Geometry* polygonal, double maxRadius);
172
173		/**
174		* Computes the maximum number of iterations allowed.
175		* Uses a heuristic based on the area of the input geometry
176		* and the tolerance distance.
177		* The number of tolerance-sized cells that cover the input geometry area
178		* is computed, times a safety factor.
179		* This prevents massive numbers of iterations and created cells
180		* for casees where the input geometry has extremely small area
181		* (e.g. is very thin).
182		*
183		* @param geom the input geometry
184		* @param toleranceDist the tolerance distance
185		* @return the maximum number of iterations allowed
186		*/
187		static std::size_t computeMaximumIterations(const geom::Geometry* geom, double toleranceDist);
188
189		private:
190
191		/* private members */
192		const geom::Geometry* inputGeom;
193		std::unique_ptr<geom::Geometry> inputGeomBoundary;
194		double tolerance;
195		IndexedFacetDistance indexedDistance;
196		IndexedPointInAreaLocator ptLocator;
197		const geom::GeometryFactory* factory;
198		bool done;
199		geom::CoordinateXY centerPt;
200		geom::CoordinateXY radiusPt;
201		double maximumRadius = -1;
202
203		/* private constant */
204		static constexpr double MAX_RADIUS_FRACTION = 0.0001;
205
206		/* private methods */
207		double distanceToBoundary(const geom::Point& pt);
208		double distanceToBoundary(double x, double y);
209		void compute();
210		void computeApproximation();
211		void createResult(const geom::CoordinateXY& c, const geom::CoordinateXY& r);
212
213		/* private class */
214		class Cell {
215		private:
216		static constexpr double SQRT2 = 1.4142135623730951;
217		double x;
218		double y;
219		double hSize;
220		double distance;
221		double maxDist;
222
223		public:
224		Cell(double p_x, double p_y, double p_hSize, double p_distanceToBoundary)
225	0	: x(p_x)
226	0	, y(p_y)
227	0	, hSize(p_hSize)
228	0	, distance(p_distanceToBoundary)
229	0	, maxDist(p_distanceToBoundary+(p_hSize*SQRT2))
230	0	{};
231
232		geom::Envelope getEnvelope() const
233	0	{
234	0	geom::Envelope env(x-hSize, x+hSize, y-hSize, y+hSize);
235	0	return env;
236	0	}
237
238		double getMaxDistance() const
239	0	{
240	0	return maxDist;
241	0	}
242		double getDistance() const
243	0	{
244	0	return distance;
245	0	}
246		double getHSize() const
247	0	{
248	0	return hSize;
249	0	}
250		double getX() const
251	0	{
252	0	return x;
253	0	}
254		double getY() const
255	0	{
256	0	return y;
257	0	}
258
259		bool operator< (const Cell& rhs) const
260	0	{
261	0	return maxDist < rhs.maxDist;
262	0	}
263
264		bool operator> (const Cell& rhs) const
265	0	{
266	0	return maxDist > rhs.maxDist;
267	0	}
268
269		bool operator==(const Cell& rhs) const
270	0	{
271	0	return maxDist == rhs.maxDist;
272	0	}
273
274		/**
275		* The Cell priority queue is sorted by the natural order of maxDistance.
276		* std::priority_queue sorts with largest first,
277		* which is what is needed for this algorithm.
278		*/
279		using CellQueue = std::priority_queue<Cell>;
280		};
281
282		void createInitialGrid(const geom::Envelope* env, Cell::CellQueue& cellQueue);
283		Cell createInteriorPointCell(const geom::Geometry* geom);
284
285		};
286
287
288		} // geos::algorithm::construct
289		} // geos::algorithm
290		} // geos