Add intersections for path ops

src/pathops/SkDCubicLineIntersection.cpp

Index: src/pathops/SkDCubicLineIntersection.cpp
===================================================================
--- src/pathops/SkDCubicLineIntersection.cpp	(revision 0)
+++ src/pathops/SkDCubicLineIntersection.cpp	(revision 0)
@@ -0,0 +1,261 @@
+/*
+ * 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 "SkIntersections.h"
+#include "SkPathOpsCubic.h"
+#include "SkPathOpsLine.h"
+
+/*
+Find the interection of a line and cubic by solving for valid t values.
+
+Analogous to line-quadratic intersection, solve line-cubic intersection by
+representing the cubic as:
+  x = a(1-t)^3 + 2b(1-t)^2t + c(1-t)t^2 + dt^3
+  y = e(1-t)^3 + 2f(1-t)^2t + g(1-t)t^2 + ht^3
+and the line as:
+  y = i*x + j  (if the line is more horizontal)
+or:
+  x = i*y + j  (if the line is more vertical)
+
+Then using Mathematica, solve for the values of t where the cubic intersects the
+line:
+
+  (in) Resultant[
+        a*(1 - t)^3 + 3*b*(1 - t)^2*t + 3*c*(1 - t)*t^2 + d*t^3 - x,
+        e*(1 - t)^3 + 3*f*(1 - t)^2*t + 3*g*(1 - t)*t^2 + h*t^3 - i*x - j, x]
+  (out) -e     +   j     +
+       3 e t   - 3 f t   -
+       3 e t^2 + 6 f t^2 - 3 g t^2 +
+         e t^3 - 3 f t^3 + 3 g t^3 - h t^3 +
+     i ( a     -
+       3 a t + 3 b t +
+       3 a t^2 - 6 b t^2 + 3 c t^2 -
+         a t^3 + 3 b t^3 - 3 c t^3 + d t^3 )
+
+if i goes to infinity, we can rewrite the line in terms of x. Mathematica:
+
+  (in) Resultant[
+        a*(1 - t)^3 + 3*b*(1 - t)^2*t + 3*c*(1 - t)*t^2 + d*t^3 - i*y - j,
+        e*(1 - t)^3 + 3*f*(1 - t)^2*t + 3*g*(1 - t)*t^2 + h*t^3 - y,       y]
+  (out)  a     -   j     -
+       3 a t   + 3 b t   +
+       3 a t^2 - 6 b t^2 + 3 c t^2 -
+         a t^3 + 3 b t^3 - 3 c t^3 + d t^3 -
+     i ( e     -
+       3 e t   + 3 f t   +
+       3 e t^2 - 6 f t^2 + 3 g t^2 -
+         e t^3 + 3 f t^3 - 3 g t^3 + h t^3 )
+
+Solving this with Mathematica produces an expression with hundreds of terms;
+instead, use Numeric Solutions recipe to solve the cubic.
+
+The near-horizontal case, in terms of:  Ax^3 + Bx^2 + Cx + D == 0
+    A =   (-(-e + 3*f - 3*g + h) + i*(-a + 3*b - 3*c + d)     )
+    B = 3*(-( e - 2*f +   g    ) + i*( a - 2*b +   c    )     )
+    C = 3*(-(-e +   f          ) + i*(-a +   b          )     )
+    D =   (-( e                ) + i*( a                ) + j )
+
+The near-vertical case, in terms of:  Ax^3 + Bx^2 + Cx + D == 0
+    A =   ( (-a + 3*b - 3*c + d) - i*(-e + 3*f - 3*g + h)     )
+    B = 3*( ( a - 2*b +   c    ) - i*( e - 2*f +   g    )     )
+    C = 3*( (-a +   b          ) - i*(-e +   f          )     )
+    D =   ( ( a                ) - i*( e                ) - j )
+
+For horizontal lines:
+(in) Resultant[
+      a*(1 - t)^3 + 3*b*(1 - t)^2*t + 3*c*(1 - t)*t^2 + d*t^3 - j,
+      e*(1 - t)^3 + 3*f*(1 - t)^2*t + 3*g*(1 - t)*t^2 + h*t^3 - y, y]
+(out)  e     -   j     -
+     3 e t   + 3 f t   +
+     3 e t^2 - 6 f t^2 + 3 g t^2 -
+       e t^3 + 3 f t^3 - 3 g t^3 + h t^3
+ */
+
+class LineCubicIntersections {
+public:
+
+LineCubicIntersections(const SkDCubic& c, const SkDLine& l, SkIntersections& i)
+    : cubic(c)
+    , line(l)
+    , intersections(i) {
+}
+
+// see parallel routine in line quadratic intersections
+int intersectRay(double roots[3]) {
+    double adj = line[1].fX - line[0].fX;
+    double opp = line[1].fY - line[0].fY;
+    SkDCubic r;
+    for (int n = 0; n < 4; ++n) {
+        r[n].fX = (cubic[n].fY - line[0].fY) * adj - (cubic[n].fX - line[0].fX) * opp;
+    }
+    double A, B, C, D;
+    SkDCubic::Coefficients(&r[0].fX, &A, &B, &C, &D);
+    return SkDCubic::RootsValidT(A, B, C, D, roots);
+}
+
+int intersect() {
+    addEndPoints();
+    double rootVals[3];
+    int roots = intersectRay(rootVals);
+    for (int index = 0; index < roots; ++index) {
+        double cubicT = rootVals[index];
+        double lineT = findLineT(cubicT);
+        if (pinTs(&cubicT, &lineT)) {
+            SkDPoint pt = line.xyAtT(lineT);
+            intersections.insert(cubicT, lineT, pt);
+        }
+    }
+    return intersections.used();
+}
+
+int horizontalIntersect(double axisIntercept, double roots[3]) {
+    double A, B, C, D;
+    SkDCubic::Coefficients(&cubic[0].fY, &A, &B, &C, &D);
+    D -= axisIntercept;
+    return SkDCubic::RootsValidT(A, B, C, D, roots);
+}
+
+int horizontalIntersect(double axisIntercept, double left, double right, bool flipped) {
+    addHorizontalEndPoints(left, right, axisIntercept);
+    double rootVals[3];
+    int roots = horizontalIntersect(axisIntercept, rootVals);
+    for (int index = 0; index < roots; ++index) {
+        double cubicT = rootVals[index];
+        SkDPoint pt = cubic.xyAtT(cubicT);
+        double lineT = (pt.fX - left) / (right - left);
+        if (pinTs(&cubicT, &lineT)) {
+            intersections.insert(cubicT, lineT, pt);
+        }
+    }
+    if (flipped) {
+        intersections.flip();
+    }
+    return intersections.used();
+}
+
+int verticalIntersect(double axisIntercept, double roots[3]) {
+    double A, B, C, D;
+    SkDCubic::Coefficients(&cubic[0].fX, &A, &B, &C, &D);
+    D -= axisIntercept;
+    return SkDCubic::RootsValidT(A, B, C, D, roots);
+}
+
+int verticalIntersect(double axisIntercept, double top, double bottom, bool flipped) {
+    addVerticalEndPoints(top, bottom, axisIntercept);
+    double rootVals[3];
+    int roots = verticalIntersect(axisIntercept, rootVals);
+    for (int index = 0; index < roots; ++index) {
+        double cubicT = rootVals[index];
+        SkDPoint pt = cubic.xyAtT(cubicT);
+        double lineT = (pt.fY - top) / (bottom - top);
+        if (pinTs(&cubicT, &lineT)) {
+            intersections.insert(cubicT, lineT, pt);
+        }
+    }
+    if (flipped) {
+        intersections.flip();
+    }
+    return intersections.used();
+}
+
+protected:
+
+void addEndPoints() {
+    for (int cIndex = 0; cIndex < 4; cIndex += 3) {
+        for (int lIndex = 0; lIndex < 2; lIndex++) {
+            if (cubic[cIndex] == line[lIndex]) {
+                intersections.insert(cIndex >> 1, lIndex, line[lIndex]);
+            }
+        }
+    }
+}
+
+void addHorizontalEndPoints(double left, double right, double y) {
+    for (int cIndex = 0; cIndex < 4; cIndex += 3) {
+        if (cubic[cIndex].fY != y) {
+            continue;
+        }
+        if (cubic[cIndex].fX == left) {
+            intersections.insert(cIndex >> 1, 0, cubic[cIndex]);
+        }
+        if (cubic[cIndex].fX == right) {
+            intersections.insert(cIndex >> 1, 1, cubic[cIndex]);
+        }
+    }
+}
+
+void addVerticalEndPoints(double top, double bottom, double x) {
+    for (int cIndex = 0; cIndex < 4; cIndex += 3) {
+        if (cubic[cIndex].fX != x) {
+            continue;
+        }
+        if (cubic[cIndex].fY == top) {
+            intersections.insert(cIndex >> 1, 0, cubic[cIndex]);
+        }
+        if (cubic[cIndex].fY == bottom) {
+            intersections.insert(cIndex >> 1, 1, cubic[cIndex]);
+        }
+    }
+}
+
+double findLineT(double t) {
+    SkDPoint xy = cubic.xyAtT(t);
+    double dx = line[1].fX - line[0].fX;
+    double dy = line[1].fY - line[0].fY;
+    if (fabs(dx) > fabs(dy)) {
+        return (xy.fX - line[0].fX) / dx;
+    }
+    return (xy.fY - line[0].fY) / dy;
+}
+
+static bool pinTs(double* cubicT, double* lineT) {
+    if (!approximately_one_or_less(*lineT)) {
+        return false;
+    }
+    if (!approximately_zero_or_more(*lineT)) {
+        return false;
+    }
+    if (precisely_less_than_zero(*cubicT)) {
+        *cubicT = 0;
+    } else if (precisely_greater_than_one(*cubicT)) {
+        *cubicT = 1;
+    }
+    if (precisely_less_than_zero(*lineT)) {
+        *lineT = 0;
+    } else if (precisely_greater_than_one(*lineT)) {
+        *lineT = 1;
+    }
+    return true;
+}
+
+private:
+
+const SkDCubic& cubic;
+const SkDLine& line;
+SkIntersections& intersections;
+};
+
+int SkIntersections::horizontal(const SkDCubic& cubic, double left, double right, double y,
+        bool flipped) {
+    LineCubicIntersections c(cubic, *(static_cast<SkDLine*>(0)), *this);
+    return c.horizontalIntersect(y, left, right, flipped);
+}
+
+int SkIntersections::vertical(const SkDCubic& cubic, double top, double bottom, double x,
+        bool flipped) {
+    LineCubicIntersections c(cubic, *(static_cast<SkDLine*>(0)), *this);
+    return c.verticalIntersect(x, top, bottom, flipped);
+}
+
+int SkIntersections::intersect(const SkDCubic& cubic, const SkDLine& line) {
+    LineCubicIntersections c(cubic, line, *this);
+    return c.intersect();
+}
+
+int SkIntersections::intersectRay(const SkDCubic& cubic, const SkDLine& line) {
+    LineCubicIntersections c(cubic, line, *this);
+    return c.intersectRay(fT[0]);
+}