Doc Classic

Class ConvexHull

  1  
  2  
  3 /*
  4  * Copyright (c) 2016 Vivid Solutions.
  5  *
  6  * All rights reserved. This program and the accompanying materials
  7  * are made available under the terms of the Eclipse Public License 2.0
  8  * and Eclipse Distribution License v. 1.0 which accompanies this distribution.
  9  * The Eclipse Public License is available at http://www.eclipse.org/legal/epl-v20.html
 10  * and the Eclipse Distribution License is available at
 11  *
 12  * http://www.eclipse.org/org/documents/edl-v10.php.
 13  */
 14 package org.locationtech.jts.algorithm;
 15  
 16 import java.util.ArrayList;
 17 import java.util.Arrays;
 18 import java.util.Comparator;
 19 import java.util.Stack;
 20 import java.util.TreeSet;
 21  
 22 import org.locationtech.jts.geom.Coordinate;
 23 import org.locationtech.jts.geom.CoordinateArrays;
 24 import org.locationtech.jts.geom.CoordinateList;
 25 import org.locationtech.jts.geom.Geometry;
 26 import org.locationtech.jts.geom.GeometryCollection;
 27 import org.locationtech.jts.geom.GeometryFactory;
 28 import org.locationtech.jts.geom.LineString;
 29 import org.locationtech.jts.geom.LinearRing;
 30 import org.locationtech.jts.geom.Point;
 31 import org.locationtech.jts.geom.Polygon;
 32 import org.locationtech.jts.util.Assert;
 33 import org.locationtech.jts.util.UniqueCoordinateArrayFilter;
 34  
 35 /**
 36  * Computes the convex hull of a {@link Geometry}.
 37  * The convex hull is the smallest convex Geometry that contains all the
 38  * points in the input Geometry.
 39  * <p>
 40  * Uses the Graham Scan algorithm.
 41  *
 42  *@version 1.7
 43  */
 44 public class ConvexHull
 45 {
 46   private GeometryFactory geomFactory;
 47   private Coordinate[] inputPts;
 48  
 49   /**
 50    * Create a new convex hull construction for the input {@link Geometry}.
 51    */
 52   public ConvexHull(Geometry geometry)
 53   {
 54     this(extractCoordinates(geometry), geometry.getFactory());
 55   }
 56   /**
 57    * Create a new convex hull construction for the input {@link Coordinate} array.
 58    */
 59   public ConvexHull(Coordinate[] pts, GeometryFactory geomFactory)
 60   {
 61     inputPts = UniqueCoordinateArrayFilter.filterCoordinates(pts);
 62     //inputPts = pts;
 63     this.geomFactory = geomFactory;
 64   }
 65  
 66   private static Coordinate[] extractCoordinates(Geometry geom)
 67   {
 68     UniqueCoordinateArrayFilter filter = new UniqueCoordinateArrayFilter();
 69     geom.apply(filter);
 70     return filter.getCoordinates();
 71   }
 72  
 73   /**
 74    * Returns a {@link Geometry} that represents the convex hull of the input
 75    * geometry.
 76    * The returned geometry contains the minimal number of points needed to
 77    * represent the convex hull.  In particular, no more than two consecutive
 78    * points will be collinear.
 79    *
 80    * @return if the convex hull contains 3 or more points, a {@link Polygon};
 81    * 2 points, a {@link LineString};
 82    * 1 point, a {@link Point};
 83    * 0 points, an empty {@link GeometryCollection}.
 84    */
 85   public Geometry getConvexHull() {
 86  
 87     if (inputPts.length == 0) {
 88       return geomFactory.createGeometryCollection();
 89     }
 90     if (inputPts.length == 1) {
 91       return geomFactory.createPoint(inputPts[0]);
 92     }
 93     if (inputPts.length == 2) {
 94       return geomFactory.createLineString(inputPts);
 95     }
 96  
 97     Coordinate[] reducedPts = inputPts;
 98     // use heuristic to reduce points, if large
 99     if (inputPts.length > 50) {
100       reducedPts = reduce(inputPts);
101     }
102     // sort points for Graham scan.
103     Coordinate[] sortedPts = preSort(reducedPts);
104  
105     // Use Graham scan to find convex hull.
106     Stack cHS = grahamScan(sortedPts);
107  
108     // Convert stack to an array.
109     Coordinate[] cH = toCoordinateArray(cHS);
110  
111     // Convert array to appropriate output geometry.
112     return lineOrPolygon(cH);
113   }
114  
115   /**
116    * An alternative to Stack.toArray, which is not present in earlier versions
117    * of Java.
118    */
119   protected Coordinate[] toCoordinateArray(Stack stack) {
120     Coordinate[] coordinates = new Coordinate[stack.size()];
121     for (int i = 0; i < stack.size(); i++) {
122       Coordinate coordinate = (Coordinate) stack.get(i);
123       coordinates[i] = coordinate;
124     }
125     return coordinates;
126   }
127  
128   /**
129    * Uses a heuristic to reduce the number of points scanned
130    * to compute the hull.
131    * The heuristic is to find a polygon guaranteed to
132    * be in (or on) the hull, and eliminate all points inside it.
133    * A quadrilateral defined by the extremal points
134    * in the four orthogonal directions
135    * can be used, but even more inclusive is
136    * to use an octilateral defined by the points in the 8 cardinal directions.
137    * <p>
138    * Note that even if the method used to determine the polygon vertices
139    * is not 100% robust, this does not affect the robustness of the convex hull.
140    * <p>
141    * To satisfy the requirements of the Graham Scan algorithm, 
142    * the returned array has at least 3 entries.
143    *
144    * @param pts the points to reduce
145    * @return the reduced list of points (at least 3)
146    */
147   private Coordinate[] reduce(Coordinate[] inputPts)
148   {
149     //Coordinate[] polyPts = computeQuad(inputPts);
150     Coordinate[] polyPts = computeOctRing(inputPts);
151     //Coordinate[] polyPts = null;
152  
153     // unable to compute interior polygon for some reason
154     if (polyPts == null)
155       return inputPts;
156  
157 //    LinearRing ring = geomFactory.createLinearRing(polyPts);
158 //    System.out.println(ring);
159  
160     // add points defining polygon
161     TreeSet reducedSet = new TreeSet();
162     for (int i = 0; i < polyPts.length; i++) {
163       reducedSet.add(polyPts[i]);
164     }
165     /**
166      * Add all unique points not in the interior poly.
167      * CGAlgorithms.isPointInRing is not defined for points actually on the ring,
168      * but this doesn't matter since the points of the interior polygon
169      * are forced to be in the reduced set.
170      */
171     for (int i = 0; i < inputPts.length; i++) {
172       if (! PointLocation.isInRing(inputPts[i], polyPts)) {
173         reducedSet.add(inputPts[i]);
174       }
175     }
176     Coordinate[] reducedPts = CoordinateArrays.toCoordinateArray(reducedSet);
177     
178     // ensure that computed array has at least 3 points (not necessarily unique)  
179     if (reducedPts.length < 3)
180       return padArray3(reducedPts); 
181     return reducedPts;
182   }
183  
184   private Coordinate[] padArray3(Coordinate[] pts)
185   {
186     Coordinate[] pad = new Coordinate[3];
187     for (int i = 0; i < pad.length; i++) {
188       if (i < pts.length) {
189         pad[i] = pts[i];
190       }
191       else
192         pad[i] = pts[0];
193     }
194     return pad;
195   }
196     
197   private Coordinate[] preSort(Coordinate[] pts) {
198     Coordinate t;
199  
200     // find the lowest point in the set. If two or more points have
201     // the same minimum y coordinate choose the one with the minimu x.
202     // This focal point is put in array location pts[0].
203     for (int i = 1; i < pts.length; i++) {
204       if ((pts[i].y < pts[0].y) || ((pts[i].y == pts[0].y) && (pts[i].x < pts[0].x))) {
205         t = pts[0];
206         pts[0] = pts[i];
207         pts[i] = t;
208       }
209     }
210  
211     // sort the points radially around the focal point.
212     Arrays.sort(pts, 1, pts.length, new RadialComparator(pts[0]));
213  
214     //radialSort(pts);
215     return pts;
216   }
217  
218   /**
219    * Uses the Graham Scan algorithm to compute the convex hull vertices.
220    * 
221    * @param c a list of points, with at least 3 entries
222    * @return a Stack containing the ordered points of the convex hull ring
223    */
224   private Stack grahamScan(Coordinate[] c) {
225     Coordinate p;
226     Stack ps = new Stack();
227     ps.push(c[0]);
228     ps.push(c[1]);
229     ps.push(c[2]);
230     for (int i = 3; i < c.length; i++) {
231       p = (Coordinate) ps.pop();
232       // check for empty stack to guard against robustness problems
233       while (
234           ! ps.empty() && 
235           Orientation.index((Coordinate) ps.peek(), p, c[i]) > 0) {
236         p = (Coordinate) ps.pop();
237       }
238       ps.push(p);
239       ps.push(c[i]);
240     }
241     ps.push(c[0]);
242     return ps;
243   }
244  
245   /**
246    *@return    whether the three coordinates are collinear and c2 lies between
247    *      c1 and c3 inclusive
248    */
249   private boolean isBetween(Coordinate c1, Coordinate c2, Coordinate c3) {
250     if (Orientation.index(c1, c2, c3) != 0) {
251       return false;
252     }
253     if (c1.x != c3.x) {
254       if (c1.x <= c2.x && c2.x <= c3.x) {
255         return true;
256       }
257       if (c3.x <= c2.x && c2.x <= c1.x) {
258         return true;
259       }
260     }
261     if (c1.y != c3.y) {
262       if (c1.y <= c2.y && c2.y <= c3.y) {
263         return true;
264       }
265       if (c3.y <= c2.y && c2.y <= c1.y) {
266         return true;
267       }
268     }
269     return false;
270   }
271  
272   private Coordinate[] computeOctRing(Coordinate[] inputPts) {
273     Coordinate[] octPts = computeOctPts(inputPts);
274     CoordinateList coordList = new CoordinateList();
275     coordList.add(octPts, false);
276  
277     // points must all lie in a line
278     if (coordList.size() < 3) {
279       return null;
280     }
281     coordList.closeRing();
282     return coordList.toCoordinateArray();
283   }
284  
285   private Coordinate[] computeOctPts(Coordinate[] inputPts)
286   {
287     Coordinate[] pts = new Coordinate[8];
288     for (int j = 0; j < pts.length; j++) {
289       pts[j] = inputPts[0];
290     }
291     for (int i = 1; i < inputPts.length; i++) {
292       if (inputPts[i].x < pts[0].x) {
293         pts[0] = inputPts[i];
294       }
295       if (inputPts[i].x - inputPts[i].y < pts[1].x - pts[1].y) {
296         pts[1] = inputPts[i];
297       }
298       if (inputPts[i].y > pts[2].y) {
299         pts[2] = inputPts[i];
300       }
301       if (inputPts[i].x + inputPts[i].y > pts[3].x + pts[3].y) {
302         pts[3] = inputPts[i];
303       }
304       if (inputPts[i].x > pts[4].x) {
305         pts[4] = inputPts[i];
306       }
307       if (inputPts[i].x - inputPts[i].y > pts[5].x - pts[5].y) {
308         pts[5] = inputPts[i];
309       }
310       if (inputPts[i].y < pts[6].y) {
311         pts[6] = inputPts[i];
312       }
313       if (inputPts[i].x + inputPts[i].y < pts[7].x + pts[7].y) {
314         pts[7] = inputPts[i];
315       }
316     }
317     return pts;
318  
319   }
320  
321 /*
322   // MD - no longer used, but keep for reference purposes
323   private Coordinate[] computeQuad(Coordinate[] inputPts) {
324     BigQuad bigQuad = bigQuad(inputPts);
325
326     // Build a linear ring defining a big poly.
327     ArrayList bigPoly = new ArrayList();
328     bigPoly.add(bigQuad.westmost);
329     if (! bigPoly.contains(bigQuad.northmost)) {
330       bigPoly.add(bigQuad.northmost);
331     }
332     if (! bigPoly.contains(bigQuad.eastmost)) {
333       bigPoly.add(bigQuad.eastmost);
334     }
335     if (! bigPoly.contains(bigQuad.southmost)) {
336       bigPoly.add(bigQuad.southmost);
337     }
338     // points must all lie in a line
339     if (bigPoly.size() < 3) {
340       return null;
341     }
342     // closing point
343     bigPoly.add(bigQuad.westmost);
344
345     Coordinate[] bigPolyArray = CoordinateArrays.toCoordinateArray(bigPoly);
346
347     return bigPolyArray;
348   }
349
350   private BigQuad bigQuad(Coordinate[] pts) {
351     BigQuad bigQuad = new BigQuad();
352     bigQuad.northmost = pts[0];
353     bigQuad.southmost = pts[0];
354     bigQuad.westmost = pts[0];
355     bigQuad.eastmost = pts[0];
356     for (int i = 1; i < pts.length; i++) {
357       if (pts[i].x < bigQuad.westmost.x) {
358         bigQuad.westmost = pts[i];
359       }
360       if (pts[i].x > bigQuad.eastmost.x) {
361         bigQuad.eastmost = pts[i];
362       }
363       if (pts[i].y < bigQuad.southmost.y) {
364         bigQuad.southmost = pts[i];
365       }
366       if (pts[i].y > bigQuad.northmost.y) {
367         bigQuad.northmost = pts[i];
368       }
369     }
370     return bigQuad;
371   }
372
373   private static class BigQuad {
374     public Coordinate northmost;
375     public Coordinate southmost;
376     public Coordinate westmost;
377     public Coordinate eastmost;
378   }
379   */
380  
381   /**
382    *@param  vertices  the vertices of a linear ring, which may or may not be
383    *      flattened (i.e. vertices collinear)
384    *@return           a 2-vertex <code>LineString</code> if the vertices are
385    *      collinear; otherwise, a <code>Polygon</code> with unnecessary
386    *      (collinear) vertices removed
387    */
388   private Geometry lineOrPolygon(Coordinate[] coordinates) {
389  
390     coordinates = cleanRing(coordinates);
391     if (coordinates.length == 3) {
392       return geomFactory.createLineString(new Coordinate[]{coordinates[0], coordinates[1]});
393 //      return new LineString(new Coordinate[]{coordinates[0], coordinates[1]},
394 //          geometry.getPrecisionModel(), geometry.getSRID());
395     }
396     LinearRing linearRing = geomFactory.createLinearRing(coordinates);
397     return geomFactory.createPolygon(linearRing);
398   }
399  
400   /**
401    *@param  vertices  the vertices of a linear ring, which may or may not be
402    *      flattened (i.e. vertices collinear)
403    *@return           the coordinates with unnecessary (collinear) vertices
404    *      removed
405    */
406   private Coordinate[] cleanRing(Coordinate[] original) {
407     Assert.equals(original[0], original[original.length - 1]);
408     ArrayList cleanedRing = new ArrayList();
409     Coordinate previousDistinctCoordinate = null;
410     for (int i = 0; i <= original.length - 2; i++) {
411       Coordinate currentCoordinate = original[i];
412       Coordinate nextCoordinate = original[i+1];
413       if (currentCoordinate.equals(nextCoordinate)) {
414         continue;
415       }
416       if (previousDistinctCoordinate != null
417           && isBetween(previousDistinctCoordinate, currentCoordinate, nextCoordinate)) {
418         continue;
419       }
420       cleanedRing.add(currentCoordinate);
421       previousDistinctCoordinate = currentCoordinate;
422     }
423     cleanedRing.add(original[original.length - 1]);
424     Coordinate[] cleanedRingCoordinates = new Coordinate[cleanedRing.size()];
425     return (Coordinate[]) cleanedRing.toArray(cleanedRingCoordinates);
426   }
427  
428  
429   /**
430    * Compares {@link Coordinate}s for their angle and distance
431    * relative to an origin.
432    *
433    * @author Martin Davis
434    * @version 1.7
435    */
436   private static class RadialComparator
437       implements Comparator
438   {
439     private Coordinate origin;
440  
441     public RadialComparator(Coordinate origin)
442     {
443       this.origin = origin;
444     }
445     public int compare(Object o1, Object o2)
446     {
447       Coordinate p1 = (Coordinate) o1;
448       Coordinate p2 = (Coordinate) o2;
449       return polarCompare(origin, p1, p2);
450     }
451  
452     /**
453      * Given two points p and q compare them with respect to their radial
454      * ordering about point o.  First checks radial ordering.
455      * If points are collinear, the comparison is based
456      * on their distance to the origin.
457      * <p>
458      * p < q iff
459      * <ul>
460      * <li>ang(o-p) < ang(o-q) (e.g. o-p-q is CCW)
461      * <li>or ang(o-p) == ang(o-q) && dist(o,p) < dist(o,q)
462      * </ul>
463      *
464      * @param o the origin
465      * @param p a point
466      * @param q another point
467      * @return -1, 0 or 1 depending on whether p is less than,
468      * equal to or greater than q
469      */
470     private static int polarCompare(Coordinate o, Coordinate p, Coordinate q)
471     {
472       double dxp = p.x - o.x;
473       double dyp = p.y - o.y;
474       double dxq = q.x - o.x;
475       double dyq = q.y - o.y;
476  
477 /*
478       // MD - non-robust
479       int result = 0;
480       double alph = Math.atan2(dxp, dyp);
481       double beta = Math.atan2(dxq, dyq);
482       if (alph < beta) {
483         result = -1;
484       }
485       if (alph > beta) {
486         result = 1;
487       }
488       if (result !=  0) return result;
489       //*/
490  
491       int orient = Orientation.index(o, p, q);
492  
493       if (orient == Orientation.COUNTERCLOCKWISE) return 1;
494       if (orient == Orientation.CLOCKWISE) return -1;
495  
496       // points are collinear - check distance
497       double op = dxp * dxp + dyp * dyp;
498       double oq = dxq * dxq + dyq * dyq;
499       if (op < oq) {
500         return -1;
501       }
502       if (op > oq) {
503         return 1;
504       }
505       return 0;
506     }
507  
508   }
509 }
510