1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/gfx/skia/trunk/src/pathops/SkDQuadIntersection.cpp Wed Dec 31 06:09:35 2014 +0100
1.3 @@ -0,0 +1,553 @@
1.4 +// Another approach is to start with the implicit form of one curve and solve
1.5 +// (seek implicit coefficients in QuadraticParameter.cpp
1.6 +// by substituting in the parametric form of the other.
1.7 +// The downside of this approach is that early rejects are difficult to come by.
1.8 +// http://planetmath.org/encyclopedia/GaloisTheoreticDerivationOfTheQuarticFormula.html#step
1.9 +
1.10 +
1.11 +#include "SkDQuadImplicit.h"
1.12 +#include "SkIntersections.h"
1.13 +#include "SkPathOpsLine.h"
1.14 +#include "SkQuarticRoot.h"
1.15 +#include "SkTArray.h"
1.16 +#include "SkTSort.h"
1.17 +
1.18 +/* given the implicit form 0 = Ax^2 + Bxy + Cy^2 + Dx + Ey + F
1.19 + * and given x = at^2 + bt + c (the parameterized form)
1.20 + * y = dt^2 + et + f
1.21 + * then
1.22 + * 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
1.23 + */
1.24 +
1.25 +static int findRoots(const SkDQuadImplicit& i, const SkDQuad& quad, double roots[4],
1.26 + bool oneHint, bool flip, int firstCubicRoot) {
1.27 + SkDQuad flipped;
1.28 + const SkDQuad& q = flip ? (flipped = quad.flip()) : quad;
1.29 + double a, b, c;
1.30 + SkDQuad::SetABC(&q[0].fX, &a, &b, &c);
1.31 + double d, e, f;
1.32 + SkDQuad::SetABC(&q[0].fY, &d, &e, &f);
1.33 + const double t4 = i.x2() * a * a
1.34 + + i.xy() * a * d
1.35 + + i.y2() * d * d;
1.36 + const double t3 = 2 * i.x2() * a * b
1.37 + + i.xy() * (a * e + b * d)
1.38 + + 2 * i.y2() * d * e;
1.39 + const double t2 = i.x2() * (b * b + 2 * a * c)
1.40 + + i.xy() * (c * d + b * e + a * f)
1.41 + + i.y2() * (e * e + 2 * d * f)
1.42 + + i.x() * a
1.43 + + i.y() * d;
1.44 + const double t1 = 2 * i.x2() * b * c
1.45 + + i.xy() * (c * e + b * f)
1.46 + + 2 * i.y2() * e * f
1.47 + + i.x() * b
1.48 + + i.y() * e;
1.49 + const double t0 = i.x2() * c * c
1.50 + + i.xy() * c * f
1.51 + + i.y2() * f * f
1.52 + + i.x() * c
1.53 + + i.y() * f
1.54 + + i.c();
1.55 + int rootCount = SkReducedQuarticRoots(t4, t3, t2, t1, t0, oneHint, roots);
1.56 + if (rootCount < 0) {
1.57 + rootCount = SkQuarticRootsReal(firstCubicRoot, t4, t3, t2, t1, t0, roots);
1.58 + }
1.59 + if (flip) {
1.60 + for (int index = 0; index < rootCount; ++index) {
1.61 + roots[index] = 1 - roots[index];
1.62 + }
1.63 + }
1.64 + return rootCount;
1.65 +}
1.66 +
1.67 +static int addValidRoots(const double roots[4], const int count, double valid[4]) {
1.68 + int result = 0;
1.69 + int index;
1.70 + for (index = 0; index < count; ++index) {
1.71 + if (!approximately_zero_or_more(roots[index]) || !approximately_one_or_less(roots[index])) {
1.72 + continue;
1.73 + }
1.74 + double t = 1 - roots[index];
1.75 + if (approximately_less_than_zero(t)) {
1.76 + t = 0;
1.77 + } else if (approximately_greater_than_one(t)) {
1.78 + t = 1;
1.79 + }
1.80 + valid[result++] = t;
1.81 + }
1.82 + return result;
1.83 +}
1.84 +
1.85 +static bool only_end_pts_in_common(const SkDQuad& q1, const SkDQuad& q2) {
1.86 +// the idea here is to see at minimum do a quick reject by rotating all points
1.87 +// to either side of the line formed by connecting the endpoints
1.88 +// if the opposite curves points are on the line or on the other side, the
1.89 +// curves at most intersect at the endpoints
1.90 + for (int oddMan = 0; oddMan < 3; ++oddMan) {
1.91 + const SkDPoint* endPt[2];
1.92 + for (int opp = 1; opp < 3; ++opp) {
1.93 + int end = oddMan ^ opp; // choose a value not equal to oddMan
1.94 + if (3 == end) { // and correct so that largest value is 1 or 2
1.95 + end = opp;
1.96 + }
1.97 + endPt[opp - 1] = &q1[end];
1.98 + }
1.99 + double origX = endPt[0]->fX;
1.100 + double origY = endPt[0]->fY;
1.101 + double adj = endPt[1]->fX - origX;
1.102 + double opp = endPt[1]->fY - origY;
1.103 + double sign = (q1[oddMan].fY - origY) * adj - (q1[oddMan].fX - origX) * opp;
1.104 + if (approximately_zero(sign)) {
1.105 + goto tryNextHalfPlane;
1.106 + }
1.107 + for (int n = 0; n < 3; ++n) {
1.108 + double test = (q2[n].fY - origY) * adj - (q2[n].fX - origX) * opp;
1.109 + if (test * sign > 0 && !precisely_zero(test)) {
1.110 + goto tryNextHalfPlane;
1.111 + }
1.112 + }
1.113 + return true;
1.114 +tryNextHalfPlane:
1.115 + ;
1.116 + }
1.117 + return false;
1.118 +}
1.119 +
1.120 +// returns false if there's more than one intercept or the intercept doesn't match the point
1.121 +// returns true if the intercept was successfully added or if the
1.122 +// original quads need to be subdivided
1.123 +static bool add_intercept(const SkDQuad& q1, const SkDQuad& q2, double tMin, double tMax,
1.124 + SkIntersections* i, bool* subDivide) {
1.125 + double tMid = (tMin + tMax) / 2;
1.126 + SkDPoint mid = q2.ptAtT(tMid);
1.127 + SkDLine line;
1.128 + line[0] = line[1] = mid;
1.129 + SkDVector dxdy = q2.dxdyAtT(tMid);
1.130 + line[0] -= dxdy;
1.131 + line[1] += dxdy;
1.132 + SkIntersections rootTs;
1.133 + rootTs.allowNear(false);
1.134 + int roots = rootTs.intersect(q1, line);
1.135 + if (roots == 0) {
1.136 + if (subDivide) {
1.137 + *subDivide = true;
1.138 + }
1.139 + return true;
1.140 + }
1.141 + if (roots == 2) {
1.142 + return false;
1.143 + }
1.144 + SkDPoint pt2 = q1.ptAtT(rootTs[0][0]);
1.145 + if (!pt2.approximatelyEqual(mid)) {
1.146 + return false;
1.147 + }
1.148 + i->insertSwap(rootTs[0][0], tMid, pt2);
1.149 + return true;
1.150 +}
1.151 +
1.152 +static bool is_linear_inner(const SkDQuad& q1, double t1s, double t1e, const SkDQuad& q2,
1.153 + double t2s, double t2e, SkIntersections* i, bool* subDivide) {
1.154 + SkDQuad hull = q1.subDivide(t1s, t1e);
1.155 + SkDLine line = {{hull[2], hull[0]}};
1.156 + const SkDLine* testLines[] = { &line, (const SkDLine*) &hull[0], (const SkDLine*) &hull[1] };
1.157 + const size_t kTestCount = SK_ARRAY_COUNT(testLines);
1.158 + SkSTArray<kTestCount * 2, double, true> tsFound;
1.159 + for (size_t index = 0; index < kTestCount; ++index) {
1.160 + SkIntersections rootTs;
1.161 + rootTs.allowNear(false);
1.162 + int roots = rootTs.intersect(q2, *testLines[index]);
1.163 + for (int idx2 = 0; idx2 < roots; ++idx2) {
1.164 + double t = rootTs[0][idx2];
1.165 +#ifdef SK_DEBUG
1.166 + SkDPoint qPt = q2.ptAtT(t);
1.167 + SkDPoint lPt = testLines[index]->ptAtT(rootTs[1][idx2]);
1.168 + SkASSERT(qPt.approximatelyPEqual(lPt));
1.169 +#endif
1.170 + if (approximately_negative(t - t2s) || approximately_positive(t - t2e)) {
1.171 + continue;
1.172 + }
1.173 + tsFound.push_back(rootTs[0][idx2]);
1.174 + }
1.175 + }
1.176 + int tCount = tsFound.count();
1.177 + if (tCount <= 0) {
1.178 + return true;
1.179 + }
1.180 + double tMin, tMax;
1.181 + if (tCount == 1) {
1.182 + tMin = tMax = tsFound[0];
1.183 + } else {
1.184 + SkASSERT(tCount > 1);
1.185 + SkTQSort<double>(tsFound.begin(), tsFound.end() - 1);
1.186 + tMin = tsFound[0];
1.187 + tMax = tsFound[tsFound.count() - 1];
1.188 + }
1.189 + SkDPoint end = q2.ptAtT(t2s);
1.190 + bool startInTriangle = hull.pointInHull(end);
1.191 + if (startInTriangle) {
1.192 + tMin = t2s;
1.193 + }
1.194 + end = q2.ptAtT(t2e);
1.195 + bool endInTriangle = hull.pointInHull(end);
1.196 + if (endInTriangle) {
1.197 + tMax = t2e;
1.198 + }
1.199 + int split = 0;
1.200 + SkDVector dxy1, dxy2;
1.201 + if (tMin != tMax || tCount > 2) {
1.202 + dxy2 = q2.dxdyAtT(tMin);
1.203 + for (int index = 1; index < tCount; ++index) {
1.204 + dxy1 = dxy2;
1.205 + dxy2 = q2.dxdyAtT(tsFound[index]);
1.206 + double dot = dxy1.dot(dxy2);
1.207 + if (dot < 0) {
1.208 + split = index - 1;
1.209 + break;
1.210 + }
1.211 + }
1.212 + }
1.213 + if (split == 0) { // there's one point
1.214 + if (add_intercept(q1, q2, tMin, tMax, i, subDivide)) {
1.215 + return true;
1.216 + }
1.217 + i->swap();
1.218 + return is_linear_inner(q2, tMin, tMax, q1, t1s, t1e, i, subDivide);
1.219 + }
1.220 + // At this point, we have two ranges of t values -- treat each separately at the split
1.221 + bool result;
1.222 + if (add_intercept(q1, q2, tMin, tsFound[split - 1], i, subDivide)) {
1.223 + result = true;
1.224 + } else {
1.225 + i->swap();
1.226 + result = is_linear_inner(q2, tMin, tsFound[split - 1], q1, t1s, t1e, i, subDivide);
1.227 + }
1.228 + if (add_intercept(q1, q2, tsFound[split], tMax, i, subDivide)) {
1.229 + result = true;
1.230 + } else {
1.231 + i->swap();
1.232 + result |= is_linear_inner(q2, tsFound[split], tMax, q1, t1s, t1e, i, subDivide);
1.233 + }
1.234 + return result;
1.235 +}
1.236 +
1.237 +static double flat_measure(const SkDQuad& q) {
1.238 + SkDVector mid = q[1] - q[0];
1.239 + SkDVector dxy = q[2] - q[0];
1.240 + double length = dxy.length(); // OPTIMIZE: get rid of sqrt
1.241 + return fabs(mid.cross(dxy) / length);
1.242 +}
1.243 +
1.244 +// FIXME ? should this measure both and then use the quad that is the flattest as the line?
1.245 +static bool is_linear(const SkDQuad& q1, const SkDQuad& q2, SkIntersections* i) {
1.246 + double measure = flat_measure(q1);
1.247 + // OPTIMIZE: (get rid of sqrt) use approximately_zero
1.248 + if (!approximately_zero_sqrt(measure)) {
1.249 + return false;
1.250 + }
1.251 + return is_linear_inner(q1, 0, 1, q2, 0, 1, i, NULL);
1.252 +}
1.253 +
1.254 +// FIXME: if flat measure is sufficiently large, then probably the quartic solution failed
1.255 +// avoid imprecision incurred with chopAt
1.256 +static void relaxed_is_linear(const SkDQuad* q1, double s1, double e1, const SkDQuad* q2,
1.257 + double s2, double e2, SkIntersections* i) {
1.258 + double m1 = flat_measure(*q1);
1.259 + double m2 = flat_measure(*q2);
1.260 + i->reset();
1.261 + const SkDQuad* rounder, *flatter;
1.262 + double sf, midf, ef, sr, er;
1.263 + if (m2 < m1) {
1.264 + rounder = q1;
1.265 + sr = s1;
1.266 + er = e1;
1.267 + flatter = q2;
1.268 + sf = s2;
1.269 + midf = (s2 + e2) / 2;
1.270 + ef = e2;
1.271 + } else {
1.272 + rounder = q2;
1.273 + sr = s2;
1.274 + er = e2;
1.275 + flatter = q1;
1.276 + sf = s1;
1.277 + midf = (s1 + e1) / 2;
1.278 + ef = e1;
1.279 + }
1.280 + bool subDivide = false;
1.281 + is_linear_inner(*flatter, sf, ef, *rounder, sr, er, i, &subDivide);
1.282 + if (subDivide) {
1.283 + relaxed_is_linear(flatter, sf, midf, rounder, sr, er, i);
1.284 + relaxed_is_linear(flatter, midf, ef, rounder, sr, er, i);
1.285 + }
1.286 + if (m2 < m1) {
1.287 + i->swapPts();
1.288 + }
1.289 +}
1.290 +
1.291 +// each time through the loop, this computes values it had from the last loop
1.292 +// if i == j == 1, the center values are still good
1.293 +// otherwise, for i != 1 or j != 1, four of the values are still good
1.294 +// and if i == 1 ^ j == 1, an additional value is good
1.295 +static bool binary_search(const SkDQuad& quad1, const SkDQuad& quad2, double* t1Seed,
1.296 + double* t2Seed, SkDPoint* pt) {
1.297 + double tStep = ROUGH_EPSILON;
1.298 + SkDPoint t1[3], t2[3];
1.299 + int calcMask = ~0;
1.300 + do {
1.301 + if (calcMask & (1 << 1)) t1[1] = quad1.ptAtT(*t1Seed);
1.302 + if (calcMask & (1 << 4)) t2[1] = quad2.ptAtT(*t2Seed);
1.303 + if (t1[1].approximatelyEqual(t2[1])) {
1.304 + *pt = t1[1];
1.305 + #if ONE_OFF_DEBUG
1.306 + SkDebugf("%s t1=%1.9g t2=%1.9g (%1.9g,%1.9g) == (%1.9g,%1.9g)\n", __FUNCTION__,
1.307 + t1Seed, t2Seed, t1[1].fX, t1[1].fY, t2[1].fX, t2[1].fY);
1.308 + #endif
1.309 + return true;
1.310 + }
1.311 + if (calcMask & (1 << 0)) t1[0] = quad1.ptAtT(*t1Seed - tStep);
1.312 + if (calcMask & (1 << 2)) t1[2] = quad1.ptAtT(*t1Seed + tStep);
1.313 + if (calcMask & (1 << 3)) t2[0] = quad2.ptAtT(*t2Seed - tStep);
1.314 + if (calcMask & (1 << 5)) t2[2] = quad2.ptAtT(*t2Seed + tStep);
1.315 + double dist[3][3];
1.316 + // OPTIMIZE: using calcMask value permits skipping some distance calcuations
1.317 + // if prior loop's results are moved to correct slot for reuse
1.318 + dist[1][1] = t1[1].distanceSquared(t2[1]);
1.319 + int best_i = 1, best_j = 1;
1.320 + for (int i = 0; i < 3; ++i) {
1.321 + for (int j = 0; j < 3; ++j) {
1.322 + if (i == 1 && j == 1) {
1.323 + continue;
1.324 + }
1.325 + dist[i][j] = t1[i].distanceSquared(t2[j]);
1.326 + if (dist[best_i][best_j] > dist[i][j]) {
1.327 + best_i = i;
1.328 + best_j = j;
1.329 + }
1.330 + }
1.331 + }
1.332 + if (best_i == 1 && best_j == 1) {
1.333 + tStep /= 2;
1.334 + if (tStep < FLT_EPSILON_HALF) {
1.335 + break;
1.336 + }
1.337 + calcMask = (1 << 0) | (1 << 2) | (1 << 3) | (1 << 5);
1.338 + continue;
1.339 + }
1.340 + if (best_i == 0) {
1.341 + *t1Seed -= tStep;
1.342 + t1[2] = t1[1];
1.343 + t1[1] = t1[0];
1.344 + calcMask = 1 << 0;
1.345 + } else if (best_i == 2) {
1.346 + *t1Seed += tStep;
1.347 + t1[0] = t1[1];
1.348 + t1[1] = t1[2];
1.349 + calcMask = 1 << 2;
1.350 + } else {
1.351 + calcMask = 0;
1.352 + }
1.353 + if (best_j == 0) {
1.354 + *t2Seed -= tStep;
1.355 + t2[2] = t2[1];
1.356 + t2[1] = t2[0];
1.357 + calcMask |= 1 << 3;
1.358 + } else if (best_j == 2) {
1.359 + *t2Seed += tStep;
1.360 + t2[0] = t2[1];
1.361 + t2[1] = t2[2];
1.362 + calcMask |= 1 << 5;
1.363 + }
1.364 + } while (true);
1.365 +#if ONE_OFF_DEBUG
1.366 + SkDebugf("%s t1=%1.9g t2=%1.9g (%1.9g,%1.9g) != (%1.9g,%1.9g) %s\n", __FUNCTION__,
1.367 + t1Seed, t2Seed, t1[1].fX, t1[1].fY, t1[2].fX, t1[2].fY);
1.368 +#endif
1.369 + return false;
1.370 +}
1.371 +
1.372 +static void lookNearEnd(const SkDQuad& q1, const SkDQuad& q2, int testT,
1.373 + const SkIntersections& orig, bool swap, SkIntersections* i) {
1.374 + if (orig.used() == 1 && orig[!swap][0] == testT) {
1.375 + return;
1.376 + }
1.377 + if (orig.used() == 2 && orig[!swap][1] == testT) {
1.378 + return;
1.379 + }
1.380 + SkDLine tmpLine;
1.381 + int testTIndex = testT << 1;
1.382 + tmpLine[0] = tmpLine[1] = q2[testTIndex];
1.383 + tmpLine[1].fX += q2[1].fY - q2[testTIndex].fY;
1.384 + tmpLine[1].fY -= q2[1].fX - q2[testTIndex].fX;
1.385 + SkIntersections impTs;
1.386 + impTs.intersectRay(q1, tmpLine);
1.387 + for (int index = 0; index < impTs.used(); ++index) {
1.388 + SkDPoint realPt = impTs.pt(index);
1.389 + if (!tmpLine[0].approximatelyEqual(realPt)) {
1.390 + continue;
1.391 + }
1.392 + if (swap) {
1.393 + i->insert(testT, impTs[0][index], tmpLine[0]);
1.394 + } else {
1.395 + i->insert(impTs[0][index], testT, tmpLine[0]);
1.396 + }
1.397 + }
1.398 +}
1.399 +
1.400 +int SkIntersections::intersect(const SkDQuad& q1, const SkDQuad& q2) {
1.401 + fMax = 4;
1.402 + // if the quads share an end point, check to see if they overlap
1.403 + for (int i1 = 0; i1 < 3; i1 += 2) {
1.404 + for (int i2 = 0; i2 < 3; i2 += 2) {
1.405 + if (q1[i1].asSkPoint() == q2[i2].asSkPoint()) {
1.406 + insert(i1 >> 1, i2 >> 1, q1[i1]);
1.407 + }
1.408 + }
1.409 + }
1.410 + SkASSERT(fUsed < 3);
1.411 + if (only_end_pts_in_common(q1, q2)) {
1.412 + return fUsed;
1.413 + }
1.414 + if (only_end_pts_in_common(q2, q1)) {
1.415 + return fUsed;
1.416 + }
1.417 + // see if either quad is really a line
1.418 + // FIXME: figure out why reduce step didn't find this earlier
1.419 + if (is_linear(q1, q2, this)) {
1.420 + return fUsed;
1.421 + }
1.422 + SkIntersections swapped;
1.423 + swapped.setMax(fMax);
1.424 + if (is_linear(q2, q1, &swapped)) {
1.425 + swapped.swapPts();
1.426 + set(swapped);
1.427 + return fUsed;
1.428 + }
1.429 + SkIntersections copyI(*this);
1.430 + lookNearEnd(q1, q2, 0, *this, false, ©I);
1.431 + lookNearEnd(q1, q2, 1, *this, false, ©I);
1.432 + lookNearEnd(q2, q1, 0, *this, true, ©I);
1.433 + lookNearEnd(q2, q1, 1, *this, true, ©I);
1.434 + int innerEqual = 0;
1.435 + if (copyI.fUsed >= 2) {
1.436 + SkASSERT(copyI.fUsed <= 4);
1.437 + double width = copyI[0][1] - copyI[0][0];
1.438 + int midEnd = 1;
1.439 + for (int index = 2; index < copyI.fUsed; ++index) {
1.440 + double testWidth = copyI[0][index] - copyI[0][index - 1];
1.441 + if (testWidth <= width) {
1.442 + continue;
1.443 + }
1.444 + midEnd = index;
1.445 + }
1.446 + for (int index = 0; index < 2; ++index) {
1.447 + double testT = (copyI[0][midEnd] * (index + 1)
1.448 + + copyI[0][midEnd - 1] * (2 - index)) / 3;
1.449 + SkDPoint testPt1 = q1.ptAtT(testT);
1.450 + testT = (copyI[1][midEnd] * (index + 1) + copyI[1][midEnd - 1] * (2 - index)) / 3;
1.451 + SkDPoint testPt2 = q2.ptAtT(testT);
1.452 + innerEqual += testPt1.approximatelyEqual(testPt2);
1.453 + }
1.454 + }
1.455 + bool expectCoincident = copyI.fUsed >= 2 && innerEqual == 2;
1.456 + if (expectCoincident) {
1.457 + reset();
1.458 + insertCoincident(copyI[0][0], copyI[1][0], copyI.fPt[0]);
1.459 + int last = copyI.fUsed - 1;
1.460 + insertCoincident(copyI[0][last], copyI[1][last], copyI.fPt[last]);
1.461 + return fUsed;
1.462 + }
1.463 + SkDQuadImplicit i1(q1);
1.464 + SkDQuadImplicit i2(q2);
1.465 + int index;
1.466 + bool flip1 = q1[2] == q2[0];
1.467 + bool flip2 = q1[0] == q2[2];
1.468 + bool useCubic = q1[0] == q2[0];
1.469 + double roots1[4];
1.470 + int rootCount = findRoots(i2, q1, roots1, useCubic, flip1, 0);
1.471 + // OPTIMIZATION: could short circuit here if all roots are < 0 or > 1
1.472 + double roots1Copy[4];
1.473 + int r1Count = addValidRoots(roots1, rootCount, roots1Copy);
1.474 + SkDPoint pts1[4];
1.475 + for (index = 0; index < r1Count; ++index) {
1.476 + pts1[index] = q1.ptAtT(roots1Copy[index]);
1.477 + }
1.478 + double roots2[4];
1.479 + int rootCount2 = findRoots(i1, q2, roots2, useCubic, flip2, 0);
1.480 + double roots2Copy[4];
1.481 + int r2Count = addValidRoots(roots2, rootCount2, roots2Copy);
1.482 + SkDPoint pts2[4];
1.483 + for (index = 0; index < r2Count; ++index) {
1.484 + pts2[index] = q2.ptAtT(roots2Copy[index]);
1.485 + }
1.486 + if (r1Count == r2Count && r1Count <= 1) {
1.487 + if (r1Count == 1 && used() == 0) {
1.488 + if (pts1[0].approximatelyEqual(pts2[0])) {
1.489 + insert(roots1Copy[0], roots2Copy[0], pts1[0]);
1.490 + } else if (pts1[0].moreRoughlyEqual(pts2[0])) {
1.491 + // experiment: try to find intersection by chasing t
1.492 + if (binary_search(q1, q2, roots1Copy, roots2Copy, pts1)) {
1.493 + insert(roots1Copy[0], roots2Copy[0], pts1[0]);
1.494 + }
1.495 + }
1.496 + }
1.497 + return fUsed;
1.498 + }
1.499 + int closest[4];
1.500 + double dist[4];
1.501 + bool foundSomething = false;
1.502 + for (index = 0; index < r1Count; ++index) {
1.503 + dist[index] = DBL_MAX;
1.504 + closest[index] = -1;
1.505 + for (int ndex2 = 0; ndex2 < r2Count; ++ndex2) {
1.506 + if (!pts2[ndex2].approximatelyEqual(pts1[index])) {
1.507 + continue;
1.508 + }
1.509 + double dx = pts2[ndex2].fX - pts1[index].fX;
1.510 + double dy = pts2[ndex2].fY - pts1[index].fY;
1.511 + double distance = dx * dx + dy * dy;
1.512 + if (dist[index] <= distance) {
1.513 + continue;
1.514 + }
1.515 + for (int outer = 0; outer < index; ++outer) {
1.516 + if (closest[outer] != ndex2) {
1.517 + continue;
1.518 + }
1.519 + if (dist[outer] < distance) {
1.520 + goto next;
1.521 + }
1.522 + closest[outer] = -1;
1.523 + }
1.524 + dist[index] = distance;
1.525 + closest[index] = ndex2;
1.526 + foundSomething = true;
1.527 + next:
1.528 + ;
1.529 + }
1.530 + }
1.531 + if (r1Count && r2Count && !foundSomething) {
1.532 + relaxed_is_linear(&q1, 0, 1, &q2, 0, 1, this);
1.533 + return fUsed;
1.534 + }
1.535 + int used = 0;
1.536 + do {
1.537 + double lowest = DBL_MAX;
1.538 + int lowestIndex = -1;
1.539 + for (index = 0; index < r1Count; ++index) {
1.540 + if (closest[index] < 0) {
1.541 + continue;
1.542 + }
1.543 + if (roots1Copy[index] < lowest) {
1.544 + lowestIndex = index;
1.545 + lowest = roots1Copy[index];
1.546 + }
1.547 + }
1.548 + if (lowestIndex < 0) {
1.549 + break;
1.550 + }
1.551 + insert(roots1Copy[lowestIndex], roots2Copy[closest[lowestIndex]],
1.552 + pts1[lowestIndex]);
1.553 + closest[lowestIndex] = -1;
1.554 + } while (++used < r1Count);
1.555 + return fUsed;
1.556 +}
