MinimumBoundingCircle.java
/*
* Copyright (c) 2016 Vivid Solutions.
*
* All rights reserved. This program and the accompanying materials
* are made available under the terms of the Eclipse Public License 2.0
* and Eclipse Distribution License v. 1.0 which accompanies this distribution.
* The Eclipse Public License is available at http://www.eclipse.org/legal/epl-v20.html
* and the Eclipse Distribution License is available at
*
* http://www.eclipse.org/org/documents/edl-v10.php.
*/
package org.locationtech.jts.algorithm;
import org.locationtech.jts.geom.Coordinate;
import org.locationtech.jts.geom.CoordinateArrays;
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.Point;
import org.locationtech.jts.geom.Triangle;
import org.locationtech.jts.util.Assert;
/**
* Computes the <b>Minimum Bounding Circle</b> (MBC)
* for the points in a {@link Geometry}.
* The MBC is the smallest circle which <tt>cover</tt>s
* all the input points
* (this is also known as the <b>Smallest Enclosing Circle</b>).
* This is equivalent to computing the Maximum Diameter
* of the input point set.
* <p>
* The computed circle can be specified in two equivalent ways,
* both of which are provide as output by this class:
* <ul>
* <li>As a centre point and a radius
* <li>By the set of points defining the circle.
* Depending on the number of points in the input
* and their relative positions, this set
* contains from 0 to 3 points.
* <ul>
* <li>0 or 1 points indicate an empty or trivial input point arrangement.
* <li>2 points define the diameter of the minimum bounding circle.
* <li>3 points define an inscribed triangle of the minimum bounding circle.
* </ul>
* </ul>
* The class can also output a {@link Geometry} which approximates the
* shape of the Minimum Bounding Circle (although as an approximation
* it is <b>not</b> guaranteed to <tt>cover</tt> all the input points.)
* <p>
* The Maximum Diameter of the input point set can
* be computed as well. The Maximum Diameter is
* defined by the pair of input points with maximum distance between them.
* The points of the maximum diameter are two of the extremal points of the Minimum Bounding Circle.
* They lie on the convex hull of the input.
* However, that the maximum diameter is not a diameter
* of the Minimum Bounding Circle in the case where the MBC is
* defined by an inscribed triangle.
*
* @author Martin Davis
*
* @see MinimumDiameter
*
*/
public class MinimumBoundingCircle
{
/*
* The algorithm used is based on the one by Jon Rokne in
* the article "An Easy Bounding Circle" in <i>Graphic Gems II</i>.
*/
private Geometry input;
private Coordinate[] extremalPts = null;
private Coordinate centre = null;
private double radius = 0.0;
/**
* Creates a new object for computing the minimum bounding circle for the
* point set defined by the vertices of the given geometry.
*
* @param geom the geometry to use to obtain the point set
*/
public MinimumBoundingCircle(Geometry geom)
{
this.input = geom;
}
/**
* Gets a geometry which represents the Minimum Bounding Circle.
* If the input is degenerate (empty or a single unique point),
* this method will return an empty geometry or a single Point geometry.
* Otherwise, a Polygon will be returned which approximates the
* Minimum Bounding Circle.
* (Note that because the computed polygon is only an approximation,
* it may not precisely contain all the input points.)
*
* @return a Geometry representing the Minimum Bounding Circle.
*/
public Geometry getCircle()
{
//TODO: ensure the output circle contains the extermal points.
//TODO: or maybe even ensure that the returned geometry contains ALL the input points?
compute();
if (centre == null)
return input.getFactory().createPolygon();
Point centrePoint = input.getFactory().createPoint(centre);
if (radius == 0.0)
return centrePoint;
return centrePoint.buffer(radius);
}
/**
* Gets a geometry representing the maximum diameter of the
* input. The maximum diameter is the longest line segment
* between any two points of the input.
* <p>
* The points are two of the extremal points of the Minimum Bounding Circle.
* They lie on the convex hull of the input.
*
* @return a LineString between the two farthest points of the input
* @return a empty LineString if the input is empty
* @return a Point if the input is a point
*/
public Geometry getMaximumDiameter() {
compute();
switch (extremalPts.length) {
case 0:
return input.getFactory().createLineString();
case 1:
return input.getFactory().createPoint(centre);
case 2:
return input.getFactory().createLineString(
new Coordinate[] { extremalPts[0], extremalPts[1] });
default: // case 3
Coordinate[] maxDiameter = farthestPoints(extremalPts);
return input.getFactory().createLineString(maxDiameter);
}
}
/**
* Gets a geometry representing a line between the two farthest points
* in the input.
* <p>
* The points are two of the extremal points of the Minimum Bounding Circle.
* They lie on the convex hull of the input.
*
* @return a LineString between the two farthest points of the input
* @return a empty LineString if the input is empty
* @return a Point if the input is a point
*
* @deprecated use #getMaximumDiameter()
*/
public Geometry getFarthestPoints() {
return getMaximumDiameter();
}
/**
* Finds the farthest pair out of 3 extremal points
* @param pts the array of extremal points
* @return the pair of farthest points
*/
private static Coordinate[] farthestPoints(Coordinate[] pts) {
double dist01 = pts[0].distance(pts[1]);
double dist12 = pts[1].distance(pts[2]);
double dist20 = pts[2].distance(pts[0]);
if (dist01 >= dist12 && dist01 >= dist20) {
return new Coordinate[] { pts[0], pts[1] };
}
if (dist12 >= dist01 && dist12 >= dist20) {
return new Coordinate[] { pts[1], pts[2] };
}
return new Coordinate[] { pts[2], pts[0] };
}
/**
* Gets a geometry representing the diameter of the computed Minimum Bounding
* Circle.
*
* @return the diameter LineString of the Minimum Bounding Circle
* @return a empty LineString if the input is empty
* @return a Point if the input is a point
*/
public Geometry getDiameter() {
compute();
switch (extremalPts.length) {
case 0:
return input.getFactory().createLineString();
case 1:
return input.getFactory().createPoint(centre);
}
// TODO: handle case of 3 extremal points, by computing a line from one of
// them through the centre point with len = 2*radius
Coordinate p0 = extremalPts[0];
Coordinate p1 = extremalPts[1];
return input.getFactory().createLineString(new Coordinate[] { p0, p1 });
}
/**
* Gets the extremal points which define the computed Minimum Bounding Circle.
* There may be zero, one, two or three of these points, depending on the number
* of points in the input and the geometry of those points.
* <ul>
* <li>0 or 1 points indicate an empty or trivial input point arrangement.
* <li>2 points define the diameter of the Minimum Bounding Circle.
* <li>3 points define an inscribed triangle of which the Minimum Bounding Circle is the circumcircle.
* The longest chords of the circle are the line segments [0-1] and [1-2]
* </ul>
*
* @return the points defining the Minimum Bounding Circle
*/
public Coordinate[] getExtremalPoints()
{
compute();
return extremalPts;
}
/**
* Gets the centre point of the computed Minimum Bounding Circle.
*
* @return the centre point of the Minimum Bounding Circle
* @return null if the input is empty
*/
public Coordinate getCentre() {
compute();
return centre;
}
/**
* Gets the radius of the computed Minimum Bounding Circle.
*
* @return the radius of the Minimum Bounding Circle
*/
public double getRadius()
{
compute();
return radius;
}
private void computeCentre()
{
switch (extremalPts.length) {
case 0:
centre = null;
break;
case 1:
centre = extremalPts[0];
break;
case 2:
centre = new Coordinate(
(extremalPts[0].x + extremalPts[1].x) / 2.0,
(extremalPts[0].y + extremalPts[1].y) / 2.0
);
break;
case 3:
centre = Triangle.circumcentre(extremalPts[0], extremalPts[1], extremalPts[2]);
break;
}
}
private void compute()
{
if (extremalPts != null) return;
computeCirclePoints();
computeCentre();
if (centre != null)
radius = centre.distance(extremalPts[0]);
}
private void computeCirclePoints()
{
// handle degenerate or trivial cases
if (input.isEmpty()) {
extremalPts = new Coordinate[0];
return;
}
if (input.getNumPoints() == 1) {
Coordinate[] pts = input.getCoordinates();
extremalPts = new Coordinate[] { new Coordinate(pts[0]) };
return;
}
/**
* The problem is simplified by reducing to the convex hull.
* Computing the convex hull also has the useful effect of eliminating duplicate points
*/
Geometry convexHull = input.convexHull();
Coordinate[] hullPts = convexHull.getCoordinates();
// strip duplicate final point, if any
Coordinate[] pts = hullPts;
if (hullPts[0].equals2D(hullPts[hullPts.length - 1])) {
pts = new Coordinate[hullPts.length - 1];
CoordinateArrays.copyDeep(hullPts, 0, pts, 0, hullPts.length - 1);
}
/**
* Optimization for the trivial case where the CH has fewer than 3 points
*/
if (pts.length <= 2) {
extremalPts = CoordinateArrays.copyDeep(pts);
return;
}
// find a point P with minimum Y ordinate
Coordinate P = lowestPoint(pts);
// find a point Q such that the angle that PQ makes with the x-axis is minimal
Coordinate Q = pointWitMinAngleWithX(pts, P);
/**
* Iterate over the remaining points to find
* a pair or triplet of points which determine the minimal circle.
* By the design of the algorithm,
* at most <tt>pts.length</tt> iterations are required to terminate
* with a correct result.
*/
for (int i = 0; i < pts.length; i++) {
Coordinate R = pointWithMinAngleWithSegment(pts, P, Q);
if (Angle.isObtuse(P, R, Q)) {
// if PRQ is obtuse, then MBC is determined by P and Q
extremalPts = new Coordinate[] { new Coordinate(P), new Coordinate(Q) };
return;
}
else if (Angle.isObtuse(R, P, Q)) {
// if RPQ is obtuse, update baseline and iterate
P = R;
continue;
}
else if (Angle.isObtuse(R, Q, P)) {
// if RQP is obtuse, update baseline and iterate
Q = R;
continue;
}
else {
// otherwise all angles are acute, and the MBC is determined by the triangle PQR
extremalPts = new Coordinate[] { new Coordinate(P), new Coordinate(Q), new Coordinate(R) };
return;
}
}
Assert.shouldNeverReachHere("Logic failure in Minimum Bounding Circle algorithm!");
}
private static Coordinate lowestPoint(Coordinate[] pts)
{
Coordinate min = pts[0];
for (int i = 1; i < pts.length; i++) {
if (pts[i].y < min.y)
min = pts[i];
}
return min;
}
private static Coordinate pointWitMinAngleWithX(Coordinate[] pts, Coordinate P)
{
double minSin = Double.MAX_VALUE;
Coordinate minAngPt = null;
for (int i = 0; i < pts.length; i++) {
Coordinate p = pts[i];
if (p == P) continue;
/**
* The sin of the angle is a simpler proxy for the angle itself
*/
double dx = p.x - P.x;
double dy = p.y - P.y;
if (dy < 0) dy = -dy;
double len = Math.hypot(dx, dy);
double sin = dy / len;
if (sin < minSin) {
minSin = sin;
minAngPt = p;
}
}
return minAngPt;
}
private static Coordinate pointWithMinAngleWithSegment(Coordinate[] pts, Coordinate P, Coordinate Q)
{
double minAng = Double.MAX_VALUE;
Coordinate minAngPt = null;
for (int i = 0; i < pts.length; i++) {
Coordinate p = pts[i];
if (p == P) continue;
if (p == Q) continue;
double ang = Angle.angleBetween(P, p, Q);
if (ang < minAng) {
minAng = ang;
minAngPt = p;
}
}
return minAngPt;
}
}