shape ops work in progress

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

25 files changed, 698 insertions, 117 deletions

diff --git a/experimental/Intersection/AddTestOutput/main.cpp b/experimental/Intersection/AddTestOutput/main.cpp
index 95955ceb78..83b315dc55 100644
--- a/experimental/Intersection/AddTestOutput/main.cpp
+++ b/experimental/Intersection/AddTestOutput/main.cpp
@@ -26,6 +26,8 @@ static bool replace(const char* fun, const char* dir, const char* filename, cons
     char* opInsert = strstr(opData.begin(), marker);
     if (!opInsert) {
         SkDebugf("%s missing marker in %s\n", fun, outFileStr.c_str());
+        opStreamOut.write(opData.begin(), opLen);
+        opStreamOut.flush();
         return false;
     }
     const char* opInsertEnd = opInsert + strlen(marker);
@@ -33,6 +35,8 @@ static bool replace(const char* fun, const char* dir, const char* filename, cons
         char* opInsert2 = strstr(opInsert, marker2);
         if (!opInsert2) {
             SkDebugf("%s missing marker second half in %s\n", fun, outFileStr.c_str());
+            opStreamOut.write(opData.begin(), opLen);
+            opStreamOut.flush();
             return false;
         }
         opInsertEnd = opInsert2 + strlen(marker2);
@@ -98,10 +102,10 @@ int main (int argc, char * const argv[]) {
         return 0;
     }
     const char simplifyMarker[] =
-            "#if 1 // set to 1 for multiple thread -- no debugging"
+            "#define FORCE_RELEASE 1  // set force release to 1 for multiple thread -- no debugging"
             ;
     const char simplifyReplace[] =
-            "#if 0 // set to 1 for multiple thread -- no debugging"
+            "#define FORCE_RELEASE 0  // set force release to 1 for multiple thread -- no debugging"
             ;
     if (!replace(argv[0], dir, "Simplify.cpp", simplifyMarker, NULL, simplifyReplace,
             sizeof(simplifyReplace) - 1)) {
diff --git a/experimental/Intersection/CubicIntersection.cpp b/experimental/Intersection/CubicIntersection.cpp
index b82b301a77..e5a2c18947 100644
--- a/experimental/Intersection/CubicIntersection.cpp
+++ b/experimental/Intersection/CubicIntersection.cpp
@@ -14,6 +14,7 @@
 
 #define DEBUG_COMPUTE_DELTA 1
 #define COMPUTE_DELTA 0
+#define DEBUG_QUAD_PART 0
 
 static const double tClipLimit = 0.8; // http://cagd.cs.byu.edu/~tom/papers/bezclip.pdf see Multiple intersections
 
@@ -256,6 +257,24 @@ static int quadPart(const Cubic& cubic, double tStart, double tEnd, Quadratic& s
     // 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);
+#if DEBUG_QUAD_PART
+    SkDebugf("%s cubic=(%1.17g,%1.17g %1.17g,%1.17g %1.17g,%1.17g %1.17g,%1.17g) t=(%1.17g,%1.17g)\n",
+            __FUNCTION__, cubic[0].x, cubic[0].y, cubic[1].x, cubic[1].y, cubic[2].x, cubic[2].y,
+            cubic[3].x, cubic[3].y, tStart, tEnd);
+    SkDebugf("%s part=(%1.17g,%1.17g %1.17g,%1.17g %1.17g,%1.17g %1.17g,%1.17g)"
+            " quad=(%1.17g,%1.17g %1.17g,%1.17g %1.17g,%1.17g)\n", __FUNCTION__, part[0].x, part[0].y,
+            part[1].x, part[1].y, part[2].x, part[2].y, part[3].x, part[3].y, quad[0].x, quad[0].y,
+            quad[1].x, quad[1].y, quad[2].x, quad[2].y);
+    SkDebugf("%s simple=(%1.17g,%1.17g", __FUNCTION__, simple[0].x, simple[0].y);
+    if (order > 1) {
+        SkDebugf(" %1.17g,%1.17g", simple[1].x, simple[1].y);
+    }
+    if (order > 2) {
+        SkDebugf(" %1.17g,%1.17g", simple[2].x, simple[2].y);
+    }
+    SkDebugf(")\n");
+    SkASSERT(order < 4 && order > 0);
+#endif
     return order;
 }
 
@@ -276,8 +295,33 @@ static void intersectWithOrder(const Quadratic& simple1, int order1, const Quadr
     }
 }
 
-static void doIntersect(const Cubic& cubic1, double t1s, double t1m, double t1e,
+static double distanceFromEnd(double t) {
+    return t > 0.5 ? 1 - t : t;
+}
+
+// OPTIMIZATION: this used to try to guess the value for delta, and that may still be worthwhile
+static void bumpForRetry(double t1, double t2, double& s1, double& e1, double& s2, double& e2) {
+    double dt1 = distanceFromEnd(t1);
+    double dt2 = distanceFromEnd(t2);
+    double delta = 1.0 / precisionUnit;
+    if (dt1 < dt2) {
+        if (t1 == dt1) {
+            s1 = SkTMax(s1 - delta, 0.);
+        } else {
+            e1 = SkTMin(e1 + delta, 1.);
+        }
+    } else {
+        if (t2 == dt2) {
+            s2 = SkTMax(s2 - delta, 0.);
+        } else {
+            e2 = SkTMin(e2 + delta, 1.);
+        }
+    }
+}
+
+static bool doIntersect(const Cubic& cubic1, double t1s, double t1m, double t1e,
         const Cubic& cubic2, double t2s, double t2m, double t2e, Intersections& i) {
+    bool result = false;
     i.upDepth();
     // divide the quadratics at the new t value and try again
     double p1s = t1s;
@@ -292,8 +336,8 @@ static void doIntersect(const Cubic& cubic1, double t1s, double t1m, double t1e,
             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) {
+            if (0.497026154 >= p1s && 0.497026535 <= p1e
+                    && 0.710440575 >= p2s && 0.710440956 <= 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) {
@@ -312,7 +356,8 @@ static void doIntersect(const Cubic& cubic1, double t1s, double t1m, double t1e,
                 }
                 Intersections xlocals;
                 intersectWithOrder(s1a, o1a, s2a, o2a, xlocals);
-            }
+                SkDebugf("xlocals.fUsed=%d\n", xlocals.used());
+            } 
         #endif
             intersectWithOrder(s1a, o1a, s2a, o2a, locals);
             for (int tIdx = 0; tIdx < locals.used(); ++tIdx) {
@@ -325,12 +370,26 @@ static void doIntersect(const Cubic& cubic1, double t1s, double t1m, double t1e,
         #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);
+                    result = true;
                 } else {
-                    doIntersect(cubic1, p1s, to1, p1e, cubic2, p2s, to2, p2e, i);
+                    result = doIntersect(cubic1, p1s, to1, p1e, cubic2, p2s, to2, p2e, i);
+                    // if both cubics curve in the same direction, the quadratic intersection
+                    // may mark a range that does not contain the cubic intersection. If no 
+                    // intersection is found, look again including the t distance of the 
+                    // of the quadratic intersection nearest a quadratic end (which in turn is
+                    // nearest the actual cubic) 
+                    if (!result) {
+                        double b1s = p1s;
+                        double b1e = p1e;
+                        double b2s = p2s;
+                        double b2e = p2e;
+                        bumpForRetry(locals.fT[0][tIdx], locals.fT[1][tIdx], b1s, b1e, b2s, b2e);
+                        result = doIntersect(cubic1, b1s, to1, b1e, cubic2, b2s, to2, b2e, i);
+                    }
                 }
             }
             p2s = p2e;
@@ -340,6 +399,7 @@ static void doIntersect(const Cubic& cubic1, double t1s, double t1m, double t1e,
         p1e = t1e;
     }
     i.downDepth();
+    return result;
 }
 
 // this flavor approximates the cubics with quads to find the intersecting ts
@@ -371,9 +431,22 @@ static bool intersect2(const Cubic& cubic1, double t1s, double t1e, const Cubic&
             const double t2 = t2s + (t2e - t2s) * tEnd2;
             Quadratic s2;
             int o2 = quadPart(cubic2, t2Start, t2, s2);
+        #if 0 && SK_DEBUG
+            if (0.497026154 >= t1Start && 0.497026535 <= t1
+                    && 0.710440575 + 0.0004 >= t2Start && 0.710440956 <= t2) {
+                Cubic cSub1, cSub2;
+                sub_divide(cubic1, t1Start, tEnd1, cSub1);
+                sub_divide(cubic2, t2Start, tEnd2, cSub2);
+                SkDebugf("t1=(%1.9g,%1.9g) t2=(%1.9g,%1.9g)\n",
+                        t1Start, t1, t2Start, t2);
+                Intersections xlocals;
+                intersectWithOrder(s1, o1, s2, o2, xlocals);
+                SkDebugf("xlocals.fUsed=%d\n", xlocals.used());
+            }
+        #endif
             Intersections locals;
             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];
@@ -404,10 +477,34 @@ static bool intersect2(const Cubic& cubic1, double t1s, double t1e, const Cubic&
                     --debugDepth;
 #endif
             #else
-                    doIntersect(cubic1, t1Start, to1, t1, cubic2, t2Start, to2, t2, i);
+                #if 0 && SK_DEBUG
+                    if (0.497026154 >= t1Start && 0.497026535 <= t1
+                            && 0.710440575 >= t2Start && 0.710440956 <= t2) {
+                        SkDebugf("t1=(%1.9g,%1.9g) t2=(%1.9g,%1.9g)\n",
+                                t1Start, t1, t2Start, t2);
+                    }
+                #endif
+                    bool found = doIntersect(cubic1, t1Start, to1, t1, cubic2, t2Start, to2, t2, i);
+                    if (!found) {
+                        double b1s = t1Start;
+                        double b1e = t1;
+                        double b2s = t2Start;
+                        double b2e = t2;
+                        bumpForRetry(locals.fT[0][tIdx], locals.fT[1][tIdx], b1s, b1e, b2s, b2e);
+                        doIntersect(cubic1, b1s, to1, b1e, cubic2, b2s, to2, b2e, i);
+                    }
             #endif
                 }
             }
+            if (locals.coincidentUsed()) {
+                SkASSERT(locals.coincidentUsed() == 2);
+                double coTs[2][2];
+                for (int tIdx = 0; tIdx < locals.coincidentUsed(); ++tIdx) {
+                    coTs[0][tIdx] = t1Start + (t1 - t1Start) * locals.fCoincidentT[0][tIdx];
+                    coTs[1][tIdx] = t2Start + (t2 - t2Start) * locals.fCoincidentT[1][tIdx];
+                }
+                i.addCoincident(coTs[0][0], coTs[0][1], coTs[1][0], coTs[1][1]);
+            }
             t2Start = t2;
         }
         t1Start = t1;
diff --git a/experimental/Intersection/CubicIntersection_Test.cpp b/experimental/Intersection/CubicIntersection_Test.cpp
index 7cfb118816..ae52a25d47 100644
--- a/experimental/Intersection/CubicIntersection_Test.cpp
+++ b/experimental/Intersection/CubicIntersection_Test.cpp
@@ -95,9 +95,35 @@ static void oneOff(const Cubic& cubic1, const Cubic& cubic2) {
 #endif
         SkASSERT(xy1.approximatelyEqual(xy2));
     }
+    for (int pt = 0; pt < intersections2.coincidentUsed(); ++pt) {
+        double tt1 = intersections2.fCoincidentT[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.fCoincidentT[1][pt2];
+        xy_at_t(cubic2, tt2, xy2.x, xy2.y);
+#if ONE_OFF_DEBUG
+        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);
+#endif
+        SkASSERT(xy1.approximatelyEqual(xy2));
+    }
 }
 
 static const Cubic testSet[] = {
+{{0,1}, {3,4}, {1,0}, {3,0}},
+{{0,1}, {0,3}, {1,0}, {4,3}},
+
+{{0, 0}, {1, 2}, {3, 4}, {4, 4}},
+{{0, 0}, {1, 2}, {3, 4}, {4, 4}},
+{{4, 4}, {3, 4}, {1, 2}, {0, 0}},
+
+{{0,1}, {2,3}, {1,0}, {1,0}},
+{{0,1}, {0,1}, {1,0}, {3,2}},
+
+{{0,2}, {0,1}, {1,0}, {1,0}},
+{{0,1}, {0,1}, {2,0}, {1,0}},
+
 {{0, 0}, {0, 1}, {1, 1}, {1, 0}},
 {{1, 0}, {0, 0}, {0, 1}, {1, 1}},
 
@@ -374,12 +400,12 @@ void CubicIntersection_RandTest() {
 }
 
 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;
+    Cubic cubic1 = {{0,1}, {3,4}, {1,0}, {3,0}};
+    Cubic cubic2 = {{0,1}, {0,3}, {1,0}, {4,3}};
+    double t1Seed = 0.496;
+    double t2Seed = 0.711;
+    double t1Step = 0.1;
+    double t2Step = 0.1;
     _Point t1[3], t2[3];
     bool toggle = true;
     do {
diff --git a/experimental/Intersection/CubicUtilities.cpp b/experimental/Intersection/CubicUtilities.cpp
index a5174c95a9..78be8df878 100644
--- a/experimental/Intersection/CubicUtilities.cpp
+++ b/experimental/Intersection/CubicUtilities.cpp
@@ -112,7 +112,9 @@ int cubicRootsReal(double A, double B, double C, double D, double s[3]) {
     bzero(str, sizeof(str));
     sprintf(str, "Solve[%1.19g x^3 + %1.19g x^2 + %1.19g x + %1.19g == 0, x]", A, B, C, D);
 #endif
-    if (approximately_zero(A)) {  // we're just a quadratic
+    if (approximately_zero_when_compared_to(A, B)
+            && approximately_zero_when_compared_to(A, C)
+            && approximately_zero_when_compared_to(A, D)) {  // we're just a quadratic
         return quadraticRootsReal(B, C, D, s);
     }
     if (approximately_zero_when_compared_to(D, A)
diff --git a/experimental/Intersection/CubicUtilities.h b/experimental/Intersection/CubicUtilities.h
index 427ce34aa8..f3c8e64bfd 100644
--- a/experimental/Intersection/CubicUtilities.h
+++ b/experimental/Intersection/CubicUtilities.h
@@ -32,4 +32,5 @@ void sub_divide(const Cubic& src, double t1, double t2, Cubic& dst);
 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 56348903d0..026d4976b7 100644
--- a/experimental/Intersection/DataTypes.h
+++ b/experimental/Intersection/DataTypes.h
@@ -12,6 +12,19 @@
 
 #include "SkTypes.h"
 
+// FIXME: move these into SkTypes.h
+template <typename T> inline T SkTMax(T a, T b) {
+    if (a < b)
+        a = b;
+    return a;
+}
+
+template <typename T> inline T SkTMin(T a, T b) {
+    if (a > b)
+        a = b;
+    return a;
+}
+
 extern bool AlmostEqualUlps(float A, float B);
 inline bool AlmostEqualUlps(double A, double B) { return AlmostEqualUlps((float) A, (float) B); }
 
@@ -55,8 +68,9 @@ inline bool approximately_zero_inverse(double x) {
     return fabs(x) > FLT_EPSILON_INVERSE;
 }
 
+// FIXME: if called multiple times with the same denom, we want to pass 1/y instead
 inline bool approximately_zero_when_compared_to(double x, double y) {
-    return fabs(x / y) < FLT_EPSILON;
+    return x == 0 || fabs(x / y) < FLT_EPSILON;
 }
 
 // Use this for comparing Ts in the range of 0 to 1. For general numbers (larger and smaller) use
@@ -144,7 +158,6 @@ inline bool approximately_zero_or_more(double x) {
 }
 
 inline bool approximately_between(double a, double b, double c) {
-    SkASSERT(a <= c);
     return a <= c ? approximately_negative(a - b) && approximately_negative(b - c)
             : approximately_negative(b - a) && approximately_negative(c - b);
 }
@@ -196,8 +209,17 @@ struct _Point {
     // return approximately_equal(a.y, y) && approximately_equal(a.x, x);
     // because that will not take the magnitude of the values
     bool approximatelyEqual(const _Point& a) const {
+#if 0
         return AlmostEqualUlps((float) x, (float) a.x)
                 && AlmostEqualUlps((float) y, (float) a.y);
+#else
+        double denom = SkTMax(fabs(x), SkTMax(fabs(y), SkTMax(fabs(a.x), fabs(a.y))));
+        if (denom == 0) {
+            return true;
+        }
+        double inv = 1 / denom;
+        return approximately_equal(x * inv, a.x * inv) && approximately_equal(y * inv, a.y * inv);
+#endif
     }
 
     bool approximatelyZero() const {
@@ -213,6 +235,11 @@ struct _Point {
         return temp.length();
     }
 
+    double distanceSquared(const _Point& a) const {
+        _Point temp = *this - a;
+        return temp.lengthSquared();
+    }
+
     double dot(const _Point& a) const {
         return x * a.x + y * a.y;
     }
@@ -229,6 +256,7 @@ struct _Point {
 
 typedef _Point _Line[2];
 typedef _Point Quadratic[3];
+typedef _Point Triangle[3];
 typedef _Point Cubic[4];
 
 struct _Rect {
@@ -294,17 +322,9 @@ struct QuadraticPair {
     _Point pts[5];
 };
 
-// FIXME: move these into SkTypes.h
-template <typename T> inline T SkTMax(T a, T b) {
-    if (a < b)
-        a = b;
-    return a;
-}
+// FIXME: move these into SkFloatingPoint.h
+#include "SkFloatingPoint.h"
 
-template <typename T> inline T SkTMin(T a, T b) {
-    if (a > b)
-        a = b;
-    return a;
-}
+#define sk_double_isnan(a) sk_float_isnan(a)
 
 #endif // __DataTypes_h__
diff --git a/experimental/Intersection/EdgeWalker_TestUtility.cpp b/experimental/Intersection/EdgeWalker_TestUtility.cpp
index ae5632e17d..cc95a0aa6f 100644
--- a/experimental/Intersection/EdgeWalker_TestUtility.cpp
+++ b/experimental/Intersection/EdgeWalker_TestUtility.cpp
@@ -58,19 +58,18 @@ static void showPathContour(SkPath::Iter& iter) {
     while ((verb = iter.next(pts)) != SkPath::kDone_Verb) {
         switch (verb) {
             case SkPath::kMove_Verb:
-                SkDebugf("path.moveTo(%1.9g, %1.9g);\n", pts[0].fX, pts[0].fY);
+                SkDebugf("path.moveTo(%1.9g,%1.9g);\n", pts[0].fX, pts[0].fY);
                 continue;
             case SkPath::kLine_Verb:
-                SkDebugf("path.lineTo(%1.9g, %1.9g);\n", pts[1].fX, pts[1].fY);
+                SkDebugf("path.lineTo(%1.9g,%1.9g);\n", pts[1].fX, pts[1].fY);
                 break;
             case SkPath::kQuad_Verb:
-                SkDebugf("path.quadTo(%1.9g, %1.9g, %1.9g, %1.9g);\n",
+                SkDebugf("path.quadTo(%1.9g,%1.9g, %1.9g,%1.9g);\n",
                     pts[1].fX, pts[1].fY, pts[2].fX, pts[2].fY);
                 break;
             case SkPath::kCubic_Verb:
-                SkDebugf("path.cubicTo(%1.9g, %1.9g, %1.9g, %1.9g);\n",
-                    pts[1].fX, pts[1].fY, pts[2].fX, pts[2].fY,
-                    pts[3].fX, pts[3].fY);
+                SkDebugf("path.cubicTo(%1.9g,%1.9g, %1.9g,%1.9g, %1.9g,%1.9g);\n",
+                    pts[1].fX, pts[1].fY, pts[2].fX, pts[2].fY, pts[3].fX, pts[3].fY);
                 break;
             case SkPath::kClose_Verb:
                 SkDebugf("path.close();\n");
diff --git a/experimental/Intersection/Intersection_Tests.cpp b/experimental/Intersection/Intersection_Tests.cpp
index a8b4416ef7..e828184e46 100644
--- a/experimental/Intersection/Intersection_Tests.cpp
+++ b/experimental/Intersection/Intersection_Tests.cpp
@@ -8,15 +8,20 @@
 #include "Intersection_Tests.h"
 
 void cubecode_test(int test);
+void parseSVG();
 
 #define TEST_QUADS_FIRST 0
 
 void Intersection_Tests() {
     int testsRun = 0;
-
+    QuadraticIntersection_OneOffTest();
     CubicIntersection_IntersectionFinder();
     CubicIntersection_OneOffTest();
+  #if 0
+    parseSVG();
+  #endif
     SimplifyNew_Test();
+    QuadraticIntersection_PointFinder();
     ShapeOps4x4CubicsThreaded_Test(testsRun);
     CubicToQuadratics_Test();
     QuadraticIntersection_Test();
diff --git a/experimental/Intersection/Intersection_Tests.h b/experimental/Intersection/Intersection_Tests.h
index 5cd37b5574..20aa80b362 100644
--- a/experimental/Intersection/Intersection_Tests.h
+++ b/experimental/Intersection/Intersection_Tests.h
@@ -27,6 +27,8 @@ void LineIntersection_Test();
 void LineParameter_Test();
 void LineQuadraticIntersection_Test();
 void MiniSimplify_Test();
+void QuadraticIntersection_OneOffTest();
+void QuadraticIntersection_PointFinder();
 void QuadLineIntersectThreaded_Test(int& );
 void QuadraticBezierClip_Test();
 void QuadraticCoincidence_Test();
diff --git a/experimental/Intersection/Intersections.cpp b/experimental/Intersection/Intersections.cpp
index cffeaec635..70e790134b 100644
--- a/experimental/Intersection/Intersections.cpp
+++ b/experimental/Intersection/Intersections.cpp
@@ -22,32 +22,32 @@ void Intersections::addCoincident(double s1, double e1, double s2, double e2) {
         if ((s1in | e1in) & (s2in | e2in)) {
             double lesser1 = SkTMin(cs1, ce1);
             index += cs1 > ce1;
-            if (s1in < lesser1) {
-                fCoincidentT[fSwap][index] = s1in;
-            } else if (e1in < lesser1) {
-                fCoincidentT[fSwap][index] = e1in;
+            if (s1 < lesser1) {
+                fCoincidentT[fSwap][index] = s1;
+            } else if (e1 < lesser1) {
+                fCoincidentT[fSwap][index] = e1;
             }
             index ^= 1;
             double greater1 = fCoincidentT[fSwap][index];
-            if (s1in > greater1) {
-                fCoincidentT[fSwap][index] = s1in;
-            } else if (e1in > greater1) {
-                fCoincidentT[fSwap][index] = e1in;
+            if (s1 > greater1) {
+                fCoincidentT[fSwap][index] = s1;
+            } else if (e1 > greater1) {
+                fCoincidentT[fSwap][index] = e1;
             }
             index &= ~1;
             double lesser2 = SkTMin(cs2, ce2);
             index += cs2 > ce2;
-            if (s2in < lesser2) {
-                fCoincidentT[fSwap ^ 1][index] = s2in;
-            } else if (e2in < lesser2) {
-                fCoincidentT[fSwap ^ 1][index] = e2in;
+            if (s2 < lesser2) {
+                fCoincidentT[fSwap ^ 1][index] = s2;
+            } else if (e2 < lesser2) {
+                fCoincidentT[fSwap ^ 1][index] = e2;
             }
             index ^= 1;
             double greater2 = fCoincidentT[fSwap ^ 1][index];
-            if (s2in > greater2) {
-                fCoincidentT[fSwap ^ 1][index] = s2in;
-            } else if (e2in > greater2) {
-                fCoincidentT[fSwap ^ 1][index] = e2in;
+            if (s2 > greater2) {
+                fCoincidentT[fSwap ^ 1][index] = s2;
+            } else if (e2 > greater2) {
+                fCoincidentT[fSwap ^ 1][index] = e2;
             }
             return;
         }
@@ -59,6 +59,9 @@ void Intersections::addCoincident(double s1, double e1, double s2, double e2) {
     fCoincidentT[fSwap][fCoincidentUsed] = e1;
     fCoincidentT[fSwap ^ 1][fCoincidentUsed] = e2;
     ++fCoincidentUsed;
+    // FIXME: assert that no entries in normal range lie between s and e
+    remove(s1, e1);
+    remove(s2, e2);
 }
 
 void Intersections::cleanUp() {
@@ -72,6 +75,12 @@ void Intersections::cleanUp() {
 void Intersections::insert(double one, double two) {
     SkASSERT(fUsed <= 1 || fT[0][0] < fT[0][1]);
     int index;
+    for (index = 0; index < fCoincidentUsed; ++index ) {
+        if (approximately_equal(fCoincidentT[0][index], one)
+                && approximately_equal(fCoincidentT[1][index], two)) {
+            return;
+        }
+    }
     for (index = 0; index < fUsed; ++index) {
         if (approximately_equal(fT[0][index], one)
                 && approximately_equal(fT[1][index], two)) {
@@ -114,3 +123,25 @@ void Intersections::insertOne(double t, int side) {
     fT[side][index] = t;
     side ? ++fUsed2 : ++fUsed;
 }
+
+void Intersections::remove(double one, double two) {
+    int index;
+    for (index = 0; index < fUsed; ++index) {
+        if (approximately_equal(fT[0][index], one)
+                && approximately_equal(fT[1][index], two)) {
+            goto foundIt;
+        }
+        if (fT[0][index] > one) {
+            return;
+        }
+    }
+    return;
+foundIt:
+    SkASSERT(fUsed > 0);
+    int remaining = fUsed - index;
+    if (remaining > 0) {
+        memmove(&fT[0][index], &fT[0][index - 1], sizeof(fT[0][0]) * remaining);
+        memmove(&fT[1][index], &fT[1][index - 1], sizeof(fT[1][0]) * remaining);
+    }
+    --fUsed;
+}
diff --git a/experimental/Intersection/Intersections.h b/experimental/Intersection/Intersections.h
index aed47bd61b..d873ba0e6e 100644
--- a/experimental/Intersection/Intersections.h
+++ b/experimental/Intersection/Intersections.h
@@ -132,6 +132,9 @@ public:
 #if SK_DEBUG
     int fDepth;
 #endif
+protected:
+    // used by addCoincident to remove ordinary intersections in range
+    void remove(double one, double two);
 private:
     int fSwap;
 };
diff --git a/experimental/Intersection/LineIntersection.cpp b/experimental/Intersection/LineIntersection.cpp
index 58d00b830b..6ab393684f 100644
--- a/experimental/Intersection/LineIntersection.cpp
+++ b/experimental/Intersection/LineIntersection.cpp
@@ -51,8 +51,9 @@ int intersect(const _Line& a, const _Line& b, double aRange[2], double bRange[2]
             || (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 ((!approximately_zero(denom) || (!approximately_zero_inverse(numerA)
+            && !approximately_zero_inverse(numerB))) && !sk_double_isnan(numerA)
+            && !sk_double_isnan(numerB)) {
         if (mayNotOverlap) {
             return 0;
         }
diff --git a/experimental/Intersection/LineUtilities.cpp b/experimental/Intersection/LineUtilities.cpp
index fa756b3d3b..6e3b966734 100644
--- a/experimental/Intersection/LineUtilities.cpp
+++ b/experimental/Intersection/LineUtilities.cpp
@@ -57,12 +57,12 @@ void sub_divide(const _Line& line, double t1, double t2, _Line& dst) {
 //            =0 for P2 on the line
 //            <0 for P2 right of the line
 //    See: the January 2001 Algorithm on Area of Triangles
-#if 0
-float isLeft( _Point P0, _Point P1, _Point P2 )
-{
-    return (float) ((P1.x - P0.x)*(P2.y - P0.y) - (P2.x - P0.x)*(P1.y - P0.y));
+//    return (float) ((P1.x - P0.x)*(P2.y - P0.y) - (P2.x - P0.x)*(P1.y - P0.y));
+double is_left(const _Line& line, const _Point& pt) {
+    _Point P0 = line[1] - line[0];
+    _Point P2 = pt - line[0];
+    return P0.cross(P2);
 }
-#endif
 
 double t_at(const _Line& line, const _Point& pt) {
     double dx = line[1].x - line[0].x;
diff --git a/experimental/Intersection/LineUtilities.h b/experimental/Intersection/LineUtilities.h
index 9d538122c0..6052875a9b 100644
--- a/experimental/Intersection/LineUtilities.h
+++ b/experimental/Intersection/LineUtilities.h
@@ -8,6 +8,7 @@
 
 bool implicitLine(const _Line& line, double& slope, double& axisIntercept);
 int reduceOrder(const _Line& line, _Line& reduced);
+double is_left(const _Line& line, const _Point& pt);
 void sub_divide(const _Line& src, double t1, double t2, _Line& dst);
 double t_at(const _Line&, const _Point& );
 void xy_at_t(const _Line& , double t, double& x, double& y);
diff --git a/experimental/Intersection/QuadraticImplicit.cpp b/experimental/Intersection/QuadraticImplicit.cpp
index 660ffe56f9..50f563356e 100644
--- a/experimental/Intersection/QuadraticImplicit.cpp
+++ b/experimental/Intersection/QuadraticImplicit.cpp
@@ -118,33 +118,6 @@ tryNextHalfPlane:
     return false;
 }
 
-// http://www.blackpawn.com/texts/pointinpoly/default.html
-static bool pointInTriangle(const _Point& pt, const _Line* testLines[]) {
-    const _Point& A = (*testLines[0])[0];
-    const _Point& B = (*testLines[1])[0];
-    const _Point& C = (*testLines[2])[0];
-
-// Compute vectors
-    _Point v0 = C - A;
-    _Point v1 = B - A;
-    _Point v2 = pt - A;
-
-// Compute dot products
-    double dot00 = v0.dot(v0);
-    double dot01 = v0.dot(v1);
-    double dot02 = v0.dot(v2);
-    double dot11 = v1.dot(v1);
-    double dot12 = v1.dot(v2);
-
-// Compute barycentric coordinates
-    double invDenom = 1 / (dot00 * dot11 - dot01 * dot01);
-    double u = (dot11 * dot02 - dot01 * dot12) * invDenom;
-    double v = (dot00 * dot12 - dot01 * dot02) * invDenom;
-
-// Check if point is in triangle
-    return (u >= 0) && (v >= 0) && (u + v < 1);
-}
-
 // returns false if there's more than one intercept or the intercept doesn't match the point
 // returns true if the intercept was successfully added or if the
 // original quads need to be subdivided
@@ -221,12 +194,12 @@ static bool isLinearInner(const Quadratic& q1, double t1s, double t1e, const Qua
     }
     _Point end;
     xy_at_t(q2, t2s, end.x, end.y);
-    bool startInTriangle = pointInTriangle(end, testLines);
+    bool startInTriangle = point_in_hull(hull, end);
     if (startInTriangle) {
         tMin = t2s;
     }
     xy_at_t(q2, t2e, end.x, end.y);
-    bool endInTriangle = pointInTriangle(end, testLines);
+    bool endInTriangle = point_in_hull(hull, end);
     if (endInTriangle) {
         tMax = t2e;
     }
diff --git a/experimental/Intersection/QuadraticIntersection_Test.cpp b/experimental/Intersection/QuadraticIntersection_Test.cpp
index 41ccb8b535..54536d87c8 100644
--- a/experimental/Intersection/QuadraticIntersection_Test.cpp
+++ b/experimental/Intersection/QuadraticIntersection_Test.cpp
@@ -8,7 +8,9 @@
 #include "CurveUtilities.h"
 #include "Intersection_Tests.h"
 #include "Intersections.h"
+#include "LineIntersection.h"
 #include "QuadraticIntersection_TestData.h"
+#include "QuadraticUtilities.h"
 #include "TestUtilities.h"
 
 const int firstQuadIntersectionTest = 9;
@@ -58,6 +60,20 @@ static void standardTestCases() {
 }
 
 static const Quadratic testSet[] = {
+  {{1.7465749139282332,1.9930452039527999}, {1.8320006564080331,1.859481345189089}, {1.8731035127758437,1.6344055934266613}},
+  {{1.8731035127758437,1.6344055934266613}, {1.89928170345231,1.5006405518943067}, {1.9223833226085514,1.3495796165215643}},
+  {{1.74657491,1.9930452}, {1.87407679,1.76762853}, {1.92238332,1.34957962}},
+  {{0.60797907,1.68776977}, {1.0447864,1.50810914}, {1.87464474,1.63655092}},
+  {{1.87464474,1.63655092}, {2.70450308,1.76499271}, {4,3}},
+
+{{1.2071879545809394,0.82163474041730045}, {1.1534203513372994,0.52790870069930229}, {1.0880000000000001,0.29599999999999982}}, //t=0.63155333662549329,0.80000000000000004
+{{0.33333333333333326,0.81481481481481488}, {0.63395173631977997,0.68744136726313931}, {1.205684411948591,0.81344322326274499}},
+{{0.33333333333333326,0.81481481481481488}, {0.63396444791444551,0.68743368362444768}, {1.205732763658403,0.81345617746834109}},//t=0.33333333333333331,0.63396444791444551
+{{1.205684411948591,0.81344322326274499}, {1.2057085875611198,0.81344969999329253}, {1.205732763658403,0.81345617746834109}},
+
+  {{1.20718795,0.82163474}, {1.15342035,0.527908701}, {1.088,0.296}},
+  {{1.20568441,0.813443223}, {1.20570859,0.8134497}, {1.20573276,0.813456177}},
+
   {{41.5072916,87.1234036}, {28.2747836,80.9545395}, {23.5780771,69.3344126}},
   {{72.9633878,95.6593007}, {42.7738746,88.4730382}, {31.1932785,80.2458029}},
 
@@ -133,7 +149,7 @@ static const Quadratic testSet[] = {
 
 const size_t testSetCount = sizeof(testSet) / sizeof(testSet[0]);
 
-#define ONE_OFF_DEBUG 0
+#define ONE_OFF_DEBUG 1
 
 static void oneOffTest1(size_t outer, size_t inner) {
     const Quadratic& quad1 = testSet[outer];
@@ -164,13 +180,19 @@ static void oneOffTest1(size_t outer, size_t inner) {
         }
 #if ONE_OFF_DEBUG
         SkDebugf("%s [%d][%d] t1=%1.9g (%1.9g, %1.9g) t2=%1.9g\n", __FUNCTION__,
-            outer, inner, tt1, tx1, tx2, tt2);
+            outer, inner, tt1, tx1, ty1, tt2);
 #endif
     }
 }
 
-static void oneOffTest() {
-//    oneOffTest1(0, 1);
+void QuadraticIntersection_OneOffTest() {
+    oneOffTest1(0, 3);
+    oneOffTest1(0, 4);
+    oneOffTest1(1, 3);
+    oneOffTest1(1, 4);
+}
+
+static void oneOffTests() {
     for (size_t outer = 0; outer < testSetCount - 1; ++outer) {
         for (size_t inner = outer + 1; inner < testSetCount; ++inner) {
             oneOffTest1(outer, inner);
@@ -204,7 +226,66 @@ static void coincidentTest() {
 }
 
 void QuadraticIntersection_Test() {
-    oneOffTest();
+    oneOffTests();
     coincidentTest();
     standardTestCases();
 }
+
+static int floatSign(double x) {
+    return x < 0 ? -1 : x > 0 ? 1 : 0;
+}
+
+static const Quadratic pointFinderTestSet[] = {
+                                                                                                                                //>=0.633974464         0.633974846 <=
+{{1.2071879545809394,0.82163474041730045}, {1.1534203513372994,0.52790870069930229}, {1.0880000000000001,0.29599999999999982}}, //t=0.63155333662549329,0.80000000000000004
+{{1.2071879545809394,0.82163474041730045}, {1.2065040319428038,0.81766753259119995}, {1.2058123269101506,0.81370135061854221}}, //t=0.63155333662549329,0.6339049773632347
+{{1.2058123269101506,0.81370135061854221}, {1.152376363978022,0.5244097415381026}, {1.0880000000000001,0.29599999999999982}},   //t=0.6339049773632347, 0.80000000000000004
+                                                                                                                                //>=0.633974083         0.633975227 <=
+{{0.33333333333333326,0.81481481481481488}, {0.63395173631977997,0.68744136726313931}, {1.205684411948591,0.81344322326274499}},//t=0.33333333333333331,0.63395173631977986
+{{0.33333333333333326,0.81481481481481488}, {0.63396444791444551,0.68743368362444768}, {1.205732763658403,0.81345617746834109}},//t=0.33333333333333331,0.63396444791444551
+{{1.205684411948591,0.81344322326274499}, {1.2057085875611198,0.81344969999329253}, {1.205732763658403,0.81345617746834109}},   //t=0.63395173631977986,0.63396444791444551
+{{1.205732763658403,0.81345617746834109}, {1.267928895828891,0.83008534558465619}, {1.3333333333333333,0.85185185185185175}},   //t=0.63396444791444551,0.66666666666666663
+};
+
+static void pointFinder(const Quadratic& q1, const Quadratic& q2) {
+    for (int index = 0; index < 3; ++index) {
+        double t = nearestT(q1, q2[index]);
+        _Point onQuad;
+        xy_at_t(q1, t, onQuad.x, onQuad.y);
+        SkDebugf("%s t=%1.9g (%1.9g,%1.9g) dist=%1.9g\n", __FUNCTION__, t, onQuad.x, onQuad.y,
+                onQuad.distance(q2[index]));
+        double left[3];
+        left[0] = is_left((const _Line&) q1[0], q2[index]);
+        left[1] = is_left((const _Line&) q1[1], q2[index]);
+        _Line diag = {q1[0], q1[2]};
+        left[2] = is_left(diag, q2[index]);
+        SkDebugf("%s left=(%d, %d, %d) inHull=%s\n", __FUNCTION__, floatSign(left[0]),
+                floatSign(left[1]), floatSign(left[2]),
+                point_in_hull(q1, q2[index]) ? "true" : "false");
+    }
+    SkDebugf("\n");
+}
+
+static void hullIntersect(const Quadratic& q1, const Quadratic& q2) {
+    SkDebugf("%s", __FUNCTION__);
+    double aRange[2], bRange[2];
+    for (int i1 = 0; i1 < 3; ++i1) {
+        _Line l1 = {q1[i1], q1[(i1 + 1) % 3]};
+        for (int i2 = 0; i2 < 3; ++i2) {
+            _Line l2 = {q2[i2], q2[(i2 + 1) % 3]};
+            if (intersect(l1, l2, aRange, bRange)) {
+                SkDebugf(" [%d,%d]", i1, i2);
+            }
+        }
+    }
+    SkDebugf("\n");
+}
+
+void QuadraticIntersection_PointFinder() {
+    pointFinder(pointFinderTestSet[0], pointFinderTestSet[4]);
+    pointFinder(pointFinderTestSet[4], pointFinderTestSet[0]);
+    pointFinder(pointFinderTestSet[0], pointFinderTestSet[6]);
+    pointFinder(pointFinderTestSet[6], pointFinderTestSet[0]);
+    hullIntersect(pointFinderTestSet[0], pointFinderTestSet[4]);
+    hullIntersect(pointFinderTestSet[0], pointFinderTestSet[6]);
+}
diff --git a/experimental/Intersection/QuadraticUtilities.cpp b/experimental/Intersection/QuadraticUtilities.cpp
index 32f95ed07f..a84009da87 100644
--- a/experimental/Intersection/QuadraticUtilities.cpp
+++ b/experimental/Intersection/QuadraticUtilities.cpp
@@ -4,23 +4,54 @@
  * Use of this source code is governed by a BSD-style license that can be
  * found in the LICENSE file.
  */
+#include "CubicUtilities.h"
 #include "QuadraticUtilities.h"
-#include <math.h>
+#include "TriangleUtilities.h"
 
-/*
+// from http://blog.gludion.com/2009/08/distance-to-quadratic-bezier-curve.html
+double nearestT(const Quadratic& quad, const _Point& pt) {
+    _Point pos = quad[0] - pt;
+    // search points P of bezier curve with PM.(dP / dt) = 0
+    // a calculus leads to a 3d degree equation :
+    _Point A = quad[1] - quad[0];
+    _Point B = quad[2] - quad[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 = cubicRootsValidT(a, b, c, d, ts);
+    double d0 = pt.distanceSquared(quad[0]);
+    double d2 = pt.distanceSquared(quad[2]);
+    double distMin = SkTMin(d0, d2);
+    int bestIndex = -1;
+    for (int index = 0; index < roots; ++index) {
+        _Point onQuad;
+        xy_at_t(quad, ts[index], onQuad.x, onQuad.y);
+        double dist = pt.distanceSquared(onQuad);
+        if (distMin > dist) {
+            distMin = dist;
+            bestIndex = index;
+        }
+    }
+    if (bestIndex >= 0) {
+        return ts[bestIndex];
+    }
+    return d0 < d2 ? 0 : 1;
+}
 
-Numeric Solutions (5.6) suggests to solve the quadratic by computing
+bool point_in_hull(const Quadratic& quad, const _Point& pt) {
+    return pointInTriangle((const Triangle&) quad, pt);
+}
 
+/*
+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
-
 */
-
-
 int add_valid_ts(double s[], int realRoots, double* t) {
     int foundRoots = 0;
     for (int index = 0; index < realRoots; ++index) {
diff --git a/experimental/Intersection/QuadraticUtilities.h b/experimental/Intersection/QuadraticUtilities.h
index f423dc6ecc..f6fb87856f 100644
--- a/experimental/Intersection/QuadraticUtilities.h
+++ b/experimental/Intersection/QuadraticUtilities.h
@@ -15,6 +15,8 @@ void chop_at(const Quadratic& src, QuadraticPair& dst, double t);
 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);
+double nearestT(const Quadratic& , const _Point& );
+bool point_in_hull(const Quadratic& , const _Point& );
 
 /* Parameterization form, given A*t*t + 2*B*t*(1-t) + C*(1-t)*(1-t)
  *
diff --git a/experimental/Intersection/QuarticRoot.cpp b/experimental/Intersection/QuarticRoot.cpp
index 389f68a0e9..50f85b59e9 100644
--- a/experimental/Intersection/QuarticRoot.cpp
+++ b/experimental/Intersection/QuarticRoot.cpp
@@ -41,8 +41,13 @@ int reducedQuarticRoots(const double t4, const double t3, const double t2, const
     sprintf(str, "Solve[%1.19g x^4 + %1.19g x^3 + %1.19g x^2 + %1.19g x + %1.19g == 0, x]",
         t4, t3, t2, t1, t0);
 #endif
-    if (approximately_zero(t4)) {
-        if (approximately_zero(t3)) {
+    if (approximately_zero_when_compared_to(t4, t0) // 0 is one root
+            && approximately_zero_when_compared_to(t4, t1)
+            && approximately_zero_when_compared_to(t4, t2)
+            && approximately_zero_when_compared_to(t4, t3)) {
+        if (approximately_zero_when_compared_to(t3, t0)
+            && approximately_zero_when_compared_to(t3, t1)
+            && approximately_zero_when_compared_to(t3, t2)) {
             return quadraticRootsReal(t2, t1, t0, roots);
         }
         return cubicRootsReal(t3, t2, t1, t0, roots);
diff --git a/experimental/Intersection/Simplify.cpp b/experimental/Intersection/Simplify.cpp
index f629a16057..f799699d10 100644
--- a/experimental/Intersection/Simplify.cpp
+++ b/experimental/Intersection/Simplify.cpp
@@ -31,12 +31,13 @@ int gDebugMaxWindValue = SK_MaxS32;
 #define APPROXIMATE_CUBICS 1
 
 #define DEBUG_UNUSED 0 // set to expose unused functions
-#define FORCE_RELEASE 0  // set force release to 1 for multiple thread -- no debugging
+#define FORCE_RELEASE 1  // set force release to 1 for multiple thread -- no debugging
 
 #if FORCE_RELEASE || defined SK_RELEASE
 
 const bool gRunTestsInOneThread = false;
 
+#define DEBUG_ACTIVE_OP 0
 #define DEBUG_ACTIVE_SPANS 0
 #define DEBUG_ACTIVE_SPANS_SHORT_FORM 0
 #define DEBUG_ADD_INTERSECTING_TS 0
@@ -59,6 +60,7 @@ const bool gRunTestsInOneThread = false;
 
 const bool gRunTestsInOneThread = true;
 
+#define DEBUG_ACTIVE_OP 1
 #define DEBUG_ACTIVE_SPANS 1
 #define DEBUG_ACTIVE_SPANS_SHORT_FORM 1
 #define DEBUG_ADD_INTERSECTING_TS 1
@@ -79,9 +81,11 @@ const bool gRunTestsInOneThread = true;
 
 #endif
 
-#define DEBUG_DUMP (DEBUG_ACTIVE_SPANS | DEBUG_CONCIDENT | DEBUG_SORT | DEBUG_PATH_CONSTRUCTION)
+#define DEBUG_DUMP (DEBUG_ACTIVE_OP | DEBUG_ACTIVE_SPANS | DEBUG_CONCIDENT | DEBUG_SORT | \
+        DEBUG_PATH_CONSTRUCTION)
 
 #if DEBUG_DUMP
+static const char* kShapeOpStr[] = {"diff", "sect", "union", "xor"};
 static const char* kLVerbStr[] = {"", "line", "quad", "cubic"};
 // static const char* kUVerbStr[] = {"", "Line", "Quad", "Cubic"};
 static int gContourID;
@@ -152,7 +156,7 @@ static int CubicIntersect(const SkPoint a[4], const SkPoint b[4], Intersections&
 #else
     intersect(aCubic, bCubic, intersections);
 #endif
-    return intersections.fUsed;
+    return intersections.fUsed ? intersections.fUsed : intersections.fCoincidentUsed;
 }
 
 static int HLineIntersect(const SkPoint a[2], SkScalar left, SkScalar right,
@@ -800,7 +804,7 @@ public:
             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__);
+ //               SkDebugf("*** %s quad is line\n", __FUNCTION__);
                 fTangent1.quadEndPoints(quad);
             }
         #endif
@@ -955,8 +959,8 @@ struct Bounds : public SkRect {
     bool isEmpty() {
         return fLeft > fRight || fTop > fBottom
                 || (fLeft == fRight && fTop == fBottom)
-                || isnan(fLeft) || isnan(fRight)
-                || isnan(fTop) || isnan(fBottom);
+                || sk_double_isnan(fLeft) || sk_double_isnan(fRight)
+                || sk_double_isnan(fTop) || sk_double_isnan(fBottom);
     }
 
     void setCubicBounds(const SkPoint a[4]) {
@@ -1316,6 +1320,10 @@ public:
             suTo = (oppSumWinding & xorSuMask) != 0;
         }
         bool result = gActiveEdge[op][miFrom][miTo][suFrom][suTo];
+#if DEBUG_ACTIVE_OP
+        SkDebugf("%s op=%s miFrom=%d miTo=%d suFrom=%d suTo=%d result=%d\n", __FUNCTION__,
+                kShapeOpStr[op], miFrom, miTo, suFrom, suTo, result);
+#endif
         SkASSERT(result != -1);
         return result;
     }
@@ -4175,8 +4183,8 @@ public:
             coincidence.fTs[swap][1] = ts.fT[0][1];
             coincidence.fTs[!swap][0] = ts.fT[1][0];
             coincidence.fTs[!swap][1] = ts.fT[1][1];
-        } else if (fSegments[index].verb() == SkPath::kQuad_Verb &&
-                other->fSegments[otherIndex].verb() == SkPath::kQuad_Verb) {
+        } else if (fSegments[index].verb() >= SkPath::kQuad_Verb &&
+                other->fSegments[otherIndex].verb() >= SkPath::kQuad_Verb) {
             coincidence.fTs[swap][0] = ts.fCoincidentT[0][0];
             coincidence.fTs[swap][1] = ts.fCoincidentT[0][1];
             coincidence.fTs[!swap][0] = ts.fCoincidentT[1][0];
@@ -5461,8 +5469,8 @@ static bool addIntersectTs(Contour* test, Contour* next) {
                     wt.addCoincident(wn, ts, swap);
                     continue;
                 }
-                if (wn.segmentType() == Work::kQuad_Segment
-                        && wt.segmentType() == Work::kQuad_Segment
+                if (wn.segmentType() >= Work::kQuad_Segment
+                        && wt.segmentType() >= Work::kQuad_Segment
                         && ts.coincidentUsed() == 2) {
                     wt.addCoincident(wn, ts, swap);
                     continue;
diff --git a/experimental/Intersection/SimplifyNew_Test.cpp b/experimental/Intersection/SimplifyNew_Test.cpp
index 97473ce984..567801bd98 100644
--- a/experimental/Intersection/SimplifyNew_Test.cpp
+++ b/experimental/Intersection/SimplifyNew_Test.cpp
@@ -3564,12 +3564,82 @@ static void cubicOp1d() {
     testShapeOp(path, pathB, kDifference_Op);
 }
 
-static void (*firstTest)() = testCubic1;
+static void cubicOp2d() {
+    SkPath path, pathB;
+    path.setFillType(SkPath::kWinding_FillType);
+    path.moveTo(0,2);
+    path.cubicTo(0,1, 1,0, 1,0);
+    path.close();
+    pathB.setFillType(SkPath::kWinding_FillType);
+    pathB.moveTo(0,1);
+    pathB.cubicTo(0,1, 2,0, 1,0);
+    pathB.close();
+    testShapeOp(path, pathB, kDifference_Op);
+}
+
+static void cubicOp3d() {
+    SkPath path, pathB;
+    path.setFillType(SkPath::kWinding_FillType);
+    path.moveTo(0,1);
+    path.cubicTo(2,3, 1,0, 1,0);
+    path.close();
+    pathB.setFillType(SkPath::kWinding_FillType);
+    pathB.moveTo(0,1);
+    pathB.cubicTo(0,1, 1,0, 3,2);
+    pathB.close();
+    testShapeOp(path, pathB, kDifference_Op);
+}
+
+static void cubicOp5d() {
+    SkPath path, pathB;
+    path.setFillType(SkPath::kWinding_FillType);
+    path.moveTo(0,1);
+    path.cubicTo(0,2, 1,0, 2,0);
+    path.close();
+    pathB.setFillType(SkPath::kWinding_FillType);
+    pathB.moveTo(0,1);
+    pathB.cubicTo(0,2, 1,0, 2,0);
+    pathB.close();
+    testShapeOp(path, pathB, kDifference_Op);
+}
+
+static void cubicOp6d() {
+    SkPath path, pathB;
+    path.setFillType(SkPath::kWinding_FillType);
+    path.moveTo(0,1);
+    path.cubicTo(0,6, 1,0, 3,0);
+    path.close();
+    pathB.setFillType(SkPath::kWinding_FillType);
+    pathB.moveTo(0,1);
+    pathB.cubicTo(0,3, 1,0, 6,0);
+    pathB.close();
+    testShapeOp(path, pathB, kDifference_Op);
+}
+
+static void cubicOp7d() {
+    SkPath path, pathB;
+    path.setFillType(SkPath::kWinding_FillType);
+    path.moveTo(0,1);
+    path.cubicTo(3,4, 1,0, 3,0);
+    path.close();
+    pathB.setFillType(SkPath::kWinding_FillType);
+    pathB.moveTo(0,1);
+    pathB.cubicTo(0,3, 1,0, 4,3);
+    pathB.close();
+    testShapeOp(path, pathB, kDifference_Op);
+}
+
+static void (*firstTest)() = 0;
 
 static struct {
     void (*fun)();
     const char* str;
 } tests[] = {
+    TEST(cubicOp7d),
+    TEST(cubicOp6d),
+    TEST(cubicOp5d),
+    TEST(cubicOp3d),
+    TEST(cubicOp2d),
     TEST(cubicOp1d),
  //   TEST(testQuadratic93),    // FIXME: gets stuck in a loop because top is unsortable
     TEST(testCubic1),
diff --git a/experimental/Intersection/TriangleUtilities.cpp b/experimental/Intersection/TriangleUtilities.cpp
new file mode 100644
index 0000000000..38a691646c
--- /dev/null
+++ b/experimental/Intersection/TriangleUtilities.cpp
@@ -0,0 +1,32 @@
+/*
+ * 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 "TriangleUtilities.h"
+
+// http://www.blackpawn.com/texts/pointinpoly/default.html
+bool pointInTriangle(const Triangle& triangle, const _Point& pt) {
+// Compute vectors
+    _Point v0 = triangle[2] - triangle[0];
+    _Point v1 = triangle[1] - triangle[0];
+    _Point v2 = pt - triangle[0];
+
+// Compute dot products
+    double dot00 = v0.dot(v0);
+    double dot01 = v0.dot(v1);
+    double dot02 = v0.dot(v2);
+    double dot11 = v1.dot(v1);
+    double dot12 = v1.dot(v2);
+
+// Compute barycentric coordinates
+    double invDenom = 1 / (dot00 * dot11 - dot01 * dot01);
+    double u = (dot11 * dot02 - dot01 * dot12) * invDenom;
+    double v = (dot00 * dot12 - dot01 * dot02) * invDenom;
+
+// Check if point is in triangle
+    return (u >= 0) && (v >= 0) && (u + v < 1);
+}
+
diff --git a/experimental/Intersection/TriangleUtilities.h b/experimental/Intersection/TriangleUtilities.h
new file mode 100644
index 0000000000..b09dbc6ef6
--- /dev/null
+++ b/experimental/Intersection/TriangleUtilities.h
@@ -0,0 +1,10 @@
+/*
+ * 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 "DataTypes.h"
+
+bool pointInTriangle(const Triangle& triangle, const _Point& pt);
diff --git a/experimental/Intersection/op.htm b/experimental/Intersection/op.htm
index 94d8adfd2b..5454fdc927 100644
--- a/experimental/Intersection/op.htm
+++ b/experimental/Intersection/op.htm
@@ -3339,11 +3339,83 @@ path.addRect(4, 13, 13, 16, SkPath::kCCW_Direction);
     pathB.close();
 </div>
 
+<div id="cubicOp2d">
+    path.setFillType(SkPath::kWinding_FillType);
+    path.moveTo(0,2);
+    path.cubicTo(0,1, 1,0, 1,0);
+    path.close();
+    pathB.setFillType(SkPath::kWinding_FillType);
+    pathB.moveTo(0,1);
+    pathB.cubicTo(0,1, 2,0, 1,0);
+    pathB.close();
+</div>
+
+<div id="cubicOp3d">
+    path.setFillType(SkPath::kWinding_FillType);
+    path.moveTo(0,1);
+    path.cubicTo(2,3, 1,0, 1,0);
+    path.close();
+    pathB.setFillType(SkPath::kWinding_FillType);
+    pathB.moveTo(0,1);
+    pathB.cubicTo(0,1, 1,0, 3,2);
+    pathB.close();
+</div>
+
+<div id="cubicOp4d">
+    path.setFillType(SkPath::kWinding_FillType);
+    path.moveTo(0,1);
+    path.cubicTo(0,2, 1,0, 2,0);
+    path.close();
+    pathB.setFillType(SkPath::kWinding_FillType);
+    pathB.moveTo(0,1);
+    pathB.cubicTo(0,2, 1,0, 2,0);
+    pathB.close();
+</div>
+
+<div id="cubicOp5d">
+    path.setFillType(SkPath::kWinding_FillType);
+    path.moveTo(0,1);
+    path.cubicTo(0,2, 1,0, 2,0);
+    path.close();
+    pathB.setFillType(SkPath::kWinding_FillType);
+    pathB.moveTo(0,1);
+    pathB.cubicTo(0,2, 1,0, 2,0);
+    pathB.close();
+</div>
+
+<div id="cubicOp6d">
+    path.setFillType(SkPath::kWinding_FillType);
+    path.moveTo(0,1);
+    path.cubicTo(0,6, 1,0, 3,0);
+    path.close();
+    pathB.setFillType(SkPath::kWinding_FillType);
+    pathB.moveTo(0,1);
+    pathB.cubicTo(0,3, 1,0, 6,0);
+    pathB.close();
+</div>
+
+<div id="cubicOp7d">
+    path.setFillType(SkPath::kWinding_FillType);
+    path.moveTo(0,1);
+    path.cubicTo(3,4, 1,0, 3,0);
+    path.close();
+    pathB.setFillType(SkPath::kWinding_FillType);
+    pathB.moveTo(0,1);
+    pathB.cubicTo(0,3, 1,0, 4,3);
+    pathB.close();
+</div>
+
 </div>
 
 <script type="text/javascript">
 
 var testDivs = [
+    cubicOp7d,
+    cubicOp6d,
+    cubicOp5d,
+    cubicOp4d,
+    cubicOp3d,
+    cubicOp2d,
     cubicOp1d,
     testQuadratic93,
     testCubic1,
diff --git a/experimental/Intersection/qc.htm b/experimental/Intersection/qc.htm
index 6bb96cc4ed..865fc13335 100644
--- a/experimental/Intersection/qc.htm
+++ b/experimental/Intersection/qc.htm
@@ -1795,11 +1795,115 @@ $7 = {{x = 24.006224853920855, y = 72.621119847810419}, {x = 29.758671200376888,
 {{x = 0, y = 1}, {x = 0, y = 2}, {x = 1, y = 0}, {x = 6, y = 1}}
 </div>
 
+<div id="cubicOp2d">
+{{0,2}, {0,1}, {1,0}, {1,0}},
+{{0,1}, {0,1}, {2,0}, {1,0}},
+  {{0,2}, {0.0185185185,1.49074074}, {0.259259259,1.03703704}},
+  {{0.259259259,1.03703704}, {0.5,0.583333333}, {0.740740741,0.296296296}},
+  {{0.740740741,0.296296296}, {0.981481481,0.00925925926}, {1,0}},
+
+  {{0,1}, {0.01953125,0.9921875}, {0.296875,0.84375}},
+  {{0.296875,0.84375}, {0.57421875,0.6953125}, {0.875,0.5}},
+  {{0.875,0.5}, {1.17578125,0.3046875}, {1.265625,0.15625}},
+  {{1.265625,0.15625}, {1.35546875,0.0078125}, {1,0}},
+</div>
+
+<div id="cubicOp2da">
+{{0.395593847,0.78966024}, {0.576287939,0.496887363}, {0.740740741,0.296296296}},
+{{0.395098308,0.789538003}, {0.634357518,0.657581172}, {0.875,0.5}},
+</div>
+
+<div id="cubicOp3d">
+{{0,1}, {2,3}, {1,0}, {1,0}},
+{{0,1}, {0,1}, {1,0}, {3,2}},
+  {{0,1}, {0.592,1.584}, {0.872,1.664}},
+  {{0.872,1.664}, {1.152,1.744}, {1.216,1.512}},
+  {{1.216,1.512}, {1.28,1.28}, {1.224,0.928}},
+  {{1.224,0.928}, {1.168,0.576}, {1.088,0.296}},
+  {{1.088,0.296}, {1.008,0.016}, {1,0}},
+
+  {{0,1}, {0,0.962962963}, {0.333333333,0.814814815}},
+  {{0.333333333,0.814814815}, {0.666666667,0.666666667}, {1.33333333,0.851851852}},
+  {{1.33333333,0.851851852}, {2,1.03703704}, {3,2}},
+</div>
+
+<div id="cubicOp3da">
+{{1.224,0.928}, {1.17328164,0.604609714}, {1.09894996,0.336624845}},
+{{1.09895195,0.33662897}, {1.22307655,0.2359436}, {1.265625,0.15625}},
+</div>
+
+<div id="cubicOp3db">
+{{x = 1.2071879545809394, y = 0.82163474041730045}, {x = 1.1534203513372994, y = 0.52790870069930229}, {x = 1.0880000000000001, y = 0.29599999999999982}}
+{{x = 1.205732763658403, y = 0.81345617746834109}, {x = 1.267928895828891, y = 0.83008534558465619}, {x = 1.3333333333333333, y = 0.85185185185185175}}
+</div>
+
+<div id="cubicOp3dc">
+part=(1.20718795,0.82163474 1.17452925,0.632190117 1.1284272,0.444233064 1.088,0.296)
+quad=(1.20718795,0.82163474 1.15342035,0.527908701 1.088,0.296)
+part=(1.20573276,0.813456177 1.24719685,0.824565605 1.28973037,0.837317532 1.33333333,0.851851852)
+quad=(1.20573276,0.813456177 1.2679289,0.830085346 1.33333333,0.851851852)
+</div>
+
+<div id="cubicOp3dd">
+{{1.20718795,0.82163474 1.17452925,0.632190117 1.1284272,0.444233064 1.088,0.296)
+{{1.20718795,0.82163474 1.15342035,0.527908701 1.088,0.296)
+{{1.20568441,0.813443223 1.20570053,0.813447541 1.20571665,0.813451859 1.20573276,0.813456177)
+{{1.20568441,0.813443223 1.20570859,0.8134497 1.20573276,0.813456177)
+{{0.33333333333333326, y = 0.81481481481481488}, {x = 0.63396444791444551, y = 0.68743368362444768}, {x = 1.205732763658403, y = 0.81345617746834109}}
+</div>
+
+<div id="cubicOp3de">
+{{1.2071879545809394,0.82163474041730045}, {1.2065040319428038,0.81766753259119995}, {1.2058123269101506,0.81370135061854221}},
+{{1.205684411948591,0.81344322326274499}, {1.2057085875611198,0.81344969999329253}, {1.205732763658403,0.81345617746834109}},
+</div>
+
+<div id="cubicOp7">
+{{0,1}, {3,4}, {1,0}, {3,0}},
+{{0,1}, {0,3}, {1,0}, {4,3}},
+
+  {{0,1}, {0.837764189,1.83435757}, {1.22841861,2.02640973}},
+  {{1.22841861,2.02640973}, {1.61907304,2.21846188}, {1.74657491,1.9930452}},
+  {{1.74657491,1.9930452}, {1.87407679,1.76762853}, {1.92238332,1.34957962}},
+  {{1.92238332,1.34957962}, {1.97220681,0.867804601}, {2.17393047,0.447689071}},
+  {{2.17393047,0.447689071}, {2.37565413,0.0275735418}, {3,0}},
+
+  {{0,1}, {-0.00234073071,1.60655471}, {0.142631845,1.70125304}},
+  {{0.142631845,1.70125304}, {0.28760442,1.79595137}, {0.60797907,1.68776977}},
+  {{0.60797907,1.68776977}, {1.0447864,1.50810914}, {1.87464474,1.63655092}},
+  {{1.87464474,1.63655092}, {2.70450308,1.76499271}, {4,3}},
+</div>
+
+<div id="cubicOp7a">
+{{x = 1.7465749139282332, y = 1.9930452039527999}, {x = 1.8417960084006277, y = 1.8552583419678612}, {x = 1.8799591210677749, y = 1.6157879692142081}, {x = 1.9223833226085514, y = 1.3495796165215643}}
+{{x = 0.6079790696638232, y = 1.6877697663020552}, {x = 0.90321659591661663, y = 1.6123550739533821}, {x = 1.3173732025571312, y = 1.5065640064343382}, {x = 1.8746447406062119, y = 1.636550924974228}}
+
+{{x = 1.7465749139282332, y = 1.9930452039527999}, {x = 1.8740767879671056, y = 1.7676285282679607}, {x = 1.9223833226085514, y = 1.3495796165215643}}
+{{x = 0.6079790696638232, y = 1.6877697663020552}, {x = 1.0447863962878021, y = 1.5081091374717195}, {x = 1.8746447406062119, y = 1.636550924974228}}
+</div>
+
+<div id="cubicOp7b">
+{{x = 1.7465749139282332, y = 1.9930452039527999}, {x = 1.8417960084006277, y = 1.8552583419678612}, {x = 1.8799591210677749, y = 1.6157879692142081}, {x = 1.9223833226085514, y = 1.3495796165215643}}
+{{x = 1.8746447406062119, y = 1.636550924974228}, {x = 2.4319162786552919, y = 1.7665378435141166}, {x = 3.1323027481129411, y = 2.1323027481129406}, {x = 4, y = 3}}
+{{x = 1.7465749139282332, y = 1.9930452039527999}, {x = 1.8740767879671056, y = 1.7676285282679607}, {x = 1.9223833226085514, y = 1.3495796165215643}}
+{{x = 1.8746447406062119, y = 1.636550924974228}, {x = 2.7045030849246219, y = 1.7649927124767357}, {x = 4, y = 3}}
+</div>
+
 </div>
 
 <script type="text/javascript">
 
 var testDivs = [
+    cubicOp7b,
+    cubicOp7a,
+    cubicOp7,
+    cubicOp3de,
+    cubicOp3dd,
+    cubicOp3dc,
+    cubicOp3db,
+    cubicOp3da,
+    cubicOp3d,
+    cubicOp2da,
+    cubicOp2d,
     cubicX,
     x1,
     x2,
@@ -2205,9 +2309,9 @@ function draw(test, title, scale) {
                     xoffset + curve[6] * unit, yoffset + curve[7] * unit);
                 break;
         }
-        if (curves == 2) ctx.strokeStyle = curves ? "red" : "blue";
+        ctx.strokeStyle = drawQuads && drawCubics && curve.length == 6 ? "red" : "black";
         ctx.stroke();
-        if (drawControlLines && curve.length >= 6) {
+        if (drawControlLines && (curve.length == 6 || curve.length == 8)) {
             ctx.strokeStyle = "rgba(0,0,0, 0.3)";
             ctx.beginPath();
             ctx.moveTo(xoffset + curve[0] * unit, yoffset + curve[1] * unit);
@@ -2215,6 +2319,7 @@ function draw(test, title, scale) {
             ctx.lineTo(xoffset + curve[4] * unit, yoffset + curve[5] * unit);
             if (curve.length == 8)
                 ctx.lineTo(xoffset + curve[6] * unit, yoffset + curve[7] * unit);
+            ctx.lineTo(xoffset + curve[0] * unit, yoffset + curve[1] * unit);
             ctx.stroke();
         }
         if (curveT >= 0 && curveT <= 1) {