4
* Copyright 2006 Nathan Hurst <njh@mail.csse.monash.edu.au>
5
* Copyright 2006 Michael G. Sloan <mgsloan@gmail.com>
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.
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
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/
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.
32
#include "convex-cover.h"
36
+ modify graham scan to work top to bottom, rather than around angles
38
+ minimum distance between convex hulls
39
+ maximum distance between convex hulls
41
+ check all degenerate cases carefully
42
+ check all algorithms meet all invariants
43
+ generalise rotating caliper algorithm (iterator/circulator?)
52
/*** SignedTriangleArea
53
* returns the area of the triangle defined by p0, p1, p2. A clockwise triangle has positive area.
56
SignedTriangleArea(Point p0, Point p1, Point p2) {
57
return cross((p1 - p0), (p2 - p0));
63
angle_cmp(Point o) : o(o) {}
66
operator()(Point a, Point b) {
73
if((da[1] == 0) && (db[1] == 0))
76
return true; // infinite tangent
78
return false; // infinite tangent
84
//assert((ata > atb) == (aa < ab));
85
double aa = atan2(da);
86
double ab = atan2(db);
91
return L2sq(da) < L2sq(db);
97
ConvexHull::find_pivot() {
100
for(unsigned i = 1; i < boundary.size(); i++)
101
if(boundary[i] <= boundary[pivot])
104
std::swap(boundary[0], boundary[pivot]);
108
ConvexHull::angle_sort() {
109
// sort points by angle (resolve ties in favor of point farther from P);
110
// we leave the first one in place as our pivot
111
std::sort(boundary.begin()+1, boundary.end(), angle_cmp(boundary[0]));
115
ConvexHull::graham_scan() {
117
for(unsigned i = 2; i < boundary.size(); i++) {
118
double o = SignedTriangleArea(boundary[stac-2],
121
if(o == 0) { // colinear - dangerous...
123
} else if(o < 0) { // anticlockwise
124
} else { // remove concavity
125
while(o >= 0 && stac > 2) {
127
o = SignedTriangleArea(boundary[stac-2],
132
boundary[stac++] = boundary[i];
134
boundary.resize(stac);
138
ConvexHull::graham() {
144
//Mathematically incorrect mod, but more useful.
145
int mod(int i, int l) {
149
//OPT: usages can often be replaced by conditions
151
/*** ConvexHull::left
152
* Tests if a point is left (outside) of a particular segment, n. */
154
ConvexHull::is_left(Point p, int n) {
155
return SignedTriangleArea((*this)[n], (*this)[n+1], p) > 0;
158
/*** ConvexHull::find_positive
159
* May return any number n where the segment n -> n + 1 (possibly looped around) in the hull such
160
* that the point is on the wrong side to be within the hull. Returns -1 if it is within the hull.*/
162
ConvexHull::find_left(Point p) {
163
int l = boundary.size(); //Who knows if C++ is smart enough to optimize this?
164
for(int i = 0; i < l; i++) {
165
if(is_left(p, i)) return i;
169
//OPT: do a spread iteration - quasi-random with no repeats and full coverage.
171
/*** ConvexHull::contains_point
172
* In order to test whether a point is inside a convex hull we can travel once around the outside making
173
* sure that each triangle made from an edge and the point has positive area. */
175
ConvexHull::contains_point(Point p) {
176
return find_left(p) == -1;
179
/*** ConvexHull::add_point
180
* to add a point we need to find whether the new point extends the boundary, and if so, what it
181
* obscures. Tarjan? Jarvis?*/
183
ConvexHull::merge(Point p) {
184
std::vector<Point> out;
186
int l = boundary.size();
189
boundary.push_back(p);
195
bool pre = is_left(p, -1);
196
for(int i = 0; i < l; i++) {
197
bool cur = is_left(p, i);
211
out.push_back(boundary[i]);
217
//OPT: quickly find an obscured point and find the bounds by extending from there. then push all points not within the bounds in order.
218
//OPT: use binary searches to find the actual starts/ends, use known rights as boundaries. may require cooperation of find_left algo.
220
/*** ConvexHull::is_clockwise
221
* We require that successive pairs of edges always turn right.
222
* proposed algorithm: walk successive edges and require triangle area is positive.
225
ConvexHull::is_clockwise() const {
228
Point first = boundary[0];
229
Point second = boundary[1];
230
for(std::vector<Point>::const_iterator it(boundary.begin()+2), e(boundary.end());
232
if(SignedTriangleArea(first, second, *it) > 0)
241
/*** ConvexHull::top_point_first
242
* 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.
243
* proposed algorithm: track lexicographic minimum while walking the list.
246
ConvexHull::top_point_first() const {
247
std::vector<Point>::const_iterator pivot = boundary.begin();
248
for(std::vector<Point>::const_iterator it(boundary.begin()+1),
251
if((*it)[1] < (*pivot)[1])
253
else if(((*it)[1] == (*pivot)[1]) &&
254
((*it)[0] < (*pivot)[0]))
257
return pivot == boundary.begin();
259
//OPT: since the Y values are orderly there should be something like a binary search to do this.
261
/*** ConvexHull::no_colinear_points
262
* We require that no three vertices are colinear.
263
proposed algorithm: We must be very careful about rounding here.
266
ConvexHull::no_colinear_points() const {
271
ConvexHull::meets_invariants() const {
272
return is_clockwise() && top_point_first() && no_colinear_points();
275
/*** ConvexHull::is_degenerate
276
* We allow three degenerate cases: empty, 1 point and 2 points. In many cases these should be handled explicitly.
279
ConvexHull::is_degenerate() const {
280
return boundary.size() < 3;
284
/* Here we really need a rotating calipers implementation. This implementation is slow and incorrect.
285
This incorrectness is a problem because it throws off the algorithms. Perhaps I will come up with
286
something better tomorrow. The incorrectness is in the order of the bridges - they must be in the
287
order of traversal around. Since the a->b and b->a bridges are seperated, they don't need to be merge
288
order, just the order of the traversal of the host hull. Currently some situations make a n->0 bridge
290
pair< map<int, int>, map<int, int> >
291
bridges(ConvexHull a, ConvexHull b) {
292
map<int, int> abridges;
293
map<int, int> bbridges;
295
for(unsigned ia = 0; ia < a.boundary.size(); ia++) {
296
for(unsigned ib = 0; ib < b.boundary.size(); ib++) {
297
Point d = b[ib] - a[ia];
298
Geom::Coord e = cross(d, a[ia - 1] - a[ia]), f = cross(d, a[ia + 1] - a[ia]);
299
Geom::Coord g = cross(d, b[ib - 1] - a[ia]), h = cross(d, b[ib + 1] - a[ia]);
300
if (e > 0 && f > 0 && g > 0 && h > 0) abridges[ia] = ib;
301
else if(e < 0 && f < 0 && g < 0 && h < 0) bbridges[ib] = ia;
305
return make_pair(abridges, bbridges);
308
std::vector<Point> bridge_points(ConvexHull a, ConvexHull b) {
310
pair< map<int, int>, map<int, int> > indices = bridges(a, b);
311
for(map<int, int>::iterator it = indices.first.begin(); it != indices.first.end(); it++) {
312
ret.push_back(a[it->first]);
313
ret.push_back(b[it->second]);
315
for(map<int, int>::iterator it = indices.second.begin(); it != indices.second.end(); it++) {
316
ret.push_back(b[it->first]);
317
ret.push_back(a[it->second]);
322
unsigned find_bottom_right(ConvexHull const &a) {
324
while(it < a.boundary.size() &&
325
a.boundary[it][Y] > a.boundary[it-1][Y])
330
/*** ConvexHull sweepline_intersection(ConvexHull a, ConvexHull b);
331
* find the intersection between two convex hulls. The intersection is also a convex hull.
332
* (Proof: take any two points both in a and in b. Any point between them is in a by convexity,
333
* and in b by convexity, thus in both. Need to prove still finite bounds.)
334
* This algorithm works by sweeping a line down both convex hulls in parallel, working out the left and right edges of the new hull.
336
ConvexHull sweepline_intersection(ConvexHull const &a, ConvexHull const &b) {
342
while(al+1 < a.boundary.size() &&
343
(a.boundary[al+1][Y] > b.boundary[bl][Y])) {
346
while(bl+1 < b.boundary.size() &&
347
(b.boundary[bl+1][Y] > a.boundary[al][Y])) {
351
// al and bl now point to the top of the first pair of edges that overlap in y value
352
// double sweep_y = std::min(a.boundary[al][Y],
353
// b.boundary[bl][Y]);
356
/*** ConvexHull intersection(ConvexHull a, ConvexHull b);
357
* find the intersection between two convex hulls. The intersection is also a convex hull.
358
* (Proof: take any two points both in a and in b. Any point between them is in a by convexity,
359
* and in b by convexity, thus in both. Need to prove still finite bounds.)
361
ConvexHull intersection(ConvexHull a, ConvexHull b) {
363
// int ai = 0, bi = 0;
364
// unsigned aj = a.boundary.size() - 1;
365
// unsigned bj = b.boundary.size() - 1;
373
/*** ConvexHull merge(ConvexHull a, ConvexHull b);
374
* find the smallest convex hull that surrounds a and b.
376
ConvexHull merge(ConvexHull a, ConvexHull b) {
379
pair< map<int, int>, map<int, int> > bpair = bridges(a, b);
380
map<int, int> ab = bpair.first;
381
map<int, int> bb = bpair.second;
388
if(a.boundary[0][1] > b.boundary[0][1]) goto start_b;
390
for(; ab.count(i) == 0; i++) {
391
ret.boundary.push_back(a[i]);
392
if(i >= static_cast<int>(a.boundary.size())) return ret;
394
if(ab[i] == 0 && i != -1) break;
398
for(; bb.count(i) == 0; i++) {
399
ret.boundary.push_back(b[i]);
400
if(i >= static_cast<int>(b.boundary.size())) return ret;
402
if(bb[i] == 0 && i != -1) break;
408
ConvexHull graham_merge(ConvexHull a, ConvexHull b) {
411
// we can avoid the find pivot step because of top_point_first
412
if(b.boundary[0] <= a.boundary[0])
415
result.boundary = a.boundary;
416
result.boundary.insert(result.boundary.end(),
417
b.boundary.begin(), b.boundary.end());
419
/** if we modified graham scan to work top to bottom as proposed in lect754.pdf we could replace the
420
angle sort with a simple merge sort type algorithm. furthermore, we could do the graham scan
421
online, avoiding a bunch of memory copies. That would probably be linear. -- njh*/
423
result.graham_scan();
428
/*ConvexCover::ConvexCover(Path const &sp) : path(&sp) {
429
cc.reserve(sp.size());
430
for(Geom::Path::const_iterator it(sp.begin()), end(sp.end()); it != end; ++it) {
431
cc.push_back(ConvexHull((*it).begin(), (*it).end()));
441
c-file-style:"stroustrup"
442
c-file-offsets:((innamespace . 0)(substatement-open . 0))
447
vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4 :
added
removed
scribus/plugins/tools/2geomtools/lib2geom/convex-cover.cpp
