| /* |
| * Copyright 2012 Google Inc. |
| * |
| * Use of this source code is governed by a BSD-style license that can be |
| * found in the LICENSE file. |
| */ |
| #include "SkGeometry.h" |
| #include "SkLineParameters.h" |
| #include "SkPathOpsConic.h" |
| #include "SkPathOpsCubic.h" |
| #include "SkPathOpsCurve.h" |
| #include "SkPathOpsLine.h" |
| #include "SkPathOpsQuad.h" |
| #include "SkPathOpsRect.h" |
| #include "SkTSort.h" |
| |
| const int SkDCubic::gPrecisionUnit = 256; // FIXME: test different values in test framework |
| |
| void SkDCubic::align(int endIndex, int ctrlIndex, SkDPoint* dstPt) const { |
| if (fPts[endIndex].fX == fPts[ctrlIndex].fX) { |
| dstPt->fX = fPts[endIndex].fX; |
| } |
| if (fPts[endIndex].fY == fPts[ctrlIndex].fY) { |
| dstPt->fY = fPts[endIndex].fY; |
| } |
| } |
| |
| // give up when changing t no longer moves point |
| // also, copy point rather than recompute it when it does change |
| double SkDCubic::binarySearch(double min, double max, double axisIntercept, |
| SearchAxis xAxis) const { |
| double t = (min + max) / 2; |
| double step = (t - min) / 2; |
| SkDPoint cubicAtT = ptAtT(t); |
| double calcPos = (&cubicAtT.fX)[xAxis]; |
| double calcDist = calcPos - axisIntercept; |
| do { |
| double priorT = t - step; |
| SkASSERT(priorT >= min); |
| SkDPoint lessPt = ptAtT(priorT); |
| if (approximately_equal_half(lessPt.fX, cubicAtT.fX) |
| && approximately_equal_half(lessPt.fY, cubicAtT.fY)) { |
| return -1; // binary search found no point at this axis intercept |
| } |
| double lessDist = (&lessPt.fX)[xAxis] - axisIntercept; |
| #if DEBUG_CUBIC_BINARY_SEARCH |
| SkDebugf("t=%1.9g calc=%1.9g dist=%1.9g step=%1.9g less=%1.9g\n", t, calcPos, calcDist, |
| step, lessDist); |
| #endif |
| double lastStep = step; |
| step /= 2; |
| if (calcDist > 0 ? calcDist > lessDist : calcDist < lessDist) { |
| t = priorT; |
| } else { |
| double nextT = t + lastStep; |
| if (nextT > max) { |
| return -1; |
| } |
| SkDPoint morePt = ptAtT(nextT); |
| if (approximately_equal_half(morePt.fX, cubicAtT.fX) |
| && approximately_equal_half(morePt.fY, cubicAtT.fY)) { |
| return -1; // binary search found no point at this axis intercept |
| } |
| double moreDist = (&morePt.fX)[xAxis] - axisIntercept; |
| if (calcDist > 0 ? calcDist <= moreDist : calcDist >= moreDist) { |
| continue; |
| } |
| t = nextT; |
| } |
| SkDPoint testAtT = ptAtT(t); |
| cubicAtT = testAtT; |
| calcPos = (&cubicAtT.fX)[xAxis]; |
| calcDist = calcPos - axisIntercept; |
| } while (!approximately_equal(calcPos, axisIntercept)); |
| return t; |
| } |
| |
| // FIXME: cache keep the bounds and/or precision with the caller? |
| double SkDCubic::calcPrecision() const { |
| SkDRect dRect; |
| dRect.setBounds(*this); // OPTIMIZATION: just use setRawBounds ? |
| double width = dRect.fRight - dRect.fLeft; |
| double height = dRect.fBottom - dRect.fTop; |
| return (width > height ? width : height) / gPrecisionUnit; |
| } |
| |
| |
| /* classic one t subdivision */ |
| static void interp_cubic_coords(const double* src, double* dst, double t) { |
| double ab = SkDInterp(src[0], src[2], t); |
| double bc = SkDInterp(src[2], src[4], t); |
| double cd = SkDInterp(src[4], src[6], t); |
| double abc = SkDInterp(ab, bc, t); |
| double bcd = SkDInterp(bc, cd, t); |
| double abcd = SkDInterp(abc, bcd, t); |
| |
| dst[0] = src[0]; |
| dst[2] = ab; |
| dst[4] = abc; |
| dst[6] = abcd; |
| dst[8] = bcd; |
| dst[10] = cd; |
| dst[12] = src[6]; |
| } |
| |
| SkDCubicPair SkDCubic::chopAt(double t) const { |
| SkDCubicPair dst; |
| if (t == 0.5) { |
| dst.pts[0] = fPts[0]; |
| dst.pts[1].fX = (fPts[0].fX + fPts[1].fX) / 2; |
| dst.pts[1].fY = (fPts[0].fY + fPts[1].fY) / 2; |
| dst.pts[2].fX = (fPts[0].fX + 2 * fPts[1].fX + fPts[2].fX) / 4; |
| dst.pts[2].fY = (fPts[0].fY + 2 * fPts[1].fY + fPts[2].fY) / 4; |
| dst.pts[3].fX = (fPts[0].fX + 3 * (fPts[1].fX + fPts[2].fX) + fPts[3].fX) / 8; |
| dst.pts[3].fY = (fPts[0].fY + 3 * (fPts[1].fY + fPts[2].fY) + fPts[3].fY) / 8; |
| dst.pts[4].fX = (fPts[1].fX + 2 * fPts[2].fX + fPts[3].fX) / 4; |
| dst.pts[4].fY = (fPts[1].fY + 2 * fPts[2].fY + fPts[3].fY) / 4; |
| dst.pts[5].fX = (fPts[2].fX + fPts[3].fX) / 2; |
| dst.pts[5].fY = (fPts[2].fY + fPts[3].fY) / 2; |
| dst.pts[6] = fPts[3]; |
| return dst; |
| } |
| interp_cubic_coords(&fPts[0].fX, &dst.pts[0].fX, t); |
| interp_cubic_coords(&fPts[0].fY, &dst.pts[0].fY, t); |
| return dst; |
| } |
| |
| void SkDCubic::Coefficients(const double* src, double* A, double* B, double* C, double* D) { |
| *A = src[6]; // d |
| *B = src[4] * 3; // 3*c |
| *C = src[2] * 3; // 3*b |
| *D = src[0]; // a |
| *A -= *D - *C + *B; // A = -a + 3*b - 3*c + d |
| *B += 3 * *D - 2 * *C; // B = 3*a - 6*b + 3*c |
| *C -= 3 * *D; // C = -3*a + 3*b |
| } |
| |
| bool SkDCubic::endsAreExtremaInXOrY() const { |
| return (between(fPts[0].fX, fPts[1].fX, fPts[3].fX) |
| && between(fPts[0].fX, fPts[2].fX, fPts[3].fX)) |
| || (between(fPts[0].fY, fPts[1].fY, fPts[3].fY) |
| && between(fPts[0].fY, fPts[2].fY, fPts[3].fY)); |
| } |
| |
| // Do a quick reject by rotating all points relative to a line formed by |
| // a pair of one cubic's points. If the 2nd cubic's points |
| // are on the line or on the opposite side from the 1st cubic's 'odd man', the |
| // curves at most intersect at the endpoints. |
| /* if returning true, check contains true if cubic's hull collapsed, making the cubic linear |
| if returning false, check contains true if the the cubic pair have only the end point in common |
| */ |
| bool SkDCubic::hullIntersects(const SkDPoint* pts, int ptCount, bool* isLinear) const { |
| bool linear = true; |
| char hullOrder[4]; |
| int hullCount = convexHull(hullOrder); |
| int end1 = hullOrder[0]; |
| int hullIndex = 0; |
| const SkDPoint* endPt[2]; |
| endPt[0] = &fPts[end1]; |
| do { |
| hullIndex = (hullIndex + 1) % hullCount; |
| int end2 = hullOrder[hullIndex]; |
| endPt[1] = &fPts[end2]; |
| double origX = endPt[0]->fX; |
| double origY = endPt[0]->fY; |
| double adj = endPt[1]->fX - origX; |
| double opp = endPt[1]->fY - origY; |
| int oddManMask = other_two(end1, end2); |
| int oddMan = end1 ^ oddManMask; |
| double sign = (fPts[oddMan].fY - origY) * adj - (fPts[oddMan].fX - origX) * opp; |
| int oddMan2 = end2 ^ oddManMask; |
| double sign2 = (fPts[oddMan2].fY - origY) * adj - (fPts[oddMan2].fX - origX) * opp; |
| if (sign * sign2 < 0) { |
| continue; |
| } |
| if (approximately_zero(sign)) { |
| sign = sign2; |
| if (approximately_zero(sign)) { |
| continue; |
| } |
| } |
| linear = false; |
| bool foundOutlier = false; |
| for (int n = 0; n < ptCount; ++n) { |
| double test = (pts[n].fY - origY) * adj - (pts[n].fX - origX) * opp; |
| if (test * sign > 0 && !precisely_zero(test)) { |
| foundOutlier = true; |
| break; |
| } |
| } |
| if (!foundOutlier) { |
| return false; |
| } |
| endPt[0] = endPt[1]; |
| end1 = end2; |
| } while (hullIndex); |
| *isLinear = linear; |
| return true; |
| } |
| |
| bool SkDCubic::hullIntersects(const SkDCubic& c2, bool* isLinear) const { |
| return hullIntersects(c2.fPts, c2.kPointCount, isLinear); |
| } |
| |
| bool SkDCubic::hullIntersects(const SkDQuad& quad, bool* isLinear) const { |
| return hullIntersects(quad.fPts, quad.kPointCount, isLinear); |
| } |
| |
| bool SkDCubic::hullIntersects(const SkDConic& conic, bool* isLinear) const { |
| |
| return hullIntersects(conic.fPts, isLinear); |
| } |
| |
| bool SkDCubic::isLinear(int startIndex, int endIndex) const { |
| if (fPts[0].approximatelyDEqual(fPts[3])) { |
| return ((const SkDQuad *) this)->isLinear(0, 2); |
| } |
| SkLineParameters lineParameters; |
| lineParameters.cubicEndPoints(*this, startIndex, endIndex); |
| // FIXME: maybe it's possible to avoid this and compare non-normalized |
| lineParameters.normalize(); |
| double tiniest = SkTMin(SkTMin(SkTMin(SkTMin(SkTMin(SkTMin(SkTMin(fPts[0].fX, fPts[0].fY), |
| fPts[1].fX), fPts[1].fY), fPts[2].fX), fPts[2].fY), fPts[3].fX), fPts[3].fY); |
| double largest = SkTMax(SkTMax(SkTMax(SkTMax(SkTMax(SkTMax(SkTMax(fPts[0].fX, fPts[0].fY), |
| fPts[1].fX), fPts[1].fY), fPts[2].fX), fPts[2].fY), fPts[3].fX), fPts[3].fY); |
| largest = SkTMax(largest, -tiniest); |
| double distance = lineParameters.controlPtDistance(*this, 1); |
| if (!approximately_zero_when_compared_to(distance, largest)) { |
| return false; |
| } |
| distance = lineParameters.controlPtDistance(*this, 2); |
| return approximately_zero_when_compared_to(distance, largest); |
| } |
| |
| bool SkDCubic::ComplexBreak(const SkPoint pointsPtr[4], SkScalar* t) { |
| SkScalar d[3]; |
| SkCubicType cubicType = SkClassifyCubic(pointsPtr, d); |
| if (cubicType == kLoop_SkCubicType) { |
| // crib code from gpu path utils that finds t values where loop self-intersects |
| // use it to find mid of t values which should be a friendly place to chop |
| SkScalar tempSqrt = SkScalarSqrt(4.f * d[0] * d[2] - 3.f * d[1] * d[1]); |
| SkScalar ls = d[1] - tempSqrt; |
| SkScalar lt = 2.f * d[0]; |
| SkScalar ms = d[1] + tempSqrt; |
| SkScalar mt = 2.f * d[0]; |
| if (roughly_between(0, ls, lt) && roughly_between(0, ms, mt)) { |
| ls = ls / lt; |
| ms = ms / mt; |
| SkASSERT(roughly_between(0, ls, 1) && roughly_between(0, ms, 1)); |
| *t = (ls + ms) / 2; |
| SkASSERT(roughly_between(0, *t, 1)); |
| return *t > 0 && *t < 1; |
| } |
| } else if (kSerpentine_SkCubicType == cubicType || kCusp_SkCubicType == cubicType) { |
| SkDCubic cubic; |
| cubic.set(pointsPtr); |
| double inflectionTs[2]; |
| int infTCount = cubic.findInflections(inflectionTs); |
| if (infTCount == 2) { |
| double maxCurvature[3]; |
| int roots = cubic.findMaxCurvature(maxCurvature); |
| #if DEBUG_CUBIC_SPLIT |
| SkDebugf("%s\n", __FUNCTION__); |
| cubic.dump(); |
| for (int index = 0; index < infTCount; ++index) { |
| SkDebugf("inflectionsTs[%d]=%1.9g ", index, inflectionTs[index]); |
| SkDPoint pt = cubic.ptAtT(inflectionTs[index]); |
| SkDVector dPt = cubic.dxdyAtT(inflectionTs[index]); |
| SkDLine perp = {{pt - dPt, pt + dPt}}; |
| perp.dump(); |
| } |
| for (int index = 0; index < roots; ++index) { |
| SkDebugf("maxCurvature[%d]=%1.9g ", index, maxCurvature[index]); |
| SkDPoint pt = cubic.ptAtT(maxCurvature[index]); |
| SkDVector dPt = cubic.dxdyAtT(maxCurvature[index]); |
| SkDLine perp = {{pt - dPt, pt + dPt}}; |
| perp.dump(); |
| } |
| #endif |
| for (int index = 0; index < roots; ++index) { |
| if (between(inflectionTs[0], maxCurvature[index], inflectionTs[1])) { |
| *t = maxCurvature[index]; |
| return *t > 0 && *t < 1; |
| } |
| } |
| } else if (infTCount == 1) { |
| *t = inflectionTs[0]; |
| return *t > 0 && *t < 1; |
| } |
| } |
| return false; |
| } |
| |
| bool SkDCubic::monotonicInX() const { |
| return precisely_between(fPts[0].fX, fPts[1].fX, fPts[3].fX) |
| && precisely_between(fPts[0].fX, fPts[2].fX, fPts[3].fX); |
| } |
| |
| bool SkDCubic::monotonicInY() const { |
| return precisely_between(fPts[0].fY, fPts[1].fY, fPts[3].fY) |
| && precisely_between(fPts[0].fY, fPts[2].fY, fPts[3].fY); |
| } |
| |
| void SkDCubic::otherPts(int index, const SkDPoint* o1Pts[kPointCount - 1]) const { |
| int offset = (int) !SkToBool(index); |
| o1Pts[0] = &fPts[offset]; |
| o1Pts[1] = &fPts[++offset]; |
| o1Pts[2] = &fPts[++offset]; |
| } |
| |
| int SkDCubic::searchRoots(double extremeTs[6], int extrema, double axisIntercept, |
| SearchAxis xAxis, double* validRoots) const { |
| extrema += findInflections(&extremeTs[extrema]); |
| extremeTs[extrema++] = 0; |
| extremeTs[extrema] = 1; |
| SkASSERT(extrema < 6); |
| SkTQSort(extremeTs, extremeTs + extrema); |
| int validCount = 0; |
| for (int index = 0; index < extrema; ) { |
| double min = extremeTs[index]; |
| double max = extremeTs[++index]; |
| if (min == max) { |
| continue; |
| } |
| double newT = binarySearch(min, max, axisIntercept, xAxis); |
| if (newT >= 0) { |
| if (validCount >= 3) { |
| return 0; |
| } |
| validRoots[validCount++] = newT; |
| } |
| } |
| return validCount; |
| } |
| |
| // cubic roots |
| |
| static const double PI = 3.141592653589793; |
| |
| // from SkGeometry.cpp (and Numeric Solutions, 5.6) |
| int SkDCubic::RootsValidT(double A, double B, double C, double D, double t[3]) { |
| double s[3]; |
| int realRoots = RootsReal(A, B, C, D, s); |
| int foundRoots = SkDQuad::AddValidTs(s, realRoots, t); |
| for (int index = 0; index < realRoots; ++index) { |
| double tValue = s[index]; |
| if (!approximately_one_or_less(tValue) && between(1, tValue, 1.00005)) { |
| for (int idx2 = 0; idx2 < foundRoots; ++idx2) { |
| if (approximately_equal(t[idx2], 1)) { |
| goto nextRoot; |
| } |
| } |
| SkASSERT(foundRoots < 3); |
| t[foundRoots++] = 1; |
| } else if (!approximately_zero_or_more(tValue) && between(-0.00005, tValue, 0)) { |
| for (int idx2 = 0; idx2 < foundRoots; ++idx2) { |
| if (approximately_equal(t[idx2], 0)) { |
| goto nextRoot; |
| } |
| } |
| SkASSERT(foundRoots < 3); |
| t[foundRoots++] = 0; |
| } |
| nextRoot: |
| ; |
| } |
| return foundRoots; |
| } |
| |
| int SkDCubic::RootsReal(double A, double B, double C, double D, double s[3]) { |
| #ifdef SK_DEBUG |
| // create a string mathematica understands |
| // GDB set print repe 15 # if repeated digits is a bother |
| // set print elements 400 # if line doesn't fit |
| char str[1024]; |
| sk_bzero(str, sizeof(str)); |
| SK_SNPRINTF(str, sizeof(str), "Solve[%1.19g x^3 + %1.19g x^2 + %1.19g x + %1.19g == 0, x]", |
| A, B, C, D); |
| SkPathOpsDebug::MathematicaIze(str, sizeof(str)); |
| #if ONE_OFF_DEBUG && ONE_OFF_DEBUG_MATHEMATICA |
| SkDebugf("%s\n", str); |
| #endif |
| #endif |
| if (approximately_zero(A) |
| && approximately_zero_when_compared_to(A, B) |
| && approximately_zero_when_compared_to(A, C) |
| && approximately_zero_when_compared_to(A, D)) { // we're just a quadratic |
| return SkDQuad::RootsReal(B, C, D, s); |
| } |
| if (approximately_zero_when_compared_to(D, A) |
| && approximately_zero_when_compared_to(D, B) |
| && approximately_zero_when_compared_to(D, C)) { // 0 is one root |
| int num = SkDQuad::RootsReal(A, B, C, s); |
| for (int i = 0; i < num; ++i) { |
| if (approximately_zero(s[i])) { |
| return num; |
| } |
| } |
| s[num++] = 0; |
| return num; |
| } |
| if (approximately_zero(A + B + C + D)) { // 1 is one root |
| int num = SkDQuad::RootsReal(A, A + B, -D, s); |
| for (int i = 0; i < num; ++i) { |
| if (AlmostDequalUlps(s[i], 1)) { |
| return num; |
| } |
| } |
| s[num++] = 1; |
| return num; |
| } |
| double a, b, c; |
| { |
| double invA = 1 / A; |
| a = B * invA; |
| b = C * invA; |
| c = D * invA; |
| } |
| double a2 = a * a; |
| double Q = (a2 - b * 3) / 9; |
| double R = (2 * a2 * a - 9 * a * b + 27 * c) / 54; |
| double R2 = R * R; |
| double Q3 = Q * Q * Q; |
| double R2MinusQ3 = R2 - Q3; |
| double adiv3 = a / 3; |
| double r; |
| double* roots = s; |
| if (R2MinusQ3 < 0) { // we have 3 real roots |
| // the divide/root can, due to finite precisions, be slightly outside of -1...1 |
| double theta = acos(SkTPin(R / sqrt(Q3), -1., 1.)); |
| double neg2RootQ = -2 * sqrt(Q); |
| |
| r = neg2RootQ * cos(theta / 3) - adiv3; |
| *roots++ = r; |
| |
| r = neg2RootQ * cos((theta + 2 * PI) / 3) - adiv3; |
| if (!AlmostDequalUlps(s[0], r)) { |
| *roots++ = r; |
| } |
| r = neg2RootQ * cos((theta - 2 * PI) / 3) - adiv3; |
| if (!AlmostDequalUlps(s[0], r) && (roots - s == 1 || !AlmostDequalUlps(s[1], r))) { |
| *roots++ = r; |
| } |
| } else { // we have 1 real root |
| double sqrtR2MinusQ3 = sqrt(R2MinusQ3); |
| double A = fabs(R) + sqrtR2MinusQ3; |
| A = SkDCubeRoot(A); |
| if (R > 0) { |
| A = -A; |
| } |
| if (A != 0) { |
| A += Q / A; |
| } |
| r = A - adiv3; |
| *roots++ = r; |
| if (AlmostDequalUlps((double) R2, (double) Q3)) { |
| r = -A / 2 - adiv3; |
| if (!AlmostDequalUlps(s[0], r)) { |
| *roots++ = r; |
| } |
| } |
| } |
| return static_cast<int>(roots - s); |
| } |
| |
| // from http://www.cs.sunysb.edu/~qin/courses/geometry/4.pdf |
| // c(t) = a(1-t)^3 + 3bt(1-t)^2 + 3c(1-t)t^2 + dt^3 |
| // c'(t) = -3a(1-t)^2 + 3b((1-t)^2 - 2t(1-t)) + 3c(2t(1-t) - t^2) + 3dt^2 |
| // = 3(b-a)(1-t)^2 + 6(c-b)t(1-t) + 3(d-c)t^2 |
| static double derivative_at_t(const double* src, double t) { |
| double one_t = 1 - t; |
| double a = src[0]; |
| double b = src[2]; |
| double c = src[4]; |
| double d = src[6]; |
| return 3 * ((b - a) * one_t * one_t + 2 * (c - b) * t * one_t + (d - c) * t * t); |
| } |
| |
| // OPTIMIZE? compute t^2, t(1-t), and (1-t)^2 and pass them to another version of derivative at t? |
| SkDVector SkDCubic::dxdyAtT(double t) const { |
| SkDVector result = { derivative_at_t(&fPts[0].fX, t), derivative_at_t(&fPts[0].fY, t) }; |
| if (result.fX == 0 && result.fY == 0) { |
| if (t == 0) { |
| result = fPts[2] - fPts[0]; |
| } else if (t == 1) { |
| result = fPts[3] - fPts[1]; |
| } else { |
| // incomplete |
| SkDebugf("!c"); |
| } |
| if (result.fX == 0 && result.fY == 0 && zero_or_one(t)) { |
| result = fPts[3] - fPts[0]; |
| } |
| } |
| return result; |
| } |
| |
| // OPTIMIZE? share code with formulate_F1DotF2 |
| int SkDCubic::findInflections(double tValues[]) const { |
| double Ax = fPts[1].fX - fPts[0].fX; |
| double Ay = fPts[1].fY - fPts[0].fY; |
| double Bx = fPts[2].fX - 2 * fPts[1].fX + fPts[0].fX; |
| double By = fPts[2].fY - 2 * fPts[1].fY + fPts[0].fY; |
| double Cx = fPts[3].fX + 3 * (fPts[1].fX - fPts[2].fX) - fPts[0].fX; |
| double Cy = fPts[3].fY + 3 * (fPts[1].fY - fPts[2].fY) - fPts[0].fY; |
| return SkDQuad::RootsValidT(Bx * Cy - By * Cx, Ax * Cy - Ay * Cx, Ax * By - Ay * Bx, tValues); |
| } |
| |
| static void formulate_F1DotF2(const double src[], double coeff[4]) { |
| double a = src[2] - src[0]; |
| double b = src[4] - 2 * src[2] + src[0]; |
| double c = src[6] + 3 * (src[2] - src[4]) - src[0]; |
| coeff[0] = c * c; |
| coeff[1] = 3 * b * c; |
| coeff[2] = 2 * b * b + c * a; |
| coeff[3] = a * b; |
| } |
| |
| /** SkDCubic'(t) = At^2 + Bt + C, where |
| A = 3(-a + 3(b - c) + d) |
| B = 6(a - 2b + c) |
| C = 3(b - a) |
| Solve for t, keeping only those that fit between 0 < t < 1 |
| */ |
| int SkDCubic::FindExtrema(const double src[], double tValues[2]) { |
| // we divide A,B,C by 3 to simplify |
| double a = src[0]; |
| double b = src[2]; |
| double c = src[4]; |
| double d = src[6]; |
| double A = d - a + 3 * (b - c); |
| double B = 2 * (a - b - b + c); |
| double C = b - a; |
| |
| return SkDQuad::RootsValidT(A, B, C, tValues); |
| } |
| |
| /* from SkGeometry.cpp |
| Looking for F' dot F'' == 0 |
| |
| A = b - a |
| B = c - 2b + a |
| C = d - 3c + 3b - a |
| |
| F' = 3Ct^2 + 6Bt + 3A |
| F'' = 6Ct + 6B |
| |
| F' dot F'' -> CCt^3 + 3BCt^2 + (2BB + CA)t + AB |
| */ |
| int SkDCubic::findMaxCurvature(double tValues[]) const { |
| double coeffX[4], coeffY[4]; |
| int i; |
| formulate_F1DotF2(&fPts[0].fX, coeffX); |
| formulate_F1DotF2(&fPts[0].fY, coeffY); |
| for (i = 0; i < 4; i++) { |
| coeffX[i] = coeffX[i] + coeffY[i]; |
| } |
| return RootsValidT(coeffX[0], coeffX[1], coeffX[2], coeffX[3], tValues); |
| } |
| |
| SkDPoint SkDCubic::ptAtT(double t) const { |
| if (0 == t) { |
| return fPts[0]; |
| } |
| if (1 == t) { |
| return fPts[3]; |
| } |
| double one_t = 1 - t; |
| double one_t2 = one_t * one_t; |
| double a = one_t2 * one_t; |
| double b = 3 * one_t2 * t; |
| double t2 = t * t; |
| double c = 3 * one_t * t2; |
| double d = t2 * t; |
| SkDPoint result = {a * fPts[0].fX + b * fPts[1].fX + c * fPts[2].fX + d * fPts[3].fX, |
| a * fPts[0].fY + b * fPts[1].fY + c * fPts[2].fY + d * fPts[3].fY}; |
| return result; |
| } |
| |
| /* |
| Given a cubic c, t1, and t2, find a small cubic segment. |
| |
| The new cubic is defined as points A, B, C, and D, where |
| s1 = 1 - t1 |
| s2 = 1 - t2 |
| A = c[0]*s1*s1*s1 + 3*c[1]*s1*s1*t1 + 3*c[2]*s1*t1*t1 + c[3]*t1*t1*t1 |
| D = c[0]*s2*s2*s2 + 3*c[1]*s2*s2*t2 + 3*c[2]*s2*t2*t2 + c[3]*t2*t2*t2 |
| |
| We don't have B or C. So We define two equations to isolate them. |
| First, compute two reference T values 1/3 and 2/3 from t1 to t2: |
| |
| c(at (2*t1 + t2)/3) == E |
| c(at (t1 + 2*t2)/3) == F |
| |
| Next, compute where those values must be if we know the values of B and C: |
| |
| _12 = A*2/3 + B*1/3 |
| 12_ = A*1/3 + B*2/3 |
| _23 = B*2/3 + C*1/3 |
| 23_ = B*1/3 + C*2/3 |
| _34 = C*2/3 + D*1/3 |
| 34_ = C*1/3 + D*2/3 |
| _123 = (A*2/3 + B*1/3)*2/3 + (B*2/3 + C*1/3)*1/3 = A*4/9 + B*4/9 + C*1/9 |
| 123_ = (A*1/3 + B*2/3)*1/3 + (B*1/3 + C*2/3)*2/3 = A*1/9 + B*4/9 + C*4/9 |
| _234 = (B*2/3 + C*1/3)*2/3 + (C*2/3 + D*1/3)*1/3 = B*4/9 + C*4/9 + D*1/9 |
| 234_ = (B*1/3 + C*2/3)*1/3 + (C*1/3 + D*2/3)*2/3 = B*1/9 + C*4/9 + D*4/9 |
| _1234 = (A*4/9 + B*4/9 + C*1/9)*2/3 + (B*4/9 + C*4/9 + D*1/9)*1/3 |
| = A*8/27 + B*12/27 + C*6/27 + D*1/27 |
| = E |
| 1234_ = (A*1/9 + B*4/9 + C*4/9)*1/3 + (B*1/9 + C*4/9 + D*4/9)*2/3 |
| = A*1/27 + B*6/27 + C*12/27 + D*8/27 |
| = F |
| E*27 = A*8 + B*12 + C*6 + D |
| F*27 = A + B*6 + C*12 + D*8 |
| |
| Group the known values on one side: |
| |
| M = E*27 - A*8 - D = B*12 + C* 6 |
| N = F*27 - A - D*8 = B* 6 + C*12 |
| M*2 - N = B*18 |
| N*2 - M = C*18 |
| B = (M*2 - N)/18 |
| C = (N*2 - M)/18 |
| */ |
| |
| static double interp_cubic_coords(const double* src, double t) { |
| double ab = SkDInterp(src[0], src[2], t); |
| double bc = SkDInterp(src[2], src[4], t); |
| double cd = SkDInterp(src[4], src[6], t); |
| double abc = SkDInterp(ab, bc, t); |
| double bcd = SkDInterp(bc, cd, t); |
| double abcd = SkDInterp(abc, bcd, t); |
| return abcd; |
| } |
| |
| SkDCubic SkDCubic::subDivide(double t1, double t2) const { |
| if (t1 == 0 || t2 == 1) { |
| if (t1 == 0 && t2 == 1) { |
| return *this; |
| } |
| SkDCubicPair pair = chopAt(t1 == 0 ? t2 : t1); |
| SkDCubic dst = t1 == 0 ? pair.first() : pair.second(); |
| return dst; |
| } |
| SkDCubic dst; |
| double ax = dst[0].fX = interp_cubic_coords(&fPts[0].fX, t1); |
| double ay = dst[0].fY = interp_cubic_coords(&fPts[0].fY, t1); |
| double ex = interp_cubic_coords(&fPts[0].fX, (t1*2+t2)/3); |
| double ey = interp_cubic_coords(&fPts[0].fY, (t1*2+t2)/3); |
| double fx = interp_cubic_coords(&fPts[0].fX, (t1+t2*2)/3); |
| double fy = interp_cubic_coords(&fPts[0].fY, (t1+t2*2)/3); |
| double dx = dst[3].fX = interp_cubic_coords(&fPts[0].fX, t2); |
| double dy = dst[3].fY = interp_cubic_coords(&fPts[0].fY, t2); |
| double mx = ex * 27 - ax * 8 - dx; |
| double my = ey * 27 - ay * 8 - dy; |
| double nx = fx * 27 - ax - dx * 8; |
| double ny = fy * 27 - ay - dy * 8; |
| /* bx = */ dst[1].fX = (mx * 2 - nx) / 18; |
| /* by = */ dst[1].fY = (my * 2 - ny) / 18; |
| /* cx = */ dst[2].fX = (nx * 2 - mx) / 18; |
| /* cy = */ dst[2].fY = (ny * 2 - my) / 18; |
| // FIXME: call align() ? |
| return dst; |
| } |
| |
| void SkDCubic::subDivide(const SkDPoint& a, const SkDPoint& d, |
| double t1, double t2, SkDPoint dst[2]) const { |
| SkASSERT(t1 != t2); |
| // this approach assumes that the control points computed directly are accurate enough |
| SkDCubic sub = subDivide(t1, t2); |
| dst[0] = sub[1] + (a - sub[0]); |
| dst[1] = sub[2] + (d - sub[3]); |
| if (t1 == 0 || t2 == 0) { |
| align(0, 1, t1 == 0 ? &dst[0] : &dst[1]); |
| } |
| if (t1 == 1 || t2 == 1) { |
| align(3, 2, t1 == 1 ? &dst[0] : &dst[1]); |
| } |
| if (AlmostBequalUlps(dst[0].fX, a.fX)) { |
| dst[0].fX = a.fX; |
| } |
| if (AlmostBequalUlps(dst[0].fY, a.fY)) { |
| dst[0].fY = a.fY; |
| } |
| if (AlmostBequalUlps(dst[1].fX, d.fX)) { |
| dst[1].fX = d.fX; |
| } |
| if (AlmostBequalUlps(dst[1].fY, d.fY)) { |
| dst[1].fY = d.fY; |
| } |
| } |
| |
| double SkDCubic::top(const SkDCubic& dCurve, double startT, double endT, SkDPoint*topPt) const { |
| double extremeTs[2]; |
| double topT = -1; |
| int roots = SkDCubic::FindExtrema(&fPts[0].fY, extremeTs); |
| for (int index = 0; index < roots; ++index) { |
| double t = startT + (endT - startT) * extremeTs[index]; |
| SkDPoint mid = dCurve.ptAtT(t); |
| if (topPt->fY > mid.fY || (topPt->fY == mid.fY && topPt->fX > mid.fX)) { |
| topT = t; |
| *topPt = mid; |
| } |
| } |
| return topT; |
| } |