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