Class MinimumBoundingCircle

1	/*
2	* Copyright (c) 2016 Vivid Solutions.
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
13	package org.locationtech.jts.algorithm;
14
15	import org.locationtech.jts.geom.Coordinate;
16	import org.locationtech.jts.geom.CoordinateArrays;
17	import org.locationtech.jts.geom.CoordinateSequence;
18	import org.locationtech.jts.geom.Geometry;
19	import org.locationtech.jts.geom.Point;
20	import org.locationtech.jts.geom.Triangle;
21	import org.locationtech.jts.util.Assert;
22
23	/**
24	* Computes the <b>Minimum Bounding Circle</b> (MBC)
25	* for the points in a {@link Geometry}.
26	* The MBC is the smallest circle which <tt>cover</tt>s
27	* all the input points
28	* (this is also known as the <b>Smallest Enclosing Circle</b>).
29	* This is equivalent to computing the Maximum Diameter
30	* of the input point set.
31	* <p>
32	* The computed circle can be specified in two equivalent ways,
33	* both of which are provide as output by this class:
34	* <ul>
35	* <li>As a centre point and a radius
36	* <li>By the set of points defining the circle.
37	* Depending on the number of points in the input
38	* and their relative positions, this set
39	* contains from 0 to 3 points.
40	* <ul>
41	* <li>0 or 1 points indicate an empty or trivial input point arrangement.
42	* <li>2 points define the diameter of the minimum bounding circle.
43	* <li>3 points define an inscribed triangle of the minimum bounding circle.
44	* </ul>
45	* </ul>
46	* The class can also output a {@link Geometry} which approximates the
47	* shape of the Minimum Bounding Circle (although as an approximation
48	* it is <b>not</b> guaranteed to <tt>cover</tt> all the input points.)
49	* <p>
50	* The Maximum Diameter of the input point set can
51	* be computed as well. The Maximum Diameter is
52	* defined by the pair of input points with maximum distance between them.
53	* The points of the maximum diameter are two of the extremal points of the Minimum Bounding Circle.
54	* They lie on the convex hull of the input.
55	* However, that the maximum diameter is not a diameter
56	* of the Minimum Bounding Circle in the case where the MBC is
57	* defined by an inscribed triangle.
58	*
59	* @author Martin Davis
60	*
61	* @see MinimumDiameter
62	*
63	*/
64	public class MinimumBoundingCircle
65	{
66	/*
67	* The algorithm used is based on the one by Jon Rokne in
68	* the article "An Easy Bounding Circle" in <i>Graphic Gems II</i>.
69	*/
70
71	private Geometry input;
72	private Coordinate[] extremalPts = null;
73	private Coordinate centre = null;
74	private double radius = 0.0;
75
76	/**
77	* Creates a new object for computing the minimum bounding circle for the
78	* point set defined by the vertices of the given geometry.
79	*
80	* @param geom the geometry to use to obtain the point set
81	*/
82	public MinimumBoundingCircle(Geometry geom)
83	{
84	this.input = geom;
85	}
86
87	/**
88	* Gets a geometry which represents the Minimum Bounding Circle.
89	* If the input is degenerate (empty or a single unique point),
90	* this method will return an empty geometry or a single Point geometry.
91	* Otherwise, a Polygon will be returned which approximates the
92	* Minimum Bounding Circle.
93	* (Note that because the computed polygon is only an approximation,
94	* it may not precisely contain all the input points.)
95	*
96	* @return a Geometry representing the Minimum Bounding Circle.
97	*/
98	public Geometry getCircle()
99	{
100	//TODO: ensure the output circle contains the extermal points.
101	//TODO: or maybe even ensure that the returned geometry contains ALL the input points?
102
103	compute();
104	if (centre == null)
105	return input.getFactory().createPolygon();
106	Point centrePoint = input.getFactory().createPoint(centre);
107	if (radius == 0.0)
108	return centrePoint;
109	return centrePoint.buffer(radius);
110	}
111
112	/**
113	* Gets a geometry representing the maximum diameter of the
114	* input. The maximum diameter is the longest line segment
115	* between any two points of the input.
116	* <p>
117	* The points are two of the extremal points of the Minimum Bounding Circle.
118	* They lie on the convex hull of the input.
119	*
120	* @return a LineString between the two farthest points of the input
121	* @return a empty LineString if the input is empty
122	* @return a Point if the input is a point
123	*/
124	public Geometry getMaximumDiameter() {
125	compute();
126	switch (extremalPts.length) {
127	case 0:
128	return input.getFactory().createLineString();
129	case 1:
130	return input.getFactory().createPoint(centre);
131	case 2:
132	return input.getFactory().createLineString(
133	new Coordinate[] { extremalPts[0], extremalPts[1] });
134	default: // case 3
135	Coordinate[] maxDiameter = farthestPoints(extremalPts);
136	return input.getFactory().createLineString(maxDiameter);
137	}
138	}
139
140	/**
141	* Gets a geometry representing a line between the two farthest points
142	* in the input.
143	* <p>
144	* The points are two of the extremal points of the Minimum Bounding Circle.
145	* They lie on the convex hull of the input.
146	*
147	* @return a LineString between the two farthest points of the input
148	* @return a empty LineString if the input is empty
149	* @return a Point if the input is a point
150	*
151	* @deprecated use #getMaximumDiameter()
152	*/
153	public Geometry getFarthestPoints() {
154	return getMaximumDiameter();
155	}
156
157	/**
158	* Finds the farthest pair out of 3 extremal points
159	* @param pts the array of extremal points
160	* @return the pair of farthest points
161	*/
162	private static Coordinate[] farthestPoints(Coordinate[] pts) {
163	double dist01 = pts[0].distance(pts[1]);
164	double dist12 = pts[1].distance(pts[2]);
165	double dist20 = pts[2].distance(pts[0]);
166	if (dist01 >= dist12 && dist01 >= dist20) {
167	return new Coordinate[] { pts[0], pts[1] };
168	}
169	if (dist12 >= dist01 && dist12 >= dist20) {
170	return new Coordinate[] { pts[1], pts[2] };
171	}
172	return new Coordinate[] { pts[2], pts[0] };
173	}
174
175	/**
176	* Gets a geometry representing the diameter of the computed Minimum Bounding
177	* Circle.
178	*
179	* @return the diameter LineString of the Minimum Bounding Circle
180	* @return a empty LineString if the input is empty
181	* @return a Point if the input is a point
182	*/
183	public Geometry getDiameter() {
184	compute();
185	switch (extremalPts.length) {
186	case 0:
187	return input.getFactory().createLineString();
188	case 1:
189	return input.getFactory().createPoint(centre);
190	}
191	// TODO: handle case of 3 extremal points, by computing a line from one of
192	// them through the centre point with len = 2*radius
193	Coordinate p0 = extremalPts[0];
194	Coordinate p1 = extremalPts[1];
195	return input.getFactory().createLineString(new Coordinate[] { p0, p1 });
196	}
197
198	/**
199	* Gets the extremal points which define the computed Minimum Bounding Circle.
200	* There may be zero, one, two or three of these points, depending on the number
201	* of points in the input and the geometry of those points.
202	* <ul>
203	* <li>0 or 1 points indicate an empty or trivial input point arrangement.
204	* <li>2 points define the diameter of the Minimum Bounding Circle.
205	* <li>3 points define an inscribed triangle of which the Minimum Bounding Circle is the circumcircle.
206	* The longest chords of the circle are the line segments [0-1] and [1-2]
207	* </ul>
208	*
209	* @return the points defining the Minimum Bounding Circle
210	*/
211	public Coordinate[] getExtremalPoints()
212	{
213	compute();
214	return extremalPts;
215	}
216
217	/**
218	* Gets the centre point of the computed Minimum Bounding Circle.
219	*
220	* @return the centre point of the Minimum Bounding Circle
221	* @return null if the input is empty
222	*/
223	public Coordinate getCentre() {
224	compute();
225	return centre;
226	}
227
228	/**
229	* Gets the radius of the computed Minimum Bounding Circle.
230	*
231	* @return the radius of the Minimum Bounding Circle
232	*/
233	public double getRadius()
234	{
235	compute();
236	return radius;
237	}
238
239	private void computeCentre()
240	{
241	switch (extremalPts.length) {
242	case 0:
243	centre = null;
244	break;
245	case 1:
246	centre = extremalPts[0];
247	break;
248	case 2:
249	centre = new Coordinate(
250	(extremalPts[0].x + extremalPts[1].x) / 2.0,
251	(extremalPts[0].y + extremalPts[1].y) / 2.0
252	);
253	break;
254	case 3:
255	centre = Triangle.circumcentre(extremalPts[0], extremalPts[1], extremalPts[2]);
256	break;
257	}
258	}
259
260	private void compute()
261	{
262	if (extremalPts != null) return;
263
264	computeCirclePoints();
265	computeCentre();
266	if (centre != null)
267	radius = centre.distance(extremalPts[0]);
268	}
269
270	private void computeCirclePoints()
271	{
272	// handle degenerate or trivial cases
273	if (input.isEmpty()) {
274	extremalPts = new Coordinate[0];
275	return;
276	}
277	if (input.getNumPoints() == 1) {
278	Coordinate[] pts = input.getCoordinates();
279	extremalPts = new Coordinate[] { new Coordinate(pts[0]) };
280	return;
281	}
282
283	/**
284	* The problem is simplified by reducing to the convex hull.
285	* Computing the convex hull also has the useful effect of eliminating duplicate points
286	*/
287	Geometry convexHull = input.convexHull();
288
289	Coordinate[] hullPts = convexHull.getCoordinates();
290
291	// strip duplicate final point, if any
292	Coordinate[] pts = hullPts;
293	if (hullPts[0].equals2D(hullPts[hullPts.length - 1])) {
294	pts = new Coordinate[hullPts.length - 1];
295	CoordinateArrays.copyDeep(hullPts, 0, pts, 0, hullPts.length - 1);
296	}
297
298	/**
299	* Optimization for the trivial case where the CH has fewer than 3 points
300	*/
301	if (pts.length <= 2) {
302	extremalPts = CoordinateArrays.copyDeep(pts);
303	return;
304	}
305
306	// find a point P with minimum Y ordinate
307	Coordinate P = lowestPoint(pts);
308
309	// find a point Q such that the angle that PQ makes with the x-axis is minimal
310	Coordinate Q = pointWitMinAngleWithX(pts, P);
311
312	/**
313	* Iterate over the remaining points to find
314	* a pair or triplet of points which determine the minimal circle.
315	* By the design of the algorithm,
316	* at most <tt>pts.length</tt> iterations are required to terminate
317	* with a correct result.
318	*/
319	for (int i = 0; i < pts.length; i++) {
320	Coordinate R = pointWithMinAngleWithSegment(pts, P, Q);
321
322	// if PRQ is obtuse, then MBC is determined by P and Q
323	if (Angle.isObtuse(P, R, Q)) {
324	extremalPts = new Coordinate[] { new Coordinate(P), new Coordinate(Q) };
325	return;
326	}
327	// if RPQ is obtuse, update baseline and iterate
328	if (Angle.isObtuse(R, P, Q)) {
329	P = R;
330	continue;
331	}
332	// if RQP is obtuse, update baseline and iterate
333	if (Angle.isObtuse(R, Q, P)) {
334	Q = R;
335	continue;
336	}
337	// otherwise all angles are acute, and the MBC is determined by the triangle PQR
338	extremalPts = new Coordinate[] { new Coordinate(P), new Coordinate(Q), new Coordinate(R) };
339	return;
340	}
341	Assert.shouldNeverReachHere("Logic failure in Minimum Bounding Circle algorithm!");
342	}
343
344	private static Coordinate lowestPoint(Coordinate[] pts)
345	{
346	Coordinate min = pts[0];
347	for (int i = 1; i < pts.length; i++) {
348	if (pts[i].y < min.y)
349	min = pts[i];
350	}
351	return min;
352	}
353
354	private static Coordinate pointWitMinAngleWithX(Coordinate[] pts, Coordinate P)
355	{
356	double minSin = Double.MAX_VALUE;
357	Coordinate minAngPt = null;
358	for (int i = 0; i < pts.length; i++) {
359
360	Coordinate p = pts[i];
361	if (p == P) continue;
362
363	/**
364	* The sin of the angle is a simpler proxy for the angle itself
365	*/
366	double dx = p.x - P.x;
367	double dy = p.y - P.y;
368	if (dy < 0) dy = -dy;
369	double len = Math.sqrt(dx * dx + dy * dy);
370	double sin = dy / len;
371
372	if (sin < minSin) {
373	minSin = sin;
374	minAngPt = p;
375	}
376	}
377	return minAngPt;
378	}
379
380	private static Coordinate pointWithMinAngleWithSegment(Coordinate[] pts, Coordinate P, Coordinate Q)
381	{
382	double minAng = Double.MAX_VALUE;
383	Coordinate minAngPt = null;
384	for (int i = 0; i < pts.length; i++) {
385
386	Coordinate p = pts[i];
387	if (p == P) continue;
388	if (p == Q) continue;
389
390	double ang = Angle.angleBetween(P, p, Q);
391	if (ang < minAng) {
392	minAng = ang;
393	minAngPt = p;
394	}
395	}
396	return minAngPt;
397
398	}
399	}
400
401
402