src/2geom/path-intersection.cpp

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