Add base types for path ops

src/pathops/SkPathOpsQuad.cpp

Index: src/pathops/SkPathOpsQuad.cpp
===================================================================
--- src/pathops/SkPathOpsQuad.cpp	(revision 0)
+++ src/pathops/SkPathOpsQuad.cpp	(revision 0)
@@ -0,0 +1,293 @@
+/*
+ * 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 "SkPathOpsQuad.h"
+#include "SkPathOpsTriangle.h"
+
+// from http://blog.gludion.com/2009/08/distance-to-quadratic-bezier-curve.html
+// (currently only used by testing)
+double SkDQuad::nearestT(const SkDPoint& pt) const {
+    SkDVector pos = fPts[0] - pt;
+    // search points P of bezier curve with PM.(dP / dt) = 0
+    // a calculus leads to a 3d degree equation :
+    SkDVector A = fPts[1] - fPts[0];
+    SkDVector B = fPts[2] - fPts[1];
+    B -= A;
+    double a = B.dot(B);
+    double b = 3 * A.dot(B);
+    double c = 2 * A.dot(A) + pos.dot(B);
+    double d = pos.dot(A);
+    double ts[3];
+    int roots = SkDCubic::RootsValidT(a, b, c, d, ts);
+    double d0 = pt.distanceSquared(fPts[0]);
+    double d2 = pt.distanceSquared(fPts[2]);
+    double distMin = SkTMin(d0, d2);
+    int bestIndex = -1;
+    for (int index = 0; index < roots; ++index) {
+        SkDPoint onQuad = xyAtT(ts[index]);
+        double dist = pt.distanceSquared(onQuad);
+        if (distMin > dist) {
+            distMin = dist;
+            bestIndex = index;
+        }
+    }
+    if (bestIndex >= 0) {
+        return ts[bestIndex];
+    }
+    return d0 < d2 ? 0 : 1;
+}
+
+bool SkDQuad::pointInHull(const SkDPoint& pt) const {
+    return ((const SkDTriangle&) fPts).contains(pt);
+}
+
+SkDPoint SkDQuad::top(double startT, double endT) const {
+    SkDQuad sub = subDivide(startT, endT);
+    SkDPoint topPt = sub[0];
+    if (topPt.fY > sub[2].fY || (topPt.fY == sub[2].fY && topPt.fX > sub[2].fX)) {
+        topPt = sub[2];
+    }
+    if (!between(sub[0].fY, sub[1].fY, sub[2].fY)) {
+        double extremeT;
+        if (FindExtrema(sub[0].fY, sub[1].fY, sub[2].fY, &extremeT)) {
+            extremeT = startT + (endT - startT) * extremeT;
+            SkDPoint test = xyAtT(extremeT);
+            if (topPt.fY > test.fY || (topPt.fY == test.fY && topPt.fX > test.fX)) {
+                topPt = test;
+            }
+        }
+    }
+    return topPt;
+}
+
+int SkDQuad::AddValidTs(double s[], int realRoots, double* t) {
+    int foundRoots = 0;
+    for (int index = 0; index < realRoots; ++index) {
+        double tValue = s[index];
+        if (approximately_zero_or_more(tValue) && approximately_one_or_less(tValue)) {
+            if (approximately_less_than_zero(tValue)) {
+                tValue = 0;
+            } else if (approximately_greater_than_one(tValue)) {
+                tValue = 1;
+            }
+            for (int idx2 = 0; idx2 < foundRoots; ++idx2) {
+                if (approximately_equal(t[idx2], tValue)) {
+                    goto nextRoot;
+                }
+            }
+            t[foundRoots++] = tValue;
+        }
+ nextRoot:
+        {}
+    }
+    return foundRoots;
+}
+
+// note: caller expects multiple results to be sorted smaller first
+// note: http://en.wikipedia.org/wiki/Loss_of_significance has an interesting
+//  analysis of the quadratic equation, suggesting why the following looks at
+//  the sign of B -- and further suggesting that the greatest loss of precision
+//  is in b squared less two a c
+int SkDQuad::RootsValidT(double A, double B, double C, double t[2]) {
+    double s[2];
+    int realRoots = RootsReal(A, B, C, s);
+    int foundRoots = AddValidTs(s, realRoots, t);
+    return foundRoots;
+}
+
+/*
+Numeric Solutions (5.6) suggests to solve the quadratic by computing
+       Q = -1/2(B + sgn(B)Sqrt(B^2 - 4 A C))
+and using the roots
+      t1 = Q / A
+      t2 = C / Q
+*/
+// this does not discard real roots <= 0 or >= 1
+int SkDQuad::RootsReal(const double A, const double B, const double C, double s[2]) {
+    const double p = B / (2 * A);
+    const double q = C / A;
+    if (approximately_zero(A) && (approximately_zero_inverse(p) || approximately_zero_inverse(q))) {
+        if (approximately_zero(B)) {
+            s[0] = 0;
+            return C == 0;
+        }
+        s[0] = -C / B;
+        return 1;
+    }
+    /* normal form: x^2 + px + q = 0 */
+    const double p2 = p * p;
+    if (!AlmostEqualUlps(p2, q) && p2 < q) {
+        return 0;
+    }
+    double sqrt_D = 0;
+    if (p2 > q) {
+        sqrt_D = sqrt(p2 - q);
+    }
+    s[0] = sqrt_D - p;
+    s[1] = -sqrt_D - p;
+    return 1 + !AlmostEqualUlps(s[0], s[1]);
+}
+
+bool SkDQuad::isLinear(int startIndex, int endIndex) const {
+    SkLineParameters lineParameters;
+    lineParameters.quadEndPoints(*this, startIndex, endIndex);
+    // FIXME: maybe it's possible to avoid this and compare non-normalized
+    lineParameters.normalize();
+    double distance = lineParameters.controlPtDistance(*this);
+    return approximately_zero(distance);
+}
+
+SkDCubic SkDQuad::toCubic() const {
+    SkDCubic cubic;
+    cubic[0] = fPts[0];
+    cubic[2] = fPts[1];
+    cubic[3] = fPts[2];
+    cubic[1].fX = (cubic[0].fX + cubic[2].fX * 2) / 3;
+    cubic[1].fY = (cubic[0].fY + cubic[2].fY * 2) / 3;
+    cubic[2].fX = (cubic[3].fX + cubic[2].fX * 2) / 3;
+    cubic[2].fY = (cubic[3].fY + cubic[2].fY * 2) / 3;
+    return cubic;
+}
+
+SkDVector SkDQuad::dxdyAtT(double t) const {
+    double a = t - 1;
+    double b = 1 - 2 * t;
+    double c = t;
+    SkDVector result = { a * fPts[0].fX + b * fPts[1].fX + c * fPts[2].fX,
+            a * fPts[0].fY + b * fPts[1].fY + c * fPts[2].fY };
+    return result;
+}
+
+SkDPoint SkDQuad::xyAtT(double t) const {
+    double one_t = 1 - t;
+    double a = one_t * one_t;
+    double b = 2 * one_t * t;
+    double c = t * t;
+    SkDPoint result = { a * fPts[0].fX + b * fPts[1].fX + c * fPts[2].fX,
+            a * fPts[0].fY + b * fPts[1].fY + c * fPts[2].fY };
+    return result;
+}
+
+/*
+Given a quadratic q, t1, and t2, find a small quadratic segment.
+
+The new quadratic is defined by A, B, and C, where
+ A = c[0]*(1 - t1)*(1 - t1) + 2*c[1]*t1*(1 - t1) + c[2]*t1*t1
+ C = c[3]*(1 - t1)*(1 - t1) + 2*c[2]*t1*(1 - t1) + c[1]*t1*t1
+
+To find B, compute the point halfway between t1 and t2:
+
+q(at (t1 + t2)/2) == D
+
+Next, compute where D must be if we know the value of B:
+
+_12 = A/2 + B/2
+12_ = B/2 + C/2
+123 = A/4 + B/2 + C/4
+    = D
+
+Group the known values on one side:
+
+B   = D*2 - A/2 - C/2
+*/
+
+static double interp_quad_coords(const double* src, double t) {
+    double ab = SkDInterp(src[0], src[2], t);
+    double bc = SkDInterp(src[2], src[4], t);
+    double abc = SkDInterp(ab, bc, t);
+    return abc;
+}
+
+bool SkDQuad::monotonicInY() const {
+    return between(fPts[0].fY, fPts[1].fY, fPts[2].fY);
+}
+
+SkDQuad SkDQuad::subDivide(double t1, double t2) const {
+    SkDQuad dst;
+    double ax = dst[0].fX = interp_quad_coords(&fPts[0].fX, t1);
+    double ay = dst[0].fY = interp_quad_coords(&fPts[0].fY, t1);
+    double dx = interp_quad_coords(&fPts[0].fX, (t1 + t2) / 2);
+    double dy = interp_quad_coords(&fPts[0].fY, (t1 + t2) / 2);
+    double cx = dst[2].fX = interp_quad_coords(&fPts[0].fX, t2);
+    double cy = dst[2].fY = interp_quad_coords(&fPts[0].fY, t2);
+    /* bx = */ dst[1].fX = 2*dx - (ax + cx)/2;
+    /* by = */ dst[1].fY = 2*dy - (ay + cy)/2;
+    return dst;
+}
+
+SkDPoint SkDQuad::subDivide(const SkDPoint& a, const SkDPoint& c, double t1, double t2) const {
+    SkDPoint b;
+    double dx = interp_quad_coords(&fPts[0].fX, (t1 + t2) / 2);
+    double dy = interp_quad_coords(&fPts[0].fY, (t1 + t2) / 2);
+    b.fX = 2 * dx - (a.fX + c.fX) / 2;
+    b.fY = 2 * dy - (a.fY + c.fY) / 2;
+    return b;
+}
+
+/* classic one t subdivision */
+static void interp_quad_coords(const double* src, double* dst, double t) {
+    double ab = SkDInterp(src[0], src[2], t);
+    double bc = SkDInterp(src[2], src[4], t);
+    dst[0] = src[0];
+    dst[2] = ab;
+    dst[4] = SkDInterp(ab, bc, t);
+    dst[6] = bc;
+    dst[8] = src[4];
+}
+
+SkDQuadPair SkDQuad::chopAt(double t) const
+{
+    SkDQuadPair dst;
+    interp_quad_coords(&fPts[0].fX, &dst.pts[0].fX, t);
+    interp_quad_coords(&fPts[0].fY, &dst.pts[0].fY, t);
+    return dst;
+}
+
+static int valid_unit_divide(double numer, double denom, double* ratio)
+{
+    if (numer < 0) {
+        numer = -numer;
+        denom = -denom;
+    }
+    if (denom == 0 || numer == 0 || numer >= denom) {
+        return 0;
+    }
+    double r = numer / denom;
+    if (r == 0) {  // catch underflow if numer <<<< denom
+        return 0;
+    }
+    *ratio = r;
+    return 1;
+}
+
+/** Quad'(t) = At + B, where
+    A = 2(a - 2b + c)
+    B = 2(b - a)
+    Solve for t, only if it fits between 0 < t < 1
+*/
+int SkDQuad::FindExtrema(double a, double b, double c, double tValue[1]) {
+    /*  At + B == 0
+        t = -B / A
+    */
+    return valid_unit_divide(a - b, a - b - b + c, tValue);
+}
+
+/* Parameterization form, given A*t*t + 2*B*t*(1-t) + C*(1-t)*(1-t)
+ *
+ * a = A - 2*B +   C
+ * b =     2*B - 2*C
+ * c =             C
+ */
+void SkDQuad::SetABC(const double* quad, double* a, double* b, double* c) {
+    *a = quad[0];      // a = A
+    *b = 2 * quad[2];  // b =     2*B
+    *c = quad[4];      // c =             C
+    *b -= *c;          // b =     2*B -   C
+    *a -= *b;          // a = A - 2*B +   C
+    *b -= *c;          // b =     2*B - 2*C
+}