ConvexHull.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 java.util.ArrayList;
import java.util.Arrays;
import java.util.Comparator;
import java.util.HashSet;
import java.util.List;
import java.util.Set;
import java.util.Stack;
import org.locationtech.jts.geom.Coordinate;
import org.locationtech.jts.geom.CoordinateArrays;
import org.locationtech.jts.geom.CoordinateList;
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.GeometryCollection;
import org.locationtech.jts.geom.GeometryFactory;
import org.locationtech.jts.geom.LineString;
import org.locationtech.jts.geom.LinearRing;
import org.locationtech.jts.geom.Point;
import org.locationtech.jts.geom.Polygon;
import org.locationtech.jts.util.Assert;
/**
* Computes the convex hull of a {@link Geometry}.
* The convex hull is the smallest convex Geometry that contains all the
* points in the input Geometry.
* <p>
* Uses the Graham Scan algorithm.
* <p>
* Incorporates heuristics to optimize checking for degenerate results,
* and to reduce the number of points processed for large inputs.
*
*@version 1.7
*/
public class ConvexHull
{
private static final int TUNING_REDUCE_SIZE = 50;
private GeometryFactory geomFactory;
private Coordinate[] inputPts;
/**
* Create a new convex hull construction for the input {@link Geometry}.
*/
public ConvexHull(Geometry geometry)
{
this(geometry.getCoordinates(), geometry.getFactory());
}
/**
* Create a new convex hull construction for the input {@link Coordinate} array.
*/
public ConvexHull(Coordinate[] pts, GeometryFactory geomFactory)
{
//-- suboptimal early uniquing - for performance testing only
//inputPts = UniqueCoordinateArrayFilter.filterCoordinates(pts);
inputPts = pts;
this.geomFactory = geomFactory;
}
/**
* Returns a {@link Geometry} that represents the convex hull of the input
* geometry.
* The returned geometry contains the minimal number of points needed to
* represent the convex hull. In particular, no more than two consecutive
* points will be collinear.
*
* @return if the convex hull contains 3 or more points, a {@link Polygon};
* 2 points, a {@link LineString};
* 1 point, a {@link Point};
* 0 points, an empty {@link GeometryCollection}.
*/
public Geometry getConvexHull() {
Geometry fewPointsGeom = createFewPointsResult();
if (fewPointsGeom != null)
return fewPointsGeom;
Coordinate[] reducedPts = inputPts;
//-- use heuristic to reduce points, if large
if (inputPts.length > TUNING_REDUCE_SIZE) {
reducedPts = reduce(inputPts);
}
else {
//-- the points must be made unique
reducedPts = extractUnique(inputPts);
}
// sort points for Graham scan.
Coordinate[] sortedPts = preSort(reducedPts);
// Use Graham scan to find convex hull.
Stack<Coordinate> cHS = grahamScan(sortedPts);
// Convert stack to an array.
Coordinate[] cH = toCoordinateArray(cHS);
// Convert array to appropriate output geometry.
// (an empty or point result will be detected earlier)
return lineOrPolygon(cH);
}
/**
* Checks if there are <= 2 unique points,
* which produce an obviously degenerate result.
* If there are more points, returns null to indicate this.
*
* This is a fast check for an obviously degenerate result.
* If the result is not obviously degenerate (at least 3 unique points found)
* the full uniquing of the entire point set is
* done only once during the reduce phase.
*
* @return a degenerate hull geometry, or null if the number of input points is large
*/
private Geometry createFewPointsResult() {
Coordinate[] uniquePts = extractUnique(inputPts, 2);
if (uniquePts == null) {
return null;
}
else if (uniquePts.length == 0) {
return geomFactory.createGeometryCollection();
}
else if (uniquePts.length == 1) {
return geomFactory.createPoint(uniquePts[0]);
}
else {
return geomFactory.createLineString(uniquePts);
}
}
private static Coordinate[] extractUnique(Coordinate[] pts) {
return extractUnique(pts, -1);
}
/**
* Extracts unique coordinates from an array of coordinates,
* up to an (optional) maximum count of values.
* If more than the given maximum of unique values are found,
* this is reported by returning <code>null</code>.
* This avoids scanning all input points if not needed.
* If the maximum points is not specified, all unique points are extracted.
*
* @param pts an array of Coordinates
* @param maxPts the maximum number of unique points to scan, or -1
* @return an array of unique values, or null
*/
private static Coordinate[] extractUnique(Coordinate[] pts, int maxPts) {
Set<Coordinate> uniquePts = new HashSet<Coordinate>();
for (Coordinate pt : pts) {
uniquePts.add(pt);
//-- if maxPts is provided, exit if more unique pts found
if (maxPts >= 0 && uniquePts.size() > maxPts) return null;
}
return CoordinateArrays.toCoordinateArray(uniquePts);
}
/**
* An alternative to Stack.toArray, which is not present in earlier versions
* of Java.
*/
protected Coordinate[] toCoordinateArray(Stack<Coordinate> stack) {
Coordinate[] coordinates = new Coordinate[stack.size()];
for (int i = 0; i < stack.size(); i++) {
Coordinate coordinate = (Coordinate) stack.get(i);
coordinates[i] = coordinate;
}
return coordinates;
}
/**
* Uses a heuristic to reduce the number of points scanned
* to compute the hull.
* The heuristic is to find a polygon guaranteed to
* be in (or on) the hull, and eliminate all points inside it.
* A quadrilateral defined by the extremal points
* in the four orthogonal directions
* can be used, but even more inclusive is
* to use an octilateral defined by the points in the 8 cardinal directions.
* <p>
* Note that even if the method used to determine the polygon vertices
* is not 100% robust, this does not affect the robustness of the convex hull.
* <p>
* To satisfy the requirements of the Graham Scan algorithm,
* the returned array has at least 3 entries.
* <p>
* This has the side effect of making the reduced points unique,
* as required by the convex hull algorithm used.
*
* @param pts the points to reduce
* @return the reduced list of points (at least 3)
*/
private Coordinate[] reduce(Coordinate[] inputPts)
{
//Coordinate[] polyPts = computeQuad(inputPts);
Coordinate[] innerPolyPts = computeInnerOctolateralRing(inputPts);
// unable to compute interior polygon for some reason
if (innerPolyPts == null)
return inputPts;
// LinearRing ring = geomFactory.createLinearRing(polyPts);
// System.out.println(ring);
// add points defining polygon
Set<Coordinate> reducedSet = new HashSet();
for (int i = 0; i < innerPolyPts.length; i++) {
reducedSet.add(innerPolyPts[i]);
}
/**
* Add all unique points not in the interior poly.
* CGAlgorithms.isPointInRing is not defined for points exactly on the ring,
* but this doesn't matter since the points of the interior polygon
* are forced to be in the reduced set.
*/
for (int i = 0; i < inputPts.length; i++) {
if (! PointLocation.isInRing(inputPts[i], innerPolyPts)) {
reducedSet.add(inputPts[i]);
}
}
Coordinate[] reducedPts = CoordinateArrays.toCoordinateArray(reducedSet);
// ensure that computed array has at least 3 points (not necessarily unique)
if (reducedPts.length < 3)
return padArray3(reducedPts);
return reducedPts;
}
private Coordinate[] padArray3(Coordinate[] pts)
{
Coordinate[] pad = new Coordinate[3];
for (int i = 0; i < pad.length; i++) {
if (i < pts.length) {
pad[i] = pts[i];
}
else
pad[i] = pts[0];
}
return pad;
}
/**
* Sorts the points radially CW around the point with minimum Y and then X.
*
* @param pts the points to sort
* @return the sorted points
*/
private Coordinate[] preSort(Coordinate[] pts) {
Coordinate t;
/**
* find the lowest point in the set. If two or more points have
* the same minimum Y coordinate choose the one with the minimum X.
* This focal point is put in array location pts[0].
*/
for (int i = 1; i < pts.length; i++) {
if ((pts[i].y < pts[0].y) || ((pts[i].y == pts[0].y) && (pts[i].x < pts[0].x))) {
t = pts[0];
pts[0] = pts[i];
pts[i] = t;
}
}
// sort the points radially around the focal point.
Arrays.sort(pts, 1, pts.length, new RadialComparator(pts[0]));
return pts;
}
/**
* Uses the Graham Scan algorithm to compute the convex hull vertices.
*
* @param c a list of points, with at least 3 entries
* @return a Stack containing the ordered points of the convex hull ring
*/
private Stack<Coordinate> grahamScan(Coordinate[] c) {
Coordinate p;
Stack<Coordinate> ps = new Stack<Coordinate>();
ps.push(c[0]);
ps.push(c[1]);
ps.push(c[2]);
for (int i = 3; i < c.length; i++) {
Coordinate cp = c[i];
p = (Coordinate) ps.pop();
// check for empty stack to guard against robustness problems
while (
! ps.empty() &&
Orientation.index((Coordinate) ps.peek(), p, cp) > 0) {
p = (Coordinate) ps.pop();
}
ps.push(p);
ps.push(cp);
}
ps.push(c[0]);
return ps;
}
/**
*@return whether the three coordinates are collinear and c2 lies between
* c1 and c3 inclusive
*/
private boolean isBetween(Coordinate c1, Coordinate c2, Coordinate c3) {
if (Orientation.index(c1, c2, c3) != 0) {
return false;
}
if (c1.x != c3.x) {
if (c1.x <= c2.x && c2.x <= c3.x) {
return true;
}
if (c3.x <= c2.x && c2.x <= c1.x) {
return true;
}
}
if (c1.y != c3.y) {
if (c1.y <= c2.y && c2.y <= c3.y) {
return true;
}
if (c3.y <= c2.y && c2.y <= c1.y) {
return true;
}
}
return false;
}
private Coordinate[] computeInnerOctolateralRing(Coordinate[] inputPts) {
Coordinate[] octPts = computeInnerOctolateralPts(inputPts);
CoordinateList coordList = new CoordinateList();
coordList.add(octPts, false);
// points must all lie in a line
if (coordList.size() < 3) {
return null;
}
coordList.closeRing();
return coordList.toCoordinateArray();
}
/**
* Computes the extremal points of an inner octolateral.
* Some points may be duplicates - these are collapsed later.
*
* @param inputPts the points to compute the octolateral for
* @return the extremal points of the octolateral
*/
private Coordinate[] computeInnerOctolateralPts(Coordinate[] inputPts)
{
Coordinate[] pts = new Coordinate[8];
for (int j = 0; j < pts.length; j++) {
pts[j] = inputPts[0];
}
for (int i = 1; i < inputPts.length; i++) {
if (inputPts[i].x < pts[0].x) {
pts[0] = inputPts[i];
}
if (inputPts[i].x - inputPts[i].y < pts[1].x - pts[1].y) {
pts[1] = inputPts[i];
}
if (inputPts[i].y > pts[2].y) {
pts[2] = inputPts[i];
}
if (inputPts[i].x + inputPts[i].y > pts[3].x + pts[3].y) {
pts[3] = inputPts[i];
}
if (inputPts[i].x > pts[4].x) {
pts[4] = inputPts[i];
}
if (inputPts[i].x - inputPts[i].y > pts[5].x - pts[5].y) {
pts[5] = inputPts[i];
}
if (inputPts[i].y < pts[6].y) {
pts[6] = inputPts[i];
}
if (inputPts[i].x + inputPts[i].y < pts[7].x + pts[7].y) {
pts[7] = inputPts[i];
}
}
return pts;
}
/**
*@param vertices the vertices of a linear ring, which may or may not be
* flattened (i.e. vertices collinear)
*@return a 2-vertex <code>LineString</code> if the vertices are
* collinear; otherwise, a <code>Polygon</code> with unnecessary
* (collinear) vertices removed
*/
private Geometry lineOrPolygon(Coordinate[] coordinates) {
coordinates = cleanRing(coordinates);
if (coordinates.length == 3) {
return geomFactory.createLineString(new Coordinate[]{coordinates[0], coordinates[1]});
}
LinearRing linearRing = geomFactory.createLinearRing(coordinates);
return geomFactory.createPolygon(linearRing);
}
/**
* Cleans a list of points by removing interior collinear vertices.
*
* @param vertices the vertices of a linear ring, which may or may not be
* flattened (i.e. vertices collinear)
* @return the coordinates with unnecessary (collinear) vertices removed
*/
private Coordinate[] cleanRing(Coordinate[] original) {
Assert.equals(original[0], original[original.length - 1]);
List<Coordinate> cleanedRing = new ArrayList<Coordinate>();
Coordinate previousDistinctCoordinate = null;
for (int i = 0; i <= original.length - 2; i++) {
Coordinate currentCoordinate = original[i];
Coordinate nextCoordinate = original[i+1];
if (currentCoordinate.equals(nextCoordinate)) {
continue;
}
if (previousDistinctCoordinate != null
&& isBetween(previousDistinctCoordinate, currentCoordinate, nextCoordinate)) {
continue;
}
cleanedRing.add(currentCoordinate);
previousDistinctCoordinate = currentCoordinate;
}
cleanedRing.add(original[original.length - 1]);
Coordinate[] cleanedRingCoordinates = new Coordinate[cleanedRing.size()];
return (Coordinate[]) cleanedRing.toArray(cleanedRingCoordinates);
}
/**
* Compares {@link Coordinate}s for their angle and distance
* relative to an origin.
* The origin is assumed to be lower in Y and then X than
* all other point inputs.
* The points are ordered CCW around the origin.
*
* @author Martin Davis
* @version 1.7
*/
private static class RadialComparator
implements Comparator<Coordinate>
{
private Coordinate origin;
/**
* Creates a new comparator using a given origin.
* The origin must be lower in Y and then X to all
* compared points.
*
* @param origin the origin of the radial comparison
*/
public RadialComparator(Coordinate origin)
{
this.origin = origin;
}
@Override
public int compare(Coordinate p1, Coordinate p2)
{
int comp = polarCompare(origin, p1, p2);
return comp;
}
/**
* Given two points p and q compare them with respect to their radial
* ordering about point o. First checks radial ordering.
* using a CCW orientation.
* If the points are collinear, the comparison is based
* on their distance to the origin.
* <p>
* p < q iff
* <ul>
* <li>ang(o-p) < ang(o-q) (e.g. o-p-q is CCW)
* <li>or ang(o-p) == ang(o-q) && dist(o,p) < dist(o,q)
* </ul>
*
* @param o the origin
* @param p a point
* @param q another point
* @return -1, 0 or 1 depending on whether p is less than,
* equal to or greater than q
*/
private static int polarCompare(Coordinate o, Coordinate p, Coordinate q)
{
int orient = Orientation.index(o, p, q);
if (orient == Orientation.COUNTERCLOCKWISE) return 1;
if (orient == Orientation.CLOCKWISE) return -1;
/**
* The points are collinear,
* so compare based on distance from the origin.
* The points p and q are >= to the origin,
* so they lie in the closed half-plane above the origin.
* If they are not in a horizontal line,
* the Y ordinate can be tested to determine distance.
* This is more robust than computing the distance explicitly.
*/
if (p.y > q.y) return 1;
if (p.y < q.y) return -1;
/**
* The points lie in a horizontal line, which should also contain the origin
* (since they are collinear).
* Also, they must be above the origin.
* Use the X ordinate to determine distance.
*/
if (p.x > q.x) return 1;
if (p.x < q.x) return -1;
// Assert: p = q
return 0;
}
}
}