Add base types for path ops

src/pathops/SkPathOpsCubic.cpp

+/*
+ * 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];
+    bzero(str, sizeof(str));
+    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;
+}