Add base types for path ops

Paths contain lines, quads, and cubics, which are
collectively curves.

To work with path intersections, intermediary curves
are constructed. For now, those intermediates use
doubles to guarantee sufficient precision.

The DVector, DPoint, DLine, DQuad, and DCubic
structs encapsulate these intermediate curves.

The DRect and DTriangle structs are created to
describe intersectable areas of interest.

The Bounds struct inherits from SkRect to create
a SkScalar-based rectangle that intersects shared
edges.

This also includes common math equalities and
debugging that the remainder of path ops builds on,
as well as a temporary top-level interface in
include/pathops/SkPathOps.h.
Review URL: https://codereview.chromium.org/12827020

git-svn-id: http://skia.googlecode.com/svn/trunk@8551 2bbb7eff-a529-9590-31e7-b0007b416f81

1 files changed, 463 insertions, 0 deletions

diff --git a/src/pathops/SkPathOpsCubic.cpp b/src/pathops/SkPathOpsCubic.cpp
new file mode 100644
index 0000000000..674213c383
--- /dev/null
+++ b/src/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;
+}