src/2geom/convex-cover.cpp

   1 /*
   2  * convex-cover.cpp
   3  *
   4  * Copyright 2006 Nathan Hurst <njh@mail.csse.monash.edu.au>
   5  * Copyright 2006 Michael G. Sloan <mgsloan@gmail.com>
   6  *
   7  * This library is free software; you can redistribute it and/or
   8  * modify it either under the terms of the GNU Lesser General Public
   9  * License version 2.1 as published by the Free Software Foundation
  10  * (the "LGPL") or, at your option, under the terms of the Mozilla
  11  * Public License Version 1.1 (the "MPL"). If you do not alter this
  12  * notice, a recipient may use your version of this file under either
  13  * the MPL or the LGPL.
  14  *
  15  * You should have received a copy of the LGPL along with this library
  16  * in the file COPYING-LGPL-2.1; if not, write to the Free Software
  17  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
  18  * You should have received a copy of the MPL along with this library
  19  * in the file COPYING-MPL-1.1
  20  *
  21  * The contents of this file are subject to the Mozilla Public License
  22  * Version 1.1 (the "License"); you may not use this file except in
  23  * compliance with the License. You may obtain a copy of the License at
  24  * http://www.mozilla.org/MPL/
  25  *
  26  * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
  27  * OF ANY KIND, either express or implied. See the LGPL or the MPL for
  28  * the specific language governing rights and limitations.
  29  *
  30  */
  31
  32 #include <2geom/convex-cover.h>
  33 #include <2geom/exception.h>
  34 #include <algorithm>
  35 #include <map>
  36
  37 /** Todo:
  38     + modify graham scan to work top to bottom, rather than around angles
  39     + intersection
  40     + minimum distance between convex hulls
  41     + maximum distance between convex hulls
  42     + hausdorf metric?
  43     + check all degenerate cases carefully
  44     + check all algorithms meet all invariants
  45     + generalise rotating caliper algorithm (iterator/circulator?)
  46 */
  47
  48 using std::vector;
  49 using std::map;
  50 using std::pair;
  51
  52 namespace Geom{
  53
  54 /*** SignedTriangleArea
  55  * returns the area of the triangle defined by p0, p1, p2.  A clockwise triangle has positive area.
  56  */
  57 double
  58 SignedTriangleArea(Point p0, Point p1, Point p2) {
  59     return cross((p1 - p0), (p2 - p0));
  60 }
  61
  62 class angle_cmp{
  63 public:
  64     Point o;
  65     angle_cmp(Point o) : o(o) {}
  66
  67 #if 0
  68     bool
  69     operator()(Point a, Point b) {
  70         // not remove this check or std::sort could crash
  71         if (a == b) return false;
  72         Point da = a - o;
  73         Point db = b - o;
  74         if (da == -db) return false;
  75
  76 #if 1
  77         double aa = da[0];
  78         double ab = db[0];
  79         if((da[1] == 0) && (db[1] == 0))
  80             return da[0] < db[0];
  81         if(da[1] == 0)
  82             return true; // infinite tangent
  83         if(db[1] == 0)
  84             return false; // infinite tangent
  85         aa = da[0] / da[1];
  86         ab = db[0] / db[1];
  87         if(aa > ab)
  88             return true;
  89 #else
  90         //assert((ata > atb) == (aa < ab));
  91         double aa = atan2(da);
  92         double ab = atan2(db);
  93         if(aa < ab)
  94             return true;
  95 #endif
  96         if(aa == ab)
  97             return L2sq(da) < L2sq(db);
  98         return false;
  99     }
 100 #else
 101     bool operator() (Point const& a, Point const&  b)
 102     {
 103         // not remove this check or std::sort could generate
 104         // a segmentation fault because it needs a strict '<'
 105         // but due to round errors a == b doesn't mean dxy == dyx
 106         if (a == b) return false;
 107         Point da = a - o;
 108         Point db = b - o;
 109         if (da == -db) return false;
 110         double dxy = da[X] * db[Y];
 111         double dyx = da[Y] * db[X];
 112         if (dxy > dyx) return true;
 113         else if (dxy < dyx) return false;
 114         return L2sq(da) < L2sq(db);
 115     }
 116 #endif
 117 };
 118
 119 void
 120 ConvexHull::find_pivot() {
 121     // Find pivot P;
 122     unsigned pivot = 0;
 123     for (unsigned i = 1; i < boundary.size(); i++)
 124         if(boundary[i] <= boundary[pivot])
 125             pivot = i;
 126
 127     std::swap(boundary[0], boundary[pivot]);
 128 }
 129
 130 void
 131 ConvexHull::angle_sort() {
 132 // sort points by angle (resolve ties in favor of point farther from P);
 133 // we leave the first one in place as our pivot
 134     std::sort(boundary.begin()+1, boundary.end(), angle_cmp(boundary[0]));
 135 }
 136
 137
 138 void
 139 ConvexHull::graham_scan() {
 140     if (boundary.size() < 4) return;
 141     unsigned stac = 2;
 142     for(unsigned int i = 2; i < boundary.size(); i++) {
 143         double o = SignedTriangleArea(boundary[stac-2],
 144                                       boundary[stac-1],
 145                                       boundary[i]);
 146         if(o == 0) { // colinear - dangerous...
 147             stac--;
 148         } else if(o < 0) { // anticlockwise
 149         } else { // remove concavity
 150             while(o >= 0 && stac > 2) {
 151                 stac--;
 152                 o = SignedTriangleArea(boundary[stac-2],
 153                                        boundary[stac-1],
 154                                        boundary[i]);
 155             }
 156         }
 157         boundary[stac++] = boundary[i];
 158     }
 159     boundary.resize(stac);
 160 }
 161
 162 void
 163 ConvexHull::graham() {
 164     if(is_degenerate()) // nothing to do
 165         return;
 166     find_pivot();
 167     angle_sort();
 168     graham_scan();
 169 }
 170
 171 //Mathematically incorrect mod, but more useful.
 172 int mod(int i, int l) {
 173     return i >= 0 ?
 174            i % l : (i % l) + l;
 175 }
 176 //OPT: usages can often be replaced by conditions
 177
 178 /*** ConvexHull::left
 179  * Tests if a point is left (outside) of a particular segment, n. */
 180 bool
 181 ConvexHull::is_left(Point p, int n) {
 182     return SignedTriangleArea((*this)[n], (*this)[n+1], p) >= 0;
 183 }
 184
 185 /*** ConvexHull::strict_left
 186  * Tests if a point is left (outside) of a particular segment, n. */
 187 bool
 188 ConvexHull::is_strict_left(Point p, int n) {
 189     return SignedTriangleArea((*this)[n], (*this)[n+1], p) > 0;
 190 }
 191
 192 /*** ConvexHull::find_positive
 193  * May return any number n where the segment n -> n + 1 (possibly looped around) in the hull such
 194  * that the point is on the wrong side to be within the hull.  Returns -1 if it is within the hull.*/
 195 int
 196 ConvexHull::find_left(Point p) {
 197     int l = boundary.size(); //Who knows if C++ is smart enough to optimize this?
 198     for(int i = 0; i < l; i++) {
 199         if(is_left(p, i)) return i;
 200     }
 201     return -1;
 202 }
 203
 204
 205 /*** ConvexHull::find_positive
 206  * May return any number n where the segment n -> n + 1 (possibly looped around) in the hull such
 207  * that the point is on the wrong side to be within the hull.  Returns -1 if it is within the hull.*/
 208 int
 209 ConvexHull::find_strict_left(Point p) {
 210     int l = boundary.size(); //Who knows if C++ is smart enough to optimize this?
 211     for(int i = 0; i < l; i++) {
 212         if(is_strict_left(p, i)) return i;
 213     }
 214     return -1;
 215 }
 216
 217 //OPT: do a spread iteration - quasi-random with no repeats and full coverage.
 218
 219 /*** ConvexHull::contains_point
 220  * In order to test whether a point is inside a convex hull we can travel once around the outside making
 221  * sure that each triangle made from an edge and the point has positive area. */
 222 bool
 223 ConvexHull::contains_point(Point p) {
 224     return find_left(p) == -1;
 225 }
 226
 227 /*** ConvexHull::strict_contains_point
 228  * In order to test whether a point is strictly inside (not on the boundary) a convex hull we can travel once around the outside making
 229  * sure that each triangle made from an edge and the point has positive area. */
 230 bool
 231 ConvexHull::strict_contains_point(Point p) {
 232     return find_strict_left(p) == -1;
 233 }
 234
 235 /*** ConvexHull::add_point
 236  * to add a point we need to find whether the new point extends the boundary, and if so, what it
 237  * obscures.  Tarjan?  Jarvis?*/
 238 void
 239 ConvexHull::merge(Point p) {
 240     std::vector<Point> out;
 241
 242     int l = boundary.size();
 243
 244     if(l < 2) {
 245         boundary.push_back(p);
 246         return;
 247     }
 248
 249     bool pushed = false;
 250
 251     bool pre = is_strict_left(p, -1);
 252     for(int i = 0; i < l; i++) {
 253         bool cur = is_strict_left(p, i);
 254         if(pre) {
 255             if(cur) {
 256                 if(!pushed) {
 257                     out.push_back(p);
 258                     pushed = true;
 259                 }
 260                 continue;
 261             }
 262             else if(!pushed) {
 263                 out.push_back(p);
 264                 pushed = true;
 265             }
 266         }
 267         out.push_back(boundary[i]);
 268         pre = cur;
 269     }
 270
 271     boundary = out;
 272 }
 273 //OPT: quickly find an obscured point and find the bounds by extending from there.  then push all points not within the bounds in order.
 274   //OPT: use binary searches to find the actual starts/ends, use known rights as boundaries.  may require cooperation of find_left algo.
 275
 276 /*** ConvexHull::is_clockwise
 277  * We require that successive pairs of edges always turn right.
 278  * proposed algorithm: walk successive edges and require triangle area is positive.
 279  */
 280 bool
 281 ConvexHull::is_clockwise() const {
 282     if(is_degenerate())
 283         return true;
 284     Point first = boundary[0];
 285     Point second = boundary[1];
 286     for(std::vector<Point>::const_iterator it(boundary.begin()+2), e(boundary.end());
 287         it != e;) {
 288         if(SignedTriangleArea(first, second, *it) > 0)
 289             return false;
 290         first = second;
 291         second = *it;
 292         ++it;
 293     }
 294     return true;
 295 }
 296
 297 /*** ConvexHull::top_point_first
 298  * We require that the first point in the convex hull has the least y coord, and that off all such points on the hull, it has the least x coord.
 299  * proposed algorithm: track lexicographic minimum while walking the list.
 300  */
 301 bool
 302 ConvexHull::top_point_first() const {
 303     std::vector<Point>::const_iterator pivot = boundary.begin();
 304     for(std::vector<Point>::const_iterator it(boundary.begin()+1),
 305             e(boundary.end());
 306         it != e; it++) {
 307         if((*it)[1] < (*pivot)[1])
 308             pivot = it;
 309         else if(((*it)[1] == (*pivot)[1]) &&
 310                 ((*it)[0] < (*pivot)[0]))
 311             pivot = it;
 312     }
 313     return pivot == boundary.begin();
 314 }
 315 //OPT: since the Y values are orderly there should be something like a binary search to do this.
 316
 317 /*** ConvexHull::no_colinear_points
 318  * We require that no three vertices are colinear.
 319 proposed algorithm:  We must be very careful about rounding here.
 320 */
 321 bool
 322 ConvexHull::no_colinear_points() const {
 323     // XXX: implement me!
 324     THROW_NOTIMPLEMENTED();
 325 }
 326
 327 bool
 328 ConvexHull::meets_invariants() const {
 329     return is_clockwise() && top_point_first() && no_colinear_points();
 330 }
 331
 332 /*** ConvexHull::is_degenerate
 333  * We allow three degenerate cases: empty, 1 point and 2 points.  In many cases these should be handled explicitly.
 334  */
 335 bool
 336 ConvexHull::is_degenerate() const {
 337     return boundary.size() < 3;
 338 }
 339
 340
 341 /* Here we really need a rotating calipers implementation.  This implementation is slow and incorrect.
 342    This incorrectness is a problem because it throws off the algorithms.  Perhaps I will come up with
 343    something better tomorrow.  The incorrectness is in the order of the bridges - they must be in the
 344    order of traversal around.  Since the a->b and b->a bridges are seperated, they don't need to be merge
 345    order, just the order of the traversal of the host hull.  Currently some situations make a n->0 bridge
 346    first.*/
 347 pair< map<int, int>, map<int, int> >
 348 bridges(ConvexHull a, ConvexHull b) {
 349     map<int, int> abridges;
 350     map<int, int> bbridges;
 351
 352     for(unsigned ia = 0; ia < a.boundary.size(); ia++) {
 353         for(unsigned ib = 0; ib < b.boundary.size(); ib++) {
 354             Point d = b[ib] - a[ia];
 355             Geom::Coord e = cross(d, a[ia - 1] - a[ia]), f = cross(d, a[ia + 1] - a[ia]);
 356             Geom::Coord g = cross(d, b[ib - 1] - a[ia]), h = cross(d, b[ib + 1] - a[ia]);
 357             if     (e > 0 && f > 0 && g > 0 && h > 0) abridges[ia] = ib;
 358             else if(e < 0 && f < 0 && g < 0 && h < 0) bbridges[ib] = ia;
 359         }
 360     }
 361
 362     return make_pair(abridges, bbridges);
 363 }
 364
 365 std::vector<Point> bridge_points(ConvexHull a, ConvexHull b) {
 366     vector<Point> ret;
 367     pair< map<int, int>, map<int, int> > indices = bridges(a, b);
 368     for(map<int, int>::iterator it = indices.first.begin(); it != indices.first.end(); it++) {
 369       ret.push_back(a[it->first]);
 370       ret.push_back(b[it->second]);
 371     }
 372     for(map<int, int>::iterator it = indices.second.begin(); it != indices.second.end(); it++) {
 373       ret.push_back(b[it->first]);
 374       ret.push_back(a[it->second]);
 375     }
 376     return ret;
 377 }
 378
 379 unsigned find_bottom_right(ConvexHull const &a) {
 380     unsigned it = 1;
 381     while(it < a.boundary.size() &&
 382           a.boundary[it][Y] > a.boundary[it-1][Y])
 383         it++;
 384     return it-1;
 385 }
 386
 387 /*** ConvexHull sweepline_intersection(ConvexHull a, ConvexHull b);
 388  * find the intersection between two convex hulls.  The intersection is also a convex hull.
 389  * (Proof: take any two points both in a and in b.  Any point between them is in a by convexity,
 390  * and in b by convexity, thus in both.  Need to prove still finite bounds.)
 391  * This algorithm works by sweeping a line down both convex hulls in parallel, working out the left and right edges of the new hull.
 392  */
 393 ConvexHull sweepline_intersection(ConvexHull const &a, ConvexHull const &b) {
 394     ConvexHull ret;
 395
 396     unsigned al = 0;
 397     unsigned bl = 0;
 398
 399     while(al+1 < a.boundary.size() &&
 400           (a.boundary[al+1][Y] > b.boundary[bl][Y])) {
 401         al++;
 402     }
 403     while(bl+1 < b.boundary.size() &&
 404           (b.boundary[bl+1][Y] > a.boundary[al][Y])) {
 405         bl++;
 406     }
 407     // al and bl now point to the top of the first pair of edges that overlap in y value
 408     //double sweep_y = std::min(a.boundary[al][Y],
 409     //                          b.boundary[bl][Y]);
 410     return ret;
 411 }
 412
 413 /*** ConvexHull intersection(ConvexHull a, ConvexHull b);
 414  * find the intersection between two convex hulls.  The intersection is also a convex hull.
 415  * (Proof: take any two points both in a and in b.  Any point between them is in a by convexity,
 416  * and in b by convexity, thus in both.  Need to prove still finite bounds.)
 417  */
 418 ConvexHull intersection(ConvexHull /*a*/, ConvexHull /*b*/) {
 419     ConvexHull ret;
 420     /*
 421     int ai = 0, bi = 0;
 422     int aj = a.boundary.size() - 1;
 423     int bj = b.boundary.size() - 1;
 424     */
 425     /*while (true) {
 426         if(a[ai]
 427     }*/
 428     return ret;
 429 }
 430
 431 /*** ConvexHull merge(ConvexHull a, ConvexHull b);
 432  * find the smallest convex hull that surrounds a and b.
 433  */
 434 ConvexHull merge(ConvexHull a, ConvexHull b) {
 435     ConvexHull ret;
 436
 437     pair< map<int, int>, map<int, int> > bpair = bridges(a, b);
 438     map<int, int> ab = bpair.first;
 439     map<int, int> bb = bpair.second;
 440
 441     ab[-1] = 0;
 442     bb[-1] = 0;
 443
 444     int i = -1; // XXX: i is int but refers to vector indices
 445
 446     if(a.boundary[0][1] > b.boundary[0][1]) goto start_b;
 447     while(true) {
 448         for(; ab.count(i) == 0; i++) {
 449             ret.boundary.push_back(a[i]);
 450             if(i >= (int)a.boundary.size()) return ret;
 451         }
 452         if(ab[i] == 0 && i != -1) break;
 453         i = ab[i];
 454         start_b:
 455
 456         for(; bb.count(i) == 0; i++) {
 457             ret.boundary.push_back(b[i]);
 458             if(i >= (int)b.boundary.size()) return ret;
 459         }
 460         if(bb[i] == 0 && i != -1) break;
 461         i = bb[i];
 462     }
 463     return ret;
 464 }
 465
 466 ConvexHull graham_merge(ConvexHull a, ConvexHull b) {
 467     ConvexHull result;
 468
 469     // we can avoid the find pivot step because of top_point_first
 470     if(b.boundary[0] <= a.boundary[0])
 471         std::swap(a, b);
 472
 473     result.boundary = a.boundary;
 474     result.boundary.insert(result.boundary.end(),
 475                            b.boundary.begin(), b.boundary.end());
 476
 477 /** if we modified graham scan to work top to bottom as proposed in lect754.pdf we could replace the
 478  angle sort with a simple merge sort type algorithm. furthermore, we could do the graham scan
 479  online, avoiding a bunch of memory copies.  That would probably be linear. -- njh*/
 480     result.angle_sort();
 481     result.graham_scan();
 482
 483     return result;
 484 }
 485 //TODO: reinstate
 486 /*ConvexCover::ConvexCover(Path const &sp) : path(&sp) {
 487     cc.reserve(sp.size());
 488     for(Geom::Path::const_iterator it(sp.begin()), end(sp.end()); it != end; ++it) {
 489         cc.push_back(ConvexHull((*it).begin(), (*it).end()));
 490     }
 491 }*/
 492
 493 double ConvexHull::centroid_and_area(Geom::Point& centroid) const {
 494     const unsigned n = boundary.size();
 495     if (n < 2)
 496         return 0;
 497     if(n < 3) {
 498         centroid = (boundary[0] + boundary[1])/2;
 499         return 0;
 500     }
 501     Geom::Point centroid_tmp(0,0);
 502     double atmp = 0;
 503     for (unsigned i = n-1, j = 0; j < n; i = j, j++) {
 504         const double ai = -cross(boundary[j], boundary[i]);
 505         atmp += ai;
 506         centroid_tmp += (boundary[j] + boundary[i])*ai; // first moment.
 507     }
 508     if (atmp != 0) {
 509         centroid = centroid_tmp / (3 * atmp);
 510     }
 511     return atmp / 2;
 512 }
 513
 514 // TODO: This can be made lg(n) using golden section/fibonacci search three starting points, say 0,
 515 // n/2, n-1 construct a new point, say (n/2 + n)/2 throw away the furthest boundary point iterate
 516 // until interval is a single value
 517 Point const * ConvexHull::furthest(Point direction) const {
 518     Point const * p = &boundary[0];
 519     double d = dot(*p, direction);
 520     for(unsigned i = 1; i < boundary.size(); i++) {
 521         double dd = dot(boundary[i], direction);
 522         if(d < dd) {
 523             p = &boundary[i];
 524             d = dd;
 525         }
 526     }
 527     return p;
 528 }
 529
 530
 531 // returns (a, (b,c)), three points which define the narrowest diameter of the hull as the pair of
 532 // lines going through b,c, and through a, parallel to b,c TODO: This can be made linear time by
 533 // moving point tc incrementally from the previous value (it can only move in one direction).  It
 534 // is currently n*O(furthest)
 535 double ConvexHull::narrowest_diameter(Point &a, Point &b, Point &c) {
 536     Point tb = boundary.back();
 537     double d = INFINITY;
 538     for(unsigned i = 0; i < boundary.size(); i++) {
 539         Point tc = boundary[i];
 540         Point n = -rot90(tb-tc);
 541         Point ta = *furthest(n);
 542         double td = dot(n, ta-tb)/dot(n,n);
 543         if(td < d) {
 544             a = ta;
 545             b = tb;
 546             c = tc;
 547             d = td;
 548         }
 549         tb = tc;
 550     }
 551     return d;
 552 }
 553
 554 };
 555
 556 /*
 557   Local Variables:
 558   mode:c++
 559   c-file-style:"stroustrup"
 560   c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
 561   indent-tabs-mode:nil
 562   fill-column:99
 563   End:
 564 */
 565 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :