src/2geom/path-intersection.cpp

   1 #include "path-intersection.h"
   2
   3 #include "ord.h"
   4
   5 //for path_direction:
   6 #include "sbasis-geometric.h"
   7
   8 namespace Geom {
   9
  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   unsigned 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];
  38
  39     if(x > bounds.right() || !bounds[Y].contains(y)) continue; //ray doesn't intersect box
  40
  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);
  45
  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     }
  97
  98     cont:(void)0;
  99   }
 100   return wind;
 101 }
 102
 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;
 131
 132     doh:
 133         //Otherwise fallback on area
 134
 135         Piecewise<D2<SBasis> > pw = p.toPwSb();
 136         double area;
 137         Point centre;
 138         Geom::centroid(pw, centre, area);
 139         return area > 0;
 140 }
 141
 142 //pair intersect code based on njh's pair-intersect
 143
 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 }
 149
 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 }
 173
 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;
 184
 185     Rect Br = B.boundsLocal(Interval(Bl, Bh));
 186     if(Br.isEmpty()) return;
 187
 188     if(!Ar.intersects(Br)) return;
 189
 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 }
 215
 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 }
 222
 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 << "]";
 233
 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;
 238
 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 }
 263
 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 }
 271
 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 }
 280
 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);
 288
 289     Curve* deriv = p[0].derivative();
 290     append(ret, curve_mono_splits(*deriv));
 291     delete deriv;
 292
 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 }
 311
 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 }
 321
 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 }
 336
 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;
 349
 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);
 352
 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]));
 356
 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;
 364
 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             }
 375
 376             for(unsigned k = 0; k < res.size(); k++) { res[k].a = i; res[k].b = jc; }
 377
 378             merge_crossings(results[i], res, i);
 379             merge_crossings(results[i], res, jc);
 380         }
 381     }
 382
 383     return results;
 384 }
 385
 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;
 391
 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]));
 396
 397     std::vector<std::vector<unsigned> > cull = sweep_bounds(rs);
 398
 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);
 401
 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;
 406
 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             }
 417
 418             for(unsigned k = 0; k < res.size(); k++) { res[k].a = i; res[k].b = j; }
 419
 420             merge_crossings(results[i], res, i);
 421             merge_crossings(results[j], res, j);
 422         }
 423     }
 424
 425     return results;
 426 }
 427 */
 428
 429
 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 }
 441
 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;
 451
 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 }
 475
 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 }
 486
 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 }
 497
 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();
 513
 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             }
 521
 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 */
 538
 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);
 550
 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 }
 566
 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 }
 571
 572 CrossingSet crossings_among(std::vector<Path> const &p) {
 573     CrossingSet results(p.size(), Crossings());
 574     if(p.empty()) return results;
 575
 576     SimpleCrosser cc;
 577
 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];
 587
 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 }
 596
 597 }