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