| // Another approach is to start with the implicit form of one curve and solve | |
| // (seek implicit coefficients in QuadraticParameter.cpp | |
| // by substituting in the parametric form of the other. | |
| // The downside of this approach is that early rejects are difficult to come by. | |
| // http://planetmath.org/encyclopedia/GaloisTheoreticDerivationOfTheQuarticFormula.html#step | |
| #include "SkDQuadImplicit.h" | |
| #include "SkIntersections.h" | |
| #include "SkPathOpsLine.h" | |
| #include "SkQuarticRoot.h" | |
| #include "SkTArray.h" | |
| #include "SkTSort.h" | |
| /* given the implicit form 0 = Ax^2 + Bxy + Cy^2 + Dx + Ey + F | |
| * and given x = at^2 + bt + c (the parameterized form) | |
| * y = dt^2 + et + f | |
| * then | |
| * 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 | |
| */ | |
| static int findRoots(const SkDQuadImplicit& i, const SkDQuad& quad, double roots[4], | |
| bool oneHint, bool flip, int firstCubicRoot) { | |
| SkDQuad flipped; | |
| const SkDQuad& q = flip ? (flipped = quad.flip()) : quad; | |
| double a, b, c; | |
| SkDQuad::SetABC(&q[0].fX, &a, &b, &c); | |
| double d, e, f; | |
| SkDQuad::SetABC(&q[0].fY, &d, &e, &f); | |
| const double t4 = i.x2() * a * a | |
| + i.xy() * a * d | |
| + i.y2() * d * d; | |
| const double t3 = 2 * i.x2() * a * b | |
| + i.xy() * (a * e + b * d) | |
| + 2 * i.y2() * d * e; | |
| const double t2 = i.x2() * (b * b + 2 * a * c) | |
| + i.xy() * (c * d + b * e + a * f) | |
| + i.y2() * (e * e + 2 * d * f) | |
| + i.x() * a | |
| + i.y() * d; | |
| const double t1 = 2 * i.x2() * b * c | |
| + i.xy() * (c * e + b * f) | |
| + 2 * i.y2() * e * f | |
| + i.x() * b | |
| + i.y() * e; | |
| const double t0 = i.x2() * c * c | |
| + i.xy() * c * f | |
| + i.y2() * f * f | |
| + i.x() * c | |
| + i.y() * f | |
| + i.c(); | |
| int rootCount = SkReducedQuarticRoots(t4, t3, t2, t1, t0, oneHint, roots); | |
| if (rootCount < 0) { | |
| rootCount = SkQuarticRootsReal(firstCubicRoot, t4, t3, t2, t1, t0, roots); | |
| } | |
| if (flip) { | |
| for (int index = 0; index < rootCount; ++index) { | |
| roots[index] = 1 - roots[index]; | |
| } | |
| } | |
| return rootCount; | |
| } | |
| static int addValidRoots(const double roots[4], const int count, double valid[4]) { | |
| int result = 0; | |
| int index; | |
| for (index = 0; index < count; ++index) { | |
| if (!approximately_zero_or_more(roots[index]) || !approximately_one_or_less(roots[index])) { | |
| continue; | |
| } | |
| double t = 1 - roots[index]; | |
| if (approximately_less_than_zero(t)) { | |
| t = 0; | |
| } else if (approximately_greater_than_one(t)) { | |
| t = 1; | |
| } | |
| valid[result++] = t; | |
| } | |
| return result; | |
| } | |
| static bool only_end_pts_in_common(const SkDQuad& q1, const SkDQuad& q2) { | |
| // the idea here is to see at minimum do a quick reject by rotating all points | |
| // to either side of the line formed by connecting the endpoints | |
| // if the opposite curves points are on the line or on the other side, the | |
| // curves at most intersect at the endpoints | |
| for (int oddMan = 0; oddMan < 3; ++oddMan) { | |
| const SkDPoint* endPt[2]; | |
| for (int opp = 1; opp < 3; ++opp) { | |
| int end = oddMan ^ opp; // choose a value not equal to oddMan | |
| if (3 == end) { // and correct so that largest value is 1 or 2 | |
| end = opp; | |
| } | |
| endPt[opp - 1] = &q1[end]; | |
| } | |
| double origX = endPt[0]->fX; | |
| double origY = endPt[0]->fY; | |
| double adj = endPt[1]->fX - origX; | |
| double opp = endPt[1]->fY - origY; | |
| double sign = (q1[oddMan].fY - origY) * adj - (q1[oddMan].fX - origX) * opp; | |
| if (approximately_zero(sign)) { | |
| goto tryNextHalfPlane; | |
| } | |
| for (int n = 0; n < 3; ++n) { | |
| double test = (q2[n].fY - origY) * adj - (q2[n].fX - origX) * opp; | |
| if (test * sign > 0 && !precisely_zero(test)) { | |
| goto tryNextHalfPlane; | |
| } | |
| } | |
| return true; | |
| tryNextHalfPlane: | |
| ; | |
| } | |
| return false; | |
| } | |
| // returns false if there's more than one intercept or the intercept doesn't match the point | |
| // returns true if the intercept was successfully added or if the | |
| // original quads need to be subdivided | |
| static bool add_intercept(const SkDQuad& q1, const SkDQuad& q2, double tMin, double tMax, | |
| SkIntersections* i, bool* subDivide) { | |
| double tMid = (tMin + tMax) / 2; | |
| SkDPoint mid = q2.ptAtT(tMid); | |
| SkDLine line; | |
| line[0] = line[1] = mid; | |
| SkDVector dxdy = q2.dxdyAtT(tMid); | |
| line[0] -= dxdy; | |
| line[1] += dxdy; | |
| SkIntersections rootTs; | |
| rootTs.allowNear(false); | |
| int roots = rootTs.intersect(q1, line); | |
| if (roots == 0) { | |
| if (subDivide) { | |
| *subDivide = true; | |
| } | |
| return true; | |
| } | |
| if (roots == 2) { | |
| return false; | |
| } | |
| SkDPoint pt2 = q1.ptAtT(rootTs[0][0]); | |
| if (!pt2.approximatelyEqual(mid)) { | |
| return false; | |
| } | |
| i->insertSwap(rootTs[0][0], tMid, pt2); | |
| return true; | |
| } | |
| static bool is_linear_inner(const SkDQuad& q1, double t1s, double t1e, const SkDQuad& q2, | |
| double t2s, double t2e, SkIntersections* i, bool* subDivide) { | |
| SkDQuad hull = q1.subDivide(t1s, t1e); | |
| SkDLine line = {{hull[2], hull[0]}}; | |
| const SkDLine* testLines[] = { &line, (const SkDLine*) &hull[0], (const SkDLine*) &hull[1] }; | |
| const size_t kTestCount = SK_ARRAY_COUNT(testLines); | |
| SkSTArray<kTestCount * 2, double, true> tsFound; | |
| for (size_t index = 0; index < kTestCount; ++index) { | |
| SkIntersections rootTs; | |
| rootTs.allowNear(false); | |
| int roots = rootTs.intersect(q2, *testLines[index]); | |
| for (int idx2 = 0; idx2 < roots; ++idx2) { | |
| double t = rootTs[0][idx2]; | |
| #if 0 // def SK_DEBUG // FIXME : accurate for error = 16, error of 17.5 seen | |
| // {{{136.08723965397621, 1648.2814535211637}, {593.49031197259478, 1190.8784277439891}, {593.49031197259478, 544.0128173828125}}} | |
| // {{{-968.181396484375, 544.0128173828125}, {592.2825927734375, 870.552490234375}, {593.435302734375, 557.8828125}}} | |
| SkDPoint qPt = q2.ptAtT(t); | |
| SkDPoint lPt = testLines[index]->ptAtT(rootTs[1][idx2]); | |
| SkASSERT(qPt.approximatelyDEqual(lPt)); | |
| #endif | |
| if (approximately_negative(t - t2s) || approximately_positive(t - t2e)) { | |
| continue; | |
| } | |
| tsFound.push_back(rootTs[0][idx2]); | |
| } | |
| } | |
| int tCount = tsFound.count(); | |
| if (tCount <= 0) { | |
| return true; | |
| } | |
| double tMin, tMax; | |
| if (tCount == 1) { | |
| tMin = tMax = tsFound[0]; | |
| } else { | |
| SkASSERT(tCount > 1); | |
| SkTQSort<double>(tsFound.begin(), tsFound.end() - 1); | |
| tMin = tsFound[0]; | |
| tMax = tsFound[tsFound.count() - 1]; | |
| } | |
| SkDPoint end = q2.ptAtT(t2s); | |
| bool startInTriangle = hull.pointInHull(end); | |
| if (startInTriangle) { | |
| tMin = t2s; | |
| } | |
| end = q2.ptAtT(t2e); | |
| bool endInTriangle = hull.pointInHull(end); | |
| if (endInTriangle) { | |
| tMax = t2e; | |
| } | |
| int split = 0; | |
| SkDVector dxy1, dxy2; | |
| if (tMin != tMax || tCount > 2) { | |
| dxy2 = q2.dxdyAtT(tMin); | |
| for (int index = 1; index < tCount; ++index) { | |
| dxy1 = dxy2; | |
| dxy2 = q2.dxdyAtT(tsFound[index]); | |
| double dot = dxy1.dot(dxy2); | |
| if (dot < 0) { | |
| split = index - 1; | |
| break; | |
| } | |
| } | |
| } | |
| if (split == 0) { // there's one point | |
| if (add_intercept(q1, q2, tMin, tMax, i, subDivide)) { | |
| return true; | |
| } | |
| i->swap(); | |
| return is_linear_inner(q2, tMin, tMax, q1, t1s, t1e, i, subDivide); | |
| } | |
| // At this point, we have two ranges of t values -- treat each separately at the split | |
| bool result; | |
| if (add_intercept(q1, q2, tMin, tsFound[split - 1], i, subDivide)) { | |
| result = true; | |
| } else { | |
| i->swap(); | |
| result = is_linear_inner(q2, tMin, tsFound[split - 1], q1, t1s, t1e, i, subDivide); | |
| } | |
| if (add_intercept(q1, q2, tsFound[split], tMax, i, subDivide)) { | |
| result = true; | |
| } else { | |
| i->swap(); | |
| result |= is_linear_inner(q2, tsFound[split], tMax, q1, t1s, t1e, i, subDivide); | |
| } | |
| return result; | |
| } | |
| static double flat_measure(const SkDQuad& q) { | |
| SkDVector mid = q[1] - q[0]; | |
| SkDVector dxy = q[2] - q[0]; | |
| double length = dxy.length(); // OPTIMIZE: get rid of sqrt | |
| return fabs(mid.cross(dxy) / length); | |
| } | |
| // FIXME ? should this measure both and then use the quad that is the flattest as the line? | |
| static bool is_linear(const SkDQuad& q1, const SkDQuad& q2, SkIntersections* i) { | |
| double measure = flat_measure(q1); | |
| // OPTIMIZE: (get rid of sqrt) use approximately_zero | |
| if (!approximately_zero_sqrt(measure)) { | |
| return false; | |
| } | |
| return is_linear_inner(q1, 0, 1, q2, 0, 1, i, NULL); | |
| } | |
| // FIXME: if flat measure is sufficiently large, then probably the quartic solution failed | |
| // avoid imprecision incurred with chopAt | |
| static void relaxed_is_linear(const SkDQuad* q1, double s1, double e1, const SkDQuad* q2, | |
| double s2, double e2, SkIntersections* i) { | |
| double m1 = flat_measure(*q1); | |
| double m2 = flat_measure(*q2); | |
| i->reset(); | |
| const SkDQuad* rounder, *flatter; | |
| double sf, midf, ef, sr, er; | |
| if (m2 < m1) { | |
| rounder = q1; | |
| sr = s1; | |
| er = e1; | |
| flatter = q2; | |
| sf = s2; | |
| midf = (s2 + e2) / 2; | |
| ef = e2; | |
| } else { | |
| rounder = q2; | |
| sr = s2; | |
| er = e2; | |
| flatter = q1; | |
| sf = s1; | |
| midf = (s1 + e1) / 2; | |
| ef = e1; | |
| } | |
| bool subDivide = false; | |
| is_linear_inner(*flatter, sf, ef, *rounder, sr, er, i, &subDivide); | |
| if (subDivide) { | |
| relaxed_is_linear(flatter, sf, midf, rounder, sr, er, i); | |
| relaxed_is_linear(flatter, midf, ef, rounder, sr, er, i); | |
| } | |
| if (m2 < m1) { | |
| i->swapPts(); | |
| } | |
| } | |
| // each time through the loop, this computes values it had from the last loop | |
| // if i == j == 1, the center values are still good | |
| // otherwise, for i != 1 or j != 1, four of the values are still good | |
| // and if i == 1 ^ j == 1, an additional value is good | |
| static bool binary_search(const SkDQuad& quad1, const SkDQuad& quad2, double* t1Seed, | |
| double* t2Seed, SkDPoint* pt) { | |
| double tStep = ROUGH_EPSILON; | |
| SkDPoint t1[3], t2[3]; | |
| int calcMask = ~0; | |
| do { | |
| if (calcMask & (1 << 1)) t1[1] = quad1.ptAtT(*t1Seed); | |
| if (calcMask & (1 << 4)) t2[1] = quad2.ptAtT(*t2Seed); | |
| if (t1[1].approximatelyEqual(t2[1])) { | |
| *pt = t1[1]; | |
| #if ONE_OFF_DEBUG | |
| SkDebugf("%s t1=%1.9g t2=%1.9g (%1.9g,%1.9g) == (%1.9g,%1.9g)\n", __FUNCTION__, | |
| t1Seed, t2Seed, t1[1].fX, t1[1].fY, t2[1].fX, t2[1].fY); | |
| #endif | |
| return true; | |
| } | |
| if (calcMask & (1 << 0)) t1[0] = quad1.ptAtT(SkTMax(0., *t1Seed - tStep)); | |
| if (calcMask & (1 << 2)) t1[2] = quad1.ptAtT(SkTMin(1., *t1Seed + tStep)); | |
| if (calcMask & (1 << 3)) t2[0] = quad2.ptAtT(SkTMax(0., *t2Seed - tStep)); | |
| if (calcMask & (1 << 5)) t2[2] = quad2.ptAtT(SkTMin(1., *t2Seed + tStep)); | |
| double dist[3][3]; | |
| // OPTIMIZE: using calcMask value permits skipping some distance calcuations | |
| // if prior loop's results are moved to correct slot for reuse | |
| dist[1][1] = t1[1].distanceSquared(t2[1]); | |
| int best_i = 1, best_j = 1; | |
| for (int i = 0; i < 3; ++i) { | |
| for (int j = 0; j < 3; ++j) { | |
| if (i == 1 && j == 1) { | |
| continue; | |
| } | |
| dist[i][j] = t1[i].distanceSquared(t2[j]); | |
| if (dist[best_i][best_j] > dist[i][j]) { | |
| best_i = i; | |
| best_j = j; | |
| } | |
| } | |
| } | |
| if (best_i == 1 && best_j == 1) { | |
| tStep /= 2; | |
| if (tStep < FLT_EPSILON_HALF) { | |
| break; | |
| } | |
| calcMask = (1 << 0) | (1 << 2) | (1 << 3) | (1 << 5); | |
| continue; | |
| } | |
| if (best_i == 0) { | |
| *t1Seed -= tStep; | |
| t1[2] = t1[1]; | |
| t1[1] = t1[0]; | |
| calcMask = 1 << 0; | |
| } else if (best_i == 2) { | |
| *t1Seed += tStep; | |
| t1[0] = t1[1]; | |
| t1[1] = t1[2]; | |
| calcMask = 1 << 2; | |
| } else { | |
| calcMask = 0; | |
| } | |
| if (best_j == 0) { | |
| *t2Seed -= tStep; | |
| t2[2] = t2[1]; | |
| t2[1] = t2[0]; | |
| calcMask |= 1 << 3; | |
| } else if (best_j == 2) { | |
| *t2Seed += tStep; | |
| t2[0] = t2[1]; | |
| t2[1] = t2[2]; | |
| calcMask |= 1 << 5; | |
| } | |
| } while (true); | |
| #if ONE_OFF_DEBUG | |
| SkDebugf("%s t1=%1.9g t2=%1.9g (%1.9g,%1.9g) != (%1.9g,%1.9g) %s\n", __FUNCTION__, | |
| t1Seed, t2Seed, t1[1].fX, t1[1].fY, t1[2].fX, t1[2].fY); | |
| #endif | |
| return false; | |
| } | |
| static void lookNearEnd(const SkDQuad& q1, const SkDQuad& q2, int testT, | |
| const SkIntersections& orig, bool swap, SkIntersections* i) { | |
| if (orig.used() == 1 && orig[!swap][0] == testT) { | |
| return; | |
| } | |
| if (orig.used() == 2 && orig[!swap][1] == testT) { | |
| return; | |
| } | |
| SkDLine tmpLine; | |
| int testTIndex = testT << 1; | |
| tmpLine[0] = tmpLine[1] = q2[testTIndex]; | |
| tmpLine[1].fX += q2[1].fY - q2[testTIndex].fY; | |
| tmpLine[1].fY -= q2[1].fX - q2[testTIndex].fX; | |
| SkIntersections impTs; | |
| impTs.intersectRay(q1, tmpLine); | |
| for (int index = 0; index < impTs.used(); ++index) { | |
| SkDPoint realPt = impTs.pt(index); | |
| if (!tmpLine[0].approximatelyPEqual(realPt)) { | |
| continue; | |
| } | |
| if (swap) { | |
| i->insert(testT, impTs[0][index], tmpLine[0]); | |
| } else { | |
| i->insert(impTs[0][index], testT, tmpLine[0]); | |
| } | |
| } | |
| } | |
| int SkIntersections::intersect(const SkDQuad& q1, const SkDQuad& q2) { | |
| fMax = 4; | |
| // if the quads share an end point, check to see if they overlap | |
| for (int i1 = 0; i1 < 3; i1 += 2) { | |
| for (int i2 = 0; i2 < 3; i2 += 2) { | |
| if (q1[i1].asSkPoint() == q2[i2].asSkPoint()) { | |
| insert(i1 >> 1, i2 >> 1, q1[i1]); | |
| } | |
| } | |
| } | |
| SkASSERT(fUsed < 3); | |
| if (only_end_pts_in_common(q1, q2)) { | |
| return fUsed; | |
| } | |
| if (only_end_pts_in_common(q2, q1)) { | |
| return fUsed; | |
| } | |
| // see if either quad is really a line | |
| // FIXME: figure out why reduce step didn't find this earlier | |
| if (is_linear(q1, q2, this)) { | |
| return fUsed; | |
| } | |
| SkIntersections swapped; | |
| swapped.setMax(fMax); | |
| if (is_linear(q2, q1, &swapped)) { | |
| swapped.swapPts(); | |
| *this = swapped; | |
| return fUsed; | |
| } | |
| SkIntersections copyI(*this); | |
| lookNearEnd(q1, q2, 0, *this, false, ©I); | |
| lookNearEnd(q1, q2, 1, *this, false, ©I); | |
| lookNearEnd(q2, q1, 0, *this, true, ©I); | |
| lookNearEnd(q2, q1, 1, *this, true, ©I); | |
| int innerEqual = 0; | |
| if (copyI.fUsed >= 2) { | |
| SkASSERT(copyI.fUsed <= 4); | |
| double width = copyI[0][1] - copyI[0][0]; | |
| int midEnd = 1; | |
| for (int index = 2; index < copyI.fUsed; ++index) { | |
| double testWidth = copyI[0][index] - copyI[0][index - 1]; | |
| if (testWidth <= width) { | |
| continue; | |
| } | |
| midEnd = index; | |
| } | |
| for (int index = 0; index < 2; ++index) { | |
| double testT = (copyI[0][midEnd] * (index + 1) | |
| + copyI[0][midEnd - 1] * (2 - index)) / 3; | |
| SkDPoint testPt1 = q1.ptAtT(testT); | |
| testT = (copyI[1][midEnd] * (index + 1) + copyI[1][midEnd - 1] * (2 - index)) / 3; | |
| SkDPoint testPt2 = q2.ptAtT(testT); | |
| innerEqual += testPt1.approximatelyEqual(testPt2); | |
| } | |
| } | |
| bool expectCoincident = copyI.fUsed >= 2 && innerEqual == 2; | |
| if (expectCoincident) { | |
| reset(); | |
| insertCoincident(copyI[0][0], copyI[1][0], copyI.fPt[0]); | |
| int last = copyI.fUsed - 1; | |
| insertCoincident(copyI[0][last], copyI[1][last], copyI.fPt[last]); | |
| return fUsed; | |
| } | |
| SkDQuadImplicit i1(q1); | |
| SkDQuadImplicit i2(q2); | |
| int index; | |
| bool flip1 = q1[2] == q2[0]; | |
| bool flip2 = q1[0] == q2[2]; | |
| bool useCubic = q1[0] == q2[0]; | |
| double roots1[4]; | |
| int rootCount = findRoots(i2, q1, roots1, useCubic, flip1, 0); | |
| // OPTIMIZATION: could short circuit here if all roots are < 0 or > 1 | |
| double roots1Copy[4]; | |
| int r1Count = addValidRoots(roots1, rootCount, roots1Copy); | |
| SkDPoint pts1[4]; | |
| for (index = 0; index < r1Count; ++index) { | |
| pts1[index] = q1.ptAtT(roots1Copy[index]); | |
| } | |
| double roots2[4]; | |
| int rootCount2 = findRoots(i1, q2, roots2, useCubic, flip2, 0); | |
| double roots2Copy[4]; | |
| int r2Count = addValidRoots(roots2, rootCount2, roots2Copy); | |
| SkDPoint pts2[4]; | |
| for (index = 0; index < r2Count; ++index) { | |
| pts2[index] = q2.ptAtT(roots2Copy[index]); | |
| } | |
| if (r1Count == r2Count && r1Count <= 1) { | |
| if (r1Count == 1 && used() == 0) { | |
| if (pts1[0].approximatelyEqual(pts2[0])) { | |
| insert(roots1Copy[0], roots2Copy[0], pts1[0]); | |
| } else { | |
| // find intersection by chasing t | |
| if (binary_search(q1, q2, roots1Copy, roots2Copy, pts1)) { | |
| insert(roots1Copy[0], roots2Copy[0], pts1[0]); | |
| } | |
| } | |
| } | |
| return fUsed; | |
| } | |
| int closest[4]; | |
| double dist[4]; | |
| bool foundSomething = false; | |
| for (index = 0; index < r1Count; ++index) { | |
| dist[index] = DBL_MAX; | |
| closest[index] = -1; | |
| for (int ndex2 = 0; ndex2 < r2Count; ++ndex2) { | |
| if (!pts2[ndex2].approximatelyEqual(pts1[index])) { | |
| continue; | |
| } | |
| double dx = pts2[ndex2].fX - pts1[index].fX; | |
| double dy = pts2[ndex2].fY - pts1[index].fY; | |
| double distance = dx * dx + dy * dy; | |
| if (dist[index] <= distance) { | |
| continue; | |
| } | |
| for (int outer = 0; outer < index; ++outer) { | |
| if (closest[outer] != ndex2) { | |
| continue; | |
| } | |
| if (dist[outer] < distance) { | |
| goto next; | |
| } | |
| closest[outer] = -1; | |
| } | |
| dist[index] = distance; | |
| closest[index] = ndex2; | |
| foundSomething = true; | |
| next: | |
| ; | |
| } | |
| } | |
| if (r1Count && r2Count && !foundSomething) { | |
| relaxed_is_linear(&q1, 0, 1, &q2, 0, 1, this); | |
| return fUsed; | |
| } | |
| int used = 0; | |
| do { | |
| double lowest = DBL_MAX; | |
| int lowestIndex = -1; | |
| for (index = 0; index < r1Count; ++index) { | |
| if (closest[index] < 0) { | |
| continue; | |
| } | |
| if (roots1Copy[index] < lowest) { | |
| lowestIndex = index; | |
| lowest = roots1Copy[index]; | |
| } | |
| } | |
| if (lowestIndex < 0) { | |
| break; | |
| } | |
| insert(roots1Copy[lowestIndex], roots2Copy[closest[lowestIndex]], | |
| pts1[lowestIndex]); | |
| closest[lowestIndex] = -1; | |
| } while (++used < r1Count); | |
| return fUsed; | |
| } |