shape ops work in progress

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

15 files changed, 292 insertions, 77 deletions

diff --git a/experimental/Intersection/CubicIntersection.cpp b/experimental/Intersection/CubicIntersection.cpp
index be3e088b25..2c6cc3f577 100644
--- a/experimental/Intersection/CubicIntersection.cpp
+++ b/experimental/Intersection/CubicIntersection.cpp
@@ -161,39 +161,64 @@ bool intersect(const Cubic& c1, const Cubic& c2, Intersections& i) {
 
 #include "CubicUtilities.h"
 
+// FIXME: ? if needed, compute the error term from the tangent length
+static double computeDelta(const Cubic& cubic, double t, double scale) {
+    double attempt = scale / precisionUnit;
+#if SK_DEBUG
+    double precision = calcPrecision(cubic);
+    _Point dxy;
+    dxdy_at_t(cubic, t, dxy);
+    _Point p1, p2;
+    xy_at_t(cubic, std::max(t - attempt, 0.), p1.x, p1.y);
+    xy_at_t(cubic, std::min(t + attempt, 1.), p2.x, p2.y);
+    double dx = p1.x - p2.x;
+    double dy = p1.y - p2.y;
+    double distSq = dx * dx + dy * dy;
+    double dist = sqrt(distSq);
+    assert(dist > precision);
+#endif
+    return attempt;
+}
+
 // this flavor approximates the cubics with quads to find the intersecting ts
 // OPTIMIZE: if this strategy proves successful, the quad approximations, or the ts used
 // to create the approximations, could be stored in the cubic segment
-// fixme: this strategy needs to add short line segments on either end, or similarly extend the
-// initial and final quadratics
-bool intersect2(const Cubic& c1, const Cubic& c2, Intersections& i) {
+// FIXME: this strategy needs to intersect the convex hull on either end with the opposite to
+// account for inset quadratics that cause the endpoint intersection to avoid detection
+// the segments can be very short -- the length of the maximum quadratic error (precision)
+// FIXME: this needs to recurse on itself, taking a range of T values and computing the new
+// t range ala is linear inner. The range can be figured by taking the dx/dy and determining
+// the fraction that matches the precision. That fraction is the change in t for the smaller cubic.
+static bool intersect2(const Cubic& cubic1, double t1s, double t1e, const Cubic& cubic2,
+        double t2s, double t2e, Intersections& i) {
+    Cubic c1, c2;
+    sub_divide(cubic1, t1s, t1e, c1);
+    sub_divide(cubic2, t2s, t2e, c2);
     SkTDArray<double> ts1;
-    double precision1 = calcPrecision(c1);
-    cubic_to_quadratics(c1, precision1, ts1);
+    cubic_to_quadratics(c1, calcPrecision(c1), ts1);
     SkTDArray<double> ts2;
-    double precision2 = calcPrecision(c2);
-    cubic_to_quadratics(c2, precision2, ts2);
-    double t1Start = 0;
+    cubic_to_quadratics(c2, calcPrecision(c2), ts2);
+    double t1Start = t1s;
     int ts1Count = ts1.count();
     for (int i1 = 0; i1 <= ts1Count; ++i1) {
-        const double t1 = i1 < ts1Count ? ts1[i1] : 1;
+        const double tEnd1 = i1 < ts1Count ? ts1[i1] : 1;
+        const double t1 = t1s + (t1e - t1s) * tEnd1;
         Cubic part1;
-        sub_divide(c1, t1Start, t1, part1);
+        sub_divide(cubic1, t1Start, t1, part1);
         Quadratic q1;
         demote_cubic_to_quad(part1, q1);
   //      start here;
         // should reduceOrder be looser in this use case if quartic is going to blow up on an
         // extremely shallow quadratic?
-        // maybe quadratics to lines need the same sort of recursive solution that I hope to find
-        // for cubics to quadratics ...
         Quadratic s1;
         int o1 = reduceOrder(q1, s1);
-        double t2Start = 0;
+        double t2Start = t2s;
         int ts2Count = ts2.count();
         for (int i2 = 0; i2 <= ts2Count; ++i2) {
-            const double t2 = i2 < ts2Count ? ts2[i2] : 1;
+            const double tEnd2 = i2 < ts2Count ? ts2[i2] : 1;
+            const double t2 = t2s + (t2e - t2s) * tEnd2;
             Cubic part2;
-            sub_divide(c2, t2Start, t2, part2);
+            sub_divide(cubic2, t2Start, t2, part2);
             Quadratic q2;
             demote_cubic_to_quad(part2, q2);
             Quadratic s2;
@@ -202,20 +227,34 @@ bool intersect2(const Cubic& c1, const Cubic& c2, Intersections& i) {
             if (o1 == 3 && o2 == 3) {
                 intersect2(q1, q2, locals);
             } else if (o1 <= 2 && o2 <= 2) {
-                i.fUsed = intersect((const _Line&) s1, (const _Line&) s2, i.fT[0], i.fT[1]);
+                locals.fUsed = intersect((const _Line&) s1, (const _Line&) s2, locals.fT[0],
+                        locals.fT[1]);
             } else if (o1 == 3 && o2 <= 2) {
-                intersect(q1, (const _Line&) s2, i);
+                intersect(q1, (const _Line&) s2, locals);
             } else {
                 SkASSERT(o1 <= 2 && o2 == 3);
-                intersect(q2, (const _Line&) s1, i);
-                for (int s = 0; s < i.fUsed; ++s) {
-                    SkTSwap(i.fT[0][s], i.fT[1][s]);
+                intersect(q2, (const _Line&) s1, locals);
+                for (int s = 0; s < locals.fUsed; ++s) {
+                    SkTSwap(locals.fT[0][s], locals.fT[1][s]);
                 }
             }
             for (int tIdx = 0; tIdx < locals.used(); ++tIdx) {
                 double to1 = t1Start + (t1 - t1Start) * locals.fT[0][tIdx];
                 double to2 = t2Start + (t2 - t2Start) * locals.fT[1][tIdx];
-                i.insert(to1, to2);
+    // if the computed t is not sufficiently precise, iterate
+                _Point p1, p2;
+                xy_at_t(cubic1, to1, p1.x, p1.y);
+                xy_at_t(cubic2, to2, p2.x, p2.y);
+                if (p1.approximatelyEqual(p2)) {
+                    i.insert(i.swapped() ? to2 : to1, i.swapped() ? to1 : to2);
+                } else {
+                    double dt1 = computeDelta(cubic1, to1, t1e - t1s);
+                    double dt2 = computeDelta(cubic2, to2, t2e - t2s);
+                    i.swap();
+                    intersect2(cubic2, std::max(to2 - dt2, 0.), std::min(to2 + dt2, 1.),
+                            cubic1, std::max(to1 - dt1, 0.), std::min(to1 + dt1, 1.), i);
+                    i.swap();
+                }
             }
             t2Start = t2;
         }
@@ -224,6 +263,69 @@ bool intersect2(const Cubic& c1, const Cubic& c2, Intersections& i) {
     return i.intersected();
 }
 
+static bool intersectEnd(const Cubic& cubic1, bool start, const Cubic& cubic2, const _Rect& bounds2,
+        Intersections& i) {
+    _Line line1;
+    line1[0] = line1[1] = cubic1[start ? 0 : 3];
+    _Point dxy1 = line1[0] - cubic1[start ? 1 : 2];
+    dxy1 /= precisionUnit;
+    line1[1] += dxy1;
+    _Rect line1Bounds;
+    line1Bounds.setBounds(line1);
+    if (!bounds2.intersects(line1Bounds)) {
+        return false;
+    }
+    _Line line2;
+    line2[0] = line2[1] = line1[0];
+    _Point dxy2 = line2[0] - cubic1[start ? 3 : 0];
+    dxy2 /= precisionUnit;
+    line2[1] += dxy2;
+#if 0 // this is so close to the first bounds test it isn't worth the short circuit test
+    _Rect line2Bounds;
+    line2Bounds.setBounds(line2);
+    if (!bounds2.intersects(line2Bounds)) {
+        return false;
+    }
+#endif
+    Intersections local1;
+    if (!intersect(cubic2, line1, local1)) {
+        return false;
+    }
+    Intersections local2;
+    if (!intersect(cubic2, line2, local2)) {
+        return false;
+    }
+    double tMin, tMax;
+    tMin = tMax = local1.fT[0][0];
+    for (int index = 1; index < local1.fUsed; ++index) {
+        tMin = std::min(tMin, local1.fT[0][index]);
+        tMax = std::max(tMax, local1.fT[0][index]);
+    }
+    for (int index = 1; index < local2.fUsed; ++index) {
+        tMin = std::min(tMin, local2.fT[0][index]);
+        tMax = std::max(tMax, local2.fT[0][index]);
+    }
+    return intersect2(cubic1, start ? 0 : 1, start ? 1.0 / precisionUnit : 1 - 1.0 / precisionUnit,
+            cubic2, tMin, tMax, i);
+}
+
+// FIXME: add intersection of convex null on cubics' ends with the opposite cubic. The hull line
+// segments can be constructed to be only as long as the calculated precision suggests. If the hull
+// line segments intersect the cubic, then use the intersections to construct a subdivision for
+// quadratic curve fitting.
+bool intersect2(const Cubic& c1, const Cubic& c2, Intersections& i) {
+    bool result = intersect2(c1, 0, 1, c2, 0, 1, i);
+    // FIXME: pass in cached bounds from caller
+    _Rect c1Bounds, c2Bounds;
+    c1Bounds.setBounds(c1); // OPTIMIZE use setRawBounds ?
+    c2Bounds.setBounds(c2);
+    result |= intersectEnd(c1, false, c2, c2Bounds, i);
+    result |= intersectEnd(c1, true, c2, c2Bounds, i);
+    result |= intersectEnd(c2, false, c1, c1Bounds, i);
+    result |= intersectEnd(c2, true, c1, c1Bounds, i);
+    return result;
+}
+
 int intersect(const Cubic& cubic, const Quadratic& quad, Intersections& i) {
     SkTDArray<double> ts;
     double precision = calcPrecision(cubic);
diff --git a/experimental/Intersection/CubicIntersection_Test.cpp b/experimental/Intersection/CubicIntersection_Test.cpp
index 61031087b3..c3ddded4be 100644
--- a/experimental/Intersection/CubicIntersection_Test.cpp
+++ b/experimental/Intersection/CubicIntersection_Test.cpp
@@ -78,14 +78,14 @@ static void oneOff(const Cubic& cubic1, const Cubic& cubic2) {
     intersect2(cubic1, cubic2, intersections2);
     for (int pt = 0; pt < intersections2.used(); ++pt) {
         double tt1 = intersections2.fT[0][pt];
-        double tx1, ty1;
-        xy_at_t(cubic1, tt1, tx1, ty1);
+        _Point xy1, xy2;
+        xy_at_t(cubic1, tt1, xy1.x, xy1.y);
         int pt2 = intersections2.fFlip ? intersections2.used() - pt - 1 : pt;
         double tt2 = intersections2.fT[1][pt2];
-        double tx2, ty2;
-        xy_at_t(cubic2, tt2, tx2, ty2);
+        xy_at_t(cubic2, tt2, xy2.x, xy2.y);
         SkDebugf("%s t1=%1.9g (%1.9g, %1.9g) (%1.9g, %1.9g) t2=%1.9g\n", __FUNCTION__,
-            tt1, tx1, ty1, tx2, ty2, tt2);
+            tt1, xy1.x, xy1.y, xy2.x, xy2.y, tt2);
+        assert(xy1.approximatelyEqual(xy2));
     }
 }
 
@@ -168,7 +168,7 @@ int depth;
 
 #define TRY_OLD 0 // old way fails on test == 1
 
-void CubicIntersection_RandTest() {
+void CubicIntersection_RandTestOld() {
     srand(0);
     const int tests = 1000000; // 10000000;
     double largestFactor = DBL_MAX;
@@ -265,3 +265,55 @@ void CubicIntersection_RandTest() {
         }
     }
 }
+
+void CubicIntersection_RandTest() {
+    srand(0);
+    const int tests = 1000000; // 10000000;
+ //   double largestFactor = DBL_MAX;
+    for (int test = 0; test < tests; ++test) {
+        Cubic cubic1, cubic2;
+        for (int i = 0; i < 4; ++i) {
+            cubic1[i].x = (double) rand() / RAND_MAX * 100;
+            cubic1[i].y = (double) rand() / RAND_MAX * 100;
+            cubic2[i].x = (double) rand() / RAND_MAX * 100;
+            cubic2[i].y = (double) rand() / RAND_MAX * 100;
+        }
+    #if DEBUG_CRASH
+        char str[1024];
+        sprintf(str, "{{%1.9g, %1.9g}, {%1.9g, %1.9g}, {%1.9g, %1.9g}, {%1.9g, %1.9g}},\n"
+            "{{%1.9g, %1.9g}, {%1.9g, %1.9g}, {%1.9g, %1.9g}, {%1.9g, %1.9g}},\n",
+                cubic1[0].x, cubic1[0].y,  cubic1[1].x, cubic1[1].y, cubic1[2].x, cubic1[2].y,
+                cubic1[3].x, cubic1[3].y,
+                cubic2[0].x, cubic2[0].y,  cubic2[1].x, cubic2[1].y, cubic2[2].x, cubic2[2].y,
+                cubic2[3].x, cubic2[3].y);
+    #endif
+        _Rect rect1, rect2;
+        rect1.setBounds(cubic1);
+        rect2.setBounds(cubic2);
+        bool boundsIntersect = rect1.left <= rect2.right && rect2.left <= rect2.right
+                && rect1.top <= rect2.bottom && rect2.top <= rect1.bottom;
+        Intersections i1, i2;
+        if (test == -1) {
+            SkDebugf("ready...\n");
+        }
+        Intersections intersections2;
+        bool newIntersects = intersect2(cubic1, cubic2, intersections2);
+        if (!boundsIntersect && newIntersects) {
+            SkDebugf("%s %d unexpected intersection boundsIntersect=%d "
+                    " newIntersects=%d\n%s %s\n", __FUNCTION__, test, boundsIntersect,
+                    newIntersects, __FUNCTION__, str);
+            assert(0);
+        }
+        for (int pt = 0; pt < intersections2.used(); ++pt) {
+            double tt1 = intersections2.fT[0][pt];
+            _Point xy1, xy2;
+            xy_at_t(cubic1, tt1, xy1.x, xy1.y);
+            int pt2 = intersections2.fFlip ? intersections2.used() - pt - 1 : pt;
+            double tt2 = intersections2.fT[1][pt2];
+            xy_at_t(cubic2, tt2, xy2.x, xy2.y);
+            SkDebugf("%s t1=%1.9g (%1.9g, %1.9g) (%1.9g, %1.9g) t2=%1.9g\n", __FUNCTION__,
+                tt1, xy1.x, xy1.y, xy2.x, xy2.y, tt2);
+            assert(xy1.approximatelyEqual(xy2));
+        }
+    }
+}
diff --git a/experimental/Intersection/CubicToQuadratics.cpp b/experimental/Intersection/CubicToQuadratics.cpp
index ff6d7a228b..66288b4316 100644
--- a/experimental/Intersection/CubicToQuadratics.cpp
+++ b/experimental/Intersection/CubicToQuadratics.cpp
@@ -137,6 +137,9 @@ static void addTs(const Cubic& cubic, double precision, double start, double end
 }
 
 // flavor that returns T values only, deferring computing the quads until they are needed
+// FIXME: when called from recursive intersect 2, this could take the original cubic
+// and do a more precise job when calling chop at and sub divide by computing the fractional ts.
+// it would still take the prechopped cubic for reduce order and find cubic inflections
 void cubic_to_quadratics(const Cubic& cubic, double precision, SkTDArray<double>& ts) {
     Cubic reduced;
     int order = reduceOrder(cubic, reduced, kReduceOrder_QuadraticsAllowed);
diff --git a/experimental/Intersection/CubicUtilities.cpp b/experimental/Intersection/CubicUtilities.cpp
index 36aa9e8fde..36fc17ee28 100644
--- a/experimental/Intersection/CubicUtilities.cpp
+++ b/experimental/Intersection/CubicUtilities.cpp
@@ -7,12 +7,15 @@
 #include "CubicUtilities.h"
 #include "QuadraticUtilities.h"
 
+const int precisionUnit = 256; // FIXME: arbitrary -- should try different values in test framework
+
+// FIXME: cache keep the bounds and/or precision with the caller?
 double calcPrecision(const Cubic& cubic) {
     _Rect dRect;
-    dRect.setBounds(cubic);
+    dRect.setBounds(cubic); // OPTIMIZATION: just use setRawBounds ?
     double width = dRect.right - dRect.left;
     double height = dRect.bottom - dRect.top;
-    return (width > height ? width : height) / 256;
+    return (width > height ? width : height) / precisionUnit;
 }
 
 void coefficients(const double* cubic, double& A, double& B, double& C, double& D) {
@@ -92,24 +95,30 @@ int cubicRoots(double A, double B, double C, double D, double t[3]) {
 // 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
-double derivativeAtT(const double* cubic, double t) {
+static double derivativeAtT(const double* cubic, double t) {
     double one_t = 1 - t;
     double a = cubic[0];
     double b = cubic[2];
     double c = cubic[4];
     double d = cubic[6];
-    return (b - a) * one_t * one_t + 2 * (c - b) * t * one_t + (d - c) * t * t;
+    return 3 * ((b - a) * one_t * one_t + 2 * (c - b) * t * one_t + (d - c) * t * t);
 }
 
-void dxdy_at_t(const Cubic& cubic, double t, double& dx, double& dy) {
-    if (&dx) {
-        dx = derivativeAtT(&cubic[0].x, t);
-    }
-    if (&dy) {
-        dy = derivativeAtT(&cubic[0].y, t);
-    }
+double dx_at_t(const Cubic& cubic, double t) {
+    return derivativeAtT(&cubic[0].x, t);
+}
+
+double dy_at_t(const Cubic& cubic, double t) {
+    return derivativeAtT(&cubic[0].y, t);
+}
+
+// OPTIMIZE? compute t^2, t(1-t), and (1-t)^2 and pass them to another version of derivative at t?
+void dxdy_at_t(const Cubic& cubic, double t, _Point& dxdy) {
+    dxdy.x = derivativeAtT(&cubic[0].x, t);
+    dxdy.y = derivativeAtT(&cubic[0].y, t);
 }
 
+
 int find_cubic_inflections(const Cubic& src, double tValues[])
 {
     double Ax = src[1].x - src[0].x;
@@ -138,6 +147,7 @@ bool rotate(const Cubic& cubic, int zero, int index, Cubic& rotPath) {
     return true;
 }
 
+#if 0 // unused for now
 double secondDerivativeAtT(const double* cubic, double t) {
     double a = cubic[0];
     double b = cubic[2];
@@ -145,6 +155,7 @@ double secondDerivativeAtT(const double* cubic, double t) {
     double d = cubic[6];
     return (c - 2 * b + a) * (1 - t) + (d - 2 * c + b) * t;
 }
+#endif
 
 void xy_at_t(const Cubic& cubic, double t, double& x, double& y) {
     double one_t = 1 - t;
diff --git a/experimental/Intersection/CubicUtilities.h b/experimental/Intersection/CubicUtilities.h
index 6159472f7e..6788b9c011 100644
--- a/experimental/Intersection/CubicUtilities.h
+++ b/experimental/Intersection/CubicUtilities.h
@@ -19,11 +19,12 @@ void coefficients(const double* cubic, double& A, double& B, double& C, double&
 int cubicRoots(double A, double B, double C, double D, double t[3]);
 int cubicRootsX(double A, double B, double C, double D, double s[3]);
 void demote_cubic_to_quad(const Cubic& cubic, Quadratic& quad);
-double derivativeAtT(const double* cubic, double t);
-void dxdy_at_t(const Cubic& , double t, double& x, double& y);
+double dx_at_t(const Cubic& , double t);
+double dy_at_t(const Cubic& , double t);
+void dxdy_at_t(const Cubic& , double t, _Point& y);
 int find_cubic_inflections(const Cubic& src, double tValues[]);
 bool rotate(const Cubic& cubic, int zero, int index, Cubic& rotPath);
-double secondDerivativeAtT(const double* cubic, double t);
 void xy_at_t(const Cubic& , double t, double& x, double& y);
 
+extern const int precisionUnit;
 #endif
diff --git a/experimental/Intersection/DataTypes.h b/experimental/Intersection/DataTypes.h
index 84c3b66201..afdb8e9182 100644
--- a/experimental/Intersection/DataTypes.h
+++ b/experimental/Intersection/DataTypes.h
@@ -135,11 +135,26 @@ struct _Point {
         return v;
     }
 
+    void operator+=(const _Point& v) {
+        x += v.x;
+        y += v.y;
+    }
+
     void operator-=(const _Point& v) {
         x -= v.x;
         y -= v.y;
     }
 
+    void operator/=(const double s) {
+        x /= s;
+        y /= s;
+    }
+
+    void operator*=(const double s) {
+        x *= s;
+        y *= s;
+    }
+
     friend bool operator==(const _Point& a, const _Point& b) {
         return a.x == b.x && a.y == b.y;
     }
@@ -187,11 +202,19 @@ struct _Rect {
     }
 
     // FIXME: used by debugging only ?
-    bool contains(const _Point& pt) {
+    bool contains(const _Point& pt) const {
         return approximately_between(left, pt.x, right)
                 && approximately_between(top, pt.y, bottom);
     }
 
+    bool intersects(_Rect& r) const {
+        assert(left < right);
+        assert(top < bottom);
+        assert(r.left < r.right);
+        assert(r.top < r.bottom);
+        return r.left < right && left < r.right && r.top < bottom && top < r.bottom;
+    }
+
     void set(const _Point& pt) {
         left = right = pt.x;
         top = bottom = pt.y;
diff --git a/experimental/Intersection/Intersection_Tests.h b/experimental/Intersection/Intersection_Tests.h
index 2ca1e056c8..6fc0c754ae 100644
--- a/experimental/Intersection/Intersection_Tests.h
+++ b/experimental/Intersection/Intersection_Tests.h
@@ -14,6 +14,7 @@ void CubicCoincidence_Test();
 void CubicIntersection_OneOffTest();
 void CubicIntersection_Test();
 void CubicIntersection_RandTest();
+void CubicIntersection_RandTestOld();
 void CubicParameterization_Test();
 void CubicReduceOrder_Test();
 void CubicsToQuadratics_RandTest();
diff --git a/experimental/Intersection/Intersections.cpp b/experimental/Intersection/Intersections.cpp
index d1e0f9bdb2..2ed6a1dc24 100644
--- a/experimental/Intersection/Intersections.cpp
+++ b/experimental/Intersection/Intersections.cpp
@@ -68,6 +68,7 @@ void Intersections::cleanUp() {
 
 }
 
+// FIXME: this doesn't respect swap, but add coincident does -- seems inconsistent
 void Intersections::insert(double one, double two) {
     assert(fUsed <= 1 || fT[0][0] < fT[0][1]);
     int index;
diff --git a/experimental/Intersection/Intersections.h b/experimental/Intersection/Intersections.h
index c8a142a19b..843957d141 100644
--- a/experimental/Intersection/Intersections.h
+++ b/experimental/Intersection/Intersections.h
@@ -12,17 +12,14 @@
 
 class Intersections {
 public:
-    Intersections()
-        : fUsed(0)
-        , fUsed2(0)
-        , fCoincidentUsed(0)
-        , fFlip(false)
-        , fUnsortable(false)
+    Intersections() 
+        : fFlip(0)
         , fSwap(0)
     {
         // OPTIMIZE: don't need to be initialized in release
         bzero(fT, sizeof(fT));
         bzero(fCoincidentT, sizeof(fCoincidentT));
+        reset();
     }
 
     void add(double one, double two) {
@@ -82,6 +79,12 @@ public:
         return fUsed == fUsed2;
     }
 
+    // leaves flip, swap alone
+    void reset() {
+        fUsed = fUsed2 = fCoincidentUsed = 0;
+        fUnsortable = false;
+    }
+
     void swap() {
         fSwap ^= true;
     }
diff --git a/experimental/Intersection/QuadraticImplicit.cpp b/experimental/Intersection/QuadraticImplicit.cpp
index 6e46557146..2d9d9b9a88 100644
--- a/experimental/Intersection/QuadraticImplicit.cpp
+++ b/experimental/Intersection/QuadraticImplicit.cpp
@@ -188,12 +188,12 @@ static bool addIntercept(const Quadratic& q1, const Quadratic& q2, double tMin,
     xy_at_t(q2, tMid, mid.x, mid.y);
     _Line line;
     line[0] = line[1] = mid;
-    double dx, dy;
-    dxdy_at_t(q2, tMid, dx, dy);
-    line[0].x -= dx;
-    line[0].y -= dy;
-    line[1].x += dx;
-    line[1].y += dy;
+    _Point dxdy;
+    dxdy_at_t(q2, tMid, dxdy);
+    line[0].x -= dxdy.x;
+    line[0].y -= dxdy.y;
+    line[1].x += dxdy.x;
+    line[1].y += dxdy.y;
     Intersections rootTs;
     int roots = intersect(q1, line, rootTs);
     assert(roots == 1);
@@ -251,10 +251,10 @@ static bool isLinearInner(const Quadratic& q1, double t1s, double t1e, const Qua
     }
     int split = 0;
     if (tMin != tMax || tCount > 2) {
-        dxdy_at_t(q2, tMin, dxy2.x, dxy2.y);
+        dxdy_at_t(q2, tMin, dxy2);
         for (int index = 1; index < tCount; ++index) {
             dxy1 = dxy2;
-            dxdy_at_t(q2, tsFound[index], dxy2.x, dxy2.y);
+            dxdy_at_t(q2, tsFound[index], dxy2);
             double dot = dxy1.dot(dxy2);
             if (dot < 0) {
                 split = index - 1;
@@ -309,6 +309,7 @@ static bool isLinear(const Quadratic& q1, const Quadratic& q2, Intersections& i)
 static bool relaxedIsLinear(const Quadratic& q1, const Quadratic& q2, Intersections& i) {
     double m1 = flatMeasure(q1);
     double m2 = flatMeasure(q2);
+    i.reset();
     if (fabs(m1) < fabs(m2)) {
         isLinearInner(q1, 0, 1, q2, 0, 1, i);
         return false;
@@ -329,7 +330,7 @@ static void unsortableExpanse(const Quadratic& q1, const Quadratic& q2, Intersec
     for (int qIdx = 0; qIdx < 2; qIdx++) {
         for (int t = 0; t < 2; t++) {
             _Point dxdy;
-            dxdy_at_t(*qs[qIdx], t, dxdy.x, dxdy.y);
+            dxdy_at_t(*qs[qIdx], t, dxdy);
             _Line perp;
             perp[0] = perp[1] = (*qs[qIdx])[t == 0 ? 0 : 2];
             perp[0].x += dxdy.y;
@@ -342,7 +343,7 @@ static void unsortableExpanse(const Quadratic& q1, const Quadratic& q2, Intersec
             if (hits) {
                 if (flip < 0) {
                     _Point dxdy2;
-                    dxdy_at_t(*qs[qIdx ^ 1], hitData.fT[0][0], dxdy2.x, dxdy2.y);
+                    dxdy_at_t(*qs[qIdx ^ 1], hitData.fT[0][0], dxdy2);
                     double dot = dxdy.dot(dxdy2);
                     flip = dot < 0;
                     i.fT[1][0] = flip;
diff --git a/experimental/Intersection/QuadraticIntersection_Test.cpp b/experimental/Intersection/QuadraticIntersection_Test.cpp
index 829dffebcc..432f614282 100644
--- a/experimental/Intersection/QuadraticIntersection_Test.cpp
+++ b/experimental/Intersection/QuadraticIntersection_Test.cpp
@@ -59,12 +59,15 @@ static void standardTestCases() {
 }
 
 static const Quadratic testSet[] = {
-{{53.774852327053594, 53.318060789841951}, {45.787877803416805, 51.393492026284981}, {46.703936967162392, 53.06860709822206}},
-{{46.703936967162392, 53.06860709822206}, {47.619996130907957, 54.74372217015916}, {53.020051653535361, 48.633140968832024}},
+{{67.426548091427676, 37.993772624988935}, {51.129513170665035, 57.542281234563646}, {44.594748190899189, 65.644267382683879}},
+{{61.336508189019057, 82.693132843213675}, {54.825078921449354, 71.663932799212432}, {47.727444217558926, 61.4049645128392}},
 
 {{67.4265481,37.9937726}, {51.1295132,57.5422812}, {44.5947482,65.6442674}},
 {{61.3365082,82.6931328}, {54.8250789,71.6639328}, {47.7274442,61.4049645}},
 
+{{53.774852327053594, 53.318060789841951}, {45.787877803416805, 51.393492026284981}, {46.703936967162392, 53.06860709822206}},
+{{46.703936967162392, 53.06860709822206}, {47.619996130907957, 54.74372217015916}, {53.020051653535361, 48.633140968832024}},
+
 {{50.934805397717923, 51.52391952648901}, {56.803308902971423, 44.246234610627596}, {69.776888596721406, 40.166645096692555}},
 {{50.230212796400401, 38.386469101526998}, {49.855620812184917, 38.818990392153609}, {56.356567496227363, 47.229909093319407}},
 
diff --git a/experimental/Intersection/QuadraticUtilities.cpp b/experimental/Intersection/QuadraticUtilities.cpp
index 95be90a460..8b02b4ce6e 100644
--- a/experimental/Intersection/QuadraticUtilities.cpp
+++ b/experimental/Intersection/QuadraticUtilities.cpp
@@ -64,16 +64,27 @@ int quadraticRoots(double A, double B, double C, double t[2]) {
     return foundRoots;
 }
 
-void dxdy_at_t(const Quadratic& quad, double t, double& x, double& y) {
+static double derivativeAtT(const double* quad, double t) {
     double a = t - 1;
     double b = 1 - 2 * t;
     double c = t;
-    if (&x) {
-        x = a * quad[0].x + b * quad[1].x + c * quad[2].x;
-    }
-    if (&y) {
-        y = a * quad[0].y + b * quad[1].y + c * quad[2].y;
-    }
+    return a * quad[0] + b * quad[2] + c * quad[4];
+}
+
+double dx_at_t(const Quadratic& quad, double t) {
+    return derivativeAtT(&quad[0].x, t);
+}
+
+double dy_at_t(const Quadratic& quad, double t) {
+    return derivativeAtT(&quad[0].y, t);
+}
+
+void dxdy_at_t(const Quadratic& quad, double t, _Point& dxy) {
+    double a = t - 1;
+    double b = 1 - 2 * t;
+    double c = t;
+    dxy.x = a * quad[0].x + b * quad[1].x + c * quad[2].x;
+    dxy.y = a * quad[0].y + b * quad[1].y + c * quad[2].y;
 }
 
 void xy_at_t(const Quadratic& quad, double t, double& x, double& y) {
diff --git a/experimental/Intersection/QuadraticUtilities.h b/experimental/Intersection/QuadraticUtilities.h
index e2f94e7f9f..13893a17e3 100644
--- a/experimental/Intersection/QuadraticUtilities.h
+++ b/experimental/Intersection/QuadraticUtilities.h
@@ -9,8 +9,9 @@
 #define QUADRATIC_UTILITIES_H
 
 #include "DataTypes.h"
-
-void dxdy_at_t(const Quadratic& , double t, double& x, double& y);
+double dx_at_t(const Quadratic& , double t);
+double dy_at_t(const Quadratic& , double t);
+void dxdy_at_t(const Quadratic& , double t, _Point& xy);
 
 /* Parameterization form, given A*t*t + 2*B*t*(1-t) + C*(1-t)*(1-t)
  *
diff --git a/experimental/Intersection/Simplify.cpp b/experimental/Intersection/Simplify.cpp
index 19390fb459..618dc30a26 100644
--- a/experimental/Intersection/Simplify.cpp
+++ b/experimental/Intersection/Simplify.cpp
@@ -305,15 +305,13 @@ static SkScalar LineDXAtT(const SkPoint a[2], double ) {
 
 static SkScalar QuadDXAtT(const SkPoint a[3], double t) {
     MAKE_CONST_QUAD(quad, a);
-    double x;
-    dxdy_at_t(quad, t, x, *(double*) 0);
+    double x = dx_at_t(quad, t);
     return SkDoubleToScalar(x);
 }
 
 static SkScalar CubicDXAtT(const SkPoint a[4], double t) {
     MAKE_CONST_CUBIC(cubic, a);
-    double x;
-    dxdy_at_t(cubic, t, x, *(double*) 0);
+    double x = dx_at_t(cubic, t);
     return SkDoubleToScalar(x);
 }
 
@@ -330,15 +328,13 @@ static SkScalar LineDYAtT(const SkPoint a[2], double ) {
 
 static SkScalar QuadDYAtT(const SkPoint a[3], double t) {
     MAKE_CONST_QUAD(quad, a);
-    double y;
-    dxdy_at_t(quad, t, *(double*) 0, y);
+    double y = dy_at_t(quad, t);
     return SkDoubleToScalar(y);
 }
 
 static SkScalar CubicDYAtT(const SkPoint a[4], double t) {
     MAKE_CONST_CUBIC(cubic, a);
-    double y;
-    dxdy_at_t(cubic, t, *(double*) 0, y);
+    double y = dy_at_t(cubic, t);
     return SkDoubleToScalar(y);
 }
 
diff --git a/experimental/Intersection/qc.htm b/experimental/Intersection/qc.htm
index 29fedab866..3fe3998b0c 100644
--- a/experimental/Intersection/qc.htm
+++ b/experimental/Intersection/qc.htm
@@ -1540,11 +1540,17 @@ $7 = {{x = 24.006224853920855, y = 72.621119847810419}, {x = 29.758671200376888,
 {{41.873704,30.8489321}, {53.3823801,24.1476487}, {77.3190123,0.63959742}},
 </div>
 
+<div id="quad14">
+{{67.4265481,37.9937726}, {51.1295132,57.5422812}, {44.5947482,65.6442674}},
+{{61.3365082,82.6931328}, {54.8250789,71.6639328}, {47.7274442,61.4049645}},
+</div>
+
 </div>
 
 <script type="text/javascript">
 
 var testDivs = [
+    quad14,
     cubic18,
     cubic17,
     cubic16,