Subversion Repositories Scribus

Rev

Rev 15347 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
12065 fschmid 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
18122 mrdocs 17
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
12065 fschmid 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
 
15347 fschmid 32
#include "convex-cover.h"
12065 fschmid 33
#include <algorithm>
34
#include <map>
35
/** Todo:
36
    + modify graham scan to work top to bottom, rather than around angles
37
    + intersection
38
    + minimum distance between convex hulls
39
    + maximum distance between convex hulls
40
    + hausdorf metric?
41
    + check all degenerate cases carefully
42
    + check all algorithms meet all invariants
43
    + generalise rotating caliper algorithm (iterator/circulator?)
44
*/
45
 
46
using std::vector;
47
using std::map;
48
using std::pair;
49
 
50
namespace Geom{
51
 
52
/*** SignedTriangleArea
53
 * returns the area of the triangle defined by p0, p1, p2.  A clockwise triangle has positive area.
54
 */
55
double
56
SignedTriangleArea(Point p0, Point p1, Point p2) {
57
    return cross((p1 - p0), (p2 - p0));
58
}
59
 
60
class angle_cmp{
61
public:
62
    Point o;
63
    angle_cmp(Point o) : o(o) {}
15347 fschmid 64
 
12065 fschmid 65
    bool
66
    operator()(Point a, Point b) {
67
        Point da = a - o;
68
        Point db = b - o;
15347 fschmid 69
 
12065 fschmid 70
#if 1
71
        double aa = da[0];
72
        double ab = db[0];
73
        if((da[1] == 0) && (db[1] == 0))
74
            return da[0] < db[0];
75
        if(da[1] == 0)
76
            return true; // infinite tangent
77
        if(db[1] == 0)
78
            return false; // infinite tangent
79
        aa = da[0] / da[1];
80
        ab = db[0] / db[1];
81
        if(aa > ab)
82
            return true;
83
#else
84
        //assert((ata > atb) == (aa < ab));
85
        double aa = atan2(da);
86
        double ab = atan2(db);
87
        if(aa < ab)
88
            return true;
89
#endif
90
        if(aa == ab)
91
            return L2sq(da) < L2sq(db);
92
        return false;
93
    }
94
};
95
 
96
void
97
ConvexHull::find_pivot() {
98
    // Find pivot P;
99
    unsigned pivot = 0;
15347 fschmid 100
    for(unsigned i = 1; i < boundary.size(); i++)
12065 fschmid 101
        if(boundary[i] <= boundary[pivot])
102
            pivot = i;
15347 fschmid 103
 
12065 fschmid 104
    std::swap(boundary[0], boundary[pivot]);
105
}
106
 
107
void
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]));
112
}
113
 
114
void
115
ConvexHull::graham_scan() {
116
    unsigned stac = 2;
15347 fschmid 117
    for(unsigned i = 2; i < boundary.size(); i++) {
118
        double o = SignedTriangleArea(boundary[stac-2],
119
                                      boundary[stac-1],
12065 fschmid 120
                                      boundary[i]);
121
        if(o == 0) { // colinear - dangerous...
122
            stac--;
123
        } else if(o < 0) { // anticlockwise
124
        } else { // remove concavity
125
            while(o >= 0 && stac > 2) {
126
                stac--;
15347 fschmid 127
                o = SignedTriangleArea(boundary[stac-2],
128
                                       boundary[stac-1],
12065 fschmid 129
                                       boundary[i]);
130
            }
131
        }
132
        boundary[stac++] = boundary[i];
133
    }
134
    boundary.resize(stac);
135
}
136
 
137
void
138
ConvexHull::graham() {
139
    find_pivot();
140
    angle_sort();
141
    graham_scan();
142
}
143
 
144
//Mathematically incorrect mod, but more useful.
145
int mod(int i, int l) {
15347 fschmid 146
    return i >= 0 ?
12065 fschmid 147
           i % l : (i % l) + l;
148
}
149
//OPT: usages can often be replaced by conditions
150
 
151
/*** ConvexHull::left
152
 * Tests if a point is left (outside) of a particular segment, n. */
153
bool
154
ConvexHull::is_left(Point p, int n) {
155
    return SignedTriangleArea((*this)[n], (*this)[n+1], p) > 0;
156
}
157
 
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.*/
161
int
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;
166
    }
