1 #include <2geom/path-intersection.h>
3 #include <2geom/ord.h>
5 //for path_direction:
6 #include <2geom/sbasis-geometric.h>
7 #include <2geom/line.h>
8 #ifdef HAVE_GSL
9 #include <gsl/gsl_vector.h>
10 #include <gsl/gsl_multiroots.h>
11 #endif
13 namespace Geom {
15 /**
16 * This function computes the winding of the path, given a reference point.
17 * Positive values correspond to counter-clockwise in the mathematical coordinate system,
18 * and clockwise in screen coordinates. This particular implementation casts a ray in
19 * the positive x direction. It iterates the path, checking for intersection with the
20 * bounding boxes. If an intersection is found, the initial/final Y value of the curve is
21 * used to derive a delta on the winding value. If the point is within the bounding box,
22 * the curve specific winding function is called.
23 */
24 int winding(Path const &path, Point p) {
25 //start on a segment which is not a horizontal line with y = p[y]
26 Path::const_iterator start;
27 for(Path::const_iterator iter = path.begin(); ; ++iter) {
28 if(iter == path.end_closed()) { return 0; }
29 if(iter->initialPoint()[Y]!=p[Y]) { start = iter; break; }
30 if(iter->finalPoint()[Y]!=p[Y]) { start = iter; break; }
31 if(iter->boundsFast()->height()!=0.){ start = iter; break; }
32 }
33 int wind = 0;
34 unsigned cnt = 0;
35 bool starting = true;
36 for (Path::const_iterator iter = start; iter != start || starting
37 ; ++iter, iter = (iter == path.end_closed()) ? path.begin() : iter )
38 {
39 cnt++;
40 if(cnt > path.size()) return wind; //some bug makes this required
41 starting = false;
42 Rect bounds = *(iter->boundsFast());
43 Coord x = p[X], y = p[Y];
45 if(x > bounds.right() || !bounds[Y].contains(y)) continue; //ray doesn't intersect box
47 Point final = iter->finalPoint();
48 Point initial = iter->initialPoint();
49 Cmp final_to_ray = cmp(final[Y], y);
50 Cmp initial_to_ray = cmp(initial[Y], y);
52 // if y is included, these will have opposite values, giving order.
53 Cmp c = cmp(final_to_ray, initial_to_ray);
54 if(x < bounds.left()) {
55 // ray goes through bbox
56 // winding delta determined by position of endpoints
57 if(final_to_ray != EQUAL_TO) {
58 wind += int(c); // GT = counter-clockwise = 1; LT = clockwise = -1; EQ = not-included = 0
59 //std::cout << int(c) << " ";
60 goto cont;
61 }
62 } else {
63 //inside bbox, use custom per-curve winding thingie
64 int delt = iter->winding(p);
65 wind += delt;
66 //std::cout << "n" << delt << " ";
67 }
68 //Handling the special case of an endpoint on the ray:
69 if(final[Y] == y) {
70 //Traverse segments until it breaks away from y
71 //99.9% of the time this will happen the first go
72 Path::const_iterator next = iter;
73 next++;
74 for(; ; next++) {
75 if(next == path.end_closed()) next = path.begin();
76 Rect bnds = *(next->boundsFast());
77 //TODO: X considerations
78 if(bnds.height() > 0) {
79 //It has diverged
80 if(bnds.contains(p)) {
81 const double fudge = 0.01;
82 if(cmp(y, next->valueAt(fudge, Y)) == initial_to_ray) {
83 wind += int(c);
84 //std::cout << "!!!!!" << int(c) << " ";
85 }
86 iter = next; // No increment, as the rest of the thing hasn't been counted.
87 } else {
88 Coord ny = next->initialPoint()[Y];
89 if(cmp(y, ny) == initial_to_ray) {
90 //Is a continuation through the ray, so counts windingwise
91 wind += int(c);
92 //std::cout << "!!!!!" << int(c) << " ";
93 }
94 iter = ++next;
95 }
96 goto cont;
97 }
98 if(next==start) return wind;
99 }
100 //Looks like it looped, which means everything's flat
101 return 0;
102 }
104 cont:(void)0;
105 }
106 return wind;
107 }
109 /**
110 * This function should only be applied to simple paths (regions), as otherwise
111 * a boolean winding direction is undefined. It returns true for fill, false for
112 * hole. Defaults to using the sign of area when it reaches funny cases.
113 */
114 bool path_direction(Path const &p) {
115 if(p.empty()) return false;
116 //could probably be more efficient, but this is a quick job
117 double y = p.initialPoint()[Y];
118 double x = p.initialPoint()[X];
119 Cmp res = cmp(p[0].finalPoint()[Y], y);
120 goto doh;
121 for(unsigned i = 1; i <= p.size(); i++) {
122 Cmp final_to_ray = cmp(p[i].finalPoint()[Y], y);
123 Cmp initial_to_ray = cmp(p[i].initialPoint()[Y], y);
124 // if y is included, these will have opposite values, giving order.
125 Cmp c = cmp(final_to_ray, initial_to_ray);
126 if(c != EQUAL_TO) {
127 std::vector<double> rs = p[i].roots(y, Y);
128 for(unsigned j = 0; j < rs.size(); j++) {
129 double nx = p[i].valueAt(rs[j], X);
130 if(nx > x) {
131 x = nx;
132 res = c;
133 }
134 }
135 } else if(final_to_ray == EQUAL_TO) goto doh;
136 }
137 return res < 0;
139 doh:
140 //Otherwise fallback on area
142 Piecewise<D2<SBasis> > pw = p.toPwSb();
143 double area;
144 Point centre;
145 Geom::centroid(pw, centre, area);
146 return area > 0;
147 }
149 //pair intersect code based on njh's pair-intersect
151 /** A little sugar for appending a list to another */
152 template<typename T>
153 void append(T &a, T const &b) {
154 a.insert(a.end(), b.begin(), b.end());
155 }
157 /**
158 * Finds the intersection between the lines defined by A0 & A1, and B0 & B1.
159 * Returns through the last 3 parameters, returning the t-values on the lines
160 * and the cross-product of the deltas (a useful byproduct). The return value
161 * indicates if the time values are within their proper range on the line segments.
162 */
163 bool
164 linear_intersect(Point A0, Point A1, Point B0, Point B1,
165 double &tA, double &tB, double &det) {
166 // kramers rule as cross products
167 Point Ad = A1 - A0,
168 Bd = B1 - B0,
169 d = B0 - A0;
170 det = cross(Ad, Bd);
171 if( 1.0 + det == 1.0 )
172 return false;
173 else
174 {
175 double detinv = 1.0 / det;
176 tA = cross(d, Bd) * detinv;
177 tB = cross(d, Ad) * detinv;
178 return tA >= 0. && tA <= 1. && tB >= 0. && tB <= 1.;
179 }
180 }
183 #if 0
184 typedef union dbl_64{
185 long long i64;
186 double d64;
187 };
189 static double EpsilonOf(double value)
190 {
191 dbl_64 s;
192 s.d64 = value;
193 if(s.i64 == 0)
194 {
195 s.i64++;
196 return s.d64 - value;
197 }
198 else if(s.i64-- < 0)
199 return s.d64 - value;
200 else
201 return value - s.d64;
202 }
203 #endif
205 #ifdef HAVE_GSL
206 struct rparams {
207 Curve const &A;
208 Curve const &B;
209 };
211 static int
212 intersect_polish_f (const gsl_vector * x, void *params,
213 gsl_vector * f)
214 {
215 const double x0 = gsl_vector_get (x, 0);
216 const double x1 = gsl_vector_get (x, 1);
218 Geom::Point dx = ((struct rparams *) params)->A(x0) -
219 ((struct rparams *) params)->B(x1);
221 gsl_vector_set (f, 0, dx[0]);
222 gsl_vector_set (f, 1, dx[1]);
224 return GSL_SUCCESS;
225 }
226 #endif
228 static void
229 intersect_polish_root (Curve const &A, double &s,
230 Curve const &B, double &t) {
231 int status;
232 size_t iter = 0;
233 std::vector<Point> as, bs;
234 as = A.pointAndDerivatives(s, 2);
235 bs = B.pointAndDerivatives(t, 2);
236 Point F = as[0] - bs[0];
237 double best = dot(F, F);
239 for(int i = 0; i < 4; i++) {
241 /**
242 we want to solve
243 J*(x1 - x0) = f(x0)
245 |dA(s)[0] -dB(t)[0]| (X1 - X0) = A(s) - B(t)
246 |dA(s)[1] -dB(t)[1]|
247 **/
249 // We're using the standard transformation matricies, which is numerically rather poor. Much better to solve the equation using elimination.
251 Matrix jack(as[1][0], as[1][1],
252 -bs[1][0], -bs[1][1],
253 0, 0);
254 Point soln = (F)*jack.inverse();
255 double ns = s - soln[0];
256 double nt = t - soln[1];
258 if (ns<0) ns=0;
259 else if (ns>1) ns=1;
260 if (nt<0) nt=0;
261 else if (nt>1) nt=1;
263 as = A.pointAndDerivatives(ns, 2);
264 bs = B.pointAndDerivatives(nt, 2);
265 F = as[0] - bs[0];
266 double trial = dot(F, F);
267 if (trial > best*0.1) // we have standards, you know
268 // At this point we could do a line search
269 break;
270 best = trial;
271 s = ns;
272 t = nt;
273 }
275 #ifdef HAVE_GSL
276 if(0) { // the GSL version is more accurate, but taints this with GPL
277 const size_t n = 2;
278 struct rparams p = {A, B};
279 gsl_multiroot_function f = {&intersect_polish_f, n, &p};
281 double x_init[2] = {s, t};
282 gsl_vector *x = gsl_vector_alloc (n);
284 gsl_vector_set (x, 0, x_init[0]);
285 gsl_vector_set (x, 1, x_init[1]);
287 const gsl_multiroot_fsolver_type *T = gsl_multiroot_fsolver_hybrids;
288 gsl_multiroot_fsolver *sol = gsl_multiroot_fsolver_alloc (T, 2);
289 gsl_multiroot_fsolver_set (sol, &f, x);
291 do
292 {
293 iter++;
294 status = gsl_multiroot_fsolver_iterate (sol);
296 if (status) /* check if solver is stuck */
297 break;
299 status =
300 gsl_multiroot_test_residual (sol->f, 1e-12);
301 }
302 while (status == GSL_CONTINUE && iter < 1000);
304 s = gsl_vector_get (sol->x, 0);
305 t = gsl_vector_get (sol->x, 1);
307 gsl_multiroot_fsolver_free (sol);
308 gsl_vector_free (x);
309 }
310 #endif
311 }
313 /**
314 * This uses the local bounds functions of curves to generically intersect two.
315 * It passes in the curves, time intervals, and keeps track of depth, while
316 * returning the results through the Crossings parameter.
317 */
318 void pair_intersect(Curve const & A, double Al, double Ah,
319 Curve const & B, double Bl, double Bh,
320 Crossings &ret, unsigned depth = 0) {
321 // std::cout << depth << "(" << Al << ", " << Ah << ")\n";
322 OptRect Ar = A.boundsLocal(Interval(Al, Ah));
323 if (!Ar) return;
325 OptRect Br = B.boundsLocal(Interval(Bl, Bh));
326 if (!Br) return;
328 if(! Ar->intersects(*Br)) 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(A.pointAt(Al), A.pointAt(Ah),
335 B.pointAt(Bl), B.pointAt(Bh),
336 tA, tB, c)) {
337 tA = tA * (Ah - Al) + Al;
338 tB = tB * (Bh - Bl) + Bl;
339 intersect_polish_root(A, tA,
340 B, tB);
341 if(depth % 2)
342 ret.push_back(Crossing(tB, tA, c < 0));
343 else
344 ret.push_back(Crossing(tA, tB, c > 0));
345 return;
346 }
347 }
348 if(depth > 12) return;
349 double mid = (Bl + Bh)/2;
350 pair_intersect(B, Bl, mid,
351 A, Al, Ah,
352 ret, depth+1);
353 pair_intersect(B, mid, Bh,
354 A, Al, Ah,
355 ret, depth+1);
356 }
358 Crossings pair_intersect(Curve const & A, Interval const &Ad,
359 Curve const & B, Interval const &Bd) {
360 Crossings ret;
361 pair_intersect(A, Ad.min(), Ad.max(), B, Bd.min(), Bd.max(), ret);
362 return ret;
363 }
365 /** A simple wrapper around pair_intersect */
366 Crossings SimpleCrosser::crossings(Curve const &a, Curve const &b) {
367 Crossings ret;
368 pair_intersect(a, 0, 1, b, 0, 1, ret);
369 return ret;
370 }
373 //same as below but curves not paths
374 void mono_intersect(Curve const &A, double Al, double Ah,
375 Curve const &B, double Bl, double Bh,
376 Crossings &ret, double tol = 0.1, unsigned depth = 0) {
377 if( Al >= Ah || Bl >= Bh) return;
378 //std::cout << " " << depth << "[" << Al << ", " << Ah << "]" << "[" << Bl << ", " << Bh << "]";
380 Point A0 = A.pointAt(Al), A1 = A.pointAt(Ah),
381 B0 = B.pointAt(Bl), B1 = B.pointAt(Bh);
382 //inline code that this implies? (without rect/interval construction)
383 Rect Ar = Rect(A0, A1), Br = Rect(B0, B1);
384 if(!Ar.intersects(Br) || A0 == A1 || B0 == B1) return;
386 if(depth > 12 || (Ar.maxExtent() < tol && Ar.maxExtent() < tol)) {
387 double tA, tB, c;
388 if(linear_intersect(A.pointAt(Al), A.pointAt(Ah),
389 B.pointAt(Bl), B.pointAt(Bh),
390 tA, tB, c)) {
391 tA = tA * (Ah - Al) + Al;
392 tB = tB * (Bh - Bl) + Bl;
393 intersect_polish_root(A, tA,
394 B, tB);
395 if(depth % 2)
396 ret.push_back(Crossing(tB, tA, c < 0));
397 else
398 ret.push_back(Crossing(tA, tB, c > 0));
399 return;
400 }
401 }
402 if(depth > 12) return;
403 double mid = (Bl + Bh)/2;
404 mono_intersect(B, Bl, mid,
405 A, Al, Ah,
406 ret, tol, depth+1);
407 mono_intersect(B, mid, Bh,
408 A, Al, Ah,
409 ret, tol, depth+1);
410 }
412 Crossings mono_intersect(Curve const & A, Interval const &Ad,
413 Curve const & B, Interval const &Bd) {
414 Crossings ret;
415 mono_intersect(A, Ad.min(), Ad.max(), B, Bd.min(), Bd.max(), ret);
416 return ret;
417 }
419 /**
420 * Takes two paths and time ranges on them, with the invariant that the
421 * paths are monotonic on the range. Splits A when the linear intersection
422 * doesn't exist or is inaccurate. Uses the fact that it is monotonic to
423 * do very fast local bounds.
424 */
425 void mono_pair(Path const &A, double Al, double Ah,
426 Path const &B, double Bl, double Bh,
427 Crossings &ret, double /*tol*/, unsigned depth = 0) {
428 if( Al >= Ah || Bl >= Bh) return;
429 std::cout << " " << depth << "[" << Al << ", " << Ah << "]" << "[" << Bl << ", " << Bh << "]";
431 Point A0 = A.pointAt(Al), A1 = A.pointAt(Ah),
432 B0 = B.pointAt(Bl), B1 = B.pointAt(Bh);
433 //inline code that this implies? (without rect/interval construction)
434 Rect Ar = Rect(A0, A1), Br = Rect(B0, B1);
435 if(!Ar.intersects(Br) || A0 == A1 || B0 == B1) return;
437 if(depth > 12 || (Ar.maxExtent() < 0.1 && Ar.maxExtent() < 0.1)) {
438 double tA, tB, c;
439 if(linear_intersect(A0, A1, B0, B1,
440 tA, tB, c)) {
441 tA = tA * (Ah - Al) + Al;
442 tB = tB * (Bh - Bl) + Bl;
443 if(depth % 2)
444 ret.push_back(Crossing(tB, tA, c < 0));
445 else
446 ret.push_back(Crossing(tA, tB, c > 0));
447 return;
448 }
449 }
450 if(depth > 12) return;
451 double mid = (Bl + Bh)/2;
452 mono_pair(B, Bl, mid,
453 A, Al, Ah,
454 ret, depth+1);
455 mono_pair(B, mid, Bh,
456 A, Al, Ah,
457 ret, depth+1);
458 }
460 /** This returns the times when the x or y derivative is 0 in the curve. */
461 std::vector<double> curve_mono_splits(Curve const &d) {
462 Curve* deriv = d.derivative();
463 std::vector<double> rs = d.roots(0, X);
464 append(rs, d.roots(0, Y));
465 delete deriv;
466 std::sort(rs.begin(), rs.end());
467 return rs;
468 }
470 /** Convenience function to add a value to each entry in a vector of doubles. */
471 std::vector<double> offset_doubles(std::vector<double> const &x, double offs) {
472 std::vector<double> ret;
473 for(unsigned i = 0; i < x.size(); i++) {
474 ret.push_back(x[i] + offs);
475 }
476 return ret;
477 }
479 /**
480 * Finds all the monotonic splits for a path. Only includes the split between
481 * curves if they switch derivative directions at that point.
482 */
483 std::vector<double> path_mono_splits(Path const &p) {
484 std::vector<double> ret;
485 if(p.empty()) return ret;
487 bool pdx=2, pdy=2; //Previous derivative direction
488 for(unsigned i = 0; i <= p.size(); i++) {
489 std::vector<double> spl = offset_doubles(curve_mono_splits(p[i]), i);
490 bool dx = p[i].initialPoint()[X] > (spl.empty()? p[i].finalPoint()[X] :
491 p.valueAt(spl.front(), X));
492 bool dy = p[i].initialPoint()[Y] > (spl.empty()? p[i].finalPoint()[Y] :
493 p.valueAt(spl.front(), Y));
494 //The direction changed, include the split time
495 if(dx != pdx || dy != pdy) {
496 ret.push_back(i);
497 pdx = dx; pdy = dy;
498 }
499 append(ret, spl);
500 }
501 return ret;
502 }
504 /**
505 * Applies path_mono_splits to multiple paths, and returns the results such that
506 * time-set i corresponds to Path i.
507 */
508 std::vector<std::vector<double> > paths_mono_splits(std::vector<Path> const &ps) {
509 std::vector<std::vector<double> > ret;
510 for(unsigned i = 0; i < ps.size(); i++)
511 ret.push_back(path_mono_splits(ps[i]));
512 return ret;
513 }
515 /**
516 * Processes the bounds for a list of paths and a list of splits on them, yielding a list of rects for each.
517 * Each entry i corresponds to path i of the input. The number of rects in each entry is guaranteed to be the
518 * number of splits for that path, subtracted by one.
519 */
520 std::vector<std::vector<Rect> > split_bounds(std::vector<Path> const &p, std::vector<std::vector<double> > splits) {
521 std::vector<std::vector<Rect> > ret;
522 for(unsigned i = 0; i < p.size(); i++) {
523 std::vector<Rect> res;
524 for(unsigned j = 1; j < splits[i].size(); j++)
525 res.push_back(Rect(p[i].pointAt(splits[i][j-1]), p[i].pointAt(splits[i][j])));
526 ret.push_back(res);
527 }
528 return ret;
529 }
531 /**
532 * This is the main routine of "MonoCrosser", and implements a monotonic strategy on multiple curves.
533 * Finds crossings between two sets of paths, yielding a CrossingSet. [0, a.size()) of the return correspond
534 * to the sorted crossings of a with paths of b. The rest of the return, [a.size(), a.size() + b.size()],
535 * corresponds to the sorted crossings of b with paths of a.
536 *
537 * This function does two sweeps, one on the bounds of each path, and after that cull, one on the curves within.
538 * This leads to a certain amount of code complexity, however, most of that is factored into the above functions
539 */
540 CrossingSet MonoCrosser::crossings(std::vector<Path> const &a, std::vector<Path> const &b) {
541 if(b.empty()) return CrossingSet(a.size(), Crossings());
542 CrossingSet results(a.size() + b.size(), Crossings());
543 if(a.empty()) return results;
545 std::vector<std::vector<double> > splits_a = paths_mono_splits(a), splits_b = paths_mono_splits(b);
546 std::vector<std::vector<Rect> > bounds_a = split_bounds(a, splits_a), bounds_b = split_bounds(b, splits_b);
548 std::vector<Rect> bounds_a_union, bounds_b_union;
549 for(unsigned i = 0; i < bounds_a.size(); i++) bounds_a_union.push_back(union_list(bounds_a[i]));
550 for(unsigned i = 0; i < bounds_b.size(); i++) bounds_b_union.push_back(union_list(bounds_b[i]));
552 std::vector<std::vector<unsigned> > cull = sweep_bounds(bounds_a_union, bounds_b_union);
553 Crossings n;
554 for(unsigned i = 0; i < cull.size(); i++) {
555 for(unsigned jx = 0; jx < cull[i].size(); jx++) {
556 unsigned j = cull[i][jx];
557 unsigned jc = j + a.size();
558 Crossings res;
560 //Sweep of the monotonic portions
561 std::vector<std::vector<unsigned> > cull2 = sweep_bounds(bounds_a[i], bounds_b[j]);
562 for(unsigned k = 0; k < cull2.size(); k++) {
563 for(unsigned lx = 0; lx < cull2[k].size(); lx++) {
564 unsigned l = cull2[k][lx];
565 mono_pair(a[i], splits_a[i][k-1], splits_a[i][k],
566 b[j], splits_b[j][l-1], splits_b[j][l],
567 res, .1);
568 }
569 }
571 for(unsigned k = 0; k < res.size(); k++) { res[k].a = i; res[k].b = jc; }
573 merge_crossings(results[i], res, i);
574 merge_crossings(results[i], res, jc);
575 }
576 }
578 return results;
579 }
581 /* This function is similar codewise to the MonoCrosser, the main difference is that it deals with
582 * only one set of paths and includes self intersection
583 CrossingSet crossings_among(std::vector<Path> const &p) {
584 CrossingSet results(p.size(), Crossings());
585 if(p.empty()) return results;
587 std::vector<std::vector<double> > splits = paths_mono_splits(p);
588 std::vector<std::vector<Rect> > prs = split_bounds(p, splits);
589 std::vector<Rect> rs;
590 for(unsigned i = 0; i < prs.size(); i++) rs.push_back(union_list(prs[i]));
592 std::vector<std::vector<unsigned> > cull = sweep_bounds(rs);
594 //we actually want to do the self-intersections, so add em in:
595 for(unsigned i = 0; i < cull.size(); i++) cull[i].push_back(i);
597 for(unsigned i = 0; i < cull.size(); i++) {
598 for(unsigned jx = 0; jx < cull[i].size(); jx++) {
599 unsigned j = cull[i][jx];
600 Crossings res;
602 //Sweep of the monotonic portions
603 std::vector<std::vector<unsigned> > cull2 = sweep_bounds(prs[i], prs[j]);
604 for(unsigned k = 0; k < cull2.size(); k++) {
605 for(unsigned lx = 0; lx < cull2[k].size(); lx++) {
606 unsigned l = cull2[k][lx];
607 mono_pair(p[i], splits[i][k-1], splits[i][k],
608 p[j], splits[j][l-1], splits[j][l],
609 res, .1);
610 }
611 }
613 for(unsigned k = 0; k < res.size(); k++) { res[k].a = i; res[k].b = j; }
615 merge_crossings(results[i], res, i);
616 merge_crossings(results[j], res, j);
617 }
618 }
620 return results;
621 }
622 */
625 Crossings curve_self_crossings(Curve const &a) {
626 Crossings res;
627 std::vector<double> spl;
628 spl.push_back(0);
629 append(spl, curve_mono_splits(a));
630 spl.push_back(1);
631 for(unsigned i = 1; i < spl.size(); i++)
632 for(unsigned j = i+1; j < spl.size(); j++)
633 pair_intersect(a, spl[i-1], spl[i], a, spl[j-1], spl[j], res);
634 return res;
635 }
637 /*
638 void mono_curve_intersect(Curve const & A, double Al, double Ah,
639 Curve const & B, double Bl, double Bh,
640 Crossings &ret, unsigned depth=0) {
641 // std::cout << depth << "(" << Al << ", " << Ah << ")\n";
642 Point A0 = A.pointAt(Al), A1 = A.pointAt(Ah),
643 B0 = B.pointAt(Bl), B1 = B.pointAt(Bh);
644 //inline code that this implies? (without rect/interval construction)
645 if(!Rect(A0, A1).intersects(Rect(B0, B1)) || A0 == A1 || B0 == B1) return;
647 //Checks the general linearity of the function
648 if((depth > 12) || (A.boundsLocal(Interval(Al, Ah), 1).maxExtent() < 0.1
649 && B.boundsLocal(Interval(Bl, Bh), 1).maxExtent() < 0.1)) {
650 double tA, tB, c;
651 if(linear_intersect(A0, A1, B0, B1, tA, tB, c)) {
652 tA = tA * (Ah - Al) + Al;
653 tB = tB * (Bh - Bl) + Bl;
654 if(depth % 2)
655 ret.push_back(Crossing(tB, tA, c < 0));
656 else
657 ret.push_back(Crossing(tA, tB, c > 0));
658 return;
659 }
660 }
661 if(depth > 12) return;
662 double mid = (Bl + Bh)/2;
663 mono_curve_intersect(B, Bl, mid,
664 A, Al, Ah,
665 ret, depth+1);
666 mono_curve_intersect(B, mid, Bh,
667 A, Al, Ah,
668 ret, depth+1);
669 }
671 std::vector<std::vector<double> > curves_mono_splits(Path const &p) {
672 std::vector<std::vector<double> > ret;
673 for(unsigned i = 0; i <= p.size(); i++) {
674 std::vector<double> spl;
675 spl.push_back(0);
676 append(spl, curve_mono_splits(p[i]));
677 spl.push_back(1);
678 ret.push_back(spl);
679 }
680 }
682 std::vector<std::vector<Rect> > curves_split_bounds(Path const &p, std::vector<std::vector<double> > splits) {
683 std::vector<std::vector<Rect> > ret;
684 for(unsigned i = 0; i < splits.size(); i++) {
685 std::vector<Rect> res;
686 for(unsigned j = 1; j < splits[i].size(); j++)
687 res.push_back(Rect(p.pointAt(splits[i][j-1]+i), p.pointAt(splits[i][j]+i)));
688 ret.push_back(res);
689 }
690 return ret;
691 }
693 Crossings path_self_crossings(Path const &p) {
694 Crossings ret;
695 std::vector<std::vector<unsigned> > cull = sweep_bounds(bounds(p));
696 std::vector<std::vector<double> > spl = curves_mono_splits(p);
697 std::vector<std::vector<Rect> > bnds = curves_split_bounds(p, spl);
698 for(unsigned i = 0; i < cull.size(); i++) {
699 Crossings res;
700 for(unsigned k = 1; k < spl[i].size(); k++)
701 for(unsigned l = k+1; l < spl[i].size(); l++)
702 mono_curve_intersect(p[i], spl[i][k-1], spl[i][k], p[i], spl[i][l-1], spl[i][l], res);
703 offset_crossings(res, i, i);
704 append(ret, res);
705 for(unsigned jx = 0; jx < cull[i].size(); jx++) {
706 unsigned j = cull[i][jx];
707 res.clear();
709 std::vector<std::vector<unsigned> > cull2 = sweep_bounds(bnds[i], bnds[j]);
710 for(unsigned k = 0; k < cull2.size(); k++) {
711 for(unsigned lx = 0; lx < cull2[k].size(); lx++) {
712 unsigned l = cull2[k][lx];
713 mono_curve_intersect(p[i], spl[i][k-1], spl[i][k], p[j], spl[j][l-1], spl[j][l], res);
714 }
715 }
717 //if(fabs(int(i)-j) == 1 || fabs(int(i)-j) == p.size()-1) {
718 Crossings res2;
719 for(unsigned k = 0; k < res.size(); k++) {
720 if(res[k].ta != 0 && res[k].ta != 1 && res[k].tb != 0 && res[k].tb != 1) {
721 res.push_back(res[k]);
722 }
723 }
724 res = res2;
725 //}
726 offset_crossings(res, i, j);
727 append(ret, res);
728 }
729 }
730 return ret;
731 }
732 */
734 Crossings self_crossings(Path const &p) {
735 Crossings ret;
736 std::vector<std::vector<unsigned> > cull = sweep_bounds(bounds(p));
737 for(unsigned i = 0; i < cull.size(); i++) {
738 Crossings res = curve_self_crossings(p[i]);
739 offset_crossings(res, i, i);
740 append(ret, res);
741 for(unsigned jx = 0; jx < cull[i].size(); jx++) {
742 unsigned j = cull[i][jx];
743 res.clear();
744 pair_intersect(p[i], 0, 1, p[j], 0, 1, res);
746 //if(fabs(int(i)-j) == 1 || fabs(int(i)-j) == p.size()-1) {
747 Crossings res2;
748 for(unsigned k = 0; k < res.size(); k++) {
749 if(res[k].ta != 0 && res[k].ta != 1 && res[k].tb != 0 && res[k].tb != 1) {
750 res2.push_back(res[k]);
751 }
752 }
753 res = res2;
754 //}
755 offset_crossings(res, i, j);
756 append(ret, res);
757 }
758 }
759 return ret;
760 }
762 void flip_crossings(Crossings &crs) {
763 for(unsigned i = 0; i < crs.size(); i++)
764 crs[i] = Crossing(crs[i].tb, crs[i].ta, crs[i].b, crs[i].a, !crs[i].dir);
765 }
767 CrossingSet crossings_among(std::vector<Path> const &p) {
768 CrossingSet results(p.size(), Crossings());
769 if(p.empty()) return results;
771 SimpleCrosser cc;
773 std::vector<std::vector<unsigned> > cull = sweep_bounds(bounds(p));
774 for(unsigned i = 0; i < cull.size(); i++) {
775 Crossings res = self_crossings(p[i]);
776 for(unsigned k = 0; k < res.size(); k++) { res[k].a = res[k].b = i; }
777 merge_crossings(results[i], res, i);
778 flip_crossings(res);
779 merge_crossings(results[i], res, i);
780 for(unsigned jx = 0; jx < cull[i].size(); jx++) {
781 unsigned j = cull[i][jx];
783 Crossings res = cc.crossings(p[i], p[j]);
784 for(unsigned k = 0; k < res.size(); k++) { res[k].a = i; res[k].b = j; }
785 merge_crossings(results[i], res, i);
786 merge_crossings(results[j], res, j);
787 }
788 }
789 return results;
790 }
792 }
794 /*
795 Local Variables:
796 mode:c++
797 c-file-style:"stroustrup"
798 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
799 indent-tabs-mode:nil
800 fill-column:99
801 End:
802 */
803 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :