src/algorithm/ConvexHull.cpp

   1 /**********************************************************************
   2  *
   3  * GEOS - Geometry Engine Open Source
   4  * http://geos.osgeo.org
   5  *
   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.
   9  *
  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.
  14  *
  15  **********************************************************************
  16  *
  17  * Last port: algorithm/ConvexHull.java r407 (JTS-1.12+)
  18  *
  19  **********************************************************************/
  20
  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>
  31
  32 #include <typeinfo>
  33 #include <algorithm>
  34
  35 #ifndef GEOS_INLINE
  36 # include "geos/algorithm/ConvexHull.inl"
  37 #endif
  38
  39 using namespace geos::geom;
  40
  41 namespace geos {
  42 namespace algorithm { // geos.algorithm
  43
  44 namespace {
  45
  46 /**
  47  * This is the class used to sort in radial order.
  48  * It's defined here as it's a static one.
  49  */
  50 class RadiallyLessThen {
  51
  52 private:
  53
  54         const Coordinate *origin;
  55
  56         int
  57         polarCompare(const Coordinate *o, const Coordinate *p,
  58                         const Coordinate *q)
  59         {
  60                 double dxp=p->x-o->x;
  61                 double dyp=p->y-o->y;
  62                 double dxq=q->x-o->x;
  63                 double dyq=q->y-o->y;
  64
  65                 int orient = CGAlgorithms::computeOrientation(*o, *p, *q);
  66
  67
  68                 if (orient == CGAlgorithms::COUNTERCLOCKWISE) return 1;
  69                 if (orient == CGAlgorithms::CLOCKWISE) return -1;
  70
  71                 // points are collinear - check distance
  72                 double op = dxp * dxp + dyp * dyp;
  73                 double oq = dxq * dxq + dyq * dyq;
  74                 if (op < oq) {
  75                         return -1;
  76                 }
  77                 if (op > oq) {
  78                         return 1;
  79                 }
  80                 return 0;
  81         }
  82
  83 public:
  84
  85         RadiallyLessThen(const Coordinate *c): origin(c) {}
  86
  87         bool operator() (const Coordinate *p1, const Coordinate *p2)
  88         {
  89                 return  ( polarCompare(origin, p1, p2) == -1 );
  90         }
  91
  92 };
  93
  94 } // unnamed namespace
  95
  96 /* private */
  97 CoordinateSequence *
  98 ConvexHull::toCoordinateSequence(Coordinate::ConstVect &cv)
  99 {
 100         const CoordinateSequenceFactory *csf =
 101                 geomFactory->getCoordinateSequenceFactory();
 102
 103         // Create a new Coordinate::Vect for feeding it to
 104         // the CoordinateSequenceFactory
 105         Coordinate::Vect *vect=new Coordinate::Vect();
 106
 107         size_t n=cv.size();
 108         vect->reserve(n); // avoid multiple reallocs
 109
 110         for (size_t i=0; i<n; ++i)
 111         {
 112                 vect->push_back(*(cv[i])); // Coordinate copy
 113         }
 114
 115         return csf->create(vect); // takes ownership of the vector
 116
 117 }
 118
 119 /* private */
 120 void
 121 ConvexHull::computeOctPts(const Coordinate::ConstVect &inputPts,
 122                 Coordinate::ConstVect &pts)
 123 {
 124         // Initialize all slots with first input coordinate
 125         pts = Coordinate::ConstVect(8, inputPts[0]);
 126
 127         for (size_t i=1, n=inputPts.size(); i<n; ++i)
 128         {
 129                 if (inputPts[i]->x < pts[0]->x) {
 130                         pts[0] = inputPts[i];
 131                 }
 132                 if (inputPts[i]->x - inputPts[i]->y < pts[1]->x - pts[1]->y) {
 133                         pts[1] = inputPts[i];
 134                 }
 135                 if (inputPts[i]->y > pts[2]->y) {
 136                         pts[2] = inputPts[i];
 137                 }
 138                 if (inputPts[i]->x + inputPts[i]->y > pts[3]->x + pts[3]->y) {
 139                         pts[3] = inputPts[i];
 140                 }
 141                 if (inputPts[i]->x > pts[4]->x) {
 142                         pts[4] = inputPts[i];
 143                 }
 144                 if (inputPts[i]->x - inputPts[i]->y > pts[5]->x - pts[5]->y) {
 145                         pts[5] = inputPts[i];
 146                 }
 147                 if (inputPts[i]->y < pts[6]->y) {
 148                         pts[6] = inputPts[i];
 149                 }
 150                 if (inputPts[i]->x + inputPts[i]->y < pts[7]->x + pts[7]->y) {
 151                         pts[7] = inputPts[i];
 152                 }
 153         }
 154
 155 }
 156
 157
 158 /* private */
 159 bool
 160 ConvexHull::computeOctRing(const Coordinate::ConstVect &inputPts,
 161         Coordinate::ConstVect &dest)
 162 {
 163         computeOctPts(inputPts, dest);
 164
 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() );
 169
 170         // points must all lie in a line
 171         if ( dest.size() < 3 ) return false;
 172
 173         // close ring
 174         dest.push_back(dest[0]);
 175
 176         return true;
 177 }
 178
 179 /* private */
 180 void
 181 ConvexHull::reduce(Coordinate::ConstVect &pts)
 182 {
 183         Coordinate::ConstVect polyPts;
 184
 185         if ( ! computeOctRing(pts, polyPts) ) {
 186         // unable to compute interior polygon for some reason
 187                 return;
 188         }
 189
 190         // add points defining polygon
 191         Coordinate::ConstSet reducedSet;
 192         reducedSet.insert(polyPts.begin(), polyPts.end());
 193
 194         /**
 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.
 200          *
 201          * @@TIP: there should be a std::algo for this
 202          */
 203         for (size_t i=0, n=pts.size(); i<n; ++i)
 204         {
 205                 if ( !CGAlgorithms::isPointInRing(*(pts[i]), polyPts) )
 206                 {
 207                         reducedSet.insert(pts[i]);
 208                 }
 209         }
 210
 211         inputPts.assign(reducedSet.begin(), reducedSet.end());
 212
 213         if ( inputPts.size() < 3 ) padArray3(inputPts);
 214 }
 215
 216 /* private */
 217 void
 218 ConvexHull::padArray3(geom::Coordinate::ConstVect &pts)
 219 {
 220   for (size_t i = pts.size(); i < 3; ++i) {
 221     pts.push_back(pts[0]);
 222   }
 223 }
 224
 225 Geometry*
 226 ConvexHull::getConvexHull()
 227 {
 228         size_t nInputPts=inputPts.size();
 229
 230         if (nInputPts==0) // Return an empty geometry
 231                 return geomFactory->createEmptyGeometry();
 232
 233         if (nInputPts==1) // Return a Point
 234         {
 235                 // Copy the Coordinate from the ConstVect
 236                 return geomFactory->createPoint(*(inputPts[0]));
 237         }
 238
 239         if (nInputPts==2) // Return a LineString
 240         {
 241                 // Copy all Coordinates from the ConstVect
 242                 CoordinateSequence *cs = toCoordinateSequence(inputPts);
 243                 return geomFactory->createLineString(cs);
 244         }
 245
 246         // use heuristic to reduce points, if large
 247         if (nInputPts > 50 )
 248         {
 249                 reduce(inputPts);
 250         }
 251
 252     GEOS_CHECK_FOR_INTERRUPTS();
 253
 254         // sort points for Graham scan.
 255         preSort(inputPts);
 256
 257     GEOS_CHECK_FOR_INTERRUPTS();
 258
 259         // Use Graham scan to find convex hull.
 260         Coordinate::ConstVect cHS;
 261         grahamScan(inputPts, cHS);
 262
 263     GEOS_CHECK_FOR_INTERRUPTS();
 264
 265         return lineOrPolygon(cHS);
 266
 267 }
 268
 269 /* private */
 270 void
 271 ConvexHull::preSort(Coordinate::ConstVect &pts)
 272 {
 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)
 277         {
 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)) )
 281                 {
 282                         const Coordinate *t=p0;
 283                         pts[0]=pi;
 284                         pts[i]=t;
 285                 }
 286         }
 287
 288         // sort the points radially around the focal point.
 289     std::sort(pts.begin(), pts.end(), RadiallyLessThen(pts[0]));
 290 }
 291
 292 /*private*/
 293 void
 294 ConvexHull::grahamScan(const Coordinate::ConstVect &c,
 295                 Coordinate::ConstVect &ps)
 296 {
 297         ps.push_back(c[0]);
 298         ps.push_back(c[1]);
 299         ps.push_back(c[2]);
 300
 301         for(size_t i=3, n=c.size(); i<n; ++i)
 302         {
 303                 const Coordinate *p = ps.back(); ps.pop_back();
 304                 while (!ps.empty() &&
 305                         CGAlgorithms::computeOrientation(
 306                     *(ps.back()), *p, *(c[i])) > 0)
 307                 {
 308                         p = ps.back(); ps.pop_back();
 309                 }
 310                 ps.push_back(p);
 311                 ps.push_back(c[i]);
 312         }
 313         ps.push_back(c[0]);
 314 }
 315
 316
 317 ///*
 318 // * Returns a pointer to input, modifying input
 319 // */
 320 //CoordinateSequence*
 321 //ConvexHull::preSort(CoordinateSequence *pts)
 322 //{
 323 //      Coordinate t;
 324 //
 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)
 330 //      {
 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)) )
 334 //              {
 335 //                      t=p0;
 336 //                      pts->setAt(pi, 0);
 337 //                      pts->setAt( t, i);
 338 //              }
 339 //      }
 340 //      // sort the points radially around the focal point.
 341 //      radialSort(pts);
 342 //      return pts;
 343 //}
 344
 345 ///*
 346 // * Returns a newly allocated CoordinateSequence object
 347 // */
 348 //CoordinateSequence*
 349 //ConvexHull::grahamScan(const CoordinateSequence *c)
 350 //{
 351 //      const Coordinate *p;
 352 //
 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));
 357 //
 358 //      p=&(c->getAt(2));
 359 //      size_t npts=c->getSize();
 360 //      for(size_t i=3; i<npts; ++i)
 361 //      {
 362 //              p=&(ps->back());
 363 //              ps->pop_back();
 364 //              while (CGAlgorithms::computeOrientation(ps->back(), *p, c->getAt(i)) > 0)
 365 //              {
 366 //                      p=&(ps->back());
 367 //                      ps->pop_back();
 368 //              }
 369 //              ps->push_back(*p);
 370 //              ps->push_back(c->getAt(i));
 371 //      }
 372 //      ps->push_back(c->getAt(0));
 373 //
 374 //      return geomFactory->getCoordinateSequenceFactory()->create(ps);
 375 //}
 376
 377 //void
 378 //ConvexHull::radialSort(CoordinateSequence *p)
 379 //{
 380 //      // A selection sort routine, assumes the pivot point is
 381 //      // the first point (i.e., p[0]).
 382 //
 383 //      const Coordinate &p0=p->getAt(0); // the pivot point
 384 //
 385 //      Coordinate t;
 386 //      size_t npts=p->getSize();
 387 //      for(size_t i=1; i<npts-1; ++i)
 388 //      {
 389 //              size_t min=i;
 390 //
 391 //              for(size_t j=i+1; j<npts; ++j)
 392 //              {
 393 //                      const Coordinate &pj=p->getAt(j);
 394 //
 395 //                      if ( polarCompare(p0, pj, p->getAt(min)) < 0 )
 396 //                      {
 397 //                              min=j;
 398 //                      }
 399 //              }
 400 //
 401 //              /*
 402 //               * Swap point[i] and point[min]
 403 //               * We can skip this if they have
 404 //               * the same value
 405 //               */
 406 //              if ( i != min )
 407 //              {
 408 //                      t=p->getAt(i);
 409 //                      p->setAt(p->getAt(min), i);
 410 //                      p->setAt(t, min);
 411 //              }
 412 //      }
 413 //}
 414
 415
 416
 417 /*private*/
 418 bool
 419 ConvexHull::isBetween(const Coordinate &c1, const Coordinate &c2, const Coordinate &c3)
 420 {
 421         if (CGAlgorithms::computeOrientation(c1, c2, c3)!=0) {
 422                 return false;
 423         }
 424         if (c1.x!=c3.x) {
 425                 if (c1.x<=c2.x && c2.x<=c3.x) {
 426                         return true;
 427                 }
 428                 if (c3.x<=c2.x && c2.x<=c1.x) {
 429                         return true;
 430                 }
 431         }
 432         if (c1.y!=c3.y) {
 433                 if (c1.y<=c2.y && c2.y<=c3.y) {
 434                         return true;
 435                 }
 436                 if (c3.y<=c2.y && c2.y<=c1.y) {
 437                         return true;
 438                 }
 439         }
 440         return false;
 441 }
 442
 443 //void
 444 //ConvexHull::makeBigQuad(const CoordinateSequence *pts, BigQuad &bigQuad)
 445 //{
 446 //      bigQuad.northmost=bigQuad.southmost=
 447 //              bigQuad.westmost=bigQuad.eastmost=pts->getAt(0);
 448 //
 449 //      size_t npts=pts->getSize();
 450 //      for (size_t i=1; i<npts; ++i)
 451 //      {
 452 //              const Coordinate &pi=pts->getAt(i);
 453 //
 454 //              if (pi.x<bigQuad.westmost.x)
 455 //                      bigQuad.westmost=pi;
 456 //
 457 //              if (pi.x>bigQuad.eastmost.x)
 458 //                      bigQuad.eastmost=pi;
 459 //
 460 //              if (pi.y<bigQuad.southmost.y)
 461 //                      bigQuad.southmost=pi;
 462 //
 463 //              if (pi.y>bigQuad.northmost.y)
 464 //                      bigQuad.northmost=pi;
 465 //      }
 466 //}
 467
 468 /* private */
 469 Geometry*
 470 ConvexHull::lineOrPolygon(const Coordinate::ConstVect &input)
 471 {
 472         Coordinate::ConstVect cleaned;
 473
 474         cleanRing(input, cleaned);
 475
 476         if (cleaned.size()==3) { // shouldn't this be 2 ??
 477                 cleaned.resize(2);
 478                 CoordinateSequence *cl1=toCoordinateSequence(cleaned);
 479                 LineString *ret = geomFactory->createLineString(cl1);
 480                 return ret;
 481         }
 482         CoordinateSequence *cl2=toCoordinateSequence(cleaned);
 483         LinearRing *linearRing=geomFactory->createLinearRing(cl2);
 484         return geomFactory->createPolygon(linearRing,NULL);
 485 }
 486
 487 /*private*/
 488 void
 489 ConvexHull::cleanRing(const Coordinate::ConstVect &original,
 490                 Coordinate::ConstVect &cleaned)
 491 {
 492         size_t npts=original.size();
 493
 494         const Coordinate *last = original[npts-1];
 495
 496         //util::Assert::equals(*(original[0]), *last);
 497         assert(last);
 498         assert(original[0]->equals2D(*last));
 499
 500         const Coordinate *prev = NULL;
 501         for (size_t i=0; i<npts-1; ++i)
 502         {
 503                 const Coordinate *curr = original[i];
 504                 const Coordinate *next = original[i+1];
 505
 506                 // skip consecutive equal coordinates
 507                 if (curr->equals2D(*next)) continue;
 508
 509                 if ( prev != NULL &&  isBetween(*prev, *curr, *next) )
 510                 {
 511                         continue;
 512                 }
 513
 514                 cleaned.push_back(curr);
 515                 prev=curr;
 516         }
 517
 518         cleaned.push_back(last);
 519
 520 }
 521
 522
 523 } // namespace geos.algorithm
 524 } // namespace geos
 525