1 /**********************************************************************
3 * GEOS - Geometry Engine Open Source
4 * http://geos.osgeo.org
6 * Copyright (C) 2011 Sandro Santilli <strk@kbt.io>
7 * Copyright (C) 2005-2006 Refractions Research Inc.
8 * Copyright (C) 2001-2002 Vivid Solutions Inc.
10 * This is free software; you can redistribute and/or modify it under
11 * the terms of the GNU Lesser General Public Licence as published
12 * by the Free Software Foundation.
13 * See the COPYING file for more information.
15 **********************************************************************
17 * Last port: algorithm/ConvexHull.java r407 (JTS-1.12+)
19 **********************************************************************/
21 #include <geos/algorithm/ConvexHull.h>
22 #include <geos/algorithm/CGAlgorithms.h>
23 #include <geos/geom/GeometryFactory.h>
24 #include <geos/geom/Coordinate.h>
25 #include <geos/geom/Point.h>
26 #include <geos/geom/Polygon.h>
27 #include <geos/geom/LineString.h>
28 #include <geos/geom/CoordinateSequence.h>
29 #include <geos/geom/CoordinateSequenceFactory.h>
30 #include <geos/util/Interrupt.h>
36 # include "geos/algorithm/ConvexHull.inl"
39 using namespace geos::geom
;
42 namespace algorithm
{ // geos.algorithm
47 * This is the class used to sort in radial order.
48 * It's defined here as it's a static one.
50 class RadiallyLessThen
{
54 const Coordinate
*origin
;
57 polarCompare(const Coordinate
*o
, const Coordinate
*p
,
65 int orient
= CGAlgorithms::computeOrientation(*o
, *p
, *q
);
68 if (orient
== CGAlgorithms::COUNTERCLOCKWISE
) return 1;
69 if (orient
== CGAlgorithms::CLOCKWISE
) return -1;
71 // points are collinear - check distance
72 double op
= dxp
* dxp
+ dyp
* dyp
;
73 double oq
= dxq
* dxq
+ dyq
* dyq
;
85 RadiallyLessThen(const Coordinate
*c
): origin(c
) {}
87 bool operator() (const Coordinate
*p1
, const Coordinate
*p2
)
89 return ( polarCompare(origin
, p1
, p2
) == -1 );
94 } // unnamed namespace
98 ConvexHull::toCoordinateSequence(Coordinate::ConstVect
&cv
)
100 const CoordinateSequenceFactory
*csf
=
101 geomFactory
->getCoordinateSequenceFactory();
103 // Create a new Coordinate::Vect for feeding it to
104 // the CoordinateSequenceFactory
105 Coordinate::Vect
*vect
=new Coordinate::Vect();
108 vect
->reserve(n
); // avoid multiple reallocs
110 for (size_t i
=0; i
<n
; ++i
)
112 vect
->push_back(*(cv
[i
])); // Coordinate copy
115 return csf
->create(vect
); // takes ownership of the vector
121 ConvexHull::computeOctPts(const Coordinate::ConstVect
&inputPts
,
122 Coordinate::ConstVect
&pts
)
124 // Initialize all slots with first input coordinate
125 pts
= Coordinate::ConstVect(8, inputPts
[0]);
127 for (size_t i
=1, n
=inputPts
.size(); i
<n
; ++i
)
129 if (inputPts
[i
]->x
< pts
[0]->x
) {
130 pts
[0] = inputPts
[i
];
132 if (inputPts
[i
]->x
- inputPts
[i
]->y
< pts
[1]->x
- pts
[1]->y
) {
133 pts
[1] = inputPts
[i
];
135 if (inputPts
[i
]->y
> pts
[2]->y
) {
136 pts
[2] = inputPts
[i
];
138 if (inputPts
[i
]->x
+ inputPts
[i
]->y
> pts
[3]->x
+ pts
[3]->y
) {
139 pts
[3] = inputPts
[i
];
141 if (inputPts
[i
]->x
> pts
[4]->x
) {
142 pts
[4] = inputPts
[i
];
144 if (inputPts
[i
]->x
- inputPts
[i
]->y
> pts
[5]->x
- pts
[5]->y
) {
145 pts
[5] = inputPts
[i
];
147 if (inputPts
[i
]->y
< pts
[6]->y
) {
148 pts
[6] = inputPts
[i
];
150 if (inputPts
[i
]->x
+ inputPts
[i
]->y
< pts
[7]->x
+ pts
[7]->y
) {
151 pts
[7] = inputPts
[i
];
160 ConvexHull::computeOctRing(const Coordinate::ConstVect
&inputPts
,
161 Coordinate::ConstVect
&dest
)
163 computeOctPts(inputPts
, dest
);
165 // Remove consecutive equal Coordinates
166 // unique() returns an iterator to the end of the resulting
167 // sequence, we erase from there to the end.
168 dest
.erase( std::unique(dest
.begin(),dest
.end()), dest
.end() );
170 // points must all lie in a line
171 if ( dest
.size() < 3 ) return false;
174 dest
.push_back(dest
[0]);
181 ConvexHull::reduce(Coordinate::ConstVect
&pts
)
183 Coordinate::ConstVect polyPts
;
185 if ( ! computeOctRing(pts
, polyPts
) ) {
186 // unable to compute interior polygon for some reason
190 // add points defining polygon
191 Coordinate::ConstSet reducedSet
;
192 reducedSet
.insert(polyPts
.begin(), polyPts
.end());
195 * Add all unique points not in the interior poly.
196 * CGAlgorithms.isPointInRing is not defined for points
197 * actually on the ring, but this doesn't matter since
198 * the points of the interior polygon are forced to be
199 * in the reduced set.
201 * @@TIP: there should be a std::algo for this
203 for (size_t i
=0, n
=pts
.size(); i
<n
; ++i
)
205 if ( !CGAlgorithms::isPointInRing(*(pts
[i
]), polyPts
) )
207 reducedSet
.insert(pts
[i
]);
211 inputPts
.assign(reducedSet
.begin(), reducedSet
.end());
213 if ( inputPts
.size() < 3 ) padArray3(inputPts
);
218 ConvexHull::padArray3(geom::Coordinate::ConstVect
&pts
)
220 for (size_t i
= pts
.size(); i
< 3; ++i
) {
221 pts
.push_back(pts
[0]);
226 ConvexHull::getConvexHull()
228 size_t nInputPts
=inputPts
.size();
230 if (nInputPts
==0) // Return an empty geometry
231 return geomFactory
->createEmptyGeometry();
233 if (nInputPts
==1) // Return a Point
235 // Copy the Coordinate from the ConstVect
236 return geomFactory
->createPoint(*(inputPts
[0]));
239 if (nInputPts
==2) // Return a LineString
241 // Copy all Coordinates from the ConstVect
242 CoordinateSequence
*cs
= toCoordinateSequence(inputPts
);
243 return geomFactory
->createLineString(cs
);
246 // use heuristic to reduce points, if large
252 GEOS_CHECK_FOR_INTERRUPTS();
254 // sort points for Graham scan.
257 GEOS_CHECK_FOR_INTERRUPTS();
259 // Use Graham scan to find convex hull.
260 Coordinate::ConstVect cHS
;
261 grahamScan(inputPts
, cHS
);
263 GEOS_CHECK_FOR_INTERRUPTS();
265 return lineOrPolygon(cHS
);
271 ConvexHull::preSort(Coordinate::ConstVect
&pts
)
273 // find the lowest point in the set. If two or more points have
274 // the same minimum y coordinate choose the one with the minimum x.
275 // This focal point is put in array location pts[0].
276 for(size_t i
=1, n
=pts
.size(); i
<n
; ++i
)
278 const Coordinate
*p0
=pts
[0]; // this will change
279 const Coordinate
*pi
=pts
[i
];
280 if ( (pi
->y
<p0
->y
) || ((pi
->y
==p0
->y
) && (pi
->x
<p0
->x
)) )
282 const Coordinate
*t
=p0
;
288 // sort the points radially around the focal point.
289 std::sort(pts
.begin(), pts
.end(), RadiallyLessThen(pts
[0]));
294 ConvexHull::grahamScan(const Coordinate::ConstVect
&c
,
295 Coordinate::ConstVect
&ps
)
301 for(size_t i
=3, n
=c
.size(); i
<n
; ++i
)
303 const Coordinate
*p
= ps
.back(); ps
.pop_back();
304 while (!ps
.empty() &&
305 CGAlgorithms::computeOrientation(
306 *(ps
.back()), *p
, *(c
[i
])) > 0)
308 p
= ps
.back(); ps
.pop_back();
318 // * Returns a pointer to input, modifying input
320 //CoordinateSequence*
321 //ConvexHull::preSort(CoordinateSequence *pts)
325 // // find the lowest point in the set. If two or more points have
326 // // the same minimum y coordinate choose the one with the minimu x.
327 // // This focal point is put in array location pts[0].
328 // size_t npts=pts->getSize();
329 // for(size_t i=1; i<npts; ++i)
331 // const Coordinate &p0=pts->getAt(0); // this will change
332 // const Coordinate &pi=pts->getAt(i);
333 // if ( (pi.y<p0.y) || ((pi.y==p0.y) && (pi.x<p0.x)) )
336 // pts->setAt(pi, 0);
337 // pts->setAt( t, i);
340 // // sort the points radially around the focal point.
346 // * Returns a newly allocated CoordinateSequence object
348 //CoordinateSequence*
349 //ConvexHull::grahamScan(const CoordinateSequence *c)
351 // const Coordinate *p;
353 // vector<Coordinate> *ps=new vector<Coordinate>();
354 // ps->push_back(c->getAt(0));
355 // ps->push_back(c->getAt(1));
356 // ps->push_back(c->getAt(2));
359 // size_t npts=c->getSize();
360 // for(size_t i=3; i<npts; ++i)
364 // while (CGAlgorithms::computeOrientation(ps->back(), *p, c->getAt(i)) > 0)
369 // ps->push_back(*p);
370 // ps->push_back(c->getAt(i));
372 // ps->push_back(c->getAt(0));
374 // return geomFactory->getCoordinateSequenceFactory()->create(ps);
378 //ConvexHull::radialSort(CoordinateSequence *p)
380 // // A selection sort routine, assumes the pivot point is
381 // // the first point (i.e., p[0]).
383 // const Coordinate &p0=p->getAt(0); // the pivot point
386 // size_t npts=p->getSize();
387 // for(size_t i=1; i<npts-1; ++i)
391 // for(size_t j=i+1; j<npts; ++j)
393 // const Coordinate &pj=p->getAt(j);
395 // if ( polarCompare(p0, pj, p->getAt(min)) < 0 )
402 // * Swap point[i] and point[min]
403 // * We can skip this if they have
409 // p->setAt(p->getAt(min), i);
419 ConvexHull::isBetween(const Coordinate
&c1
, const Coordinate
&c2
, const Coordinate
&c3
)
421 if (CGAlgorithms::computeOrientation(c1
, c2
, c3
)!=0) {
425 if (c1
.x
<=c2
.x
&& c2
.x
<=c3
.x
) {
428 if (c3
.x
<=c2
.x
&& c2
.x
<=c1
.x
) {
433 if (c1
.y
<=c2
.y
&& c2
.y
<=c3
.y
) {
436 if (c3
.y
<=c2
.y
&& c2
.y
<=c1
.y
) {
444 //ConvexHull::makeBigQuad(const CoordinateSequence *pts, BigQuad &bigQuad)
446 // bigQuad.northmost=bigQuad.southmost=
447 // bigQuad.westmost=bigQuad.eastmost=pts->getAt(0);
449 // size_t npts=pts->getSize();
450 // for (size_t i=1; i<npts; ++i)
452 // const Coordinate &pi=pts->getAt(i);
454 // if (pi.x<bigQuad.westmost.x)
455 // bigQuad.westmost=pi;
457 // if (pi.x>bigQuad.eastmost.x)
458 // bigQuad.eastmost=pi;
460 // if (pi.y<bigQuad.southmost.y)
461 // bigQuad.southmost=pi;
463 // if (pi.y>bigQuad.northmost.y)
464 // bigQuad.northmost=pi;
470 ConvexHull::lineOrPolygon(const Coordinate::ConstVect
&input
)
472 Coordinate::ConstVect cleaned
;
474 cleanRing(input
, cleaned
);
476 if (cleaned
.size()==3) { // shouldn't this be 2 ??
478 CoordinateSequence
*cl1
=toCoordinateSequence(cleaned
);
479 LineString
*ret
= geomFactory
->createLineString(cl1
);
482 CoordinateSequence
*cl2
=toCoordinateSequence(cleaned
);
483 LinearRing
*linearRing
=geomFactory
->createLinearRing(cl2
);
484 return geomFactory
->createPolygon(linearRing
,NULL
);
489 ConvexHull::cleanRing(const Coordinate::ConstVect
&original
,
490 Coordinate::ConstVect
&cleaned
)
492 size_t npts
=original
.size();
494 const Coordinate
*last
= original
[npts
-1];
496 //util::Assert::equals(*(original[0]), *last);
498 assert(original
[0]->equals2D(*last
));
500 const Coordinate
*prev
= NULL
;
501 for (size_t i
=0; i
<npts
-1; ++i
)
503 const Coordinate
*curr
= original
[i
];
504 const Coordinate
*next
= original
[i
+1];
506 // skip consecutive equal coordinates
507 if (curr
->equals2D(*next
)) continue;
509 if ( prev
!= NULL
&& isBetween(*prev
, *curr
, *next
) )
514 cleaned
.push_back(curr
);
518 cleaned
.push_back(last
);
523 } // namespace geos.algorithm