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