GCC Code Coverage Report

Directory: ./
File: src/2geom/path-intersection.cpp
Date: 2024-03-18 17:01:34
Exec Total Coverage
Lines: 8 275 2.9%
Functions: 1 23 4.3%
Branches: 4 490 0.8%
Line Branch Exec Source
1 #include <2geom/path-intersection.h>
2
3 #include <2geom/ord.h>
4
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
12
13 namespace Geom {
14
15 /// Compute winding number of the path at the specified point
16 ✗ int winding(Path const &path, Point const &p) {
17 ✗ return path.winding(p);
18 }
19
20 /**
21 * This function should only be applied to simple paths (regions), as otherwise
22 * a boolean winding direction is undefined. It returns true for fill, false for
23 * hole. Defaults to using the sign of area when it reaches funny cases.
24 */
25 14 bool path_direction(Path const &p) {
26
2/4
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 14 times.
 14 if(p.empty()) return false;
27
28 /*goto doh;
29 //could probably be more efficient, but this is a quick job
30 double y = p.initialPoint()[Y];
31 double x = p.initialPoint()[X];
32 Cmp res = cmp(p[0].finalPoint()[Y], y);
33 for(unsigned i = 1; i < p.size(); i++) {
34 Cmp final_to_ray = cmp(p[i].finalPoint()[Y], y);
35 Cmp initial_to_ray = cmp(p[i].initialPoint()[Y], y);
36 // if y is included, these will have opposite values, giving order.
37 Cmp c = cmp(final_to_ray, initial_to_ray);
38 if(c != EQUAL_TO) {
39 std::vector<double> rs = p[i].roots(y, Y);
40 for(unsigned j = 0; j < rs.size(); j++) {
41 double nx = p[i].valueAt(rs[j], X);
42 if(nx > x) {
43 x = nx;
44 res = c;
45 }
46 }
47 } else if(final_to_ray == EQUAL_TO) goto doh;
48 }
49 return res < 0;
50
51 doh:*/
52 //Otherwise fallback on area
53
54
1/2
✓ Branch 2 taken 14 times.
✗ Branch 3 not taken.
 14 Piecewise<D2<SBasis> > pw = p.toPwSb();
55 14 double area;
56 14 Point centre;
57
1/2
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
 14 Geom::centroid(pw, centre, area);
58 14 return area > 0;
59 14 }
60
61 //pair intersect code based on njh's pair-intersect
62
63 /** A little sugar for appending a list to another */
64 template<typename T>
65 ✗ void append(T &a, T const &b) {
66 ✗ a.insert(a.end(), b.begin(), b.end());
67 ✗ }
68
69 /**
70 * Finds the intersection between the lines defined by A0 & A1, and B0 & B1.
71 * Returns through the last 3 parameters, returning the t-values on the lines
72 * and the cross-product of the deltas (a useful byproduct). The return value
73 * indicates if the time values are within their proper range on the line segments.
74 */
75 bool
76 ✗ linear_intersect(Point const &A0, Point const &A1, Point const &B0, Point const &B1,
77 double &tA, double &tB, double &det) {
78 ✗ bool both_lines_non_zero = (!are_near(A0, A1)) && (!are_near(B0, B1));
79
80 // Cramer's rule as cross products
81 ✗ Point Ad = A1 - A0,
82 ✗ Bd = B1 - B0,
83 ✗ d = B0 - A0;
84 ✗ det = cross(Ad, Bd);
85
86 ✗ double det_rel = det; // Calculate the determinant of the normalized vectors
87 ✗ if (both_lines_non_zero) {
88 ✗ det_rel /= Ad.length();
89 ✗ det_rel /= Bd.length();
90 }
91
92 ✗ if( fabs(det_rel) < 1e-12 ) { // If the cross product is NEARLY zero,
93 // Then one of the linesegments might have length zero
94 ✗ if (both_lines_non_zero) {
95 // If that's not the case, then we must have either:
96 // - parallel lines, with no intersections, or
97 // - coincident lines, with an infinite number of intersections
98 // Either is quite useless, so we'll just bail out
99 ✗ return false;
100 } // Else, one of the linesegments is zero, and we might still be able to calculate a single intersection point
101 } // Else we haven't bailed out, and we'll try to calculate the intersections
102
103 ✗ double detinv = 1.0 / det;
104 ✗ tA = cross(d, Bd) * detinv;
105 ✗ tB = cross(d, Ad) * detinv;
106 ✗ return (tA >= 0.) && (tA <= 1.) && (tB >= 0.) && (tB <= 1.);
107 }
108
109
110 #if 0
111 typedef union dbl_64{
112 long long i64;
113 double d64;
114 };
115
116 static double EpsilonOf(double value)
117 {
118 dbl_64 s;
119 s.d64 = value;
120 if(s.i64 == 0)
121 {
122 s.i64++;
123 return s.d64 - value;
124 }
125 else if(s.i64-- < 0)
126 return s.d64 - value;
127 else
128 return value - s.d64;
129 }
130 #endif
131
132 #ifdef HAVE_GSL
133 struct rparams {
134 Curve const &A;
135 Curve const &B;
136 };
137
138 static int
139 ✗ intersect_polish_f (const gsl_vector * x, void *params,
140 gsl_vector * f)
141 {
142 ✗ const double x0 = gsl_vector_get (x, 0);
143 ✗ const double x1 = gsl_vector_get (x, 1);
144
145 ✗ Geom::Point dx = ((struct rparams *) params)->A(x0) -
146 ✗ ((struct rparams *) params)->B(x1);
147
148 ✗ gsl_vector_set (f, 0, dx[0]);
149 ✗ gsl_vector_set (f, 1, dx[1]);
150
151 ✗ return GSL_SUCCESS;
152 }
153 #endif
154
155 static void
156 ✗ intersect_polish_root (Curve const &A, double &s, Curve const &B, double &t)
157 {
158 ✗ std::vector<Point> as, bs;
159 ✗ as = A.pointAndDerivatives(s, 2);
160 ✗ bs = B.pointAndDerivatives(t, 2);
161 ✗ Point F = as[0] - bs[0];
162 ✗ double best = dot(F, F);
163
164 ✗ for(int i = 0; i < 4; i++) {
165
166 /**
167 we want to solve
168 J*(x1 - x0) = f(x0)
169
170 |dA(s)[0] -dB(t)[0]| (X1 - X0) = A(s) - B(t)
171 |dA(s)[1] -dB(t)[1]|
172 **/
173
174 // We're using the standard transformation matricies, which is numerically rather poor. Much better to solve the equation using elimination.
175
176 ✗ Affine jack(as[1][0], as[1][1],
177 ✗ -bs[1][0], -bs[1][1],
178 ✗ 0, 0);
179 ✗ Point soln = (F)*jack.inverse();
180 ✗ double ns = s - soln[0];
181 ✗ double nt = t - soln[1];
182
183 ✗ if (ns<0) ns=0;
184 ✗ else if (ns>1) ns=1;
185 ✗ if (nt<0) nt=0;
186 ✗ else if (nt>1) nt=1;
187
188 ✗ as = A.pointAndDerivatives(ns, 2);
189 ✗ bs = B.pointAndDerivatives(nt, 2);
190 ✗ F = as[0] - bs[0];
191 ✗ double trial = dot(F, F);
192 ✗ if (trial > best*0.1) // we have standards, you know
193 // At this point we could do a line search
194 ✗ break;
195 ✗ best = trial;
196 ✗ s = ns;
197 ✗ t = nt;
198 }
199
200 #ifdef HAVE_GSL
201 if(0) { // the GSL version is more accurate, but taints this with GPL
202 int status;
203 size_t iter = 0;
204 const size_t n = 2;
205 struct rparams p = {A, B};
206 gsl_multiroot_function f = {&intersect_polish_f, n, &p};
207
208 double x_init[2] = {s, t};
209 gsl_vector *x = gsl_vector_alloc (n);
210
211 gsl_vector_set (x, 0, x_init[0]);
212 gsl_vector_set (x, 1, x_init[1]);
213
214 const gsl_multiroot_fsolver_type *T = gsl_multiroot_fsolver_hybrids;
215 gsl_multiroot_fsolver *sol = gsl_multiroot_fsolver_alloc (T, 2);
216 gsl_multiroot_fsolver_set (sol, &f, x);
217
218 do
219 {
220 iter++;
221 status = gsl_multiroot_fsolver_iterate (sol);
222
223 if (status) /* check if solver is stuck */
224 break;
225
226 status =
227 gsl_multiroot_test_residual (sol->f, 1e-12);
228 }
229 while (status == GSL_CONTINUE && iter < 1000);
230
231 s = gsl_vector_get (sol->x, 0);
232 t = gsl_vector_get (sol->x, 1);
233
234 gsl_multiroot_fsolver_free (sol);
235 gsl_vector_free (x);
236 }
237 #endif
238 ✗ }
239
240 /**
241 * This uses the local bounds functions of curves to generically intersect two.
242 * It passes in the curves, time intervals, and keeps track of depth, while
243 * returning the results through the Crossings parameter.
244 */
245 ✗ void pair_intersect(Curve const & A, double Al, double Ah,
246 Curve const & B, double Bl, double Bh,
247 Crossings &ret, unsigned depth = 0) {
248 // std::cout << depth << "(" << Al << ", " << Ah << ")\n";
249 ✗ OptRect Ar = A.boundsLocal(Interval(Al, Ah));
250 ✗ if (!Ar) return;
251
252 ✗ OptRect Br = B.boundsLocal(Interval(Bl, Bh));
253 ✗ if (!Br) return;
254
255 ✗ if(! Ar->intersects(*Br)) return;
256
257 //Checks the general linearity of the function
258 ✗ if((depth > 12)) { // || (A.boundsLocal(Interval(Al, Ah), 1).maxExtent() < 0.1
259 //&& B.boundsLocal(Interval(Bl, Bh), 1).maxExtent() < 0.1)) {
260 ✗ double tA, tB, c;
261 ✗ if(linear_intersect(A.pointAt(Al), A.pointAt(Ah),
262 ✗ B.pointAt(Bl), B.pointAt(Bh),
263 tA, tB, c)) {
264 ✗ tA = tA * (Ah - Al) + Al;
265 ✗ tB = tB * (Bh - Bl) + Bl;
266 ✗ intersect_polish_root(A, tA,
267 B, tB);
268 ✗ if(depth % 2)
269 ✗ ret.push_back(Crossing(tB, tA, c < 0));
270 else
271 ✗ ret.push_back(Crossing(tA, tB, c > 0));
272 ✗ return;
273 }
274 }
275 ✗ if(depth > 12) return;
276 ✗ double mid = (Bl + Bh)/2;
277 ✗ pair_intersect(B, Bl, mid,
278 A, Al, Ah,
279 ret, depth+1);
280 ✗ pair_intersect(B, mid, Bh,
281 A, Al, Ah,
282 ret, depth+1);
283 }
284
285 ✗ Crossings pair_intersect(Curve const & A, Interval const &Ad,
286 Curve const & B, Interval const &Bd) {
287 ✗ Crossings ret;
288 ✗ pair_intersect(A, Ad.min(), Ad.max(), B, Bd.min(), Bd.max(), ret);
289 ✗ return ret;
290 ✗ }
291
292 /** A simple wrapper around pair_intersect */
293 ✗ Crossings SimpleCrosser::crossings(Curve const &a, Curve const &b) {
294 ✗ Crossings ret;
295 ✗ pair_intersect(a, 0, 1, b, 0, 1, ret);
296 ✗ return ret;
297 ✗ }
298
299
300 //same as below but curves not paths
301 ✗ void mono_intersect(Curve const &A, double Al, double Ah,
302 Curve const &B, double Bl, double Bh,
303 Crossings &ret, double tol = 0.1, unsigned depth = 0) {
304 ✗ if( Al >= Ah || Bl >= Bh) return;
305 //std::cout << " " << depth << "[" << Al << ", " << Ah << "]" << "[" << Bl << ", " << Bh << "]";
306
307 ✗ Point A0 = A.pointAt(Al), A1 = A.pointAt(Ah),
308 ✗ B0 = B.pointAt(Bl), B1 = B.pointAt(Bh);
309 //inline code that this implies? (without rect/interval construction)
310 ✗ Rect Ar = Rect(A0, A1), Br = Rect(B0, B1);
311 ✗ if(!Ar.intersects(Br) || A0 == A1 || B0 == B1) return;
312
313 ✗ if(depth > 12 || (Ar.maxExtent() < tol && Ar.maxExtent() < tol)) {
314 ✗ double tA, tB, c;
315 ✗ if(linear_intersect(A.pointAt(Al), A.pointAt(Ah),
316 ✗ B.pointAt(Bl), B.pointAt(Bh),
317 tA, tB, c)) {
318 ✗ tA = tA * (Ah - Al) + Al;
319 ✗ tB = tB * (Bh - Bl) + Bl;
320 ✗ intersect_polish_root(A, tA,
321 B, tB);
322 ✗ if(depth % 2)
323 ✗ ret.push_back(Crossing(tB, tA, c < 0));
324 else
325 ✗ ret.push_back(Crossing(tA, tB, c > 0));
326 ✗ return;
327 }
328 }
329 ✗ if(depth > 12) return;
330 ✗ double mid = (Bl + Bh)/2;
331 ✗ mono_intersect(B, Bl, mid,
332 A, Al, Ah,
333 ret, tol, depth+1);
334 ✗ mono_intersect(B, mid, Bh,
335 A, Al, Ah,
336 ret, tol, depth+1);
337 }
338
339 ✗ Crossings mono_intersect(Curve const & A, Interval const &Ad,
340 Curve const & B, Interval const &Bd) {
341 ✗ Crossings ret;
342 ✗ mono_intersect(A, Ad.min(), Ad.max(), B, Bd.min(), Bd.max(), ret);
343 ✗ return ret;
344 ✗ }
345
346 /**
347 * Takes two paths and time ranges on them, with the invariant that the
348 * paths are monotonic on the range. Splits A when the linear intersection
349 * doesn't exist or is inaccurate. Uses the fact that it is monotonic to
350 * do very fast local bounds.
351 */
352 ✗ void mono_pair(Path const &A, double Al, double Ah,
353 Path const &B, double Bl, double Bh,
354 Crossings &ret, double /*tol*/, unsigned depth = 0) {
355 ✗ if( Al >= Ah || Bl >= Bh) return;
356 ✗ std::cout << " " << depth << "[" << Al << ", " << Ah << "]" << "[" << Bl << ", " << Bh << "]";
357
358 ✗ Point A0 = A.pointAt(Al), A1 = A.pointAt(Ah),
359 ✗ B0 = B.pointAt(Bl), B1 = B.pointAt(Bh);
360 //inline code that this implies? (without rect/interval construction)
361 ✗ Rect Ar = Rect(A0, A1), Br = Rect(B0, B1);
362 ✗ if(!Ar.intersects(Br) || A0 == A1 || B0 == B1) return;
363
364 ✗ if(depth > 12 || (Ar.maxExtent() < 0.1 && Ar.maxExtent() < 0.1)) {
365 ✗ double tA, tB, c;
366 ✗ if(linear_intersect(A0, A1, B0, B1,
367 tA, tB, c)) {
368 ✗ tA = tA * (Ah - Al) + Al;
369 ✗ tB = tB * (Bh - Bl) + Bl;
370 ✗ if(depth % 2)
371 ✗ ret.push_back(Crossing(tB, tA, c < 0));
372 else
373 ✗ ret.push_back(Crossing(tA, tB, c > 0));
374 ✗ return;
375 }
376 }
377 ✗ if(depth > 12) return;
378 ✗ double mid = (Bl + Bh)/2;
379 ✗ mono_pair(B, Bl, mid,
380 A, Al, Ah,
381 ✗ ret, depth+1);
382 ✗ mono_pair(B, mid, Bh,
383 A, Al, Ah,
384 ✗ ret, depth+1);
385 }
386
387 /** This returns the times when the x or y derivative is 0 in the curve. */
388 ✗ std::vector<double> curve_mono_splits(Curve const &d) {
389 ✗ Curve* deriv = d.derivative();
390 ✗ std::vector<double> rs = deriv->roots(0, X);
391 ✗ append(rs, deriv->roots(0, Y));
392 ✗ delete deriv;
393 ✗ std::sort(rs.begin(), rs.end());
394 ✗ return rs;
395 ✗ }
396
397 /** Convenience function to add a value to each entry in a vector of doubles. */
398 ✗ std::vector<double> offset_doubles(std::vector<double> const &x, double offs) {
399 ✗ std::vector<double> ret;
400 ✗ for(double i : x) {
401 ✗ ret.push_back(i + offs);
402 }
403 ✗ return ret;
404 ✗ }
405
406 /**
407 * Finds all the monotonic splits for a path. Only includes the split between
408 * curves if they switch derivative directions at that point.
409 */
410 ✗ std::vector<double> path_mono_splits(Path const &p) {
411 ✗ std::vector<double> ret;
412 ✗ if(p.empty()) return ret;
413
414 ✗ int pdx = 2, pdy = 2; // Previous derivative direction
415 ✗ for(unsigned i = 0; i < p.size(); i++) {
416 ✗ std::vector<double> spl = offset_doubles(curve_mono_splits(p[i]), i);
417 ✗ int dx = p[i].initialPoint()[X] > (spl.empty() ? p[i].finalPoint()[X] : p.valueAt(spl.front(), X)) ? 1 : 0;
418 ✗ int dy = p[i].initialPoint()[Y] > (spl.empty() ? p[i].finalPoint()[Y] : p.valueAt(spl.front(), Y)) ? 1 : 0;
419 //The direction changed, include the split time
420 ✗ if(dx != pdx || dy != pdy) {
421 ✗ ret.push_back(i);
422 ✗ pdx = dx; pdy = dy;
423 }
424 ✗ append(ret, spl);
425 ✗ }
426 ✗ return ret;
427 ✗ }
428
429 /**
430 * Applies path_mono_splits to multiple paths, and returns the results such that
431 * time-set i corresponds to Path i.
432 */
433 ✗ std::vector<std::vector<double> > paths_mono_splits(PathVector const &ps) {
434 ✗ std::vector<std::vector<double> > ret;
435 ✗ for(const auto & p : ps)
436 ✗ ret.push_back(path_mono_splits(p));
437 ✗ return ret;
438 ✗ }
439
440 /**
441 * Processes the bounds for a list of paths and a list of splits on them, yielding a list of rects for each.
442 * Each entry i corresponds to path i of the input. The number of rects in each entry is guaranteed to be the
443 * number of splits for that path, subtracted by one.
444 */
445 ✗ std::vector<std::vector<Rect> > split_bounds(PathVector const &p, std::vector<std::vector<double> > splits) {
446 ✗ std::vector<std::vector<Rect> > ret;
447 ✗ for(unsigned i = 0; i < p.size(); i++) {
448 ✗ std::vector<Rect> res;
449 ✗ for(unsigned j = 1; j < splits[i].size(); j++)
450 ✗ res.emplace_back(p[i].pointAt(splits[i][j-1]), p[i].pointAt(splits[i][j]));
451 ✗ ret.push_back(res);
452 ✗ }
453 ✗ return ret;
454 ✗ }
455
456 /**
457 * This is the main routine of "MonoCrosser", and implements a monotonic strategy on multiple curves.
458 * Finds crossings between two sets of paths, yielding a CrossingSet. [0, a.size()) of the return correspond
459 * to the sorted crossings of a with paths of b. The rest of the return, [a.size(), a.size() + b.size()],
460 * corresponds to the sorted crossings of b with paths of a.
461 *
462 * This function does two sweeps, one on the bounds of each path, and after that cull, one on the curves within.
463 * This leads to a certain amount of code complexity, however, most of that is factored into the above functions
464 */
465 ✗ CrossingSet MonoCrosser::crossings(PathVector const &a, PathVector const &b) {
466 ✗ if(b.empty()) return CrossingSet(a.size(), Crossings());
467 ✗ CrossingSet results(a.size() + b.size(), Crossings());
468 ✗ if(a.empty()) return results;
469
470 ✗ std::vector<std::vector<double> > splits_a = paths_mono_splits(a), splits_b = paths_mono_splits(b);
471 ✗ std::vector<std::vector<Rect> > bounds_a = split_bounds(a, splits_a), bounds_b = split_bounds(b, splits_b);
472
473 ✗ std::vector<Rect> bounds_a_union, bounds_b_union;
474 ✗ for(auto & i : bounds_a) bounds_a_union.push_back(union_list(i));
475 ✗ for(auto & i : bounds_b) bounds_b_union.push_back(union_list(i));
476
477 ✗ std::vector<std::vector<unsigned> > cull = sweep_bounds(bounds_a_union, bounds_b_union);
478 ✗ Crossings n;
479 ✗ for(unsigned i = 0; i < cull.size(); i++) {
480 ✗ for(unsigned jx = 0; jx < cull[i].size(); jx++) {
481 ✗ unsigned j = cull[i][jx];
482 ✗ unsigned jc = j + a.size();
483 ✗ Crossings res;
484
485 //Sweep of the monotonic portions
486 ✗ std::vector<std::vector<unsigned> > cull2 = sweep_bounds(bounds_a[i], bounds_b[j]);
487 ✗ for(unsigned k = 0; k < cull2.size(); k++) {
488 ✗ for(unsigned lx = 0; lx < cull2[k].size(); lx++) {
489 ✗ unsigned l = cull2[k][lx];
490 ✗ mono_pair(a[i], splits_a[i][k-1], splits_a[i][k],
491 ✗ b[j], splits_b[j][l-1], splits_b[j][l],
492 res, .1);
493 }
494 }
495
496 ✗ for(auto & re : res) { re.a = i; re.b = jc; }
497
498 ✗ merge_crossings(results[i], res, i);
499 ✗ merge_crossings(results[i], res, jc);
500 ✗ }
501 }
502
503 ✗ return results;
504 ✗ }
505
506 /* This function is similar codewise to the MonoCrosser, the main difference is that it deals with
507 * only one set of paths and includes self intersection
508 CrossingSet crossings_among(PathVector const &p) {
509 CrossingSet results(p.size(), Crossings());
510 if(p.empty()) return results;
511
512 std::vector<std::vector<double> > splits = paths_mono_splits(p);
513 std::vector<std::vector<Rect> > prs = split_bounds(p, splits);
514 std::vector<Rect> rs;
515 for(unsigned i = 0; i < prs.size(); i++) rs.push_back(union_list(prs[i]));
516
517 std::vector<std::vector<unsigned> > cull = sweep_bounds(rs);
518
519 //we actually want to do the self-intersections, so add em in:
520 for(unsigned i = 0; i < cull.size(); i++) cull[i].push_back(i);
521
522 for(unsigned i = 0; i < cull.size(); i++) {
523 for(unsigned jx = 0; jx < cull[i].size(); jx++) {
524 unsigned j = cull[i][jx];
525 Crossings res;
526
527 //Sweep of the monotonic portions
528 std::vector<std::vector<unsigned> > cull2 = sweep_bounds(prs[i], prs[j]);
529 for(unsigned k = 0; k < cull2.size(); k++) {
530 for(unsigned lx = 0; lx < cull2[k].size(); lx++) {
531 unsigned l = cull2[k][lx];
532 mono_pair(p[i], splits[i][k-1], splits[i][k],
533 p[j], splits[j][l-1], splits[j][l],
534 res, .1);
535 }
536 }
537
538 for(unsigned k = 0; k < res.size(); k++) { res[k].a = i; res[k].b = j; }
539
540 merge_crossings(results[i], res, i);
541 merge_crossings(results[j], res, j);
542 }
543 }
544
545 return results;
546 }
547 */
548
549
550 ✗ Crossings curve_self_crossings(Curve const &a) {
551 ✗ Crossings res;
552 ✗ std::vector<double> spl;
553 ✗ spl.push_back(0);
554 ✗ append(spl, curve_mono_splits(a));
555 ✗ spl.push_back(1);
556 ✗ for(unsigned i = 1; i < spl.size(); i++)
557 ✗ for(unsigned j = i+1; j < spl.size(); j++)
558 ✗ pair_intersect(a, spl[i-1], spl[i], a, spl[j-1], spl[j], res);
559 ✗ return res;
560 ✗ }
561
562 /*
563 void mono_curve_intersect(Curve const & A, double Al, double Ah,
564 Curve const & B, double Bl, double Bh,
565 Crossings &ret, unsigned depth=0) {
566 // std::cout << depth << "(" << Al << ", " << Ah << ")\n";
567 Point A0 = A.pointAt(Al), A1 = A.pointAt(Ah),
568 B0 = B.pointAt(Bl), B1 = B.pointAt(Bh);
569 //inline code that this implies? (without rect/interval construction)
570 if(!Rect(A0, A1).intersects(Rect(B0, B1)) || A0 == A1 || B0 == B1) return;
571
572 //Checks the general linearity of the function
573 if((depth > 12) || (A.boundsLocal(Interval(Al, Ah), 1).maxExtent() < 0.1
574 && B.boundsLocal(Interval(Bl, Bh), 1).maxExtent() < 0.1)) {
575 double tA, tB, c;
576 if(linear_intersect(A0, A1, B0, B1, tA, tB, c)) {
577 tA = tA * (Ah - Al) + Al;
578 tB = tB * (Bh - Bl) + Bl;
579 if(depth % 2)
580 ret.push_back(Crossing(tB, tA, c < 0));
581 else
582 ret.push_back(Crossing(tA, tB, c > 0));
583 return;
584 }
585 }
586 if(depth > 12) return;
587 double mid = (Bl + Bh)/2;
588 mono_curve_intersect(B, Bl, mid,
589 A, Al, Ah,
590 ret, depth+1);
591 mono_curve_intersect(B, mid, Bh,
592 A, Al, Ah,
593 ret, depth+1);
594 }
595
596 std::vector<std::vector<double> > curves_mono_splits(Path const &p) {
597 std::vector<std::vector<double> > ret;
598 for(unsigned i = 0; i <= p.size(); i++) {
599 std::vector<double> spl;
600 spl.push_back(0);
601 append(spl, curve_mono_splits(p[i]));
602 spl.push_back(1);
603 ret.push_back(spl);
604 }
605 }
606
607 std::vector<std::vector<Rect> > curves_split_bounds(Path const &p, std::vector<std::vector<double> > splits) {
608 std::vector<std::vector<Rect> > ret;
609 for(unsigned i = 0; i < splits.size(); i++) {
610 std::vector<Rect> res;
611 for(unsigned j = 1; j < splits[i].size(); j++)
612 res.push_back(Rect(p.pointAt(splits[i][j-1]+i), p.pointAt(splits[i][j]+i)));
613 ret.push_back(res);
614 }
615 return ret;
616 }
617
618 Crossings path_self_crossings(Path const &p) {
619 Crossings ret;
620 std::vector<std::vector<unsigned> > cull = sweep_bounds(bounds(p));
621 std::vector<std::vector<double> > spl = curves_mono_splits(p);
622 std::vector<std::vector<Rect> > bnds = curves_split_bounds(p, spl);
623 for(unsigned i = 0; i < cull.size(); i++) {
624 Crossings res;
625 for(unsigned k = 1; k < spl[i].size(); k++)
626 for(unsigned l = k+1; l < spl[i].size(); l++)
627 mono_curve_intersect(p[i], spl[i][k-1], spl[i][k], p[i], spl[i][l-1], spl[i][l], res);
628 offset_crossings(res, i, i);
629 append(ret, res);
630 for(unsigned jx = 0; jx < cull[i].size(); jx++) {
631 unsigned j = cull[i][jx];
632 res.clear();
633
634 std::vector<std::vector<unsigned> > cull2 = sweep_bounds(bnds[i], bnds[j]);
635 for(unsigned k = 0; k < cull2.size(); k++) {
636 for(unsigned lx = 0; lx < cull2[k].size(); lx++) {
637 unsigned l = cull2[k][lx];
638 mono_curve_intersect(p[i], spl[i][k-1], spl[i][k], p[j], spl[j][l-1], spl[j][l], res);
639 }
640 }
641
642 //if(fabs(int(i)-j) == 1 || fabs(int(i)-j) == p.size()-1) {
643 Crossings res2;
644 for(unsigned k = 0; k < res.size(); k++) {
645 if(res[k].ta != 0 && res[k].ta != 1 && res[k].tb != 0 && res[k].tb != 1) {
646 res.push_back(res[k]);
647 }
648 }
649 res = res2;
650 //}
651 offset_crossings(res, i, j);
652 append(ret, res);
653 }
654 }
655 return ret;
656 }
657 */
658
659 ✗ Crossings self_crossings(Path const &p) {
660 ✗ Crossings ret;
661 ✗ std::vector<std::vector<unsigned> > cull = sweep_bounds(bounds(p));
662 ✗ for(unsigned i = 0; i < cull.size(); i++) {
663 ✗ Crossings res = curve_self_crossings(p[i]);
664 ✗ offset_crossings(res, i, i);
665 ✗ append(ret, res);
666 ✗ for(unsigned jx = 0; jx < cull[i].size(); jx++) {
667 ✗ unsigned j = cull[i][jx];
668 ✗ res.clear();
669 ✗ pair_intersect(p[i], 0, 1, p[j], 0, 1, res);
670
671 //if(fabs(int(i)-j) == 1 || fabs(int(i)-j) == p.size()-1) {
672 ✗ Crossings res2;
673 ✗ for(auto & re : res) {
674 ✗ if(re.ta != 0 && re.ta != 1 && re.tb != 0 && re.tb != 1) {
675 ✗ res2.push_back(re);
676 }
677 }
678 ✗ res = res2;
679 //}
680 ✗ offset_crossings(res, i, j);
681 ✗ append(ret, res);
682 ✗ }
683 ✗ }
684 ✗ return ret;
685 ✗ }
686
687 ✗ void flip_crossings(Crossings &crs) {
688 ✗ for(auto & cr : crs)
689 ✗ cr = Crossing(cr.tb, cr.ta, cr.b, cr.a, !cr.dir);
690 ✗ }
691
692 ✗ CrossingSet crossings_among(PathVector const &p) {
693 ✗ CrossingSet results(p.size(), Crossings());
694 ✗ if(p.empty()) return results;
695
696 ✗ SimpleCrosser cc;
697
698 ✗ std::vector<std::vector<unsigned> > cull = sweep_bounds(bounds(p));
699 ✗ for(unsigned i = 0; i < cull.size(); i++) {
700 ✗ Crossings res = self_crossings(p[i]);
701 ✗ for(auto & re : res) { re.a = re.b = i; }
702 ✗ merge_crossings(results[i], res, i);
703 ✗ flip_crossings(res);
704 ✗ merge_crossings(results[i], res, i);
705 ✗ for(unsigned jx = 0; jx < cull[i].size(); jx++) {
706 ✗ unsigned j = cull[i][jx];
707
708 ✗ Crossings res = cc.crossings(p[i], p[j]);
709 ✗ for(auto & re : res) { re.a = i; re.b = j; }
710 ✗ merge_crossings(results[i], res, i);
711 ✗ merge_crossings(results[j], res, j);
712 ✗ }
713 ✗ }
714 ✗ return results;
715 ✗ }
716
717 }
718
719 /*
720 Local Variables:
721 mode:c++
722 c-file-style:"stroustrup"
723 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
724 indent-tabs-mode:nil
725 fill-column:99
726 End:
727 */
728 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
729