1 #include <2geom/path-intersection.h>
3 #include <2geom/ord.h>
5 //for path_direction:
6 #include <2geom/sbasis-geometric.h>
7 #include <gsl/gsl_vector.h>
8 #include <gsl/gsl_multiroots.h>
10 namespace Geom {
12 /* This function computes the winding of the path, given a reference point.
13 * Positive values correspond to counter-clockwise in the mathematical coordinate system,
14 * and clockwise in screen coordinates. This particular implementation casts a ray in
15 * the positive x direction. It iterates the path, checking for intersection with the
16 * bounding boxes. If an intersection is found, the initial/final Y value of the curve is
17 * used to derive a delta on the winding value. If the point is within the bounding box,
18 * the curve specific winding function is called.
19 */
20 int winding(Path const &path, Point p) {
21 //start on a segment which is not a horizontal line with y = p[y]
22 Path::const_iterator start;
23 for(Path::const_iterator iter = path.begin(); ; ++iter) {
24 if(iter == path.end_closed()) { return 0; }
25 if(iter->initialPoint()[Y]!=p[Y]) { start = iter; break; }
26 if(iter->finalPoint()[Y]!=p[Y]) { start = iter; break; }
27 if(iter->boundsFast().height()!=0.){ start = iter; break; }
28 }
29 int wind = 0;
30 unsigned cnt = 0;
31 bool starting = true;
32 for (Path::const_iterator iter = start; iter != start || starting
33 ; ++iter, iter = (iter == path.end_closed()) ? path.begin() : iter )
34 {
35 cnt++;
36 if(cnt > path.size()) return wind; //some bug makes this required
37 starting = false;
38 Rect bounds = iter->boundsFast();
39 Coord x = p[X], y = p[Y];
41 if(x > bounds.right() || !bounds[Y].contains(y)) continue; //ray doesn't intersect box
43 Point final = iter->finalPoint();
44 Point initial = iter->initialPoint();
45 Cmp final_to_ray = cmp(final[Y], y);
46 Cmp initial_to_ray = cmp(initial[Y], y);
48 // if y is included, these will have opposite values, giving order.
49 Cmp c = cmp(final_to_ray, initial_to_ray);
50 if(x < bounds.left()) {
51 // ray goes through bbox
52 // winding delta determined by position of endpoints
53 if(final_to_ray != EQUAL_TO) {
54 wind += int(c); // GT = counter-clockwise = 1; LT = clockwise = -1; EQ = not-included = 0
55 //std::cout << int(c) << " ";
56 goto cont;
57 }
58 } else {
59 //inside bbox, use custom per-curve winding thingie
60 int delt = iter->winding(p);
61 wind += delt;
62 //std::cout << "n" << delt << " ";
63 }
64 //Handling the special case of an endpoint on the ray:
65 if(final[Y] == y) {
66 //Traverse segments until it breaks away from y
67 //99.9% of the time this will happen the first go
68 Path::const_iterator next = iter;
69 next++;
70 for(; ; next++) {
71 if(next == path.end_closed()) next = path.begin();
72 Rect bnds = next->boundsFast();
73 //TODO: X considerations
74 if(bnds.height() > 0) {
75 //It has diverged
76 if(bnds.contains(p)) {
77 const double fudge = 0.01;
78 if(cmp(y, next->valueAt(fudge, Y)) == initial_to_ray) {
79 wind += int(c);
80 std::cout << "!!!!!" << int(c) << " ";
81 }
82 iter = next; // No increment, as the rest of the thing hasn't been counted.
83 } else {
84 Coord ny = next->initialPoint()[Y];
85 if(cmp(y, ny) == initial_to_ray) {
86 //Is a continuation through the ray, so counts windingwise
87 wind += int(c);
88 std::cout << "!!!!!" << int(c) << " ";
89 }
90 iter = ++next;
91 }
92 goto cont;
93 }
94 if(next==start) return wind;
95 }
96 //Looks like it looped, which means everything's flat
97 return 0;
98 }
100 cont:(void)0;
101 }
102 return wind;
103 }
105 /* This function should only be applied to simple paths (regions), as otherwise
106 * a boolean winding direction is undefined. It returns true for fill, false for
107 * hole. Defaults to using the sign of area when it reaches funny cases.
108 */
109 bool path_direction(Path const &p) {
110 if(p.empty()) return false;
111 //could probably be more efficient, but this is a quick job
112 double y = p.initialPoint()[Y];
113 double x = p.initialPoint()[X];
114 Cmp res = cmp(p[0].finalPoint()[Y], y);
115 goto doh;
116 for(unsigned i = 1; i <= p.size(); i++) {
117 Cmp final_to_ray = cmp(p[i].finalPoint()[Y], y);
118 Cmp initial_to_ray = cmp(p[i].initialPoint()[Y], y);
119 // if y is included, these will have opposite values, giving order.
120 Cmp c = cmp(final_to_ray, initial_to_ray);
121 if(c != EQUAL_TO) {
122 std::vector<double> rs = p[i].roots(y, Y);
123 for(unsigned j = 0; j < rs.size(); j++) {
124 double nx = p[i].valueAt(rs[j], X);
125 if(nx > x) {
126 x = nx;
127 res = c;
128 }
129 }
130 } else if(final_to_ray == EQUAL_TO) goto doh;
131 }
132 return res < 0;
134 doh:
135 //Otherwise fallback on area
137 Piecewise<D2<SBasis> > pw = p.toPwSb();
138 double area;
139 Point centre;
140 Geom::centroid(pw, centre, area);
141 return area > 0;
142 }
144 //pair intersect code based on njh's pair-intersect
146 // A little sugar for appending a list to another
147 template<typename T>
148 void append(T &a, T const &b) {
149 a.insert(a.end(), b.begin(), b.end());
150 }
152 /* Finds the intersection between the lines defined by A0 & A1, and B0 & B1.
153 * Returns through the last 3 parameters, returning the t-values on the lines
154 * and the cross-product of the deltas (a useful byproduct). The return value
155 * indicates if the time values are within their proper range on the line segments.
156 */
157 bool
158 linear_intersect(Point A0, Point A1, Point B0, Point B1,
159 double &tA, double &tB, double &det) {
160 // kramers rule as cross products
161 Point Ad = A1 - A0,
162 Bd = B1 - B0,
163 d = B0 - A0;
164 det = cross(Ad, Bd);
165 if( 1.0 + det == 1.0 )
166 return false;
167 else
168 {
169 double detinv = 1.0 / det;
170 tA = cross(d, Bd) * detinv;
171 tB = cross(d, Ad) * detinv;
172 return tA >= 0. && tA <= 1. && tB >= 0. && tB <= 1.;
173 }
174 }
177 typedef union dbl_64{
178 long long i64;
179 double d64;
180 };
182 static double EpsilonOf(double value)
183 {
184 dbl_64 s;
185 s.d64 = value;
186 if(s.i64 == 0)
187 {
188 s.i64++;
189 return s.d64 - value;
190 }
191 else if(s.i64-- < 0)
192 return s.d64 - value;
193 else
194 return value - s.d64;
195 }
197 struct rparams {
198 Curve const &A;
199 Curve const &B;
200 };
202 static int
203 intersect_polish_f (const gsl_vector * x, void *params,
204 gsl_vector * f)
205 {
206 const double x0 = gsl_vector_get (x, 0);
207 const double x1 = gsl_vector_get (x, 1);
209 Geom::Point dx = ((struct rparams *) params)->A(x0) -
210 ((struct rparams *) params)->B(x1);
212 gsl_vector_set (f, 0, dx[0]);
213 gsl_vector_set (f, 1, dx[1]);
215 return GSL_SUCCESS;
216 }
218 static void
219 intersect_polish_root (Curve const &A, double &s,
220 Curve const &B, double &t) {
221 int status;
222 size_t iter = 0;
224 const size_t n = 2;
225 struct rparams p = {A, B};
226 gsl_multiroot_function f = {&intersect_polish_f, n, &p};
228 double x_init[2] = {s, t};
229 gsl_vector *x = gsl_vector_alloc (n);
231 gsl_vector_set (x, 0, x_init[0]);
232 gsl_vector_set (x, 1, x_init[1]);
234 const gsl_multiroot_fsolver_type *T = gsl_multiroot_fsolver_hybrids;
235 gsl_multiroot_fsolver *sol = gsl_multiroot_fsolver_alloc (T, 2);
236 gsl_multiroot_fsolver_set (sol, &f, x);
238 do
239 {
240 iter++;
241 status = gsl_multiroot_fsolver_iterate (sol);
243 if (status) /* check if solver is stuck */
244 break;
246 status =
247 gsl_multiroot_test_residual (sol->f, 1e-12);
248 }
249 while (status == GSL_CONTINUE && iter < 1000);
251 s = gsl_vector_get (sol->x, 0);
252 t = gsl_vector_get (sol->x, 1);
254 gsl_multiroot_fsolver_free (sol);
255 gsl_vector_free (x);
256 }
258 /* This uses the local bounds functions of curves to generically intersect two.
259 * It passes in the curves, time intervals, and keeps track of depth, while
260 * returning the results through the Crossings parameter.
261 */
262 void pair_intersect(Curve const & A, double Al, double Ah,
263 Curve const & B, double Bl, double Bh,
264 Crossings &ret, unsigned depth=0) {
265 // std::cout << depth << "(" << Al << ", " << Ah << ")\n";
266 Rect Ar = A.boundsLocal(Interval(Al, Ah));
267 if(Ar.isEmpty()) return;
269 Rect Br = B.boundsLocal(Interval(Bl, Bh));
270 if(Br.isEmpty()) return;
272 if(!Ar.intersects(Br)) return;
274 //Checks the general linearity of the function
275 if((depth > 12)) { // || (A.boundsLocal(Interval(Al, Ah), 1).maxExtent() < 0.1
276 //&& B.boundsLocal(Interval(Bl, Bh), 1).maxExtent() < 0.1)) {
277 double tA, tB, c;
278 if(linear_intersect(A.pointAt(Al), A.pointAt(Ah),
279 B.pointAt(Bl), B.pointAt(Bh),
280 tA, tB, c)) {
281 tA = tA * (Ah - Al) + Al;
282 tB = tB * (Bh - Bl) + Bl;
283 intersect_polish_root(A, tA,
284 B, tB);
285 if(depth % 2)
286 ret.push_back(Crossing(tB, tA, c < 0));
287 else
288 ret.push_back(Crossing(tA, tB, c > 0));
289 return;
290 }
291 }
292 if(depth > 12) return;
293 double mid = (Bl + Bh)/2;
294 pair_intersect(B, Bl, mid,
295 A, Al, Ah,
296 ret, depth+1);
297 pair_intersect(B, mid, Bh,
298 A, Al, Ah,
299 ret, depth+1);
300 }
301 // A simple wrapper around pair_intersect
302 Crossings SimpleCrosser::crossings(Curve const &a, Curve const &b) {
303 Crossings ret;
304 pair_intersect(a, 0, 1, b, 0, 1, ret);
305 return ret;
306 }
308 /* Takes two paths and time ranges on them, with the invariant that the
309 * paths are monotonic on the range. Splits A when the linear intersection
310 * doesn't exist or is inaccurate. Uses the fact that it is monotonic to
311 * do very fast local bounds.
312 */
313 void mono_pair(Path const &A, double Al, double Ah,
314 Path const &B, double Bl, double Bh,
315 Crossings &ret, double /*tol*/, unsigned depth = 0) {
316 if( Al >= Ah || Bl >= Bh) return;
317 std::cout << " " << depth << "[" << Al << ", " << Ah << "]" << "[" << Bl << ", " << Bh << "]";
319 Point A0 = A.pointAt(Al), A1 = A.pointAt(Ah),
320 B0 = B.pointAt(Bl), B1 = B.pointAt(Bh);
321 //inline code that this implies? (without rect/interval construction)
322 if(!Rect(A0, A1).intersects(Rect(B0, B1)) || A0 == A1 || B0 == B1) return;
324 //Checks the general linearity of the function
325 //if((depth > 12) || (A.boundsLocal(Interval(Al, Ah), 1).maxExtent() < 0.1
326 // && B.boundsLocal(Interval(Bl, Bh), 1).maxExtent() < 0.1)) {
327 double tA, tB, c;
328 if(linear_intersect(A0, A1, B0, B1,
329 tA, tB, c)) {
330 tA = tA * (Ah - Al) + Al;
331 tB = tB * (Bh - Bl) + Bl;
332 if(depth % 2)
333 ret.push_back(Crossing(tB, tA, c < 0));
334 else
335 ret.push_back(Crossing(tA, tB, c > 0));
336 return;
337 }
338 //}
339 if(depth > 12) return;
340 double mid = (Bl + Bh)/2;
341 mono_pair(B, Bl, mid,
342 A, Al, Ah,
343 ret, depth+1);
344 mono_pair(B, mid, Bh,
345 A, Al, Ah,
346 ret, depth+1);
347 }
349 // This returns the times when the x or y derivative is 0 in the curve.
350 std::vector<double> curve_mono_splits(Curve const &d) {
351 std::vector<double> rs = d.roots(0, X);
352 append(rs, d.roots(0, Y));
353 std::sort(rs.begin(), rs.end());
354 return rs;
355 }
357 // Convenience function to add a value to each entry in a vector of doubles.
358 std::vector<double> offset_doubles(std::vector<double> const &x, double offs) {
359 std::vector<double> ret;
360 for(unsigned i = 0; i < x.size(); i++) {
361 ret.push_back(x[i] + offs);
362 }
363 return ret;
364 }
366 /* Finds all the monotonic splits for a path. Only includes the split between
367 * curves if they switch derivative directions at that point.
368 */
369 std::vector<double> path_mono_splits(Path const &p) {
370 std::vector<double> ret;
371 if(p.empty()) return ret;
372 ret.push_back(0);
374 Curve* deriv = p[0].derivative();
375 append(ret, curve_mono_splits(*deriv));
376 delete deriv;
378 bool pdx=2, pdy=2; //Previous derivative direction
379 for(unsigned i = 0; i <= p.size(); i++) {
380 deriv = p[i].derivative();
381 std::vector<double> spl = offset_doubles(curve_mono_splits(*deriv), i);
382 delete deriv;
383 bool dx = p[i].initialPoint()[X] > (spl.empty()? p[i].finalPoint()[X] :
384 p.valueAt(spl.front(), X));
385 bool dy = p[i].initialPoint()[Y] > (spl.empty()? p[i].finalPoint()[Y] :
386 p.valueAt(spl.front(), Y));
387 //The direction changed, include the split time
388 if(dx != pdx || dy != pdy) {
389 ret.push_back(i);
390 pdx = dx; pdy = dy;
391 }
392 append(ret, spl);
393 }
394 return ret;
395 }
397 /* Applies path_mono_splits to multiple paths, and returns the results such that
398 * time-set i corresponds to Path i.
399 */
400 std::vector<std::vector<double> > paths_mono_splits(std::vector<Path> const &ps) {
401 std::vector<std::vector<double> > ret;
402 for(unsigned i = 0; i < ps.size(); i++)
403 ret.push_back(path_mono_splits(ps[i]));
404 return ret;
405 }
407 /* Processes the bounds for a list of paths and a list of splits on them, yielding a list of rects for each.
408 * Each entry i corresponds to path i of the input. The number of rects in each entry is guaranteed to be the
409 * number of splits for that path, subtracted by one.
410 */
411 std::vector<std::vector<Rect> > split_bounds(std::vector<Path> const &p, std::vector<std::vector<double> > splits) {
412 std::vector<std::vector<Rect> > ret;
413 for(unsigned i = 0; i < p.size(); i++) {
414 std::vector<Rect> res;
415 for(unsigned j = 1; j < splits[i].size(); j++)
416 res.push_back(Rect(p[i].pointAt(splits[i][j-1]), p[i].pointAt(splits[i][j])));
417 ret.push_back(res);
418 }
419 return ret;
420 }
422 /* This is the main routine of "MonoCrosser", and implements a monotonic strategy on multiple curves.
423 * Finds crossings between two sets of paths, yielding a CrossingSet. [0, a.size()) of the return correspond
424 * to the sorted crossings of a with paths of b. The rest of the return, [a.size(), a.size() + b.size()],
425 * corresponds to the sorted crossings of b with paths of a.
426 *
427 * This function does two sweeps, one on the bounds of each path, and after that cull, one on the curves within.
428 * This leads to a certain amount of code complexity, however, most of that is factored into the above functions
429 */
430 CrossingSet MonoCrosser::crossings(std::vector<Path> const &a, std::vector<Path> const &b) {
431 if(b.empty()) return CrossingSet(a.size(), Crossings());
432 CrossingSet results(a.size() + b.size(), Crossings());
433 if(a.empty()) return results;
435 std::vector<std::vector<double> > splits_a = paths_mono_splits(a), splits_b = paths_mono_splits(b);
436 std::vector<std::vector<Rect> > bounds_a = split_bounds(a, splits_a), bounds_b = split_bounds(b, splits_b);
438 std::vector<Rect> bounds_a_union, bounds_b_union;
439 for(unsigned i = 0; i < bounds_a.size(); i++) bounds_a_union.push_back(union_list(bounds_a[i]));
440 for(unsigned i = 0; i < bounds_b.size(); i++) bounds_b_union.push_back(union_list(bounds_b[i]));
442 std::vector<std::vector<unsigned> > cull = sweep_bounds(bounds_a_union, bounds_b_union);
443 Crossings n;
444 for(unsigned i = 0; i < cull.size(); i++) {
445 for(unsigned jx = 0; jx < cull[i].size(); jx++) {
446 unsigned j = cull[i][jx];
447 unsigned jc = j + a.size();
448 Crossings res;
450 //Sweep of the monotonic portions
451 std::vector<std::vector<unsigned> > cull2 = sweep_bounds(bounds_a[i], bounds_b[j]);
452 for(unsigned k = 0; k < cull2.size(); k++) {
453 for(unsigned lx = 0; lx < cull2[k].size(); lx++) {
454 unsigned l = cull2[k][lx];
455 mono_pair(a[i], splits_a[i][k-1], splits_a[i][k],
456 b[j], splits_b[j][l-1], splits_b[j][l],
457 res, .1);
458 }
459 }
461 for(unsigned k = 0; k < res.size(); k++) { res[k].a = i; res[k].b = jc; }
463 merge_crossings(results[i], res, i);
464 merge_crossings(results[i], res, jc);
465 }
466 }
468 return results;
469 }
471 /* This function is similar codewise to the MonoCrosser, the main difference is that it deals with
472 * only one set of paths and includes self intersection
473 CrossingSet crossings_among(std::vector<Path> const &p) {
474 CrossingSet results(p.size(), Crossings());
475 if(p.empty()) return results;
477 std::vector<std::vector<double> > splits = paths_mono_splits(p);
478 std::vector<std::vector<Rect> > prs = split_bounds(p, splits);
479 std::vector<Rect> rs;
480 for(unsigned i = 0; i < prs.size(); i++) rs.push_back(union_list(prs[i]));
482 std::vector<std::vector<unsigned> > cull = sweep_bounds(rs);
484 //we actually want to do the self-intersections, so add em in:
485 for(unsigned i = 0; i < cull.size(); i++) cull[i].push_back(i);
487 for(unsigned i = 0; i < cull.size(); i++) {
488 for(unsigned jx = 0; jx < cull[i].size(); jx++) {
489 unsigned j = cull[i][jx];
490 Crossings res;
492 //Sweep of the monotonic portions
493 std::vector<std::vector<unsigned> > cull2 = sweep_bounds(prs[i], prs[j]);
494 for(unsigned k = 0; k < cull2.size(); k++) {
495 for(unsigned lx = 0; lx < cull2[k].size(); lx++) {
496 unsigned l = cull2[k][lx];
497 mono_pair(p[i], splits[i][k-1], splits[i][k],
498 p[j], splits[j][l-1], splits[j][l],
499 res, .1);
500 }
501 }
503 for(unsigned k = 0; k < res.size(); k++) { res[k].a = i; res[k].b = j; }
505 merge_crossings(results[i], res, i);
506 merge_crossings(results[j], res, j);
507 }
508 }
510 return results;
511 }
512 */
515 Crossings curve_self_crossings(Curve const &a) {
516 Crossings res;
517 std::vector<double> spl;
518 spl.push_back(0);
519 append(spl, curve_mono_splits(a));
520 spl.push_back(1);
521 for(unsigned i = 1; i < spl.size(); i++)
522 for(unsigned j = i+1; j < spl.size(); j++)
523 pair_intersect(a, spl[i-1], spl[i], a, spl[j-1], spl[j], res);
524 return res;
525 }
527 /*
528 void mono_curve_intersect(Curve const & A, double Al, double Ah,
529 Curve const & B, double Bl, double Bh,
530 Crossings &ret, unsigned depth=0) {
531 // std::cout << depth << "(" << Al << ", " << Ah << ")\n";
532 Point A0 = A.pointAt(Al), A1 = A.pointAt(Ah),
533 B0 = B.pointAt(Bl), B1 = B.pointAt(Bh);
534 //inline code that this implies? (without rect/interval construction)
535 if(!Rect(A0, A1).intersects(Rect(B0, B1)) || A0 == A1 || B0 == B1) return;
537 //Checks the general linearity of the function
538 if((depth > 12) || (A.boundsLocal(Interval(Al, Ah), 1).maxExtent() < 0.1
539 && B.boundsLocal(Interval(Bl, Bh), 1).maxExtent() < 0.1)) {
540 double tA, tB, c;
541 if(linear_intersect(A0, A1, B0, B1, tA, tB, c)) {
542 tA = tA * (Ah - Al) + Al;
543 tB = tB * (Bh - Bl) + Bl;
544 if(depth % 2)
545 ret.push_back(Crossing(tB, tA, c < 0));
546 else
547 ret.push_back(Crossing(tA, tB, c > 0));
548 return;
549 }
550 }
551 if(depth > 12) return;
552 double mid = (Bl + Bh)/2;
553 mono_curve_intersect(B, Bl, mid,
554 A, Al, Ah,
555 ret, depth+1);
556 mono_curve_intersect(B, mid, Bh,
557 A, Al, Ah,
558 ret, depth+1);
559 }
561 std::vector<std::vector<double> > curves_mono_splits(Path const &p) {
562 std::vector<std::vector<double> > ret;
563 for(unsigned i = 0; i <= p.size(); i++) {
564 std::vector<double> spl;
565 spl.push_back(0);
566 append(spl, curve_mono_splits(p[i]));
567 spl.push_back(1);
568 ret.push_back(spl);
569 }
570 }
572 std::vector<std::vector<Rect> > curves_split_bounds(Path const &p, std::vector<std::vector<double> > splits) {
573 std::vector<std::vector<Rect> > ret;
574 for(unsigned i = 0; i < splits.size(); i++) {
575 std::vector<Rect> res;
576 for(unsigned j = 1; j < splits[i].size(); j++)
577 res.push_back(Rect(p.pointAt(splits[i][j-1]+i), p.pointAt(splits[i][j]+i)));
578 ret.push_back(res);
579 }
580 return ret;
581 }
583 Crossings path_self_crossings(Path const &p) {
584 Crossings ret;
585 std::vector<std::vector<unsigned> > cull = sweep_bounds(bounds(p));
586 std::vector<std::vector<double> > spl = curves_mono_splits(p);
587 std::vector<std::vector<Rect> > bnds = curves_split_bounds(p, spl);
588 for(unsigned i = 0; i < cull.size(); i++) {
589 Crossings res;
590 for(unsigned k = 1; k < spl[i].size(); k++)
591 for(unsigned l = k+1; l < spl[i].size(); l++)
592 mono_curve_intersect(p[i], spl[i][k-1], spl[i][k], p[i], spl[i][l-1], spl[i][l], res);
593 offset_crossings(res, i, i);
594 append(ret, res);
595 for(unsigned jx = 0; jx < cull[i].size(); jx++) {
596 unsigned j = cull[i][jx];
597 res.clear();
599 std::vector<std::vector<unsigned> > cull2 = sweep_bounds(bnds[i], bnds[j]);
600 for(unsigned k = 0; k < cull2.size(); k++) {
601 for(unsigned lx = 0; lx < cull2[k].size(); lx++) {
602 unsigned l = cull2[k][lx];
603 mono_curve_intersect(p[i], spl[i][k-1], spl[i][k], p[j], spl[j][l-1], spl[j][l], res);
604 }
605 }
607 //if(fabs(int(i)-j) == 1 || fabs(int(i)-j) == p.size()-1) {
608 Crossings res2;
609 for(unsigned k = 0; k < res.size(); k++) {
610 if(res[k].ta != 0 && res[k].ta != 1 && res[k].tb != 0 && res[k].tb != 1) {
611 res.push_back(res[k]);
612 }
613 }
614 res = res2;
615 //}
616 offset_crossings(res, i, j);
617 append(ret, res);
618 }
619 }
620 return ret;
621 }
622 */
624 Crossings self_crossings(Path const &p) {
625 Crossings ret;
626 std::vector<std::vector<unsigned> > cull = sweep_bounds(bounds(p));
627 for(unsigned i = 0; i < cull.size(); i++) {
628 Crossings res = curve_self_crossings(p[i]);
629 offset_crossings(res, i, i);
630 append(ret, res);
631 for(unsigned jx = 0; jx < cull[i].size(); jx++) {
632 unsigned j = cull[i][jx];
633 res.clear();
634 pair_intersect(p[i], 0, 1, p[j], 0, 1, res);
636 //if(fabs(int(i)-j) == 1 || fabs(int(i)-j) == p.size()-1) {
637 Crossings res2;
638 for(unsigned k = 0; k < res.size(); k++) {
639 if(res[k].ta != 0 && res[k].ta != 1 && res[k].tb != 0 && res[k].tb != 1) {
640 res2.push_back(res[k]);
641 }
642 }
643 res = res2;
644 //}
645 offset_crossings(res, i, j);
646 append(ret, res);
647 }
648 }
649 return ret;
650 }
652 void flip_crossings(Crossings &crs) {
653 for(unsigned i = 0; i < crs.size(); i++)
654 crs[i] = Crossing(crs[i].tb, crs[i].ta, crs[i].b, crs[i].a, !crs[i].dir);
655 }
657 CrossingSet crossings_among(std::vector<Path> const &p) {
658 CrossingSet results(p.size(), Crossings());
659 if(p.empty()) return results;
661 SimpleCrosser cc;
663 std::vector<std::vector<unsigned> > cull = sweep_bounds(bounds(p));
664 for(unsigned i = 0; i < cull.size(); i++) {
665 Crossings res = self_crossings(p[i]);
666 for(unsigned k = 0; k < res.size(); k++) { res[k].a = res[k].b = i; }
667 merge_crossings(results[i], res, i);
668 flip_crossings(res);
669 merge_crossings(results[i], res, i);
670 for(unsigned jx = 0; jx < cull[i].size(); jx++) {
671 unsigned j = cull[i][jx];
673 Crossings res = cc.crossings(p[i], p[j]);
674 for(unsigned k = 0; k < res.size(); k++) { res[k].a = i; res[k].b = j; }
675 merge_crossings(results[i], res, i);
676 merge_crossings(results[j], res, j);
677 }
678 }
679 return results;
680 }
682 }
684 /*
685 Local Variables:
686 mode:c++
687 c-file-style:"stroustrup"
688 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
689 indent-tabs-mode:nil
690 fill-column:99
691 End:
692 */
693 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :