Class MaximumInscribedCircle

1	/*
2	* Copyright (c) 2020 Martin Davis.
3	*
4	* All rights reserved. This program and the accompanying materials
5	* are made available under the terms of the Eclipse Public License 2.0
6	* and Eclipse Distribution License v. 1.0 which accompanies this distribution.
7	* The Eclipse Public License is available at http://www.eclipse.org/legal/epl-v20.html
8	* and the Eclipse Distribution License is available at
9	*
10	* http://www.eclipse.org/org/documents/edl-v10.php.
11	*/
12	package org.locationtech.jts.algorithm.construct;
13
14	import java.util.PriorityQueue;
15
16	import org.locationtech.jts.algorithm.locate.IndexedPointInAreaLocator;
17	import org.locationtech.jts.geom.Coordinate;
18	import org.locationtech.jts.geom.Envelope;
19	import org.locationtech.jts.geom.Geometry;
20	import org.locationtech.jts.geom.GeometryFactory;
21	import org.locationtech.jts.geom.LineString;
22	import org.locationtech.jts.geom.Location;
23	import org.locationtech.jts.geom.MultiPolygon;
24	import org.locationtech.jts.geom.Point;
25	import org.locationtech.jts.geom.Polygon;
26	import org.locationtech.jts.operation.distance.IndexedFacetDistance;
27
28	/**
29	* Constructs the Maximum Inscribed Circle for a
30	* polygonal {@link Geometry}, up to a specified tolerance.
31	* The Maximum Inscribed Circle is determined by a point in the interior of the area
32	* which has the farthest distance from the area boundary,
33	* along with a boundary point at that distance.
34	* <p>
35	* In the context of geography the center of the Maximum Inscribed Circle
36	* is known as the <b>Pole of Inaccessibility</b>.
37	* A cartographic use case is to determine a suitable point
38	* to place a map label within a polygon.
39	* <p>
40	* The radius length of the Maximum Inscribed Circle is a
41	* measure of how "narrow" a polygon is. It is the
42	* distance at which the negative buffer becomes empty.
43	* <p>
44	* The class supports polygons with holes and multipolygons.
45	* <p>
46	* The implementation uses a successive-approximation technique
47	* over a grid of square cells covering the area geometry.
48	* The grid is refined using a branch-and-bound algorithm.
49	* Point containment and distance are computed in a performant
50	* way by using spatial indexes.
51	*
52	* <h3>Future Enhancements</h3>
53	* <ul>
54	* <li>Support a polygonal constraint on placement of center
55	* </ul>
56	*
57	* @author Martin Davis
58	*
59	*/
60	public class MaximumInscribedCircle {
61
62	/**
63	* Computes the center point of the Maximum Inscribed Circle
64	* of a polygonal geometry, up to a given tolerance distance.
65	*
66	* @param polygonal a polygonal geometry
67	* @param tolerance the distance tolerance for computing the center point
68	* @return the center point of the maximum inscribed circle
69	*/
70	public static Point getCenter(Geometry polygonal, double tolerance) {
71	MaximumInscribedCircle mic = new MaximumInscribedCircle(polygonal, tolerance);
72	return mic.getCenter();
73	}
74
75	/**
76	* Computes a radius line of the Maximum Inscribed Circle
77	* of a polygonal geometry, up to a given tolerance distance.
78	*
79	* @param polygonal a polygonal geometry
80	* @param tolerance the distance tolerance for computing the center point
81	* @return a line from the center to a point on the circle
82	*/
83	public static LineString getRadiusLine(Geometry polygonal, double tolerance) {
84	MaximumInscribedCircle mic = new MaximumInscribedCircle(polygonal, tolerance);
85	return mic.getRadiusLine();
86	}
87
88	private Geometry inputGeom;
89	private double tolerance;
90
91	private GeometryFactory factory;
92	private IndexedPointInAreaLocator ptLocater;
93	private IndexedFacetDistance indexedDistance;
94	private Cell centerCell = null;
95	private Coordinate centerPt = null;
96	private Coordinate radiusPt;
97	private Point centerPoint;
98	private Point radiusPoint;
99
100	/**
101	* Creates a new instance of a Maximum Inscribed Circle computation.
102	*
103	* @param polygonal an areal geometry
104	* @param tolerance the distance tolerance for computing the centre point
105	*/
106	public MaximumInscribedCircle(Geometry polygonal, double tolerance) {
107	if (! (polygonal instanceof Polygon \|\| polygonal instanceof MultiPolygon)) {
108	throw new IllegalArgumentException("Input geometry must be a Polygon or MultiPolygon");
109	}
110	if (polygonal.isEmpty()) {
111	throw new IllegalArgumentException("Empty input geometry is not supported");
112	}
113
114	this.inputGeom = polygonal;
115	this.factory = polygonal.getFactory();
116	this.tolerance = tolerance;
117	ptLocater = new IndexedPointInAreaLocator(polygonal);
118	indexedDistance = new IndexedFacetDistance( polygonal.getBoundary() );
119	}
120
121	/**
122	* Gets the center point of the maximum inscribed circle
123	* (up to the tolerance distance).
124	*
125	* @return the center point of the maximum inscribed circle
126	*/
127	public Point getCenter() {
128	compute();
129	return centerPoint;
130	}
131
132	/**
133	* Gets a point defining the radius of the Maximum Inscribed Circle.
134	* This is a point on the boundary which is
135	* nearest to the computed center of the Maximum Inscribed Circle.
136	* The line segment from the center to this point
137	* is a radius of the constructed circle, and this point
138	* lies on the boundary of the circle.
139	*
140	* @return a point defining the radius of the Maximum Inscribed Circle
141	*/
142	public Point getRadiusPoint() {
143	compute();
144	return radiusPoint;
145	}
146
147	/**
148	* Gets a line representing a radius of the Largest Empty Circle.
149	*
150	* @return a line from the center of the circle to a point on the edge
151	*/
152	public LineString getRadiusLine() {
153	compute();
154	LineString radiusLine = factory.createLineString(
155	new Coordinate[] { centerPt.copy(), radiusPt.copy() });
156	return radiusLine;
157	}
158
159	/**
160	* Computes the signed distance from a point to the area boundary.
161	* Points outside the polygon are assigned a negative distance.
162	* Their containing cells will be last in the priority queue
163	* (but may still end up being tested since they may need to be refined).
164	*
165	* @param p the point to compute the distance for
166	* @return the signed distance to the area boundary (negative indicates outside the area)
167	*/
168	private double distanceToBoundary(Point p) {
169	double dist = indexedDistance.distance(p);
170	boolean isOutide = Location.EXTERIOR == ptLocater.locate(p.getCoordinate());
171	if (isOutide) return -dist;
172	return dist;
173	}
174
175	private double distanceToBoundary(double x, double y) {
176	Coordinate coord = new Coordinate(x, y);
177	Point pt = factory.createPoint(coord);
178	return distanceToBoundary(pt);
179	}
180
181	private void compute() {
182	// check if already computed
183	if (centerCell != null) return;
184
185	// Priority queue of cells, ordered by maximum distance from boundary
186	PriorityQueue<Cell> cellQueue = new PriorityQueue<>();
187
188	createInitialGrid(inputGeom.getEnvelopeInternal(), cellQueue);
189
190	// use the area centroid as the initial candidate center point
191	Cell farthestCell = createCentroidCell(inputGeom);
192	//int totalCells = cellQueue.size();
193
194	/**
195	* Carry out the branch-and-bound search
196	* of the cell space
197	*/
198	while (! cellQueue.isEmpty()) {
199	// pick the most promising cell from the queue
200	Cell cell = cellQueue.remove();
201	//System.out.println(factory.toGeometry(cell.getEnvelope()));
202
203	// update the center cell if the candidate is further from the boundary
204	if (cell.getDistance() > farthestCell.getDistance()) {
205	farthestCell = cell;
206	}
207	/**
208	* Refine this cell if the potential distance improvement
209	* is greater than the required tolerance.
210	* Otherwise the cell is pruned (not investigated further),
211	* since no point in it is further than
212	* the current farthest distance.
213	*/
214	double potentialIncrease = cell.getMaxDistance() - farthestCell.getDistance();
215	if (potentialIncrease > tolerance) {
216	// split the cell into four sub-cells
217	double h2 = cell.getHSide() / 2;
218	cellQueue.add( createCell( cell.getX() - h2, cell.getY() - h2, h2));
219	cellQueue.add( createCell( cell.getX() + h2, cell.getY() - h2, h2));
220	cellQueue.add( createCell( cell.getX() - h2, cell.getY() + h2, h2));
221	cellQueue.add( createCell( cell.getX() + h2, cell.getY() + h2, h2));
222	//totalCells += 4;
223	}
224	}
225	// the farthest cell is the best approximation to the MIC center
226	centerCell = farthestCell;
227	centerPt = new Coordinate(centerCell.getX(), centerCell.getY());
228	centerPoint = factory.createPoint(centerPt);
229	// compute radius point
230	Coordinate[] nearestPts = indexedDistance.nearestPoints(centerPoint);
231	radiusPt = nearestPts[0].copy();
232	radiusPoint = factory.createPoint(radiusPt);
233	}
234
235	/**
236	* Initializes the queue with a grid of cells covering
237	* the extent of the area.
238	*
239	* @param env the area extent to cover
240	* @param cellQueue the queue to initialize
241	*/
242	private void createInitialGrid(Envelope env, PriorityQueue<Cell> cellQueue) {
243	double minX = env.getMinX();
244	double maxX = env.getMaxX();
245	double minY = env.getMinY();
246	double maxY = env.getMaxY();
247	double width = env.getWidth();
248	double height = env.getHeight();
249	double cellSize = Math.min(width, height);
250	double hSide = cellSize / 2.0;
251
252	// compute initial grid of cells to cover area
253	for (double x = minX; x < maxX; x += cellSize) {
254	for (double y = minY; y < maxY; y += cellSize) {
255	cellQueue.add(createCell(x + hSide, y + hSide, hSide));
256	}
257	}
258	}
259
260	private Cell createCell(double x, double y, double hSide) {
261	return new Cell(x, y, hSide, distanceToBoundary(x, y));
262	}
263
264	// create a cell centered on area centroid
265	private Cell createCentroidCell(Geometry geom) {
266	Point p = geom.getCentroid();
267	return new Cell(p.getX(), p.getY(), 0, distanceToBoundary(p));
268	}
269
270	/**
271	* A square grid cell centered on a given point,
272	* with a given half-side size, and having a given distance
273	* to the area boundary.
274	* The maximum possible distance from any point in the cell to the
275	* boundary can be computed, and is used
276	* as the ordering and upper-bound function in
277	* the branch-and-bound algorithm.
278	*
279	*/
280	private static class Cell implements Comparable<Cell> {
281
282	private static final double SQRT2 = 1.4142135623730951;
283
284	private double x;
285	private double y;
286	private double hSide;
287	private double distance;
288	private double maxDist;
289
290	Cell(double x, double y, double hSide, double distanceToBoundary) {
291	this.x = x; // cell center x
292	this.y = y; // cell center y
293	this.hSide = hSide; // half the cell size
294
295	// the distance from cell center to area boundary
296	distance = distanceToBoundary;
297
298	// the maximum possible distance to area boundary for points in this cell
299	this.maxDist = distance + hSide * SQRT2;
300	}
301
302	public Envelope getEnvelope() {
303	return new Envelope(x - hSide, x + hSide, y - hSide, y + hSide);
304	}
305
306	public double getMaxDistance() {
307	return maxDist;
308	}
309
310	public double getDistance() {
311	return distance;
312	}
313
314	public double getHSide() {
315	return hSide;
316	}
317
318	public double getX() {
319	return x;
320	}
321
322	public double getY() {
323	return y;
324	}
325
326	/**
327	* A cell is greater iff its maximum possible distance is larger.
328	*/
329	public int compareTo(Cell o) {
330	return (int) (o.maxDist - this.maxDist);
331	}
332
333	}
334
335	}
336