1 #include "path-intersection.h"
3 #include "ord.h"
5 //for path_direction:
6 #include "sbasis-geometric.h"
8 namespace Geom {
10 /* This function computes the winding of the path, given a reference point.
11 * Positive values correspond to counter-clockwise in the mathematical coordinate system,
12 * and clockwise in screen coordinates. This particular implementation casts a ray in
13 * the positive x direction. It iterates the path, checking for intersection with the
14 * bounding boxes. If an intersection is found, the initial/final Y value of the curve is
15 * used to derive a delta on the winding value. If the point is within the bounding box,
16 * the curve specific winding function is called.
17 */
18 int winding(Path const &path, Point p) {
19 //start on a segment which is not a horizontal line with y = p[y]
20 Path::const_iterator start;
21 for(Path::const_iterator iter = path.begin(); ; ++iter) {
22 if(iter == path.end_closed()) { return 0; }
23 if(iter->initialPoint()[Y]!=p[Y]) { start = iter; break; }
24 if(iter->finalPoint()[Y]!=p[Y]) { start = iter; break; }
25 if(iter->boundsFast().height()!=0.){ start = iter; break; }
26 }
27 int wind = 0;
28 int cnt = 0;
29 bool starting = true;
30 for (Path::const_iterator iter = start; iter != start || starting
31 ; ++iter, iter = (iter == path.end_closed()) ? path.begin() : iter )
32 {
33 cnt++;
34 if(cnt > path.size()) return wind; //some bug makes this required
35 starting = false;
36 Rect bounds = iter->boundsFast();
37 Coord x = p[X], y = p[Y];
39 if(x > bounds.right() || !bounds[Y].contains(y)) continue; //ray doesn't intersect box
41 Point final = iter->finalPoint();
42 Point initial = iter->initialPoint();
43 Cmp final_to_ray = cmp(final[Y], y);
44 Cmp initial_to_ray = cmp(initial[Y], y);
46 // if y is included, these will have opposite values, giving order.
47 Cmp c = cmp(final_to_ray, initial_to_ray);
48 if(x < bounds.left()) {
49 // ray goes through bbox
50 // winding delta determined by position of endpoints
51 if(final_to_ray != EQUAL_TO) {
52 wind += int(c); // GT = counter-clockwise = 1; LT = clockwise = -1; EQ = not-included = 0
53 //std::cout << int(c) << " ";
54 goto cont;
55 }
56 } else {
57 //inside bbox, use custom per-curve winding thingie
58 int delt = iter->winding(p);
59 wind += delt;
60 //std::cout << "n" << delt << " ";
61 }
62 //Handling the special case of an endpoint on the ray:
63 if(final[Y] == y) {
64 //Traverse segments until it breaks away from y
65 //99.9% of the time this will happen the first go
66 Path::const_iterator next = iter;
67 next++;
68 for(; ; next++) {
69 if(next == path.end_closed()) next = path.begin();
70 Rect bnds = next->boundsFast();
71 //TODO: X considerations
72 if(bnds.height() > 0) {
73 //It has diverged
74 if(bnds.contains(p)) {
75 const double fudge = 0.01;
76 if(cmp(y, next->valueAt(fudge, Y)) == initial_to_ray) {
77 wind += int(c);
78 std::cout << "!!!!!" << int(c) << " ";
79 }
80 iter = next; // No increment, as the rest of the thing hasn't been counted.
81 } else {
82 Coord ny = next->initialPoint()[Y];
83 if(cmp(y, ny) == initial_to_ray) {
84 //Is a continuation through the ray, so counts windingwise
85 wind += int(c);
86 std::cout << "!!!!!" << int(c) << " ";
87 }
88 iter = ++next;
89 }
90 goto cont;
91 }
92 if(next==start) return wind;
93 }
94 //Looks like it looped, which means everything's flat
95 return 0;
96 }
98 cont:(void)0;
99 }
100 return wind;
101 }
103 /* This function should only be applied to simple paths (regions), as otherwise
104 * a boolean winding direction is undefined. It returns true for fill, false for
105 * hole. Defaults to using the sign of area when it reaches funny cases.
106 */
107 bool path_direction(Path const &p) {
108 if(p.empty()) return false;
109 //could probably be more efficient, but this is a quick job
110 double y = p.initialPoint()[Y];
111 double x = p.initialPoint()[X];
112 Cmp res = cmp(p[0].finalPoint()[Y], y);
113 goto doh;
114 for(unsigned i = 1; i <= p.size(); i++) {
115 Cmp final_to_ray = cmp(p[i].finalPoint()[Y], y);
116 Cmp initial_to_ray = cmp(p[i].initialPoint()[Y], y);
117 // if y is included, these will have opposite values, giving order.
118 Cmp c = cmp(final_to_ray, initial_to_ray);
119 if(c != EQUAL_TO) {
120 std::vector<double> rs = p[i].roots(y, Y);
121 for(unsigned j = 0; j < rs.size(); j++) {
122 double nx = p[i].valueAt(rs[j], X);
123 if(nx > x) {
124 x = nx;
125 res = c;
126 }
127 }
128 } else if(final_to_ray == EQUAL_TO) goto doh;
129 }
130 return res < 0;
132 doh:
133 //Otherwise fallback on area
135 Piecewise<D2<SBasis> > pw = p.toPwSb();
136 double area;
137 Point centre;
138 Geom::centroid(pw, centre, area);
139 return area > 0;
140 }
142 //pair intersect code based on njh's pair-intersect
144 // A little sugar for appending a list to another
145 template<typename T>
146 void append(T &a, T const &b) {
147 a.insert(a.end(), b.begin(), b.end());
148 }
150 /* Finds the intersection between the lines defined by A0 & A1, and B0 & B1.
151 * Returns through the last 3 parameters, returning the t-values on the lines
152 * and the cross-product of the deltas (a useful byproduct). The return value
153 * indicates if the time values are within their proper range on the line segments.
154 */
155 bool
156 linear_intersect(Point A0, Point A1, Point B0, Point B1,
157 double &tA, double &tB, double &det) {
158 // kramers rule as cross products
159 Point Ad = A1 - A0,
160 Bd = B1 - B0,
161 d = B0 - A0;
162 det = cross(Ad, Bd);
163 if( 1.0 + det == 1.0 )
164 return false;
165 else
166 {
167 double detinv = 1.0 / det;
168 tA = cross(d, Bd) * detinv;
169 tB = cross(d, Ad) * detinv;
170 return tA >= 0. && tA <= 1. && tB >= 0. && tB <= 1.;
171 }
172 }
174 /* This uses the local bounds functions of curves to generically intersect two.
175 * It passes in the curves, time intervals, and keeps track of depth, while
176 * returning the results through the Crossings parameter.
177 */
178 void pair_intersect(Curve const & A, double Al, double Ah,
179 Curve const & B, double Bl, double Bh,
180 Crossings &ret, unsigned depth=0) {
181 // std::cout << depth << "(" << Al << ", " << Ah << ")\n";
182 Rect Ar = A.boundsLocal(Interval(Al, Ah));
183 if(Ar.isEmpty()) return;
185 Rect Br = B.boundsLocal(Interval(Bl, Bh));
186 if(Br.isEmpty()) return;
188 if(!Ar.intersects(Br)) return;
190 //Checks the general linearity of the function
191 if((depth > 12)) { // || (A.boundsLocal(Interval(Al, Ah), 1).maxExtent() < 0.1
192 //&& B.boundsLocal(Interval(Bl, Bh), 1).maxExtent() < 0.1)) {
193 double tA, tB, c;
194 if(linear_intersect(A.pointAt(Al), A.pointAt(Ah),
195 B.pointAt(Bl), B.pointAt(Bh),
196 tA, tB, c)) {
197 tA = tA * (Ah - Al) + Al;
198 tB = tB * (Bh - Bl) + Bl;
199 if(depth % 2)
200 ret.push_back(Crossing(tB, tA, c < 0));
201 else
202 ret.push_back(Crossing(tA, tB, c > 0));
203 return;
204 }
205 }
206 if(depth > 12) return;
207 double mid = (Bl + Bh)/2;
208 pair_intersect(B, Bl, mid,
209 A, Al, Ah,
210 ret, depth+1);
211 pair_intersect(B, mid, Bh,
212 A, Al, Ah,
213 ret, depth+1);
214 }
216 // A simple wrapper around pair_intersect
217 Crossings SimpleCrosser::crossings(Curve const &a, Curve const &b) {
218 Crossings ret;
219 pair_intersect(a, 0, 1, b, 0, 1, ret);
220 return ret;
221 }
223 /* Takes two paths and time ranges on them, with the invariant that the
224 * paths are monotonic on the range. Splits A when the linear intersection
225 * doesn't exist or is inaccurate. Uses the fact that it is monotonic to
226 * do very fast local bounds.
227 */
228 void mono_pair(Path const &A, double Al, double Ah,
229 Path const &B, double Bl, double Bh,
230 Crossings &ret, double tol, unsigned depth = 0) {
231 if( Al >= Ah || Bl >= Bh) return;
232 std::cout << " " << depth << "[" << Al << ", " << Ah << "]" << "[" << Bl << ", " << Bh << "]";
234 Point A0 = A.pointAt(Al), A1 = A.pointAt(Ah),
235 B0 = B.pointAt(Bl), B1 = B.pointAt(Bh);
236 //inline code that this implies? (without rect/interval construction)
237 if(!Rect(A0, A1).intersects(Rect(B0, B1)) || A0 == A1 || B0 == B1) return;
239 //Checks the general linearity of the function
240 //if((depth > 12) || (A.boundsLocal(Interval(Al, Ah), 1).maxExtent() < 0.1
241 // && B.boundsLocal(Interval(Bl, Bh), 1).maxExtent() < 0.1)) {
242 double tA, tB, c;
243 if(linear_intersect(A0, A1, B0, B1,
244 tA, tB, c)) {
245 tA = tA * (Ah - Al) + Al;
246 tB = tB * (Bh - Bl) + Bl;
247 if(depth % 2)
248 ret.push_back(Crossing(tB, tA, c < 0));
249 else
250 ret.push_back(Crossing(tA, tB, c > 0));
251 return;
252 }
253 //}
254 if(depth > 12) return;
255 double mid = (Bl + Bh)/2;
256 mono_pair(B, Bl, mid,
257 A, Al, Ah,
258 ret, depth+1);
259 mono_pair(B, mid, Bh,
260 A, Al, Ah,
261 ret, depth+1);
262 }
264 // This returns the times when the x or y derivative is 0 in the curve.
265 std::vector<double> curve_mono_splits(Curve const &d) {
266 std::vector<double> rs = d.roots(0, X);
267 append(rs, d.roots(0, Y));
268 std::sort(rs.begin(), rs.end());
269 return rs;
270 }
272 // Convenience function to add a value to each entry in a vector of doubles.
273 std::vector<double> offset_doubles(std::vector<double> const &x, double offs) {
274 std::vector<double> ret;
275 for(unsigned i = 0; i < x.size(); i++) {
276 ret.push_back(x[i] + offs);
277 }
278 return ret;
279 }
281 /* Finds all the monotonic splits for a path. Only includes the split between
282 * curves if they switch derivative directions at that point.
283 */
284 std::vector<double> path_mono_splits(Path const &p) {
285 std::vector<double> ret;
286 if(p.empty()) return ret;
287 ret.push_back(0);
289 Curve* deriv = p[0].derivative();
290 append(ret, curve_mono_splits(*deriv));
291 delete deriv;
293 bool pdx=2, pdy=2; //Previous derivative direction
294 for(unsigned i = 0; i <= p.size(); i++) {
295 deriv = p[i].derivative();
296 std::vector<double> spl = offset_doubles(curve_mono_splits(*deriv), i);
297 delete deriv;
298 bool dx = p[i].initialPoint()[X] > (spl.empty()? p[i].finalPoint()[X] :
299 p.valueAt(spl.front(), X));
300 bool dy = p[i].initialPoint()[Y] > (spl.empty()? p[i].finalPoint()[Y] :
301 p.valueAt(spl.front(), Y));
302 //The direction changed, include the split time
303 if(dx != pdx || dy != pdy) {
304 ret.push_back(i);
305 pdx = dx; pdy = dy;
306 }
307 append(ret, spl);
308 }
309 return ret;
310 }
312 /* Applies path_mono_splits to multiple paths, and returns the results such that
313 * time-set i corresponds to Path i.
314 */
315 std::vector<std::vector<double> > paths_mono_splits(std::vector<Path> const &ps) {
316 std::vector<std::vector<double> > ret;
317 for(unsigned i = 0; i < ps.size(); i++)
318 ret.push_back(path_mono_splits(ps[i]));
319 return ret;
320 }
322 /* Processes the bounds for a list of paths and a list of splits on them, yielding a list of rects for each.
323 * Each entry i corresponds to path i of the input. The number of rects in each entry is guaranteed to be the
324 * number of splits for that path, subtracted by one.
325 */
326 std::vector<std::vector<Rect> > split_bounds(std::vector<Path> const &p, std::vector<std::vector<double> > splits) {
327 std::vector<std::vector<Rect> > ret;
328 for(unsigned i = 0; i < p.size(); i++) {
329 std::vector<Rect> res;
330 for(unsigned j = 1; j < splits[i].size(); j++)
331 res.push_back(Rect(p[i].pointAt(splits[i][j-1]), p[i].pointAt(splits[i][j])));
332 ret.push_back(res);
333 }
334 return ret;
335 }
337 /* This is the main routine of "MonoCrosser", and implements a monotonic strategy on multiple curves.
338 * Finds crossings between two sets of paths, yielding a CrossingSet. [0, a.size()) of the return correspond
339 * to the sorted crossings of a with paths of b. The rest of the return, [a.size(), a.size() + b.size()],
340 * corresponds to the sorted crossings of b with paths of a.
341 *
342 * This function does two sweeps, one on the bounds of each path, and after that cull, one on the curves within.
343 * This leads to a certain amount of code complexity, however, most of that is factored into the above functions
344 */
345 CrossingSet MonoCrosser::crossings(std::vector<Path> const &a, std::vector<Path> const &b) {
346 if(b.empty()) return CrossingSet(a.size(), Crossings());
347 CrossingSet results(a.size() + b.size(), Crossings());
348 if(a.empty()) return results;
350 std::vector<std::vector<double> > splits_a = paths_mono_splits(a), splits_b = paths_mono_splits(b);
351 std::vector<std::vector<Rect> > bounds_a = split_bounds(a, splits_a), bounds_b = split_bounds(b, splits_b);
353 std::vector<Rect> bounds_a_union, bounds_b_union;
354 for(unsigned i = 0; i < bounds_a.size(); i++) bounds_a_union.push_back(union_list(bounds_a[i]));
355 for(unsigned i = 0; i < bounds_b.size(); i++) bounds_b_union.push_back(union_list(bounds_b[i]));
357 std::vector<std::vector<unsigned> > cull = sweep_bounds(bounds_a_union, bounds_b_union);
358 Crossings n;
359 for(unsigned i = 0; i < cull.size(); i++) {
360 for(unsigned jx = 0; jx < cull[i].size(); jx++) {
361 unsigned j = cull[i][jx];
362 unsigned jc = j + a.size();
363 Crossings res;
365 //Sweep of the monotonic portions
366 std::vector<std::vector<unsigned> > cull2 = sweep_bounds(bounds_a[i], bounds_b[j]);
367 for(unsigned k = 0; k < cull2.size(); k++) {
368 for(unsigned lx = 0; lx < cull2[k].size(); lx++) {
369 unsigned l = cull2[k][lx];
370 mono_pair(a[i], splits_a[i][k-1], splits_a[i][k],
371 b[j], splits_b[j][l-1], splits_b[j][l],
372 res, .1);
373 }
374 }
376 for(unsigned k = 0; k < res.size(); k++) { res[k].a = i; res[k].b = jc; }
378 merge_crossings(results[i], res, i);
379 merge_crossings(results[i], res, jc);
380 }
381 }
383 return results;
384 }
386 /* This function is similar codewise to the MonoCrosser, the main difference is that it deals with
387 * only one set of paths and includes self intersection
388 CrossingSet crossings_among(std::vector<Path> const &p) {
389 CrossingSet results(p.size(), Crossings());
390 if(p.empty()) return results;
392 std::vector<std::vector<double> > splits = paths_mono_splits(p);
393 std::vector<std::vector<Rect> > prs = split_bounds(p, splits);
394 std::vector<Rect> rs;
395 for(unsigned i = 0; i < prs.size(); i++) rs.push_back(union_list(prs[i]));
397 std::vector<std::vector<unsigned> > cull = sweep_bounds(rs);
399 //we actually want to do the self-intersections, so add em in:
400 for(unsigned i = 0; i < cull.size(); i++) cull[i].push_back(i);
402 for(unsigned i = 0; i < cull.size(); i++) {
403 for(unsigned jx = 0; jx < cull[i].size(); jx++) {
404 unsigned j = cull[i][jx];
405 Crossings res;
407 //Sweep of the monotonic portions
408 std::vector<std::vector<unsigned> > cull2 = sweep_bounds(prs[i], prs[j]);
409 for(unsigned k = 0; k < cull2.size(); k++) {
410 for(unsigned lx = 0; lx < cull2[k].size(); lx++) {
411 unsigned l = cull2[k][lx];
412 mono_pair(p[i], splits[i][k-1], splits[i][k],
413 p[j], splits[j][l-1], splits[j][l],
414 res, .1);
415 }
416 }
418 for(unsigned k = 0; k < res.size(); k++) { res[k].a = i; res[k].b = j; }
420 merge_crossings(results[i], res, i);
421 merge_crossings(results[j], res, j);
422 }
423 }
425 return results;
426 }
427 */
430 Crossings curve_self_crossings(Curve const &a) {
431 Crossings res;
432 std::vector<double> spl;
433 spl.push_back(0);
434 append(spl, curve_mono_splits(a));
435 spl.push_back(1);
436 for(unsigned i = 1; i < spl.size(); i++)
437 for(unsigned j = i+1; j < spl.size(); j++)
438 pair_intersect(a, spl[i-1], spl[i], a, spl[j-1], spl[j], res);
439 return res;
440 }
442 /*
443 void mono_curve_intersect(Curve const & A, double Al, double Ah,
444 Curve const & B, double Bl, double Bh,
445 Crossings &ret, unsigned depth=0) {
446 // std::cout << depth << "(" << Al << ", " << Ah << ")\n";
447 Point A0 = A.pointAt(Al), A1 = A.pointAt(Ah),
448 B0 = B.pointAt(Bl), B1 = B.pointAt(Bh);
449 //inline code that this implies? (without rect/interval construction)
450 if(!Rect(A0, A1).intersects(Rect(B0, B1)) || A0 == A1 || B0 == B1) return;
452 //Checks the general linearity of the function
453 if((depth > 12) || (A.boundsLocal(Interval(Al, Ah), 1).maxExtent() < 0.1
454 && B.boundsLocal(Interval(Bl, Bh), 1).maxExtent() < 0.1)) {
455 double tA, tB, c;
456 if(linear_intersect(A0, A1, B0, B1, tA, tB, c)) {
457 tA = tA * (Ah - Al) + Al;
458 tB = tB * (Bh - Bl) + Bl;
459 if(depth % 2)
460 ret.push_back(Crossing(tB, tA, c < 0));
461 else
462 ret.push_back(Crossing(tA, tB, c > 0));
463 return;
464 }
465 }
466 if(depth > 12) return;
467 double mid = (Bl + Bh)/2;
468 mono_curve_intersect(B, Bl, mid,
469 A, Al, Ah,
470 ret, depth+1);
471 mono_curve_intersect(B, mid, Bh,
472 A, Al, Ah,
473 ret, depth+1);
474 }
476 std::vector<std::vector<double> > curves_mono_splits(Path const &p) {
477 std::vector<std::vector<double> > ret;
478 for(unsigned i = 0; i <= p.size(); i++) {
479 std::vector<double> spl;
480 spl.push_back(0);
481 append(spl, curve_mono_splits(p[i]));
482 spl.push_back(1);
483 ret.push_back(spl);
484 }
485 }
487 std::vector<std::vector<Rect> > curves_split_bounds(Path const &p, std::vector<std::vector<double> > splits) {
488 std::vector<std::vector<Rect> > ret;
489 for(unsigned i = 0; i < splits.size(); i++) {
490 std::vector<Rect> res;
491 for(unsigned j = 1; j < splits[i].size(); j++)
492 res.push_back(Rect(p.pointAt(splits[i][j-1]+i), p.pointAt(splits[i][j]+i)));
493 ret.push_back(res);
494 }
495 return ret;
496 }
498 Crossings path_self_crossings(Path const &p) {
499 Crossings ret;
500 std::vector<std::vector<unsigned> > cull = sweep_bounds(bounds(p));
501 std::vector<std::vector<double> > spl = curves_mono_splits(p);
502 std::vector<std::vector<Rect> > bnds = curves_split_bounds(p, spl);
503 for(unsigned i = 0; i < cull.size(); i++) {
504 Crossings res;
505 for(unsigned k = 1; k < spl[i].size(); k++)
506 for(unsigned l = k+1; l < spl[i].size(); l++)
507 mono_curve_intersect(p[i], spl[i][k-1], spl[i][k], p[i], spl[i][l-1], spl[i][l], res);
508 offset_crossings(res, i, i);
509 append(ret, res);
510 for(unsigned jx = 0; jx < cull[i].size(); jx++) {
511 unsigned j = cull[i][jx];
512 res.clear();
514 std::vector<std::vector<unsigned> > cull2 = sweep_bounds(bnds[i], bnds[j]);
515 for(unsigned k = 0; k < cull2.size(); k++) {
516 for(unsigned lx = 0; lx < cull2[k].size(); lx++) {
517 unsigned l = cull2[k][lx];
518 mono_curve_intersect(p[i], spl[i][k-1], spl[i][k], p[j], spl[j][l-1], spl[j][l], res);
519 }
520 }
522 //if(fabs(int(i)-j) == 1 || fabs(int(i)-j) == p.size()-1) {
523 Crossings res2;
524 for(unsigned k = 0; k < res.size(); k++) {
525 if(res[k].ta != 0 && res[k].ta != 1 && res[k].tb != 0 && res[k].tb != 1) {
526 res.push_back(res[k]);
527 }
528 }
529 res = res2;
530 //}
531 offset_crossings(res, i, j);
532 append(ret, res);
533 }
534 }
535 return ret;
536 }
537 */
539 Crossings self_crossings(Path const &p) {
540 Crossings ret;
541 std::vector<std::vector<unsigned> > cull = sweep_bounds(bounds(p));
542 for(unsigned i = 0; i < cull.size(); i++) {
543 Crossings res = curve_self_crossings(p[i]);
544 offset_crossings(res, i, i);
545 append(ret, res);
546 for(unsigned jx = 0; jx < cull[i].size(); jx++) {
547 unsigned j = cull[i][jx];
548 res.clear();
549 pair_intersect(p[i], 0, 1, p[j], 0, 1, res);
551 //if(fabs(int(i)-j) == 1 || fabs(int(i)-j) == p.size()-1) {
552 Crossings res2;
553 for(unsigned k = 0; k < res.size(); k++) {
554 if(res[k].ta != 0 && res[k].ta != 1 && res[k].tb != 0 && res[k].tb != 1) {
555 res2.push_back(res[k]);
556 }
557 }
558 res = res2;
559 //}
560 offset_crossings(res, i, j);
561 append(ret, res);
562 }
563 }
564 return ret;
565 }
567 void flip_crossings(Crossings &crs) {
568 for(unsigned i = 0; i < crs.size(); i++)
569 crs[i] = Crossing(crs[i].tb, crs[i].ta, crs[i].b, crs[i].a, !crs[i].dir);
570 }
572 CrossingSet crossings_among(std::vector<Path> const &p) {
573 CrossingSet results(p.size(), Crossings());
574 if(p.empty()) return results;
576 SimpleCrosser cc;
578 std::vector<std::vector<unsigned> > cull = sweep_bounds(bounds(p));
579 for(unsigned i = 0; i < cull.size(); i++) {
580 Crossings res = self_crossings(p[i]);
581 for(unsigned k = 0; k < res.size(); k++) { res[k].a = res[k].b = i; }
582 merge_crossings(results[i], res, i);
583 flip_crossings(res);
584 merge_crossings(results[i], res, i);
585 for(unsigned jx = 0; jx < cull[i].size(); jx++) {
586 unsigned j = cull[i][jx];
588 Crossings res = cc.crossings(p[i], p[j]);
589 for(unsigned k = 0; k < res.size(); k++) { res[k].a = i; res[k].b = j; }
590 merge_crossings(results[i], res, i);
591 merge_crossings(results[j], res, j);
592 }
593 }
594 return results;
595 }
597 }