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 */
32 #include <2geom/convex-cover.h>
33 #include <2geom/exception.h>
34 #include <algorithm>
35 #include <map>
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 */
48 using std::vector;
49 using std::map;
50 using std::pair;
52 namespace Geom{
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 }
62 class angle_cmp{
63 public:
64 Point o;
65 angle_cmp(Point o) : o(o) {}
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;
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 };
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;
127 std::swap(boundary[0], boundary[pivot]);
128 }
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 }
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 }
162 void
163 ConvexHull::graham() {
164 if(is_degenerate()) // nothing to do
165 return;
166 find_pivot();
167 angle_sort();
168 graham_scan();
169 }
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
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 }
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 }
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 }
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 }
217 //OPT: do a spread iteration - quasi-random with no repeats and full coverage.
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 }
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 }
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;
242 int l = boundary.size();
244 if(l < 2) {
245 boundary.push_back(p);
246 return;
247 }
249 bool pushed = false;
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 }
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.
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 }
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.
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 }
327 bool
328 ConvexHull::meets_invariants() const {
329 return is_clockwise() && top_point_first() && no_colinear_points();
330 }
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 }
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;
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 }
362 return make_pair(abridges, bbridges);
363 }
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 }
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 }
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;
396 unsigned al = 0;
397 unsigned bl = 0;
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 }
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 }
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;
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;
441 ab[-1] = 0;
442 bb[-1] = 0;
444 int i = -1; // XXX: i is int but refers to vector indices
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:
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 }
466 ConvexHull graham_merge(ConvexHull a, ConvexHull b) {
467 ConvexHull result;
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);
473 result.boundary = a.boundary;
474 result.boundary.insert(result.boundary.end(),
475 b.boundary.begin(), b.boundary.end());
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();
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 }*/
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 }
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 }
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 }
554 };
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:encoding=utf-8:textwidth=99 :