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