167
    return -1;
168
}
15347 fschmid 169
//OPT: do a spread iteration - quasi-random with no repeats and full coverage. 
12065 fschmid 170
 
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. */
174
bool
175
ConvexHull::contains_point(Point p) {
176
    return find_left(p) == -1;
177
}
178
 
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?*/
182
void
183
ConvexHull::merge(Point p) {
184
    std::vector<Point> out;
185
 
186
    int l = boundary.size();
187
 
188
    if(l < 2) {
189
        boundary.push_back(p);
190
        return;
191
    }
192
 
193
    bool pushed = false;
194
 
15347 fschmid 195
    bool pre = is_left(p, -1);
12065 fschmid 196
    for(int i = 0; i < l; i++) {
15347 fschmid 197
        bool cur = is_left(p, i);
12065 fschmid 198
        if(pre) {
199
            if(cur) {
200
                if(!pushed) {
201
                    out.push_back(p);
202
                    pushed = true;
203
                }
204
                continue;
205
            }
206
            else if(!pushed) {
207
                out.push_back(p);
208
                pushed = true;
209
            }
210
        }
211
        out.push_back(boundary[i]);
212
        pre = cur;
213
    }
15347 fschmid 214
 
12065 fschmid 215
    boundary = out;
216
}
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.
219
 
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.
223
 */
224
bool
225
ConvexHull::is_clockwise() const {
226
    if(is_degenerate())
227
        return true;
228
    Point first = boundary[0];
229
    Point second = boundary[1];
230
    for(std::vector<Point>::const_iterator it(boundary.begin()+2), e(boundary.end());
231
        it != e;) {
232
        if(SignedTriangleArea(first, second, *it) > 0)
233
            return false;
234
        first = second;
235
        second = *it;
236
        ++it;
237
    }
238
    return true;
239
}
240
 
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.
244
 */
245
bool
246
ConvexHull::top_point_first() const {
247
    std::vector<Point>::const_iterator pivot = boundary.begin();
15347 fschmid 248
    for(std::vector<Point>::const_iterator it(boundary.begin()+1),
12065 fschmid 249
            e(boundary.end());
250
        it != e; it++) {
251
        if((*it)[1] < (*pivot)[1])
252
            pivot = it;
15347 fschmid 253
        else if(((*it)[1] == (*pivot)[1]) &&
12065 fschmid 254
                ((*it)[0] < (*pivot)[0]))
255
            pivot = it;
256
    }
257
    return pivot == boundary.begin();
258
}
259
//OPT: since the Y values are orderly there should be something like a binary search to do this.
260
 
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.
264
*/
265
bool
266
ConvexHull::no_colinear_points() const {
15347 fschmid 267
        return true;
12065 fschmid 268
}
269
 
270
bool
271
ConvexHull::meets_invariants() const {
272
    return is_clockwise() && top_point_first() && no_colinear_points();
273
}
274
 
275
/*** ConvexHull::is_degenerate
276
 * We allow three degenerate cases: empty, 1 point and 2 points.  In many cases these should be handled explicitly.
277
 */
278
bool
279
ConvexHull::is_degenerate() const {
280
    return boundary.size() < 3;
281
}
282
 
283
 
15347 fschmid 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
289
   first.*/
290
pair< map<int, int>, map<int, int> >
291
bridges(ConvexHull a, ConvexHull b) {
292
    map<int, int> abridges;
293
    map<int, int> bbridges;
12065 fschmid 294
 
15347 fschmid 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;
302
        }
12065 fschmid 303
    }
15347 fschmid 304
 
305
    return make_pair(abridges, bbridges);
12065 fschmid 306
}
307
 
15347 fschmid 308
std::vector<Point> bridge_points(ConvexHull a, ConvexHull b) {
309
    vector<Point> ret;
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]);
12065 fschmid 314
    }
15347 fschmid 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]);
318
    }
12065 fschmid 319
    return ret;
320
}
321
 
322
unsigned find_bottom_right(ConvexHull const &a) {
323
    unsigned it = 1;
15347 fschmid 324
    while(it < a.boundary.size() &&
12065 fschmid 325
          a.boundary[it][Y] > a.boundary[it-1][Y])
326
        it++;
327
    return it-1;
328
}
329
 
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.
335
 */
336
ConvexHull sweepline_intersection(ConvexHull const &a, ConvexHull const &b) {
337
    ConvexHull ret;
15347 fschmid 338
 
12065 fschmid 339
    unsigned al = 0;
340
    unsigned bl = 0;
15347 fschmid 341
 
12065 fschmid 342
    while(al+1 < a.boundary.size() &&
343
          (a.boundary[al+1][Y] > b.boundary[bl][Y])) {
344
        al++;
345
    }
346
    while(bl+1 < b.boundary.size() &&
347
          (b.boundary[bl+1][Y] > a.boundary[al][Y])) {
348
        bl++;
349
    }
15347 fschmid 350
    return ret;
15344 fschmid 351
    // al and bl now point to the top of the first pair of edges that overlap in y value
15347 fschmid 352
//    double sweep_y = std::min(a.boundary[al][Y],
353
//                              b.boundary[bl][Y]);
12065 fschmid 354
}
355
 
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.)
360
 */
15347 fschmid 361
ConvexHull intersection(ConvexHull a, ConvexHull b) {
12065 fschmid 362
    ConvexHull ret;
15347 fschmid 363
//    int ai = 0, bi = 0;
364
//    unsigned aj = a.boundary.size() - 1;
365
//    unsigned bj = b.boundary.size() - 1;
366
 
12065 fschmid 367
    /*while (true) {
368
        if(a[ai]
369
    }*/
370
    return ret;
371
}
372
 
373
/*** ConvexHull merge(ConvexHull a, ConvexHull b);
374
 * find the smallest convex hull that surrounds a and b.
375
 */
376
ConvexHull merge(ConvexHull a, ConvexHull b) {
377
    ConvexHull ret;
378
 
15347 fschmid 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;
12065 fschmid 382
 
383
    ab[-1] = 0;
384
    bb[-1] = 0;
385
 
15347 fschmid 386
    int i = -1;
12065 fschmid 387
 
388
    if(a.boundary[0][1] > b.boundary[0][1]) goto start_b;
389
    while(true) {
390
        for(; ab.count(i) == 0; i++) {
391
            ret.boundary.push_back(a[i]);
15347 fschmid 392
            if(i >= static_cast<int>(a.boundary.size())) return ret;
12065 fschmid 393
        }
394
        if(ab[i] == 0 && i != -1) break;
395
        i = ab[i];
396
        start_b:
15347 fschmid 397
 
12065 fschmid 398
        for(; bb.count(i) == 0; i++) {
399
            ret.boundary.push_back(b[i]);
15347 fschmid 400
            if(i >= static_cast<int>(b.boundary.size())) return ret;
12065 fschmid 401
        }
402
        if(bb[i] == 0 && i != -1) break;
403
        i = bb[i];
404
    }
405
    return ret;
406
}
407
 
408
ConvexHull graham_merge(ConvexHull a, ConvexHull b) {
409
    ConvexHull result;
15347 fschmid 410
 
12065 fschmid 411
    // we can avoid the find pivot step because of top_point_first
412
    if(b.boundary[0] <= a.boundary[0])
413
        std::swap(a, b);
15347 fschmid 414
 
12065 fschmid 415
    result.boundary = a.boundary;
15347 fschmid 416
    result.boundary.insert(result.boundary.end(),
12065 fschmid 417
                           b.boundary.begin(), b.boundary.end());
15347 fschmid 418
 
12065 fschmid 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*/
422
    result.angle_sort();
423
    result.graham_scan();
15347 fschmid 424
 
12065 fschmid 425
    return result;
426
}
427
//TODO: reinstate
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()));
432
    }
433
}*/
434
 
435
 
436
};
437
 
438
/*
439
  Local Variables:
440
  mode:c++
441
  c-file-style:"stroustrup"
15347 fschmid 442
  c-file-offsets:((innamespace . 0)(substatement-open . 0))
12065 fschmid 443
  indent-tabs-mode:nil
15347 fschmid 444
  c-brace-offset:0
12065 fschmid 445
  fill-column:99
446
  End:
15347 fschmid 447
  vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4 :
12065 fschmid 448
*/
15347 fschmid 449
 
450