src/third_party/skia/experimental/Intersection/QuadraticImplicit.cpp

   1 // Another approach is to start with the implicit form of one curve and solve
   2 // (seek implicit coefficients in QuadraticParameter.cpp
   3 // by substituting in the parametric form of the other.
   4 // The downside of this approach is that early rejects are difficult to come by.
   5 // http://planetmath.org/encyclopedia/GaloisTheoreticDerivationOfTheQuarticFormula.html#step
   6
   7
   8 #include "CubicUtilities.h"
   9 #include "CurveIntersection.h"
  10 #include "Intersections.h"
  11 #include "QuadraticParameterization.h"
  12 #include "QuarticRoot.h"
  13 #include "QuadraticUtilities.h"
  14 #include "TSearch.h"
  15
  16 #ifdef SK_DEBUG
  17 #include "LineUtilities.h"
  18 #endif
  19
  20 /* given the implicit form 0 = Ax^2 + Bxy + Cy^2 + Dx + Ey + F
  21  * and given x = at^2 + bt + c  (the parameterized form)
  22  *           y = dt^2 + et + f
  23  * then
  24  * 0 = A(at^2+bt+c)(at^2+bt+c)+B(at^2+bt+c)(dt^2+et+f)+C(dt^2+et+f)(dt^2+et+f)+D(at^2+bt+c)+E(dt^2+et+f)+F
  25  */
  26
  27 static int findRoots(const QuadImplicitForm& i, const Quadratic& q2, double roots[4],
  28         bool oneHint, int firstCubicRoot) {
  29     double a, b, c;
  30     set_abc(&q2[0].x, a, b, c);
  31     double d, e, f;
  32     set_abc(&q2[0].y, d, e, f);
  33     const double t4 =     i.x2() *  a * a
  34                     +     i.xy() *  a * d
  35                     +     i.y2() *  d * d;
  36     const double t3 = 2 * i.x2() *  a * b
  37                     +     i.xy() * (a * e +     b * d)
  38                     + 2 * i.y2() *  d * e;
  39     const double t2 =     i.x2() * (b * b + 2 * a * c)
  40                     +     i.xy() * (c * d +     b * e + a * f)
  41                     +     i.y2() * (e * e + 2 * d * f)
  42                     +     i.x()  *  a
  43                     +     i.y()  *  d;
  44     const double t1 = 2 * i.x2() *  b * c
  45                     +     i.xy() * (c * e + b * f)
  46                     + 2 * i.y2() *  e * f
  47                     +     i.x()  *  b
  48                     +     i.y()  *  e;
  49     const double t0 =     i.x2() *  c * c
  50                     +     i.xy() *  c * f
  51                     +     i.y2() *  f * f
  52                     +     i.x()  *  c
  53                     +     i.y()  *  f
  54                     +     i.c();
  55     int rootCount = reducedQuarticRoots(t4, t3, t2, t1, t0, oneHint, roots);
  56     if (rootCount >= 0) {
  57         return rootCount;
  58     }
  59     return quarticRootsReal(firstCubicRoot, t4, t3, t2, t1, t0, roots);
  60 }
  61
  62 static int addValidRoots(const double roots[4], const int count, double valid[4]) {
  63     int result = 0;
  64     int index;
  65     for (index = 0; index < count; ++index) {
  66         if (!approximately_zero_or_more(roots[index]) || !approximately_one_or_less(roots[index])) {
  67             continue;
  68         }
  69         double t = 1 - roots[index];
  70         if (approximately_less_than_zero(t)) {
  71             t = 0;
  72         } else if (approximately_greater_than_one(t)) {
  73             t = 1;
  74         }
  75         valid[result++] = t;
  76     }
  77     return result;
  78 }
  79
  80 static bool onlyEndPtsInCommon(const Quadratic& q1, const Quadratic& q2, Intersections& i) {
  81 // the idea here is to see at minimum do a quick reject by rotating all points
  82 // to either side of the line formed by connecting the endpoints
  83 // if the opposite curves points are on the line or on the other side, the
  84 // curves at most intersect at the endpoints
  85     for (int oddMan = 0; oddMan < 3; ++oddMan) {
  86         const _Point* endPt[2];
  87         for (int opp = 1; opp < 3; ++opp) {
  88             int end = oddMan ^ opp;
  89             if (end == 3) {
  90                 end = opp;
  91             }
  92             endPt[opp - 1] = &q1[end];
  93         }
  94         double origX = endPt[0]->x;
  95         double origY = endPt[0]->y;
  96         double adj = endPt[1]->x - origX;
  97         double opp = endPt[1]->y - origY;
  98         double sign = (q1[oddMan].y - origY) * adj - (q1[oddMan].x - origX) * opp;
  99         if (approximately_zero(sign)) {
 100             goto tryNextHalfPlane;
 101         }
 102         for (int n = 0; n < 3; ++n) {
 103             double test = (q2[n].y - origY) * adj - (q2[n].x - origX) * opp;
 104             if (test * sign > 0) {
 105                 goto tryNextHalfPlane;
 106             }
 107         }
 108         for (int i1 = 0; i1 < 3; i1 += 2) {
 109             for (int i2 = 0; i2 < 3; i2 += 2) {
 110                 if (q1[i1] == q2[i2]) {
 111                     i.insert(i1 >> 1, i2 >> 1, q1[i1]);
 112                 }
 113             }
 114         }
 115         SkASSERT(i.fUsed < 3);
 116         return true;
 117 tryNextHalfPlane:
 118         ;
 119     }
 120     return false;
 121 }
 122
 123 // returns false if there's more than one intercept or the intercept doesn't match the point
 124 // returns true if the intercept was successfully added or if the
 125 // original quads need to be subdivided
 126 static bool addIntercept(const Quadratic& q1, const Quadratic& q2, double tMin, double tMax,
 127         Intersections& i, bool* subDivide) {
 128     double tMid = (tMin + tMax) / 2;
 129     _Point mid;
 130     xy_at_t(q2, tMid, mid.x, mid.y);
 131     _Line line;
 132     line[0] = line[1] = mid;
 133     _Vector dxdy = dxdy_at_t(q2, tMid);
 134     line[0] -= dxdy;
 135     line[1] += dxdy;
 136     Intersections rootTs;
 137     int roots = intersect(q1, line, rootTs);
 138     if (roots == 0) {
 139         if (subDivide) {
 140             *subDivide = true;
 141         }
 142         return true;
 143     }
 144     if (roots == 2) {
 145         return false;
 146     }
 147     _Point pt2;
 148     xy_at_t(q1, rootTs.fT[0][0], pt2.x, pt2.y);
 149     if (!pt2.approximatelyEqualHalf(mid)) {
 150         return false;
 151     }
 152     i.insertSwap(rootTs.fT[0][0], tMid, pt2);
 153     return true;
 154 }
 155
 156 static bool isLinearInner(const Quadratic& q1, double t1s, double t1e, const Quadratic& q2,
 157         double t2s, double t2e, Intersections& i, bool* subDivide) {
 158     Quadratic hull;
 159     sub_divide(q1, t1s, t1e, hull);
 160     _Line line = {hull[2], hull[0]};
 161     const _Line* testLines[] = { &line, (const _Line*) &hull[0], (const _Line*) &hull[1] };
 162     size_t testCount = sizeof(testLines) / sizeof(testLines[0]);
 163     SkTDArray<double> tsFound;
 164     for (size_t index = 0; index < testCount; ++index) {
 165         Intersections rootTs;
 166         int roots = intersect(q2, *testLines[index], rootTs);
 167         for (int idx2 = 0; idx2 < roots; ++idx2) {
 168             double t = rootTs.fT[0][idx2];
 169 #ifdef SK_DEBUG
 170         _Point qPt, lPt;
 171         xy_at_t(q2, t, qPt.x, qPt.y);
 172         xy_at_t(*testLines[index], rootTs.fT[1][idx2], lPt.x, lPt.y);
 173         SkASSERT(qPt.approximatelyEqual(lPt));
 174 #endif
 175             if (approximately_negative(t - t2s) || approximately_positive(t - t2e)) {
 176                 continue;
 177             }
 178             *tsFound.append() = rootTs.fT[0][idx2];
 179         }
 180     }
 181     int tCount = tsFound.count();
 182     if (!tCount) {
 183         return true;
 184     }
 185     double tMin, tMax;
 186     if (tCount == 1) {
 187         tMin = tMax = tsFound[0];
 188     } else if (tCount > 1) {
 189         QSort<double>(tsFound.begin(), tsFound.end() - 1);
 190         tMin = tsFound[0];
 191         tMax = tsFound[tsFound.count() - 1];
 192     }
 193     _Point end;
 194     xy_at_t(q2, t2s, end.x, end.y);
 195     bool startInTriangle = point_in_hull(hull, end);
 196     if (startInTriangle) {
 197         tMin = t2s;
 198     }
 199     xy_at_t(q2, t2e, end.x, end.y);
 200     bool endInTriangle = point_in_hull(hull, end);
 201     if (endInTriangle) {
 202         tMax = t2e;
 203     }
 204     int split = 0;
 205     _Vector dxy1, dxy2;
 206     if (tMin != tMax || tCount > 2) {
 207         dxy2 = dxdy_at_t(q2, tMin);
 208         for (int index = 1; index < tCount; ++index) {
 209             dxy1 = dxy2;
 210             dxy2 = dxdy_at_t(q2, tsFound[index]);
 211             double dot = dxy1.dot(dxy2);
 212             if (dot < 0) {
 213                 split = index - 1;
 214                 break;
 215             }
 216         }
 217
 218     }
 219     if (split == 0) { // there's one point
 220         if (addIntercept(q1, q2, tMin, tMax, i, subDivide)) {
 221             return true;
 222         }
 223         i.swap();
 224         return isLinearInner(q2, tMin, tMax, q1, t1s, t1e, i, subDivide);
 225     }
 226     // At this point, we have two ranges of t values -- treat each separately at the split
 227     bool result;
 228     if (addIntercept(q1, q2, tMin, tsFound[split - 1], i, subDivide)) {
 229         result = true;
 230     } else {
 231         i.swap();
 232         result = isLinearInner(q2, tMin, tsFound[split - 1], q1, t1s, t1e, i, subDivide);
 233     }
 234     if (addIntercept(q1, q2, tsFound[split], tMax, i, subDivide)) {
 235         result = true;
 236     } else {
 237         i.swap();
 238         result |= isLinearInner(q2, tsFound[split], tMax, q1, t1s, t1e, i, subDivide);
 239     }
 240     return result;
 241 }
 242
 243 static double flatMeasure(const Quadratic& q) {
 244     _Vector mid = q[1] - q[0];
 245     _Vector dxy = q[2] - q[0];
 246     double length = dxy.length(); // OPTIMIZE: get rid of sqrt
 247     return fabs(mid.cross(dxy) / length);
 248 }
 249
 250 // FIXME ? should this measure both and then use the quad that is the flattest as the line?
 251 static bool isLinear(const Quadratic& q1, const Quadratic& q2, Intersections& i) {
 252     double measure = flatMeasure(q1);
 253     // OPTIMIZE: (get rid of sqrt) use approximately_zero
 254     if (!approximately_zero_sqrt(measure)) {
 255         return false;
 256     }
 257     return isLinearInner(q1, 0, 1, q2, 0, 1, i, NULL);
 258 }
 259
 260 // FIXME: if flat measure is sufficiently large, then probably the quartic solution failed
 261 static void relaxedIsLinear(const Quadratic& q1, const Quadratic& q2, Intersections& i) {
 262     double m1 = flatMeasure(q1);
 263     double m2 = flatMeasure(q2);
 264 #ifdef SK_DEBUG
 265     double min = SkTMin(m1, m2);
 266     if (min > 5) {
 267         SkDebugf("%s maybe not flat enough.. %1.9g\n", __FUNCTION__, min);
 268     }
 269 #endif
 270     i.reset();
 271     const Quadratic& rounder = m2 < m1 ? q1 : q2;
 272     const Quadratic& flatter = m2 < m1 ? q2 : q1;
 273     bool subDivide = false;
 274     isLinearInner(flatter, 0, 1, rounder, 0, 1, i, &subDivide);
 275     if (subDivide) {
 276         QuadraticPair pair;
 277         chop_at(flatter, pair, 0.5);
 278         Intersections firstI, secondI;
 279         relaxedIsLinear(pair.first(), rounder, firstI);
 280         for (int index = 0; index < firstI.used(); ++index) {
 281             i.insert(firstI.fT[0][index] * 0.5, firstI.fT[1][index], firstI.fPt[index]);
 282         }
 283         relaxedIsLinear(pair.second(), rounder, secondI);
 284         for (int index = 0; index < secondI.used(); ++index) {
 285             i.insert(0.5 + secondI.fT[0][index] * 0.5, secondI.fT[1][index], secondI.fPt[index]);
 286         }
 287     }
 288     if (m2 < m1) {
 289         i.swapPts();
 290     }
 291 }
 292
 293 #if 0
 294 static void unsortableExpanse(const Quadratic& q1, const Quadratic& q2, Intersections& i) {
 295     const Quadratic* qs[2] = { &q1, &q2 };
 296     // need t values for start and end of unsortable expanse on both curves
 297     // try projecting lines parallel to the end points
 298     i.fT[0][0] = 0;
 299     i.fT[0][1] = 1;
 300     int flip = -1; // undecided
 301     for (int qIdx = 0; qIdx < 2; qIdx++) {
 302         for (int t = 0; t < 2; t++) {
 303             _Point dxdy;
 304             dxdy_at_t(*qs[qIdx], t, dxdy);
 305             _Line perp;
 306             perp[0] = perp[1] = (*qs[qIdx])[t == 0 ? 0 : 2];
 307             perp[0].x += dxdy.y;
 308             perp[0].y -= dxdy.x;
 309             perp[1].x -= dxdy.y;
 310             perp[1].y += dxdy.x;
 311             Intersections hitData;
 312             int hits = intersectRay(*qs[qIdx ^ 1], perp, hitData);
 313             SkASSERT(hits <= 1);
 314             if (hits) {
 315                 if (flip < 0) {
 316                     _Point dxdy2;
 317                     dxdy_at_t(*qs[qIdx ^ 1], hitData.fT[0][0], dxdy2);
 318                     double dot = dxdy.dot(dxdy2);
 319                     flip = dot < 0;
 320                     i.fT[1][0] = flip;
 321                     i.fT[1][1] = !flip;
 322                 }
 323                 i.fT[qIdx ^ 1][t ^ flip] = hitData.fT[0][0];
 324             }
 325         }
 326     }
 327     i.fUnsortable = true; // failed, probably coincident or near-coincident
 328     i.fUsed = 2;
 329 }
 330 #endif
 331
 332 // each time through the loop, this computes values it had from the last loop
 333 // if i == j == 1, the center values are still good
 334 // otherwise, for i != 1 or j != 1, four of the values are still good
 335 // and if i == 1 ^ j == 1, an additional value is good
 336 static bool binarySearch(const Quadratic& quad1, const Quadratic& quad2, double& t1Seed,
 337         double& t2Seed, _Point& pt) {
 338     double tStep = ROUGH_EPSILON;
 339     _Point t1[3], t2[3];
 340     int calcMask = ~0;
 341     do {
 342         if (calcMask & (1 << 1)) t1[1] = xy_at_t(quad1, t1Seed);
 343         if (calcMask & (1 << 4)) t2[1] = xy_at_t(quad2, t2Seed);
 344         if (t1[1].approximatelyEqual(t2[1])) {
 345             pt = t1[1];
 346     #if ONE_OFF_DEBUG
 347             SkDebugf("%s t1=%1.9g t2=%1.9g (%1.9g,%1.9g) == (%1.9g,%1.9g)\n", __FUNCTION__,
 348                     t1Seed, t2Seed, t1[1].x, t1[1].y, t1[2].x, t1[2].y);
 349     #endif
 350             return true;
 351         }
 352         if (calcMask & (1 << 0)) t1[0] = xy_at_t(quad1, t1Seed - tStep);
 353         if (calcMask & (1 << 2)) t1[2] = xy_at_t(quad1, t1Seed + tStep);
 354         if (calcMask & (1 << 3)) t2[0] = xy_at_t(quad2, t2Seed - tStep);
 355         if (calcMask & (1 << 5)) t2[2] = xy_at_t(quad2, t2Seed + tStep);
 356         double dist[3][3];
 357         // OPTIMIZE: using calcMask value permits skipping some distance calcuations
 358         //   if prior loop's results are moved to correct slot for reuse
 359         dist[1][1] = t1[1].distanceSquared(t2[1]);
 360         int best_i = 1, best_j = 1;
 361         for (int i = 0; i < 3; ++i) {
 362             for (int j = 0; j < 3; ++j) {
 363                 if (i == 1 && j == 1) {
 364                     continue;
 365                 }
 366                 dist[i][j] = t1[i].distanceSquared(t2[j]);
 367                 if (dist[best_i][best_j] > dist[i][j]) {
 368                     best_i = i;
 369                     best_j = j;
 370                 }
 371             }
 372         }
 373         if (best_i == 1 && best_j == 1) {
 374             tStep /= 2;
 375             if (tStep < FLT_EPSILON_HALF) {
 376                 break;
 377             }
 378             calcMask = (1 << 0) | (1 << 2) | (1 << 3) | (1 << 5);
 379             continue;
 380         }
 381         if (best_i == 0) {
 382             t1Seed -= tStep;
 383             t1[2] = t1[1];
 384             t1[1] = t1[0];
 385             calcMask = 1 << 0;
 386         } else if (best_i == 2) {
 387             t1Seed += tStep;
 388             t1[0] = t1[1];
 389             t1[1] = t1[2];
 390             calcMask = 1 << 2;
 391         } else {
 392             calcMask = 0;
 393         }
 394         if (best_j == 0) {
 395             t2Seed -= tStep;
 396             t2[2] = t2[1];
 397             t2[1] = t2[0];
 398             calcMask |= 1 << 3;
 399         } else if (best_j == 2) {
 400             t2Seed += tStep;
 401             t2[0] = t2[1];
 402             t2[1] = t2[2];
 403             calcMask |= 1 << 5;
 404         }
 405     } while (true);
 406 #if ONE_OFF_DEBUG
 407     SkDebugf("%s t1=%1.9g t2=%1.9g (%1.9g,%1.9g) != (%1.9g,%1.9g) %s\n", __FUNCTION__,
 408         t1Seed, t2Seed, t1[1].x, t1[1].y, t1[2].x, t1[2].y);
 409 #endif
 410     return false;
 411 }
 412
 413 bool intersect2(const Quadratic& q1, const Quadratic& q2, Intersections& i) {
 414     // if the quads share an end point, check to see if they overlap
 415
 416     if (onlyEndPtsInCommon(q1, q2, i)) {
 417         return i.intersected();
 418     }
 419     if (onlyEndPtsInCommon(q2, q1, i)) {
 420         i.swapPts();
 421         return i.intersected();
 422     }
 423     // see if either quad is really a line
 424     if (isLinear(q1, q2, i)) {
 425         return i.intersected();
 426     }
 427     if (isLinear(q2, q1, i)) {
 428         i.swapPts();
 429         return i.intersected();
 430     }
 431     QuadImplicitForm i1(q1);
 432     QuadImplicitForm i2(q2);
 433     if (i1.implicit_match(i2)) {
 434         // FIXME: compute T values
 435         // compute the intersections of the ends to find the coincident span
 436         bool useVertical = fabs(q1[0].x - q1[2].x) < fabs(q1[0].y - q1[2].y);
 437         double t;
 438         if ((t = axialIntersect(q1, q2[0], useVertical)) >= 0) {
 439             i.insertCoincident(t, 0, q2[0]);
 440         }
 441         if ((t = axialIntersect(q1, q2[2], useVertical)) >= 0) {
 442             i.insertCoincident(t, 1, q2[2]);
 443         }
 444         useVertical = fabs(q2[0].x - q2[2].x) < fabs(q2[0].y - q2[2].y);
 445         if ((t = axialIntersect(q2, q1[0], useVertical)) >= 0) {
 446             i.insertCoincident(0, t, q1[0]);
 447         }
 448         if ((t = axialIntersect(q2, q1[2], useVertical)) >= 0) {
 449             i.insertCoincident(1, t, q1[2]);
 450         }
 451         SkASSERT(i.coincidentUsed() <= 2);
 452         return i.coincidentUsed() > 0;
 453     }
 454     int index;
 455     bool useCubic = q1[0] == q2[0] || q1[0] == q2[2] || q1[2] == q2[0];
 456     double roots1[4];
 457     int rootCount = findRoots(i2, q1, roots1, useCubic, 0);
 458     // OPTIMIZATION: could short circuit here if all roots are < 0 or > 1
 459     double roots1Copy[4];
 460     int r1Count = addValidRoots(roots1, rootCount, roots1Copy);
 461     _Point pts1[4];
 462     for (index = 0; index < r1Count; ++index) {
 463         xy_at_t(q1, roots1Copy[index], pts1[index].x, pts1[index].y);
 464     }
 465     double roots2[4];
 466     int rootCount2 = findRoots(i1, q2, roots2, useCubic, 0);
 467     double roots2Copy[4];
 468     int r2Count = addValidRoots(roots2, rootCount2, roots2Copy);
 469     _Point pts2[4];
 470     for (index = 0; index < r2Count; ++index) {
 471         xy_at_t(q2, roots2Copy[index], pts2[index].x, pts2[index].y);
 472     }
 473     if (r1Count == r2Count && r1Count <= 1) {
 474         if (r1Count == 1) {
 475             if (pts1[0].approximatelyEqualHalf(pts2[0])) {
 476                 i.insert(roots1Copy[0], roots2Copy[0], pts1[0]);
 477             } else if (pts1[0].moreRoughlyEqual(pts2[0])) {
 478                 // experiment: see if a different cubic solution provides the correct quartic answer
 479             #if 0
 480                 for (int cu1 = 0; cu1 < 3; ++cu1) {
 481                     rootCount = findRoots(i2, q1, roots1, useCubic, cu1);
 482                     r1Count = addValidRoots(roots1, rootCount, roots1Copy);
 483                     if (r1Count == 0) {
 484                         continue;
 485                     }
 486                     for (int cu2 = 0; cu2 < 3; ++cu2) {
 487                         if (cu1 == 0 && cu2 == 0) {
 488                             continue;
 489                         }
 490                         rootCount2 = findRoots(i1, q2, roots2, useCubic, cu2);
 491                         r2Count = addValidRoots(roots2, rootCount2, roots2Copy);
 492                         if (r2Count == 0) {
 493                             continue;
 494                         }
 495                         SkASSERT(r1Count == 1 && r2Count == 1);
 496                         SkDebugf("*** [%d,%d] (%1.9g,%1.9g) %s (%1.9g,%1.9g)\n", cu1, cu2,
 497                                 pts1[0].x, pts1[0].y, pts1[0].approximatelyEqualHalf(pts2[0])
 498                                 ? "==" : "!=", pts2[0].x, pts2[0].y);
 499                     }
 500                 }
 501             #endif
 502                 // experiment: try to find intersection by chasing t
 503                 rootCount = findRoots(i2, q1, roots1, useCubic, 0);
 504                 r1Count = addValidRoots(roots1, rootCount, roots1Copy);
 505                 rootCount2 = findRoots(i1, q2, roots2, useCubic, 0);
 506                 r2Count = addValidRoots(roots2, rootCount2, roots2Copy);
 507                 if (binarySearch(q1, q2, roots1Copy[0], roots2Copy[0], pts1[0])) {
 508                     i.insert(roots1Copy[0], roots2Copy[0], pts1[0]);
 509                 }
 510             }
 511         }
 512         return i.intersected();
 513     }
 514     int closest[4];
 515     double dist[4];
 516     bool foundSomething = false;
 517     for (index = 0; index < r1Count; ++index) {
 518         dist[index] = DBL_MAX;
 519         closest[index] = -1;
 520         for (int ndex2 = 0; ndex2 < r2Count; ++ndex2) {
 521             if (!pts2[ndex2].approximatelyEqualHalf(pts1[index])) {
 522                 continue;
 523             }
 524             double dx = pts2[ndex2].x - pts1[index].x;
 525             double dy = pts2[ndex2].y - pts1[index].y;
 526             double distance = dx * dx + dy * dy;
 527             if (dist[index] <= distance) {
 528                 continue;
 529             }
 530             for (int outer = 0; outer < index; ++outer) {
 531                 if (closest[outer] != ndex2) {
 532                     continue;
 533                 }
 534                 if (dist[outer] < distance) {
 535                     goto next;
 536                 }
 537                 closest[outer] = -1;
 538             }
 539             dist[index] = distance;
 540             closest[index] = ndex2;
 541             foundSomething = true;
 542         next:
 543             ;
 544         }
 545     }
 546     if (r1Count && r2Count && !foundSomething) {
 547         relaxedIsLinear(q1, q2, i);
 548         return i.intersected();
 549     }
 550     int used = 0;
 551     do {
 552         double lowest = DBL_MAX;
 553         int lowestIndex = -1;
 554         for (index = 0; index < r1Count; ++index) {
 555             if (closest[index] < 0) {
 556                 continue;
 557             }
 558             if (roots1Copy[index] < lowest) {
 559                 lowestIndex = index;
 560                 lowest = roots1Copy[index];
 561             }
 562         }
 563         if (lowestIndex < 0) {
 564             break;
 565         }
 566         i.insert(roots1Copy[lowestIndex], roots2Copy[closest[lowestIndex]],
 567                 pts1[lowestIndex]);
 568         closest[lowestIndex] = -1;
 569     } while (++used < r1Count);
 570     i.fFlip = false;
 571     return i.intersected();
 572 }