shape ops work in progress

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

16 files changed, 541 insertions, 227 deletions

diff --git a/experimental/Intersection/CubicIntersection.cpp b/experimental/Intersection/CubicIntersection.cpp
index bedbfa100c..ae165f9b09 100644
--- a/experimental/Intersection/CubicIntersection.cpp
+++ b/experimental/Intersection/CubicIntersection.cpp
@@ -12,6 +12,9 @@
 #include "LineIntersection.h"
 #include "LineUtilities.h"
 
+#define DEBUG_COMPUTE_DELTA 1
+#define COMPUTE_DELTA 0
+
 static const double tClipLimit = 0.8; // http://cagd.cs.byu.edu/~tom/papers/bezclip.pdf see Multiple intersections
 
 class CubicIntersections : public Intersections {
@@ -162,27 +165,61 @@ bool intersect(const Cubic& c1, const Cubic& c2, Intersections& i) {
     return c.intersect();
 }
 
-#include "CubicUtilities.h"
-
+#if COMPUTE_DELTA
 static void cubicTangent(const Cubic& cubic, double t, _Line& tangent, _Point& pt, _Point& dxy) {
     xy_at_t(cubic, t, tangent[0].x, tangent[0].y);
     pt = tangent[1] = tangent[0];
     dxdy_at_t(cubic, t, dxy);
+    if (dxy.approximatelyZero()) {
+        if (approximately_zero(t)) {
+            SkASSERT(cubic[0].approximatelyEqual(cubic[1]));
+            dxy = cubic[2];
+            dxy -= cubic[0];
+        } else {
+            SkASSERT(approximately_equal(t, 1));
+            SkASSERT(cubic[3].approximatelyEqual(cubic[2]));
+            dxy = cubic[3];
+            dxy -= cubic[1];
+        }
+        SkASSERT(!dxy.approximatelyZero());
+    }
     tangent[0] -= dxy;
     tangent[1] += dxy;
+#if DEBUG_COMPUTE_DELTA
+    SkDebugf("%s t=%1.9g tangent=(%1.9g,%1.9g %1.9g,%1.9g)"
+            " pt=(%1.9g %1.9g) dxy=(%1.9g %1.9g)\n", __FUNCTION__, t,
+            tangent[0].x, tangent[0].y, tangent[1].x, tangent[1].y, pt.x, pt.y,
+            dxy.x, dxy.y);
+#endif
 }
+#endif
 
+#if COMPUTE_DELTA
 static double cubicDelta(const _Point& dxy, _Line& tangent, double scale)  {
     double tangentLen = dxy.length();
     tangent[0] -= tangent[1];
     double intersectLen = tangent[0].length();
     double result = intersectLen / tangentLen + scale;
+#if DEBUG_COMPUTE_DELTA
+    SkDebugf("%s tangent=(%1.9g,%1.9g %1.9g,%1.9g) intersectLen=%1.9g tangentLen=%1.9g scale=%1.9g"
+            " result=%1.9g\n", __FUNCTION__, tangent[0].x, tangent[0].y, tangent[1].x, tangent[1].y,
+            intersectLen, tangentLen, scale, result);
+#endif
     return result;
 }
+#endif
 
+#if COMPUTE_DELTA
 // FIXME: after testing, make this static
-void computeDelta(const Cubic& c1, double t1, double scale1, const Cubic& c2, double t2,
+static void computeDelta(const Cubic& c1, double t1, double scale1, const Cubic& c2, double t2,
         double scale2, double& delta1, double& delta2) {
+#if DEBUG_COMPUTE_DELTA
+    SkDebugf("%s c1=(%1.9g,%1.9g %1.9g,%1.9g %1.9g,%1.9g %1.9g,%1.9g) t1=%1.9g scale1=%1.9g"
+            " c2=(%1.9g,%1.9g %1.9g,%1.9g %1.9g,%1.9g %1.9g,%1.9g) t2=%1.9g scale2=%1.9g\n",
+            __FUNCTION__,
+            c1[0].x, c1[0].y, c1[1].x, c1[1].y, c1[2].x, c1[2].y, c1[3].x, c1[3].y, t1, scale1,
+            c2[0].x, c2[0].y, c2[1].x, c2[1].y, c2[2].x, c2[2].y, c2[3].x, c2[3].y, t2, scale2);
+#endif
     _Line tangent1, tangent2, line1, line2;
     _Point dxy1, dxy2;
     cubicTangent(c1, t1, line1, tangent1[0], dxy1);
@@ -209,6 +246,101 @@ void computeDelta(const Cubic& c1, double t1, double scale1, const Cubic& c2, do
 #if SK_DEBUG
 int debugDepth;
 #endif
+#endif
+
+static int quadPart(const Cubic& cubic, double tStart, double tEnd, Quadratic& simple) {
+    Cubic part;
+    sub_divide(cubic, tStart, tEnd, part);
+    Quadratic quad;
+    demote_cubic_to_quad(part, quad);
+    // FIXME: should reduceOrder be looser in this use case if quartic is going to blow up on an
+    // extremely shallow quadratic?
+    int order = reduceOrder(quad, simple);
+    return order;
+}
+
+static void intersectWithOrder(const Quadratic& simple1, int order1, const Quadratic& simple2,
+        int order2, Intersections& i) {
+    if (order1 == 3 && order2 == 3) {
+        intersect2(simple1, simple2, i);
+    } else if (order1 <= 2 && order2 <= 2) {
+        i.fUsed = intersect((const _Line&) simple1, (const _Line&) simple2, i.fT[0], i.fT[1]);
+    } else if (order1 == 3 && order2 <= 2) {
+        intersect(simple1, (const _Line&) simple2, i);
+    } else {
+        SkASSERT(order1 <= 2 && order2 == 3);
+        intersect(simple2, (const _Line&) simple1, i);
+        for (int s = 0; s < i.fUsed; ++s) {
+            SkTSwap(i.fT[0][s], i.fT[1][s]);
+        }
+    }
+}
+
+static void doIntersect(const Cubic& cubic1, double t1s, double t1m, double t1e,
+        const Cubic& cubic2, double t2s, double t2m, double t2e, Intersections& i) {
+    i.upDepth();
+    // divide the quadratics at the new t value and try again
+    double p1s = t1s;
+    double p1e = t1m;
+    for (int p1 = 0; p1 < 2; ++p1) {
+        Quadratic s1a;
+        int o1a = quadPart(cubic1, p1s, p1e, s1a);
+        double p2s = t2s;
+        double p2e = t2m;
+        for (int p2 = 0; p2 < 2; ++p2) {
+            Quadratic s2a;
+            int o2a = quadPart(cubic2, p2s, p2e, s2a);
+            Intersections locals;
+        #if 0 && SK_DEBUG
+            if (0.366025505 >= p1s && 0.366025887 <= p1e
+                    && 0.633974342 >= p2s && 0.633975487 <= p2e) {
+                SkDebugf("t1=(%1.9g,%1.9g) o1=%d t2=(%1.9g,%1.9g) o2=%d\n",
+                    p1s, p1e, o1a, p2s, p2e, o2a);
+                if (o1a == 2) {
+                    SkDebugf("{{%1.9g,%1.9g}, {%1.9g,%1.9g}},\n",
+                            s1a[0].x, s1a[0].y, s1a[1].x, s1a[1].y);
+                } else {
+                    SkDebugf("{{%1.9g,%1.9g}, {%1.9g,%1.9g}, {%1.9g,%1.9g}},\n",
+                            s1a[0].x, s1a[0].y, s1a[1].x, s1a[1].y, s1a[2].x, s1a[2].y);
+                }
+                if (o2a == 2) {
+                    SkDebugf("{{%1.9g,%1.9g}, {%1.9g,%1.9g}},\n",
+                            s2a[0].x, s2a[0].y, s2a[1].x, s2a[1].y);
+                } else {
+                    SkDebugf("{{%1.9g,%1.9g}, {%1.9g,%1.9g}, {%1.9g,%1.9g}},\n",
+                            s2a[0].x, s2a[0].y, s2a[1].x, s2a[1].y, s2a[2].x, s2a[2].y);
+                }
+                Intersections xlocals;
+                intersectWithOrder(s1a, o1a, s2a, o2a, xlocals);
+            } 
+        #endif
+            intersectWithOrder(s1a, o1a, s2a, o2a, locals);
+            for (int tIdx = 0; tIdx < locals.used(); ++tIdx) {
+                double to1 = p1s + (p1e - p1s) * locals.fT[0][tIdx];
+                double to2 = p2s + (p2e - p2s) * locals.fT[1][tIdx];
+    // 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 0 && SK_DEBUG
+                SkDebugf("to1=%1.9g p1=(%1.9g,%1.9g) to2=%1.9g p2=(%1.9g,%1.9g) d=%1.9g\n",
+                    to1, p1.x, p1.y, to2, p2.x, p2.y, p1.distance(p2));
+                    
+        #endif
+                if (p1.approximatelyEqual(p2)) {
+                    i.insert(i.swapped() ? to2 : to1, i.swapped() ? to1 : to2);
+                } else {
+                    doIntersect(cubic1, p1s, to1, p1e, cubic2, p2s, to2, p2e, i);
+                }
+            }
+            p2s = p2e;
+            p2e = t2e;
+        }
+        p1s = p1e;
+        p1e = t1e;
+    }
+    i.downDepth();
+}
 
 // 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
@@ -216,9 +348,6 @@ int debugDepth;
 // 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, double precisionScale, Intersections& i) {
     Cubic c1, c2;
@@ -233,41 +362,18 @@ static bool intersect2(const Cubic& cubic1, double t1s, double t1e, const Cubic&
     for (int i1 = 0; i1 <= ts1Count; ++i1) {
         const double tEnd1 = i1 < ts1Count ? ts1[i1] : 1;
         const double t1 = t1s + (t1e - t1s) * tEnd1;
-        Cubic 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?
         Quadratic s1;
-        int o1 = reduceOrder(q1, s1);
+        int o1 = quadPart(cubic1, t1Start, t1, s1);
         double t2Start = t2s;
         int ts2Count = ts2.count();
         for (int i2 = 0; i2 <= ts2Count; ++i2) {
             const double tEnd2 = i2 < ts2Count ? ts2[i2] : 1;
             const double t2 = t2s + (t2e - t2s) * tEnd2;
-            Cubic part2;
-            sub_divide(cubic2, t2Start, t2, part2);
-            Quadratic q2;
-            demote_cubic_to_quad(part2, q2);
             Quadratic s2;
-            double o2 = reduceOrder(q2, s2);
+            int o2 = quadPart(cubic2, t2Start, t2, s2);
             Intersections locals;
-            if (o1 == 3 && o2 == 3) {
-                intersect2(q1, q2, locals);
-            } else if (o1 <= 2 && o2 <= 2) {
-                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, locals);
-            } else {
-                SkASSERT(o1 <= 2 && o2 == 3);
-                intersect(q2, (const _Line&) s1, locals);
-                for (int s = 0; s < locals.fUsed; ++s) {
-                    SkTSwap(locals.fT[0][s], locals.fT[1][s]);
-                }
-            }
+            intersectWithOrder(s1, o1, s2, o2, locals);
+            
             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];
@@ -278,6 +384,7 @@ static bool intersect2(const Cubic& cubic1, double t1s, double t1e, const Cubic&
                 if (p1.approximatelyEqual(p2)) {
                     i.insert(i.swapped() ? to2 : to1, i.swapped() ? to1 : to2);
                 } else {
+            #if COMPUTE_DELTA
                     double dt1, dt2;
                     computeDelta(cubic1, to1, (t1e - t1s), cubic2, to2, (t2e - t2s), dt1, dt2);
                     double scale = precisionScale;
@@ -296,6 +403,9 @@ static bool intersect2(const Cubic& cubic1, double t1s, double t1e, const Cubic&
 #if SK_DEBUG
                     --debugDepth;
 #endif
+            #else
+                    doIntersect(cubic1, t1Start, to1, t1, cubic2, t2Start, to2, t2, i);
+            #endif
                 }
             }
             t2Start = t2;
@@ -308,8 +418,15 @@ static bool intersect2(const Cubic& cubic1, double t1s, double t1e, const Cubic&
 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];
+    line1[1] = cubic1[start ? 0 : 3];
+    if (line1[1].approximatelyEqual(cubic2[0]) || line1[1].approximatelyEqual(cubic2[3])) {
+        return false;
+    }
+    line1[0] = line1[1];
     _Point dxy1 = line1[0] - cubic1[start ? 1 : 2];
+    if (dxy1.approximatelyZero()) {
+        dxy1 = line1[0] - cubic1[start ? 2 : 1];
+    }
     dxy1 /= precisionUnit;
     line1[1] += dxy1;
     _Rect line1Bounds;
@@ -320,6 +437,7 @@ static bool intersectEnd(const Cubic& cubic1, bool start, const Cubic& cubic2, c
     _Line line2;
     line2[0] = line2[1] = line1[0];
     _Point dxy2 = line2[0] - cubic1[start ? 3 : 0];
+    SkASSERT(!dxy2.approximatelyZero());
     dxy2 /= precisionUnit;
     line2[1] += dxy2;
 #if 0 // this is so close to the first bounds test it isn't worth the short circuit test
@@ -347,19 +465,19 @@ static bool intersectEnd(const Cubic& cubic1, bool start, const Cubic& cubic2, c
         tMin = SkTMin(tMin, local2.fT[0][index]);
         tMax = SkTMax(tMax, local2.fT[0][index]);
     }
-#if SK_DEBUG
+#if SK_DEBUG && COMPUTE_DELTA
     debugDepth = 0;
 #endif
     return intersect2(cubic1, start ? 0 : 1, start ? 1.0 / precisionUnit : 1 - 1.0 / precisionUnit,
             cubic2, tMin, tMax, 1, i);
 }
 
-// FIXME: add intersection of convex null on cubics' ends with the opposite cubic. The hull line
+// FIXME: add intersection of convex hull 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) {
-#if SK_DEBUG
+#if SK_DEBUG && COMPUTE_DELTA
     debugDepth = 0;
 #endif
     bool result = intersect2(c1, 0, 1, c2, 0, 1, 1, i);
diff --git a/experimental/Intersection/CubicIntersection_Test.cpp b/experimental/Intersection/CubicIntersection_Test.cpp
index d6816e3400..66c4d17dd7 100644
--- a/experimental/Intersection/CubicIntersection_Test.cpp
+++ b/experimental/Intersection/CubicIntersection_Test.cpp
@@ -98,11 +98,17 @@ static void oneOff(const Cubic& cubic1, const Cubic& cubic2) {
 }
 
 static const Cubic testSet[] = {
+{{0, 0}, {0, 1}, {1, 1}, {1, 0}},
+{{1, 0}, {0, 0}, {0, 1}, {1, 1}},
+
 {{0, 1}, {0, 2}, {1, 0}, {1, 0}},
 {{0, 1}, {0, 1}, {1, 0}, {2, 0}},
 
-{{0, 0}, {0, 1}, {1, 1}, {1, 0}},
-{{1, 0}, {0, 0}, {0, 1}, {1, 1}},
+{{0, 1}, {1, 6}, {1, 0}, {2, 0}},
+{{0, 1}, {0, 2}, {1, 0}, {6, 1}},
+
+{{0, 1}, {5, 6}, {1, 0}, {1, 0}},
+{{0, 1}, {0, 1}, {1, 0}, {6, 5}},
 
 {{95.837747722788592, 45.025976907939643}, {16.564570095652982, 0.72959763963222402}, {63.209855865319199, 68.047528419665767}, {57.640240647662544, 59.524565264361243}},
 {{51.593891741518817, 38.53849970667553}, {62.34752929878772, 74.924924725166022}, {74.810149322641152, 34.17966562983564}, {29.368398119401373, 94.66719277886078}},
@@ -153,7 +159,7 @@ void CubicIntersection_OneOffTest() {
     }
 }
 
-#define DEBUG_CRASH 1
+#define DEBUG_CRASH 0
 
 class CubicChopper {
 public:
@@ -250,14 +256,18 @@ void CubicIntersection_RandTestOld() {
         }
         bool newIntersects = intersect2(cubic1, cubic2, i2);
         if (!boundsIntersect && (oldIntersects || newIntersects)) {
+    #if DEBUG_CRASH
             SkDebugf("%s %d unexpected intersection boundsIntersect=%d oldIntersects=%d"
                     " newIntersects=%d\n%s %s\n", __FUNCTION__, test, boundsIntersect,
                     oldIntersects, newIntersects, __FUNCTION__, str);
+    #endif
             SkASSERT(0);
         }
         if (oldIntersects && !newIntersects) {
+    #if DEBUG_CRASH
             SkDebugf("%s %d missing intersection oldIntersects=%d newIntersects=%d\n%s %s\n",
                     __FUNCTION__, test, oldIntersects, newIntersects, __FUNCTION__, str);
+    #endif
             SkASSERT(0);
         }
         if (!oldIntersects && !newIntersects) {
@@ -292,13 +302,17 @@ void CubicIntersection_RandTestOld() {
                 calc1, delta1, factor1, calc2, delta2, factor2);
         if (factor1 < largestFactor) {
             SkDebugf("WE HAVE A WINNER! %1.9g\n", factor1);
+    #if DEBUG_CRASH
             SkDebugf("%s\n", str);
+    #endif
             oneOff(cubic1, cubic2);
             largestFactor = factor1;
         }
         if (factor2 < largestFactor) {
             SkDebugf("WE HAVE A WINNER! %1.9g\n", factor2);
+    #if DEBUG_CRASH
             SkDebugf("%s\n", str);
+    #endif
             oneOff(cubic1, cubic2);
             largestFactor = factor2;
         }
@@ -336,9 +350,11 @@ void CubicIntersection_RandTest() {
         Intersections intersections2;
         bool newIntersects = intersect2(cubic1, cubic2, intersections2);
         if (!boundsIntersect && newIntersects) {
+    #if DEBUG_CRASH
             SkDebugf("%s %d unexpected intersection boundsIntersect=%d "
                     " newIntersects=%d\n%s %s\n", __FUNCTION__, test, boundsIntersect,
                     newIntersects, __FUNCTION__, str);
+    #endif
             SkASSERT(0);
         }
         for (int pt = 0; pt < intersections2.used(); ++pt) {
@@ -357,47 +373,91 @@ void CubicIntersection_RandTest() {
     }
 }
 
-static Cubic deltaTestSet[] = {
-  {{1, 4}, {1, 4.*2/3}, {1, 4.*1/3}, {1, 0}},
-  {{0, 3}, {1, 2}, {2, 1}, {3, 0}},
-  {{1, 4}, {1, 4.*2/3}, {1, 4.*1/3}, {1, 0}},
-  {{3.5, 1}, {2.5, 2}, {1.5, 3}, {0.5, 4}}
-};
-
-size_t deltaTestSetLen = sizeof(deltaTestSet) / sizeof(deltaTestSet[0]);
-
-static double deltaTestSetT[] = {
-    3./8,
-    5./12,
-    6./8,
-    9./12
-};
-
-size_t deltaTestSetTLen = sizeof(deltaTestSetT) / sizeof(deltaTestSetT[0]);
-
-static double expectedT[] = {
-    0.5,
-    1./3,
-    1./8,
-    5./6
-};
-
-size_t expectedTLen = sizeof(expectedT) / sizeof(expectedT[0]);
-
-// FIXME: this test no longer valid -- does not take minimum scale contribution into account
-void CubicIntersection_ComputeDeltaTest() {
-    SkASSERT(deltaTestSetLen == deltaTestSetTLen);
-    SkASSERT(expectedTLen == deltaTestSetTLen);
-    for (size_t index = 0; index < deltaTestSetLen; index += 2) {
-        const Cubic& c1 = deltaTestSet[index];
-        const Cubic& c2 = deltaTestSet[index + 1];
-        double t1 = deltaTestSetT[index];
-        double t2 = deltaTestSetT[index + 1];
-        double d1, d2;
-        computeDelta(c1, t1, 1, c2, t2, 1, d1, d2);
-        SkASSERT(approximately_equal(t1 + d1, expectedT[index])
-            || approximately_equal(t1 - d1, expectedT[index]));
-        SkASSERT(approximately_equal(t2 + d2, expectedT[index + 1])
-            || approximately_equal(t2 - d2, expectedT[index + 1]));
+void CubicIntersection_IntersectionFinder() {
+    Cubic cubic1 = {{0, 1}, {0, 2}, {1, 0}, {1, 0}};
+    Cubic cubic2 = {{0, 1}, {0, 2}, {1, 0}, {6, 1}};
+    double t1Seed = 0.19923954998177532;
+    double t2Seed = 0.17140596934291233;
+    double t1Step = 0.00035501449125494022;
+    double t2Step = 0.00020896171344569892;
+    _Point t1[3], t2[3];
+    bool toggle = true;
+    do {
+        xy_at_t(cubic1, t1Seed - t1Step, t1[0].x, t1[0].y);
+        xy_at_t(cubic1, t1Seed,          t1[1].x, t1[1].y);
+        xy_at_t(cubic1, t1Seed + t1Step, t1[2].x, t1[2].y);
+        xy_at_t(cubic2, t2Seed - t2Step, t2[0].x, t2[0].y);
+        xy_at_t(cubic2, t2Seed,          t2[1].x, t2[1].y);
+        xy_at_t(cubic2, t2Seed + t2Step, t2[2].x, t2[2].y);
+        double dist[3][3];
+        dist[1][1] = t1[1].distance(t2[1]);
+        int best_i = 1, best_j = 1;
+        for (int i = 0; i < 3; ++i) {
+            for (int j = 0; j < 3; ++j) {
+                if (i == 1 && j == 1) {
+                    continue;
+                }
+                dist[i][j] = t1[i].distance(t2[j]);
+                if (dist[best_i][best_j] > dist[i][j]) {
+                    best_i = i;
+                    best_j = j;
+                }
+            }
+        }
+        if (best_i == 0) {
+            t1Seed -= t1Step;
+        } else if (best_i == 2) {
+            t1Seed += t1Step;
+        }
+        if (best_j == 0) {
+            t2Seed -= t2Step;
+        } else if (best_j == 2) {
+            t2Seed += t2Step;
+        }
+        if (best_i == 1 && best_j == 1) {
+            if ((toggle ^= true)) {
+                t1Step /= 2;
+            } else {
+                t2Step /= 2;
+            }
+        }
+    } while (!t1[1].approximatelyEqual(t2[1]));
+    t1Step = t2Step = 0.1;
+    double t10 = t1Seed - t1Step * 2;
+    double t12 = t1Seed + t1Step * 2;
+    double t20 = t2Seed - t2Step * 2;
+    double t22 = t2Seed + t2Step * 2;
+    _Point test;
+    while (!approximately_zero(t1Step)) {
+        xy_at_t(cubic1, t10, test.x, test.y);
+        t10 += t1[1].approximatelyEqual(test) ? -t1Step : t1Step;
+        t1Step /= 2;
+    }
+    t1Step = 0.1;
+    while (!approximately_zero(t1Step)) {
+        xy_at_t(cubic1, t12, test.x, test.y);
+        t12 -= t1[1].approximatelyEqual(test) ? -t1Step : t1Step;
+        t1Step /= 2;
+    }
+    while (!approximately_zero(t2Step)) {
+        xy_at_t(cubic2, t20, test.x, test.y);
+        t20 += t2[1].approximatelyEqual(test) ? -t2Step : t2Step;
+        t2Step /= 2;
     }
+    t2Step = 0.1;
+    while (!approximately_zero(t2Step)) {
+        xy_at_t(cubic2, t22, test.x, test.y);
+        t22 -= t2[1].approximatelyEqual(test) ? -t2Step : t2Step;
+        t2Step /= 2;
+    }
+    SkDebugf("%s t1=(%1.9g<%1.9g<%1.9g) t2=(%1.9g<%1.9g<%1.9g)\n", __FUNCTION__,
+        t10, t1Seed, t12, t20, t2Seed, t22);
+}
+
+#if 0
+void CubicIntersection_CoincidentTest() {
+    Cubic cubic1 = {{0, 1}, {0, 2}, {1, 0}, {1, 0}};
+    Cubic cubic2 = {{0, 1}, {0, 2}, {1, 0}, {6, 1}};
+    
 }
+#endif
diff --git a/experimental/Intersection/CubicUtilities.cpp b/experimental/Intersection/CubicUtilities.cpp
index 70f10c17a4..a5174c95a9 100644
--- a/experimental/Intersection/CubicUtilities.cpp
+++ b/experimental/Intersection/CubicUtilities.cpp
@@ -115,7 +115,9 @@ int cubicRootsReal(double A, double B, double C, double D, double s[3]) {
     if (approximately_zero(A)) {  // we're just a quadratic
         return quadraticRootsReal(B, C, D, s);
     }
-    if (approximately_zero(D)) { // 0 is one root
+    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 = quadraticRootsReal(A, B, C, s);
         for (int i = 0; i < num; ++i) {
             if (approximately_zero(s[i])) {
@@ -231,7 +233,6 @@ void dxdy_at_t(const Cubic& cubic, double t, _Point& dxdy) {
     dxdy.y = derivativeAtT(&cubic[0].y, t);
 }
 
-
 int find_cubic_inflections(const Cubic& src, double tValues[])
 {
     double Ax = src[1].x - src[0].x;
diff --git a/experimental/Intersection/CubicUtilities.h b/experimental/Intersection/CubicUtilities.h
index 186ed43b6a..427ce34aa8 100644
--- a/experimental/Intersection/CubicUtilities.h
+++ b/experimental/Intersection/CubicUtilities.h
@@ -15,9 +15,6 @@ double calcPrecision(const Cubic& cubic);
 double calcPrecision(const Cubic& cubic, double t, double scale);
 #endif
 void chop_at(const Cubic& src, CubicPair& dst, double t);
-// FIXME: should be private static but public here for testing
-void computeDelta(const Cubic& c1, double t1, double scale1, const Cubic& c2, double t2,
-    double scale2, double& delta1, double& delta2);
 double cube_root(double x);
 int cubic_to_quadratics(const Cubic& cubic, double precision,
         SkTDArray<Quadratic>& quadratics);
diff --git a/experimental/Intersection/CurveIntersection.h b/experimental/Intersection/CurveIntersection.h
index 5f2745633f..34cd6fb6c2 100644
--- a/experimental/Intersection/CurveIntersection.h
+++ b/experimental/Intersection/CurveIntersection.h
@@ -50,6 +50,7 @@ bool intersect2(const Cubic& cubic1, const Cubic& cubic2, Intersections& );
 bool intersect(const Cubic& cubic, Intersections& i); // return true if cubic self-intersects
 int intersect(const Cubic& cubic, const Quadratic& quad, Intersections& );
 int intersect(const Cubic& cubic, const _Line& line, Intersections& );
+int intersectRay(const Cubic& quad, const _Line& line, Intersections& i);
 bool intersect(const Quadratic& q1, const Quadratic& q2, Intersections& );
 int intersect(const Quadratic& quad, const _Line& line, Intersections& );
 // the following flavor uses the implicit form instead of convex hulls
diff --git a/experimental/Intersection/DataTypes.cpp b/experimental/Intersection/DataTypes.cpp
index 2be8938e07..dc4fea1b82 100644
--- a/experimental/Intersection/DataTypes.cpp
+++ b/experimental/Intersection/DataTypes.cpp
@@ -9,17 +9,6 @@
 #include <sys/types.h>
 #include <stdlib.h>
 
-#ifdef __cplusplus
-extern "C" {
-#endif
-
-void    *memcpy(void *, const void *, size_t);
-
-#ifdef __cplusplus
-}
-#endif
-
-
 #if USE_EPSILON
 const double PointEpsilon = 0.000001;
 const double SquaredEpsilon = PointEpsilon * PointEpsilon;
@@ -32,13 +21,14 @@ union Float_t
 {
     Float_t(float num = 0.0f) : f(num) {}
     // Portable extraction of components.
-    bool Negative() const { return (i >> 31) != 0; }
+    bool negative() const { return (i >> 31) != 0; }
+#if 0 // unused
     int32_t RawMantissa() const { return i & ((1 << 23) - 1); }
     int32_t RawExponent() const { return (i >> 23) & 0xFF; }
-
+#endif
     int32_t i;
     float f;
-#ifdef _DEBUG
+#if SK_DEBUG
     struct
     {   // Bitfields for exploration. Do not use in production code.
         uint32_t mantissa : 23;
@@ -54,7 +44,7 @@ bool AlmostEqualUlps(float A, float B)
     Float_t uB(B);
 
     // Different signs means they do not match.
-    if (uA.Negative() != uB.Negative())
+    if (uA.negative() != uB.negative())
     {
         // Check for equality to make sure +0==-0
         return A == B;
diff --git a/experimental/Intersection/DataTypes.h b/experimental/Intersection/DataTypes.h
index e3d18afadd..ecf3d3676d 100644
--- a/experimental/Intersection/DataTypes.h
+++ b/experimental/Intersection/DataTypes.h
@@ -55,6 +55,10 @@ inline bool approximately_zero_inverse(double x) {
     return fabs(x) > FLT_EPSILON_INVERSE;
 }
 
+inline bool approximately_zero_when_compared_to(double x, double y) {
+    return fabs(x / y) < FLT_EPSILON;
+}
+
 // Use this for comparing Ts in the range of 0 to 1. For general numbers (larger and smaller) use
 // AlmostEqualUlps instead.
 inline bool approximately_equal(double x, double y) {
@@ -195,11 +199,20 @@ struct _Point {
         return AlmostEqualUlps((float) x, (float) a.x)
                 && AlmostEqualUlps((float) y, (float) a.y);
     }
+    
+    bool approximatelyZero() const {
+        return approximately_zero(x) && approximately_zero(y);
+    }
 
     double cross(const _Point& a) const {
         return x * a.y - y * a.x;
     }
 
+    double distance(const _Point& a) const {
+        _Point temp = *this - a;
+        return temp.length();
+    }
+
     double dot(const _Point& a) const {
         return x * a.x + y * a.y;
     }
diff --git a/experimental/Intersection/Intersection_Tests.cpp b/experimental/Intersection/Intersection_Tests.cpp
index fc68da46b4..cbbf4cb79a 100644
--- a/experimental/Intersection/Intersection_Tests.cpp
+++ b/experimental/Intersection/Intersection_Tests.cpp
@@ -13,9 +13,10 @@ void cubecode_test(int test);
 
 void Intersection_Tests() {
     int testsRun = 0;
-
-    SimplifyNew_Test();
+    
+    CubicIntersection_IntersectionFinder();
     CubicIntersection_OneOffTest();
+    SimplifyNew_Test();
     ShapeOps4x4CubicsThreaded_Test(testsRun);
     CubicToQuadratics_Test();
     QuadraticIntersection_Test();
diff --git a/experimental/Intersection/Intersection_Tests.h b/experimental/Intersection/Intersection_Tests.h
index c5a6d33f5a..5cd37b5574 100644
--- a/experimental/Intersection/Intersection_Tests.h
+++ b/experimental/Intersection/Intersection_Tests.h
@@ -11,7 +11,7 @@ void ConvexHull_Test();
 void ConvexHull_X_Test();
 void CubicBezierClip_Test();
 void CubicCoincidence_Test();
-void CubicIntersection_ComputeDeltaTest();
+void CubicIntersection_IntersectionFinder();
 void CubicIntersection_OneOffTest();
 void CubicIntersection_Test();
 void CubicIntersection_RandTest();
@@ -27,6 +27,13 @@ void LineIntersection_Test();
 void LineParameter_Test();
 void LineQuadraticIntersection_Test();
 void MiniSimplify_Test();
+void QuadLineIntersectThreaded_Test(int& );
+void QuadraticBezierClip_Test();
+void QuadraticCoincidence_Test();
+void QuadraticIntersection_Test();
+void QuadraticParameterization_Test();
+void QuadraticReduceOrder_Test();
+void QuarticRoot_Test();
 void SimplifyAddIntersectingTs_Test();
 void SimplifyAngle_Test();
 void SimplifyDegenerate4x4TrianglesThreaded_Test(int& );
@@ -43,10 +50,3 @@ void Simplify4x4RectsThreaded_Test(int& );
 void SimplifyRectangularPaths_Test();
 void ShapeOps4x4CubicsThreaded_Test(int& );
 void ShapeOps4x4RectsThreaded_Test(int& );
-void QuadLineIntersectThreaded_Test(int& );
-void QuadraticBezierClip_Test();
-void QuadraticCoincidence_Test();
-void QuadraticIntersection_Test();
-void QuadraticParameterization_Test();
-void QuadraticReduceOrder_Test();
-void QuarticRoot_Test();
diff --git a/experimental/Intersection/Intersections.h b/experimental/Intersection/Intersections.h
index 779ff33b3f..6dd2777136 100644
--- a/experimental/Intersection/Intersections.h
+++ b/experimental/Intersection/Intersections.h
@@ -12,6 +12,9 @@ public:
     Intersections()
         : fFlip(0)
         , fSwap(0)
+#if SK_DEBUG
+        , fDepth(0)
+#endif
     {
         // OPTIMIZE: don't need to be initialized in release
         bzero(fT, sizeof(fT));
@@ -52,9 +55,15 @@ public:
     // remove once curve/curve intersections are improved
     void cleanUp();
 
-    int coincidentUsed() const{
+    int coincidentUsed() const {
         return fCoincidentUsed;
     }
+    
+#if SK_DEBUG
+    int depth() const {
+        return fDepth;
+    }
+#endif
 
     void offset(int base, double start, double end) {
         for (int index = base; index < fUsed; ++index) {
@@ -105,6 +114,14 @@ public:
         return fUsed;
     }
 
+    void downDepth() {
+        SkASSERT(--fDepth >= 0);
+    }
+
+    void upDepth() {
+        SkASSERT(++fDepth < 16);
+    }
+
     double fT[2][9];
     double fCoincidentT[2][9];
     int fUsed;
@@ -112,6 +129,9 @@ public:
     int fCoincidentUsed;
     bool fFlip;
     bool fUnsortable;
+#if SK_DEBUG
+    int fDepth;
+#endif
 private:
     int fSwap;
 };
diff --git a/experimental/Intersection/LineCubicIntersection.cpp b/experimental/Intersection/LineCubicIntersection.cpp
index cc63099fd1..faed89fe91 100644
--- a/experimental/Intersection/LineCubicIntersection.cpp
+++ b/experimental/Intersection/LineCubicIntersection.cpp
@@ -287,3 +287,8 @@ int intersect(const Cubic& cubic, const _Line& line, Intersections& i) {
     LineCubicIntersections c(cubic, line, i);
     return c.intersect();
 }
+
+int intersectRay(const Cubic& cubic, const _Line& line, Intersections& i) {
+    LineCubicIntersections c(cubic, line, i);
+    return c.intersectRay(i.fT[0]);
+}
diff --git a/experimental/Intersection/LineIntersection.cpp b/experimental/Intersection/LineIntersection.cpp
index f5501db11a..58d00b830b 100644
--- a/experimental/Intersection/LineIntersection.cpp
+++ b/experimental/Intersection/LineIntersection.cpp
@@ -43,97 +43,81 @@ int intersect(const _Line& a, const _Line& b, double aRange[2], double bRange[2]
              byLen  * axLen         -   ayLen          * bxLen == 0 ( == denom )
      */
     double denom = byLen * axLen - ayLen * bxLen;
-    if (approximately_zero(denom)) {
-       /* See if the axis intercepts match:
-                  ay - ax * ayLen / axLen  ==          by - bx * ayLen / axLen
-         axLen * (ay - ax * ayLen / axLen) == axLen * (by - bx * ayLen / axLen)
-         axLen *  ay - ax * ayLen          == axLen *  by - bx * ayLen
-        */
-        // FIXME: need to use AlmostEqualUlps variant instead
-        if (approximately_equal_squared(axLen * a[0].y - ayLen * a[0].x,
-                axLen * b[0].y - ayLen * b[0].x)) {
-            const double* aPtr;
-            const double* bPtr;
-            if (fabs(axLen) > fabs(ayLen) || fabs(bxLen) > fabs(byLen)) {
-                aPtr = &a[0].x;
-                bPtr = &b[0].x;
-            } else {
-                aPtr = &a[0].y;
-                bPtr = &b[0].y;
-            }
-        #if 0 // sorting edges fails to preserve original direction
-            double aMin = aPtr[0];
-            double aMax = aPtr[2];
-            double bMin = bPtr[0];
-            double bMax = bPtr[2];
-            if (aMin > aMax) {
-                SkTSwap(aMin, aMax);
-            }
-            if (bMin > bMax) {
-                SkTSwap(bMin, bMax);
-            }
-            if (aMax < bMin || bMax < aMin) {
-                return 0;
-            }
-            if (aRange) {
-                aRange[0] = bMin <= aMin ? 0 : (bMin - aMin) / (aMax - aMin);
-                aRange[1] = bMax >= aMax ? 1 : (bMax - aMin) / (aMax - aMin);
-            }
-            int bIn = (aPtr[0] - aPtr[2]) * (bPtr[0] - bPtr[2]) < 0;
-            if (bRange) {
-                bRange[bIn] = aMin <= bMin ? 0 : (aMin - bMin) / (bMax - bMin);
-                bRange[!bIn] = aMax >= bMax ? 1 : (aMax - bMin) / (bMax - bMin);
-            }
-            return 1 + ((aRange[0] != aRange[1]) || (bRange[0] != bRange[1]));
-        #else
-            SkASSERT(aRange);
-            SkASSERT(bRange);
-            double a0 = aPtr[0];
-            double a1 = aPtr[2];
-            double b0 = bPtr[0];
-            double b1 = bPtr[2];
-            // OPTIMIZATION: restructure to reject before the divide
-            // e.g., if ((a0 - b0) * (a0 - a1) < 0 || abs(a0 - b0) > abs(a0 - a1))
-            // (except efficient)
-            double at0 = (a0 - b0) / (a0 - a1);
-            double at1 = (a0 - b1) / (a0 - a1);
-            if ((at0 < 0 && at1 < 0) || (at0 > 1 && at1 > 1)) {
-                return 0;
-            }
-            aRange[0] = SkTMax(SkTMin(at0, 1.0), 0.0);
-            aRange[1] = SkTMax(SkTMin(at1, 1.0), 0.0);
-            int bIn = (a0 - a1) * (b0 - b1) < 0;
-            double bDenom = b0 - b1;
-            if (approximately_zero(bDenom)) {
-                bRange[0] = bRange[1] = 0;
-            } else {
-                bRange[bIn] = SkTMax(SkTMin((b0 - a0) / bDenom, 1.0), 0.0);
-                bRange[!bIn] = SkTMax(SkTMin((b0 - a1) / bDenom, 1.0), 0.0);
-            }
-            bool second = fabs(aRange[0] - aRange[1]) > FLT_EPSILON;
-            SkASSERT((fabs(bRange[0] - bRange[1]) <= FLT_EPSILON) ^ second);
-            return 1 + second;
-        #endif
-        }
-    }
     double ab0y = a[0].y - b[0].y;
     double ab0x = a[0].x - b[0].x;
     double numerA = ab0y * bxLen - byLen * ab0x;
-    if ((numerA < 0 && denom > numerA) || (numerA > 0 && denom < numerA)) {
-        return 0;
-    }
     double numerB = ab0y * axLen - ayLen * ab0x;
-    if ((numerB < 0 && denom > numerB) || (numerB > 0 && denom < numerB)) {
+    bool mayNotOverlap = (numerA < 0 && denom > numerA) || (numerA > 0 && denom < numerA)
+            || (numerB < 0 && denom > numerB) || (numerB > 0 && denom < numerB);
+    numerA /= denom;
+    numerB /= denom;
+    if (!approximately_zero(denom) || (!approximately_zero_inverse(numerA) &&
+            !approximately_zero_inverse(numerB))) {
+        if (mayNotOverlap) {
+            return 0;
+        }
+        if (aRange) {
+            aRange[0] = numerA;
+        }
+        if (bRange) {
+            bRange[0] = numerB;
+        }
+        return 1;
+    }
+   /* See if the axis intercepts match:
+              ay - ax * ayLen / axLen  ==          by - bx * ayLen / axLen
+     axLen * (ay - ax * ayLen / axLen) == axLen * (by - bx * ayLen / axLen)
+     axLen *  ay - ax * ayLen          == axLen *  by - bx * ayLen
+    */
+    // FIXME: need to use AlmostEqualUlps variant instead
+    if (!approximately_equal_squared(axLen * a[0].y - ayLen * a[0].x,
+            axLen * b[0].y - ayLen * b[0].x)) {
         return 0;
     }
-    /* Is the intersection along the the segments */
-    if (aRange) {
-        aRange[0] = numerA / denom;
+    const double* aPtr;
+    const double* bPtr;
+    if (fabs(axLen) > fabs(ayLen) || fabs(bxLen) > fabs(byLen)) {
+        aPtr = &a[0].x;
+        bPtr = &b[0].x;
+    } else {
+        aPtr = &a[0].y;
+        bPtr = &b[0].y;
     }
-    if (bRange) {
-        bRange[0] = numerB / denom;
+    SkASSERT(aRange);
+    SkASSERT(bRange);
+    double a0 = aPtr[0];
+    double a1 = aPtr[2];
+    double b0 = bPtr[0];
+    double b1 = bPtr[2];
+    // OPTIMIZATION: restructure to reject before the divide
+    // e.g., if ((a0 - b0) * (a0 - a1) < 0 || abs(a0 - b0) > abs(a0 - a1))
+    // (except efficient)
+    double aDenom = a0 - a1;
+    if (approximately_zero(aDenom)) {
+        if (!between(b0, a0, b1)) {
+            return 0;
+        }
+        aRange[0] = aRange[1] = 0;
+    } else {
+        double at0 = (a0 - b0) / aDenom;
+        double at1 = (a0 - b1) / aDenom;
+        if ((at0 < 0 && at1 < 0) || (at0 > 1 && at1 > 1)) {
+            return 0;
+        }
+        aRange[0] = SkTMax(SkTMin(at0, 1.0), 0.0);
+        aRange[1] = SkTMax(SkTMin(at1, 1.0), 0.0);
     }
-    return 1;
+    double bDenom = b0 - b1;
+    if (approximately_zero(bDenom)) {
+        bRange[0] = bRange[1] = 0;
+    } else {
+        int bIn = aDenom * bDenom < 0;
+        bRange[bIn] = SkTMax(SkTMin((b0 - a0) / bDenom, 1.0), 0.0);
+        bRange[!bIn] = SkTMax(SkTMin((b0 - a1) / bDenom, 1.0), 0.0);
+    }
+    bool second = fabs(aRange[0] - aRange[1]) > FLT_EPSILON;
+    SkASSERT((fabs(bRange[0] - bRange[1]) <= FLT_EPSILON) ^ second);
+    return 1 + second;
 }
 
 int horizontalIntersect(const _Line& line, double y, double tRange[2]) {
diff --git a/experimental/Intersection/QuarticRoot.cpp b/experimental/Intersection/QuarticRoot.cpp
index 759e209569..389f68a0e9 100644
--- a/experimental/Intersection/QuarticRoot.cpp
+++ b/experimental/Intersection/QuarticRoot.cpp
@@ -47,7 +47,10 @@ int reducedQuarticRoots(const double t4, const double t3, const double t2, const
         }
         return cubicRootsReal(t3, t2, t1, t0, roots);
     }
-    if (approximately_zero(t0)) { // 0 is one root
+    if (approximately_zero_when_compared_to(t0, t1) // 0 is one root
+            && approximately_zero_when_compared_to(t0, t2)
+            && approximately_zero_when_compared_to(t0, t3)
+            && approximately_zero_when_compared_to(t0, t4)) {
         int num = cubicRootsReal(t4, t3, t2, t1, roots);
         for (int i = 0; i < num; ++i) {
             if (approximately_zero(roots[i])) {
diff --git a/experimental/Intersection/Simplify.cpp b/experimental/Intersection/Simplify.cpp
index 9b4ebc4481..f629a16057 100644
--- a/experimental/Intersection/Simplify.cpp
+++ b/experimental/Intersection/Simplify.cpp
@@ -31,7 +31,7 @@ int gDebugMaxWindValue = SK_MaxS32;
 #define APPROXIMATE_CUBICS 1
 
 #define DEBUG_UNUSED 0 // set to expose unused functions
-#define FORCE_RELEASE 1  // set force release to 1 for multiple thread -- no debugging
+#define FORCE_RELEASE 0  // set force release to 1 for multiple thread -- no debugging
 
 #if FORCE_RELEASE || defined SK_RELEASE
 
@@ -215,6 +215,11 @@ static void LineXYAtT(const SkPoint a[2], double t, SkPoint* out) {
     out->fY = SkDoubleToScalar(y);
 }
 
+static void LineXYAtT(const SkPoint a[2], double t, _Point* out) {
+    MAKE_CONST_LINE(line, a);
+    xy_at_t(line, t, out->x, out->y);
+}
+
 static void QuadXYAtT(const SkPoint a[3], double t, SkPoint* out) {
     MAKE_CONST_QUAD(quad, a);
     double x, y;
@@ -236,6 +241,11 @@ static void CubicXYAtT(const SkPoint a[4], double t, SkPoint* out) {
     out->fY = SkDoubleToScalar(y);
 }
 
+static void CubicXYAtT(const SkPoint a[4], double t, _Point* out) {
+    MAKE_CONST_CUBIC(cubic, a);
+    xy_at_t(cubic, t, out->x, out->y);
+}
+
 static void (* const SegmentXYAtT[])(const SkPoint [], double , SkPoint* ) = {
     NULL,
     LineXYAtT,
@@ -243,6 +253,13 @@ static void (* const SegmentXYAtT[])(const SkPoint [], double , SkPoint* ) = {
     CubicXYAtT
 };
 
+static void (* const SegmentXYAtT2[])(const SkPoint [], double , _Point* ) = {
+    NULL,
+    LineXYAtT,
+    QuadXYAtT,
+    CubicXYAtT
+};
+
 static SkScalar LineXAtT(const SkPoint a[2], double t) {
     MAKE_CONST_LINE(aLine, a);
     double x;
@@ -505,12 +522,25 @@ static int QuadRayIntersect(const SkPoint a[3], const SkPoint b[2],
 }
 #endif
 
-static int QuadRayIntersect(const SkPoint a[3], const _Line& bLine,
-        Intersections& intersections) {
+static int QuadRayIntersect(const SkPoint a[3], const _Line& bLine, Intersections& intersections) {
     MAKE_CONST_QUAD(aQuad, a);
     return intersectRay(aQuad, bLine, intersections);
 }
 
+static int CubicRayIntersect(const SkPoint a[3], const _Line& bLine, Intersections& intersections) {
+    MAKE_CONST_CUBIC(aCubic, a);
+    return intersectRay(aCubic, bLine, intersections);
+}
+
+static int (* const SegmentRayIntersect[])(const SkPoint [], const _Line& , Intersections&) = {
+    NULL,
+    NULL,
+    QuadRayIntersect,
+    CubicRayIntersect
+};
+
+
+
 static bool LineVertical(const SkPoint a[2], double startT, double endT) {
     MAKE_CONST_LINE(aLine, a);
     double x[2];
@@ -642,8 +672,8 @@ public:
             rh.fUnsortable = true;
             return this < &rh; // even with no solution, return a stable sort
         }
-        SkASSERT(fVerb == SkPath::kQuad_Verb); // worry about cubics later
-        SkASSERT(rh.fVerb == SkPath::kQuad_Verb);
+        SkASSERT(fVerb >= SkPath::kQuad_Verb);
+        SkASSERT(rh.fVerb >= SkPath::kQuad_Verb);
         // FIXME: until I can think of something better, project a ray from the
         // end of the shorter tangent to midway between the end points
         // through both curves and use the resulting angle to sort
@@ -655,15 +685,16 @@ public:
         int roots, rroots;
         bool flip = false;
         do {
-            const Quadratic& q = (len < rlen) ^ flip ? fQ : rh.fQ;
-            double midX = (q[0].x + q[2].x) / 2;
-            double midY = (q[0].y + q[2].y) / 2;
-            ray[0] = q[1];
-            ray[1].x = midX;
-            ray[1].y = midY;
+            bool useThis = (len < rlen) ^ flip;
+            const Cubic& part = useThis ? fCurvePart : rh.fCurvePart;
+            SkPath::Verb partVerb = useThis ? fVerb : rh.fVerb;
+            ray[0] = partVerb == SkPath::kCubic_Verb && part[0].approximatelyEqual(part[1]) ?
+                part[2] : part[1];
+            ray[1].x = (part[0].x + part[partVerb].x) / 2;
+            ray[1].y = (part[0].y + part[partVerb].y) / 2;
             SkASSERT(ray[0] != ray[1]);
-            roots = QuadRayIntersect(fPts, ray, i);
-            rroots = QuadRayIntersect(rh.fPts, ray, ri);
+            roots = (*SegmentRayIntersect[fVerb])(fPts, ray, i);
+            rroots = (*SegmentRayIntersect[rh.fVerb])(rh.fPts, ray, ri);
         } while ((roots == 0 || rroots == 0) && (flip ^= true));
         if (roots == 0 || rroots == 0) {
             // FIXME: we don't have a solution in this case. The interim solution
@@ -678,7 +709,7 @@ public:
         double dx, dy, dist;
         int index;
         for (index = 0; index < roots; ++index) {
-            QuadXYAtT(fPts, i.fT[0][index], &loc);
+            (*SegmentXYAtT2[fVerb])(fPts, i.fT[0][index], &loc);
             dx = loc.x - ray[0].x;
             dy = loc.y - ray[0].y;
             dist = dx * dx + dy * dy;
@@ -687,7 +718,7 @@ public:
             }
         }
         for (index = 0; index < rroots; ++index) {
-            QuadXYAtT(rh.fPts, ri.fT[0][index], &loc);
+            (*SegmentXYAtT2[rh.fVerb])(rh.fPts, ri.fT[0][index], &loc);
             dx = loc.x - ray[0].x;
             dy = loc.y - ray[0].y;
             dist = dx * dx + dy * dy;
@@ -763,35 +794,35 @@ public:
             fTangent1.lineEndPoints(l);
             fSide = 0;
             break;
-        case SkPath::kQuad_Verb:
-            QuadSubDivideHD(fPts, startT, endT, fQ);
-            fTangent1.quadEndPoints(fQ, 0, 1);
+        case SkPath::kQuad_Verb: {
+            Quadratic& quad = (Quadratic&)fCurvePart;
+            QuadSubDivideHD(fPts, startT, endT, quad);
+            fTangent1.quadEndPoints(quad, 0, 1);
         #if 1 // FIXME: try enabling this and see if a) it's called and b) does it break anything
             if (dx() == 0 && dy() == 0) {
                 SkDebugf("*** %s quad is line\n", __FUNCTION__);
-                fTangent1.quadEndPoints(fQ);
+                fTangent1.quadEndPoints(quad);
             }
         #endif
-            fSide = -fTangent1.pointDistance(fQ[2]); // not normalized -- compare sign only
-            break;
+            fSide = -fTangent1.pointDistance(fCurvePart[2]); // not normalized -- compare sign only
+            } break;
         case SkPath::kCubic_Verb: {
-            Cubic c;
             int nextC = 2;
-            CubicSubDivideHD(fPts, startT, endT, c);
-            fTangent1.cubicEndPoints(c, 0, 1);
+            CubicSubDivideHD(fPts, startT, endT, fCurvePart);
+            fTangent1.cubicEndPoints(fCurvePart, 0, 1);
             if (dx() == 0 && dy() == 0) {
-                fTangent1.cubicEndPoints(c, 0, 2);
+                fTangent1.cubicEndPoints(fCurvePart, 0, 2);
                 nextC = 3;
         #if 1 // FIXME: try enabling this and see if a) it's called and b) does it break anything
                 if (dx() == 0 && dy() == 0) {
                     SkDebugf("*** %s cubic is line\n");
-                    fTangent1.cubicEndPoints(c, 0, 3);
+                    fTangent1.cubicEndPoints(fCurvePart, 0, 3);
                 }
         #endif
             }
-            fSide = -fTangent1.pointDistance(c[nextC]); // not normalized -- compare sign only
+            fSide = -fTangent1.pointDistance(fCurvePart[nextC]); // compare sign only
             if (nextC == 2 && approximately_zero(fSide)) {
-                fSide = -fTangent1.pointDistance(c[3]);
+                fSide = -fTangent1.pointDistance(fCurvePart[3]);
             }
             } break;
         default:
@@ -876,7 +907,7 @@ public:
 
 private:
     const SkPoint* fPts;
-    Quadratic fQ;
+    Cubic fCurvePart;
     SkPath::Verb fVerb;
     double fSide;
     LineParameters fTangent1;
diff --git a/experimental/Intersection/SimplifyNew_Test.cpp b/experimental/Intersection/SimplifyNew_Test.cpp
index aa0a45bb62..97473ce984 100644
--- a/experimental/Intersection/SimplifyNew_Test.cpp
+++ b/experimental/Intersection/SimplifyNew_Test.cpp
@@ -3564,7 +3564,7 @@ static void cubicOp1d() {
     testShapeOp(path, pathB, kDifference_Op);
 }
 
-static void (*firstTest)() = 0;
+static void (*firstTest)() = testCubic1;
 
 static struct {
     void (*fun)();
diff --git a/experimental/Intersection/qc.htm b/experimental/Intersection/qc.htm
index 2053b9b9de..6bb96cc4ed 100644
--- a/experimental/Intersection/qc.htm
+++ b/experimental/Intersection/qc.htm
@@ -1717,11 +1717,101 @@ $7 = {{x = 24.006224853920855, y = 72.621119847810419}, {x = 29.758671200376888,
   {{1.03703704,0.259259259}, {1.49074074,0.0185185185}, {2,0}},
 </div>
 
+<div id="cubicTest1">
+{{0, 1}, {5, 6}, {1, 0}, {1, 0}},
+{{0, 1}, {0, 1}, {1, 0}, {6, 5}},
+  {{0,1}, {1.474,2.466}, {2.024,2.816}},
+  {{2.024,2.816}, {2.574,3.166}, {2.512,2.808}},
+  {{2.512,2.808}, {2.45,2.45}, {2.088,1.792}},
+  {{2.088,1.792}, {1.726,1.134}, {1.376,0.584}},
+  {{1.376,0.584}, {1.026,0.034}, {1,0}},
+
+  {{0,1}, {-0.0277777778,0.935185185}, {0.444444444,0.925925926}},
+  {{0.444444444,0.925925926}, {0.916666667,0.916666667}, {2.22222222,1.74074074}},
+  {{2.22222222,1.74074074}, {3.52777778,2.56481481}, {6,5}},
+</div>
+
+<div id="cubicTest2">
+{{fX = 0, fY = 1}, {fX = 0, fY = 2}, {fX = 1, fY = 0}, {fX = 5, fY = 0}}
+{{fX = 0, fY = 1}, {fX = 0, fY = 5}, {fX = 1, fY = 0}, {fX = 2, fY = 0}}
+</div>
+
+<div id="cubicTest3">
+{{x = 0, y = 1}, {x = 1, y = 6}, {x = 1, y = 0}, {x = 2, y = 0}}
+{{x = 0, y = 1}, {x = 0, y = 2}, {x = 1, y = 0}, {x = 6, y = 1}}
+  {{0,1}, {0.296,2.466}, {0.496,2.816}},
+  {{0.496,2.816}, {0.696,3.166}, {0.848,2.808}},
+  {{0.848,2.808}, {1,2.45}, {1.152,1.792}},
+  {{1.152,1.792}, {1.304,1.134}, {1.504,0.584}},
+  {{1.504,0.584}, {1.704,0.034}, {2,0}},
+
+  {{0,1}, {-0.040150997,1.4850292}, {0.586736465,1.17347293}},
+  {{0.586736465,1.17347293}, {1.11546634,0.936942311}, {2.40073788,0.7574854}},
+  {{2.40073788,0.7574854}, {3.68600943,0.578028489}, {6,1}},
+</div>
+
+<div id="lineTest1">
+{{x = 1.3834888994942065, y = 0.93137912059586503}, {x = 1.3835184247369658, y = 0.93128386633972826}}
+{{x = 1.3834889487833637, y = 0.93138119876148517}, {x = 1.3835243637966514, y = 0.93137251518247632}}
+</div>
+
+<div id="lineQuad1">
+{{x = 0.5, y = 0.875}, {x = 0.65728065326954233, y = 0.63483143515742313}, {x = 0.78907579690950191, y = 0.3959256577144401}}
+{{x = 0.78906720197447666, y = 0.39595053259806418}, {x = 0.78966023971986588, y = 0.39559384654294}}}
+</div>
+
+<div id="lineQuad2">
+{{x = 0.78907579690950191, y = 0.3959256577144401}, {x = 0.78953800250469452, y = 0.39509830801181633}}
+{{x = 0.29629629629629628, y = 0.74074074074074081}, {x = 0.49667506598603178, y = 0.57646999653254238}, {x = 0.78906720197447666, y = 0.39595053259806418}}
+</div>
+
+<div id="x1">
+{{0.5,0.875}, {0.657581172,0.634357518}, {0.789538003,0.395098308}},
+{{0.296296296,0.740740741}, {0.496887363,0.576287939}, {0.78966024,0.395593847}},
+</div>
+
+<div id="x2">
+{{0.5,0.875}, {0.657280653,0.634831435}, {0.789075797,0.395925658}},
+{{0.296296296,0.740740741}, {0.496675066,0.576469997}, {0.789067202,0.395950533}},
+</div>
+
+<div id="x3">
+{{0.789075797,0.395925658}, {0.789538003,0.395098308}},
+{{0.296296296,0.740740741}, {0.496675066,0.576469997}, {0.789067202,0.395950533}},
+</div>
+
+<div id="x4">
+{{0.789538003,0.395098308}, {0.818060471,0.344004048}, {0.84375,0.296875}},
+{{0.296296296,0.740740741}, {0.496887363,0.576287939}, {0.78966024,0.395593847}},
+</div>
+
+<div id="x5">
+{{0.789538003,0.395098308}, {0.818060471,0.344004048}, {0.84375,0.296875}},
+{{0.78966024,0.395593847}, {0.908006316,0.324112697}, {1.03703704,0.259259259}},
+</div>
+
+<div id="cubicX">
+{{x = 0, y = 1}, {x = 0, y = 2}, {x = 1, y = 0}, {x = 1, y = 0}}
+{{x = 0, y = 1}, {x = 0, y = 2}, {x = 1, y = 0}, {x = 6, y = 1}}
+</div>
+
 </div>
 
 <script type="text/javascript">
 
 var testDivs = [
+    cubicX,
+    x1,
+    x2,
+    x3,
+    x4,
+    x5,
+    lineQuad2,
+    lineQuad1,
+    lineTest1,
+    cubicTest3,
+    cubicTest2,
+    cubicTest1,
     cubicOp1d,
     testCubic1b,
     testCubic1a,