Improve accuracy of cubic classification

- Updates the logic to reflect the Loop-Blinn paper instead of the GPU
  gems website.
- Removes the threshold for detecting local cusps. The serpentine
  codepath works for these cusps anyway, so what we really want to know
  is whether the discriminant is negative.
- Makes sure to not scale the inflection function by 1/0.
- Shifts the inflection function coefficients in d[] so they match the
  paper.
- Stores the cubic discriminant in d[0].

Bug: skia:
Change-Id: I909a522a0fd27c9c8dfbc27d968bc43eeb7a416f
Reviewed-on: https://skia-review.googlesource.com/13304
Reviewed-by: Greg Daniel <egdaniel@google.com>
Commit-Queue: Chris Dalton <csmartdalton@google.com>

6 files changed, 97 insertions, 85 deletions

diff --git a/bench/CubicKLMBench.cpp b/bench/CubicKLMBench.cpp
index 3c8f740bc7..1cdb068c08 100644
--- a/bench/CubicKLMBench.cpp
+++ b/bench/CubicKLMBench.cpp
@@ -22,12 +22,12 @@ public:
         fPoints[3].set(x3, y3);
 
         fName = "cubic_klm_";
-        SkScalar d[3];
+        SkScalar d[4];
         switch (SkClassifyCubic(fPoints, d)) {
-            case kSerpentine_SkCubicType:
+            case SkCubicType::kSerpentine:
                 fName.append("serp");
                 break;
-            case kLoop_SkCubicType:
+            case SkCubicType::kLoop:
                 fName.append("loop");
                 break;
             default:
diff --git a/src/core/SkGeometry.cpp b/src/core/SkGeometry.cpp
index c2e9b60fb8..17ff43cb8b 100644
--- a/src/core/SkGeometry.cpp
+++ b/src/core/SkGeometry.cpp
@@ -531,30 +531,34 @@ int SkChopCubicAtInflections(const SkPoint src[], SkPoint dst[10]) {
     return count + 1;
 }
 
-// See http://http.developer.nvidia.com/GPUGems3/gpugems3_ch25.html (from the book GPU Gems 3)
-// discr(I) = d0^2 * (3*d1^2 - 4*d0*d2)
+// See "Resolution Independent Curve Rendering using Programmable Graphics Hardware"
+// https://www.microsoft.com/en-us/research/wp-content/uploads/2005/01/p1000-loop.pdf
+// discr(I) = 3*d2^2 - 4*d1*d3
 // Classification:
-// discr(I) > 0        Serpentine
-// discr(I) = 0        Cusp
-// discr(I) < 0        Loop
-// d0 = d1 = 0         Quadratic
-// d0 = d1 = d2 = 0    Line
-// p0 = p1 = p2 = p3   Point
-static SkCubicType classify_cubic(const SkPoint p[4], const SkScalar d[3]) {
-    if (p[0] == p[1] && p[0] == p[2] && p[0] == p[3]) {
-        return kPoint_SkCubicType;
-    }
-    const SkScalar discr = d[0] * d[0] * (3.f * d[1] * d[1] - 4.f * d[0] * d[2]);
-    if (discr > SK_ScalarNearlyZero) {
-        return kSerpentine_SkCubicType;
-    } else if (discr < -SK_ScalarNearlyZero) {
-        return kLoop_SkCubicType;
+// d1 != 0, discr(I) > 0        Serpentine
+// d1 != 0, discr(I) < 0        Loop
+// d1 != 0, discr(I) = 0        Cusp (with inflection at infinity)
+// d1 = 0, d2 != 0              Cusp (with cusp at infinity)
+// d1 = d2 = 0, d3 != 0         Quadratic
+// d1 = d2 = d3 = 0             Line or Point
+static SkCubicType classify_cubic(SkScalar d[4]) {
+    if (!SkScalarNearlyZero(d[1])) {
+        d[0] = 3 * d[2] * d[2] - 4 * d[1] * d[3];
+        if (d[0] > 0) {
+            return SkCubicType::kSerpentine;
+        } else if (d[0] < 0) {
+            return SkCubicType::kLoop;
+        } else {
+            SkASSERT(0 == d[0]); // Detect NaN.
+            return SkCubicType::kLocalCusp;
+        }
     } else {
-        if (SkScalarAbs(d[0]) < SK_ScalarNearlyZero && SkScalarAbs(d[1]) < SK_ScalarNearlyZero) {
-            return ((SkScalarAbs(d[2]) < SK_ScalarNearlyZero) ? kLine_SkCubicType
-                                                              : kQuadratic_SkCubicType);
+        if (!SkScalarNearlyZero(d[2])) {
+            return SkCubicType::kInfiniteCusp;
+        } else if (!SkScalarNearlyZero(d[3])) {
+            return SkCubicType::kQuadratic;
         } else {
-            return kCusp_SkCubicType;
+            return SkCubicType::kLineOrPoint;
         }
     }
 }
@@ -577,7 +581,7 @@ static SkScalar calc_dot_cross_cubic(const SkPoint& p0, const SkPoint& p1, const
 // a2 = p1 . (p0 x p3)
 // a3 = p2 . (p1 x p0)
 // Places the values of d1, d2, d3 in array d passed in
-static void calc_cubic_inflection_func(const SkPoint p[4], SkScalar d[3]) {
+static void calc_cubic_inflection_func(const SkPoint p[4], SkScalar d[4]) {
     SkScalar a1 = calc_dot_cross_cubic(p[0], p[3], p[2]);
     SkScalar a2 = calc_dot_cross_cubic(p[1], p[0], p[3]);
     SkScalar a3 = calc_dot_cross_cubic(p[2], p[1], p[0]);
@@ -586,19 +590,22 @@ static void calc_cubic_inflection_func(const SkPoint p[4], SkScalar d[3]) {
     SkScalar max = SkScalarAbs(a1);
     max = SkMaxScalar(max, SkScalarAbs(a2));
     max = SkMaxScalar(max, SkScalarAbs(a3));
-    max = 1.f/max;
-    a1 = a1 * max;
-    a2 = a2 * max;
-    a3 = a3 * max;
+    if (0 != max) {
+        max = 1.f/max;
+        a1 = a1 * max;
+        a2 = a2 * max;
+        a3 = a3 * max;
+    }
 
-    d[2] = 3.f * a3;
-    d[1] = d[2] - a2;
-    d[0] = d[1] - a2 + a1;
+    d[3] = 3.f * a3;
+    d[2] = d[3] - a2;
+    d[1] = d[2] - a2 + a1;
+    d[0] = 0;
 }
 
-SkCubicType SkClassifyCubic(const SkPoint src[4], SkScalar d[3]) {
+SkCubicType SkClassifyCubic(const SkPoint src[4], SkScalar d[4]) {
     calc_cubic_inflection_func(src, d);
-    return classify_cubic(src, d);
+    return classify_cubic(d);
 }
 
 template <typename T> void bubble_sort(T array[], int count) {
diff --git a/src/core/SkGeometry.h b/src/core/SkGeometry.h
index 55d763b967..91b4d2d895 100644
--- a/src/core/SkGeometry.h
+++ b/src/core/SkGeometry.h
@@ -158,19 +158,26 @@ int SkChopCubicAtMaxCurvature(const SkPoint src[4], SkPoint dst[13],
 bool SkChopMonoCubicAtX(SkPoint src[4], SkScalar y, SkPoint dst[7]);
 bool SkChopMonoCubicAtY(SkPoint src[4], SkScalar x, SkPoint dst[7]);
 
-enum SkCubicType {
-    kSerpentine_SkCubicType,
-    kCusp_SkCubicType,
-    kLoop_SkCubicType,
-    kQuadratic_SkCubicType,
-    kLine_SkCubicType,
-    kPoint_SkCubicType
+enum class SkCubicType {
+    kSerpentine,
+    kLoop,
+    kLocalCusp,     // Cusp at a non-infinite parameter value with an inflection at t=infinity.
+    kInfiniteCusp,  // Cusp with a cusp at t=infinity and a local inflection.
+    kQuadratic,
+    kLineOrPoint
 };
 
-/** Returns the cubic classification. Pass scratch storage for computing inflection data,
-    which can be used with additional work to find the loop intersections and so on.
+/** Returns the cubic classification.
+
+    d[] is filled with the cubic inflection function coefficients. Furthermore, since d0 is always
+    zero for integral curves, if the cubic type is kSerpentine, kLoop, or kLocalCusp then d[0] will
+    instead contain the cubic discriminant: 3*d2^2 - 4*d1*d3.
+
+    See "Resolution Independent Curve Rendering using Programmable Graphics Hardware",
+    4.2 Curve Categorization
+    https://www.microsoft.com/en-us/research/wp-content/uploads/2005/01/p1000-loop.pdf
 */
-SkCubicType SkClassifyCubic(const SkPoint p[4], SkScalar inflection[3]);
+SkCubicType SkClassifyCubic(const SkPoint p[4], SkScalar d[4]);
 
 ///////////////////////////////////////////////////////////////////////////////
 
diff --git a/src/gpu/GrPathUtils.cpp b/src/gpu/GrPathUtils.cpp
index 536f3c3acb..57bb86fc5b 100644
--- a/src/gpu/GrPathUtils.cpp
+++ b/src/gpu/GrPathUtils.cpp
@@ -660,16 +660,17 @@ static void negate_kl(SkMatrix* klm) {
     }
 }
 
-static void calc_serp_klm(const SkPoint pts[4], const SkScalar d[3], SkMatrix* klm) {
+static void calc_serp_klm(const SkPoint pts[4], const SkScalar d[4], SkMatrix* klm) {
     SkMatrix CIT;
     int skipCol = calc_inverse_transpose_power_basis_matrix(pts, &CIT);
 
-    const SkScalar root = SkScalarSqrt(9 * d[1] * d[1] - 12 * d[0] * d[2]);
+    SkASSERT(d[0] >= 0);
+    const SkScalar root = SkScalarSqrt(3 * d[0]);
 
-    const SkScalar tl = 3 * d[1] + root;
-    const SkScalar sl = 6 * d[0];
-    const SkScalar tm = 3 * d[1] - root;
-    const SkScalar sm = 6 * d[0];
+    const SkScalar tl = 3 * d[2] + root;
+    const SkScalar sl = 6 * d[1];
+    const SkScalar tm = 3 * d[2] - root;
+    const SkScalar sm = 6 * d[1];
 
     SkMatrix klmCoeffs;
     int col = 0;
@@ -699,10 +700,10 @@ static void calc_serp_klm(const SkPoint pts[4], const SkScalar d[3], SkMatrix* k
 
     klm->setConcat(klmCoeffs, CIT);
 
-    // If d0 > 0 we need to flip the orientation of our curve
+    // If d1 > 0 we need to flip the orientation of our curve
     // This is done by negating the k and l values
     // We want negative distance values to be on the inside
-    if (d[0] > 0) {
+    if (d[1] > 0) {
         negate_kl(klm);
     }
 }
@@ -844,16 +845,17 @@ int GrPathUtils::chopCubicAtLoopIntersection(const SkPoint src[4], SkPoint dst[1
     // Homogeneous parametric values at the loop double point.
     SkScalar td, sd, te, se;
 
-    SkScalar d[3];
+    SkScalar d[4];
     SkCubicType cType = SkClassifyCubic(src, d);
 
     int chop_count = 0;
-    if (kLoop_SkCubicType == cType) {
-        SkScalar tempSqrt = SkScalarSqrt(4.f * d[0] * d[2] - 3.f * d[1] * d[1]);
-        td = d[1] + tempSqrt;
-        sd = 2.f * d[0];
-        te = d[1] - tempSqrt;
-        se = 2.f * d[0];
+    if (SkCubicType::kLoop == cType) {
+        SkASSERT(d[0] < 0);
+        const SkScalar tempSqrt = SkScalarSqrt(-d[0]);
+        td = d[2] + tempSqrt;
+        sd = 2.f * d[1];
+        te = d[2] - tempSqrt;
+        se = 2.f * d[1];
 
         t1 = td / sd;
         t2 = te / se;
@@ -898,26 +900,20 @@ int GrPathUtils::chopCubicAtLoopIntersection(const SkPoint src[4], SkPoint dst[1
 
     if (klm) {
         switch (cType) {
-            case kSerpentine_SkCubicType:
+            case SkCubicType::kSerpentine:
+            case SkCubicType::kLocalCusp:
                 calc_serp_klm(src, d, klm);
                 break;
-            case kLoop_SkCubicType:
-                calc_loop_klm(src, d[0], td, sd, te, se, klm);
+            case SkCubicType::kLoop:
+                calc_loop_klm(src, d[1], td, sd, te, se, klm);
                 break;
-            case kCusp_SkCubicType:
-                if (0 != d[0]) {
-                    // FIXME: SkClassifyCubic has a tolerance, but we need an exact classification
-                    // here to be sure we won't get a negative in the square root.
-                    calc_serp_klm(src, d, klm);
-                } else {
-                    calc_inf_cusp_klm(src, d[1], d[2], klm);
-                }
+            case SkCubicType::kInfiniteCusp:
+                calc_inf_cusp_klm(src, d[2], d[3], klm);
                 break;
-            case kQuadratic_SkCubicType:
-                calc_quadratic_klm(src, d[2], klm);
+            case SkCubicType::kQuadratic:
+                calc_quadratic_klm(src, d[3], klm);
                 break;
-            case kLine_SkCubicType:
-            case kPoint_SkCubicType:
+            case SkCubicType::kLineOrPoint:
                 calc_line_klm(src, klm);
                 break;
         };
diff --git a/src/pathops/SkPathOpsCubic.cpp b/src/pathops/SkPathOpsCubic.cpp
index d842e2cc0c..1e74eb0afd 100644
--- a/src/pathops/SkPathOpsCubic.cpp
+++ b/src/pathops/SkPathOpsCubic.cpp
@@ -248,17 +248,18 @@ int SkDCubic::ComplexBreak(const SkPoint pointsPtr[4], SkScalar* t) {
     if (cubic.monotonicInX() && cubic.monotonicInY()) {
         return 0;
     }
-    SkScalar d[3];
+    SkScalar d[4];
     SkCubicType cubicType = SkClassifyCubic(pointsPtr, d);
     switch (cubicType) {
-        case kLoop_SkCubicType: {
+        case SkCubicType::kLoop: {
             // crib code from gpu path utils that finds t values where loop self-intersects
             // use it to find mid of t values which should be a friendly place to chop
-            SkScalar tempSqrt = SkScalarSqrt(4.f * d[0] * d[2] - 3.f * d[1] * d[1]);
-            SkScalar ls = d[1] - tempSqrt;
-            SkScalar lt = 2.f * d[0];
-            SkScalar ms = d[1] + tempSqrt;
-            SkScalar mt = 2.f * d[0];
+            SkASSERT(d[0] < 0);
+            SkScalar tempSqrt = SkScalarSqrt(-d[0]);
+            SkScalar ls = d[2] - tempSqrt;
+            SkScalar lt = 2.f * d[1];
+            SkScalar ms = d[2] + tempSqrt;
+            SkScalar mt = 2.f * d[1];
             if (roughly_between(0, ls, lt) && roughly_between(0, ms, mt)) {
                 ls = ls / lt;
                 ms = ms / mt;
@@ -269,8 +270,9 @@ int SkDCubic::ComplexBreak(const SkPoint pointsPtr[4], SkScalar* t) {
             }
         }
         // fall through if no t value found
-        case kSerpentine_SkCubicType:
-        case kCusp_SkCubicType: {
+        case SkCubicType::kSerpentine:
+        case SkCubicType::kLocalCusp:
+        case SkCubicType::kInfiniteCusp: {
             double inflectionTs[2];
             int infTCount = cubic.findInflections(inflectionTs);
             double maxCurvature[3];
diff --git a/tests/PathOpsCubicIntersectionTest.cpp b/tests/PathOpsCubicIntersectionTest.cpp
index 638ecb008c..66becf304a 100644
--- a/tests/PathOpsCubicIntersectionTest.cpp
+++ b/tests/PathOpsCubicIntersectionTest.cpp
@@ -646,11 +646,11 @@ static void selfOneOff(skiatest::Reporter* reporter, int index) {
         c[i] = cubic.fPts[i].asSkPoint();
     }
     SkScalar loopT[3];
-    SkScalar d[3];
+    SkScalar d[4];
     SkCubicType cubicType = SkClassifyCubic(c, d);
     int breaks = SkDCubic::ComplexBreak(c, loopT);
     SkASSERT(breaks < 2);
-    if (breaks && cubicType == SkCubicType::kLoop_SkCubicType) {
+    if (breaks && cubicType == SkCubicType::kLoop) {
         SkIntersections i;
         SkPoint twoCubics[7];
         SkChopCubicAt(c, twoCubics, loopT[0]);