diff options
Diffstat (limited to 'pathops/SkPathOpsCubic.cpp')
-rw-r--r-- | pathops/SkPathOpsCubic.cpp | 463 |
1 files changed, 463 insertions, 0 deletions
diff --git a/pathops/SkPathOpsCubic.cpp b/pathops/SkPathOpsCubic.cpp new file mode 100644 index 00000000..674213c3 --- /dev/null +++ b/pathops/SkPathOpsCubic.cpp @@ -0,0 +1,463 @@ +/* + * 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 "SkLineParameters.h" +#include "SkPathOpsCubic.h" +#include "SkPathOpsLine.h" +#include "SkPathOpsQuad.h" +#include "SkPathOpsRect.h" + +const int SkDCubic::gPrecisionUnit = 256; // FIXME: test different values in test framework + +// 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; +} + +bool SkDCubic::clockwise() const { + double sum = (fPts[0].fX - fPts[3].fX) * (fPts[0].fY + fPts[3].fY); + for (int idx = 0; idx < 3; ++idx) { + sum += (fPts[idx + 1].fX - fPts[idx].fX) * (fPts[idx + 1].fY + fPts[idx].fY); + } + return sum <= 0; +} + +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::controlsContainedByEnds() const { + SkDVector startTan = fPts[1] - fPts[0]; + if (startTan.fX == 0 && startTan.fY == 0) { + startTan = fPts[2] - fPts[0]; + } + SkDVector endTan = fPts[2] - fPts[3]; + if (endTan.fX == 0 && endTan.fY == 0) { + endTan = fPts[1] - fPts[3]; + } + if (startTan.dot(endTan) >= 0) { + return false; + } + SkDLine startEdge = {{fPts[0], fPts[0]}}; + startEdge[1].fX -= startTan.fY; + startEdge[1].fY += startTan.fX; + SkDLine endEdge = {{fPts[3], fPts[3]}}; + endEdge[1].fX -= endTan.fY; + endEdge[1].fY += endTan.fX; + double leftStart1 = startEdge.isLeft(fPts[1]); + if (leftStart1 * startEdge.isLeft(fPts[2]) < 0) { + return false; + } + double leftEnd1 = endEdge.isLeft(fPts[1]); + if (leftEnd1 * endEdge.isLeft(fPts[2]) < 0) { + return false; + } + return leftStart1 * leftEnd1 >= 0; +} + +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)); +} + +bool SkDCubic::isLinear(int startIndex, int endIndex) const { + SkLineParameters lineParameters; + lineParameters.cubicEndPoints(*this, startIndex, endIndex); + // FIXME: maybe it's possible to avoid this and compare non-normalized + lineParameters.normalize(); + double distance = lineParameters.controlPtDistance(*this, 1); + if (!approximately_zero(distance)) { + return false; + } + distance = lineParameters.controlPtDistance(*this, 2); + return approximately_zero(distance); +} + +bool SkDCubic::monotonicInY() const { + return between(fPts[0].fY, fPts[1].fY, fPts[3].fY) + && between(fPts[0].fY, fPts[2].fY, fPts[3].fY); +} + +bool SkDCubic::serpentine() const { + if (!controlsContainedByEnds()) { + return false; + } + double wiggle = (fPts[0].fX - fPts[2].fX) * (fPts[0].fY + fPts[2].fY); + for (int idx = 0; idx < 2; ++idx) { + wiggle += (fPts[idx + 1].fX - fPts[idx].fX) * (fPts[idx + 1].fY + fPts[idx].fY); + } + double waggle = (fPts[1].fX - fPts[3].fX) * (fPts[1].fY + fPts[3].fY); + for (int idx = 1; idx < 3; ++idx) { + waggle += (fPts[idx + 1].fX - fPts[idx].fX) * (fPts[idx + 1].fY + fPts[idx].fY); + } + return wiggle * waggle < 0; +} + +// 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); + 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); + mathematica_ize(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 (AlmostEqualUlps(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 + double theta = acos(R / sqrt(Q3)); + double neg2RootQ = -2 * sqrt(Q); + + r = neg2RootQ * cos(theta / 3) - adiv3; + *roots++ = r; + + r = neg2RootQ * cos((theta + 2 * PI) / 3) - adiv3; + if (!AlmostEqualUlps(s[0], r)) { + *roots++ = r; + } + r = neg2RootQ * cos((theta - 2 * PI) / 3) - adiv3; + if (!AlmostEqualUlps(s[0], r) && (roots - s == 1 || !AlmostEqualUlps(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 (AlmostEqualUlps(R2, Q3)) { + r = -A / 2 - adiv3; + if (!AlmostEqualUlps(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) }; + 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(double a, double b, double c, double d, double tValues[2]) { + // we divide A,B,C by 3 to simplify + 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::top(double startT, double endT) const { + SkDCubic sub = subDivide(startT, endT); + SkDPoint topPt = sub[0]; + if (topPt.fY > sub[3].fY || (topPt.fY == sub[3].fY && topPt.fX > sub[3].fX)) { + topPt = sub[3]; + } + double extremeTs[2]; + if (!sub.monotonicInY()) { + int roots = FindExtrema(sub[0].fY, sub[1].fY, sub[2].fY, sub[3].fY, extremeTs); + for (int index = 0; index < roots; ++index) { + double t = startT + (endT - startT) * extremeTs[index]; + SkDPoint mid = xyAtT(t); + if (topPt.fY > mid.fY || (topPt.fY == mid.fY && topPt.fX > mid.fX)) { + topPt = mid; + } + } + } + return topPt; +} + +SkDPoint SkDCubic::xyAtT(double t) const { + 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) { + return *this; + } + 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; + return dst; +} + +void SkDCubic::subDivide(const SkDPoint& a, const SkDPoint& d, + double t1, double t2, SkDPoint dst[2]) const { + 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 mx = ex * 27 - a.fX * 8 - d.fX; + double my = ey * 27 - a.fY * 8 - d.fY; + double nx = fx * 27 - a.fX - d.fX * 8; + double ny = fy * 27 - a.fY - d.fY * 8; + /* bx = */ dst[0].fX = (mx * 2 - nx) / 18; + /* by = */ dst[0].fY = (my * 2 - ny) / 18; + /* cx = */ dst[1].fX = (nx * 2 - mx) / 18; + /* cy = */ dst[1].fY = (ny * 2 - my) / 18; +} + +/* 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; +} |