ccpr: Normalize the cubic inflection function instead of its roots

When solving for KLM, switches back to normalizing the cubic's
inflection function rather than both individual roots. Also performs
some general code clean up for SkClassifyCubic.

Bug: skia:
Change-Id: Id513e7e02c50a8709f3eccf92fad9e5134d73d83
Reviewed-on: https://skia-review.googlesource.com/121201
Reviewed-by: Cary Clark <caryclark@skia.org>
Commit-Queue: Chris Dalton <csmartdalton@google.com>

3 files changed, 101 insertions, 146 deletions

diff --git a/src/core/SkGeometry.cpp b/src/core/SkGeometry.cpp
index d42a3b1961..6ed8911320 100644
--- a/src/core/SkGeometry.cpp
+++ b/src/core/SkGeometry.cpp
@@ -542,63 +542,28 @@ static double calc_dot_cross_cubic(const SkPoint& p0, const SkPoint& p1, const S
     return (xComp + yComp + wComp);
 }
 
-// Calc coefficients of I(s,t) where roots of I are inflection points of curve
-// I(s,t) = t*(3*d0*s^2 - 3*d1*s*t + d2*t^2)
-// d0 = a1 - 2*a2+3*a3
-// d1 = -a2 + 3*a3
-// d2 = 3*a3
-// a1 = p0 . (p3 x p2)
-// 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], double d[4]) {
-    const double a1 = calc_dot_cross_cubic(p[0], p[3], p[2]);
-    const double a2 = calc_dot_cross_cubic(p[1], p[0], p[3]);
-    const double a3 = calc_dot_cross_cubic(p[2], p[1], p[0]);
-
-    d[3] = 3 * a3;
-    d[2] = d[3] - a2;
-    d[1] = d[2] - a2 + a1;
-    d[0] = 0;
-}
-
-static void normalize_t_s(double t[], double s[], int count) {
-    // Keep the exponents at or below zero to avoid overflow down the road.
-    for (int i = 0; i < count; ++i) {
-        SkASSERT(0 != s[i]); // classify_cubic should not call this method when s[i] is 0 or NaN.
-
-        uint64_t bitsT, bitsS;
-        memcpy(&bitsT, &t[i], sizeof(double));
-        memcpy(&bitsS, &s[i], sizeof(double));
-
-        uint64_t maxExponent = SkTMax(bitsT & 0x7ff0000000000000, bitsS & 0x7ff0000000000000);
+// Returns a positive power of 2 that, when multiplied by n, and excepting the two edge cases listed
+// below, shifts the exponent of n to yield a magnitude somewhere inside [1..2).
+// Returns 2^1023 if abs(n) < 2^-1022 (including 0).
+// Returns NaN if n is Inf or NaN.
+inline static double previous_inverse_pow2(double n) {
+    uint64_t bits;
+    memcpy(&bits, &n, sizeof(double));
+    bits = ((1023llu*2 << 52) + ((1llu << 52) - 1)) - bits; // exp=-exp
+    bits &= (0x7ffllu) << 52; // mantissa=1.0, sign=0
+    memcpy(&n, &bits, sizeof(double));
+    return n;
+}
+
+inline static void write_cubic_inflection_roots(double t0, double s0, double t1, double s1,
+                                                double* t, double* s) {
+    t[0] = t0;
+    s[0] = s0;
 
-#ifdef SK_DEBUG
-        uint64_t maxExponentValue = maxExponent >> 52;
-        // Ensure max(absT,absS) is NOT in denormalized form. SkClassifyCubic is given fp32 points,
-        // and does not call this method when s==0, so this should never happen.
-        SkASSERT(0 != maxExponentValue);
-        // Ensure 1/max(absT,absS) will NOT be in denormalized form. SkClassifyCubic is given fp32
-        // points, so this should never happen.
-        SkASSERT(2046 != maxExponentValue);
-#endif
-
-        // Pick a normalizer that scales the larger exponent to 1 (aka 1023 in biased form), but
-        // does NOT change the mantissa (thus preserving accuracy).
-        double normalizer;
-        uint64_t normalizerExponent = (uint64_t(1023 * 2) << 52) - maxExponent;
-        memcpy(&normalizer, &normalizerExponent, sizeof(double));
-
-        t[i] *= normalizer;
-        s[i] *= normalizer;
-    }
-}
-
-static void sort_and_orient_t_s(double t[2], double s[2]) {
     // This copysign/abs business orients the implicit function so positive values are always on the
     // "left" side of the curve.
-    t[1] = -copysign(t[1], t[1] * s[1]);
-    s[1] = -fabs(s[1]);
+    t[1] = -copysign(t1, t1 * s1);
+    s[1] = -fabs(s1);
 
     // Ensure t[0]/s[0] <= t[1]/s[1] (s[1] is negative from above).
     if (copysign(s[1], s[0]) * t[0] > -fabs(s[0]) * t[1]) {
@@ -607,78 +572,78 @@ static void sort_and_orient_t_s(double t[2], double s[2]) {
     }
 }
 
-// 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:
-// 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(const double d[4], double t[2], double s[2]) {
-    if (0 == d[1]) {
-        if (0 == d[2]) {
-            if (t && s) {
-                t[0] = t[1] = 1;
-                s[0] = s[1] = 0; // infinity
-            }
-            return 0 == d[3] ? SkCubicType::kLineOrPoint : SkCubicType::kQuadratic;
-        }
-        if (t && s) {
-            t[0] = d[3];
-            s[0] = 3 * d[2];
-            normalize_t_s(t, s, 1);
-            t[1] = 1;
-            s[1] = 0; // infinity
-        }
-        return SkCubicType::kCuspAtInfinity;
+SkCubicType SkClassifyCubic(const SkPoint P[4], double t[2], double s[2], double d[4]) {
+    // Find the cubic's inflection function, I = [T^3  -3T^2  3T  -1] dot D. (D0 will always be 0
+    // for integral cubics.)
+    //
+    // 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
+    double A1 = calc_dot_cross_cubic(P[0], P[3], P[2]);
+    double A2 = calc_dot_cross_cubic(P[1], P[0], P[3]);
+    double A3 = calc_dot_cross_cubic(P[2], P[1], P[0]);
+
+    double D3 = 3 * A3;
+    double D2 = D3 - A2;
+    double D1 = D2 - A2 + A1;
+
+    // Shift the exponents in D so the largest magnitude falls somewhere in 1..2. This protects us
+    // from overflow down the road while solving for roots and KLM functionals.
+    double Dmax = std::max(std::max(fabs(D1), fabs(D2)), fabs(D3));
+    double norm = previous_inverse_pow2(Dmax);
+    D1 *= norm;
+    D2 *= norm;
+    D3 *= norm;
+
+    if (d) {
+        d[3] = D3;
+        d[2] = D2;
+        d[1] = D1;
+        d[0] = 0;
     }
 
-    const double discr = 3 * d[2] * d[2] - 4 * d[1] * d[3];
-    if (discr > 0) {
-        if (t && s) {
-            const double q = 3 * d[2] + copysign(sqrt(3 * discr), d[2]);
-            t[0] = q;
-            s[0] = 6 * d[1];
-            t[1] = 2 * d[3];
-            s[1] = q;
-            normalize_t_s(t, s, 2);
-            sort_and_orient_t_s(t, s);
-        }
-        return SkCubicType::kSerpentine;
-    } else if (discr < 0) {
-        if (t && s) {
-            const double q = d[2] + copysign(sqrt(-discr), d[2]);
-            t[0] = q;
-            s[0] = 2 * d[1];
-            t[1] = 2 * (d[2] * d[2] - d[3] * d[1]);
-            s[1] = d[1] * q;
-            normalize_t_s(t, s, 2);
-            sort_and_orient_t_s(t, s);
+    // Now use the inflection function to classify the cubic.
+    //
+    // See "Resolution Independent Curve Rendering using Programmable Graphics Hardware",
+    // 4.4 Integral Cubics:
+    //
+    // https://www.microsoft.com/en-us/research/wp-content/uploads/2005/01/p1000-loop.pdf
+    if (0 != D1) {
+        double discr = 3*D2*D2 - 4*D1*D3;
+        if (discr > 0) { // Serpentine.
+            if (t && s) {
+                double q = 3*D2 + copysign(sqrt(3*discr), D2);
+                write_cubic_inflection_roots(q, 6*D1, 2*D3, q, t, s);
+            }
+            return SkCubicType::kSerpentine;
+        } else if (discr < 0) { // Loop.
+            if (t && s) {
+                double q = D2 + copysign(sqrt(-discr), D2);
+                write_cubic_inflection_roots(q, 2*D1, 2*(D2*D2 - D3*D1), D1*q, t, s);
+            }
+            return SkCubicType::kLoop;
+        } else { // Cusp.
+            if (t && s) {
+                write_cubic_inflection_roots(D2, 2*D1, D2, 2*D1, t, s);
+            }
+            return SkCubicType::kLocalCusp;
         }
-        return SkCubicType::kLoop;
     } else {
-        if (t && s) {
-            t[0] = d[2];
-            s[0] = 2 * d[1];
-            normalize_t_s(t, s, 1);
-            t[1] = t[0];
-            s[1] = s[0];
-            sort_and_orient_t_s(t, s);
+        if (0 != D2) { // Cusp at T=infinity.
+            if (t && s) {
+                write_cubic_inflection_roots(D3, 3*D2, 1, 0, t, s); // T1=infinity.
+            }
+            return SkCubicType::kCuspAtInfinity;
+        } else { // Degenerate.
+            if (t && s) {
+                write_cubic_inflection_roots(1, 0, 1, 0, t, s); // T0=T1=infinity.
+            }
+            return 0 != D3 ? SkCubicType::kQuadratic : SkCubicType::kLineOrPoint;
         }
-        return SkCubicType::kLocalCusp;
     }
 }
 
-SkCubicType SkClassifyCubic(const SkPoint src[4], double t[2], double s[2], double d[4]) {
-    double localD[4];
-    double* dd = d ? d : localD;
-    calc_cubic_inflection_func(src, dd);
-    return classify_cubic(dd, t, s);
-}
-
 template <typename T> void bubble_sort(T array[], int count) {
     for (int i = count - 1; i > 0; --i)
         for (int j = i; j > 0; --j)
diff --git a/src/gpu/ccpr/GrCCCubicShader.cpp b/src/gpu/ccpr/GrCCCubicShader.cpp
index f8a8352383..9a57c1cd44 100644
--- a/src/gpu/ccpr/GrCCCubicShader.cpp
+++ b/src/gpu/ccpr/GrCCCubicShader.cpp
@@ -15,28 +15,6 @@ using Shader = GrCCCoverageProcessor::Shader;
 
 void GrCCCubicShader::emitSetupCode(GrGLSLVertexGeoBuilder* s, const char* pts,
                                     const char* wind, const char** /*outHull4*/) const {
-    // Define a function that normalizes the homogeneous coordinates T=t/s in order to avoid
-    // exponent overflow.
-    SkString normalizeHomogCoordFn;
-    GrShaderVar coord("coord", kFloat2_GrSLType);
-    s->emitFunction(kFloat2_GrSLType, "normalize_homogeneous_coord", 1, &coord,
-                    s->getProgramBuilder()->shaderCaps()->fpManipulationSupport()
-                            // Exponent manipulation version: Scale the exponents so the larger
-                            // component has a magnitude in 1..2.
-                            // (Neither component should be infinity because ccpr crops big paths.)
-                            ? "int exp;"
-                              "frexp(max(abs(coord.t), abs(coord.s)), exp);"
-                              "return coord * ldexp(1, 1 - exp);"
-
-                            // Division version: Divide by the component with the larger magnitude.
-                            // (Both should not be 0 because ccpr catches degenerate cubics.)
-                            : "bool swap = abs(coord.t) > abs(coord.s);"
-                              "coord = swap ? coord.ts : coord;"
-                              "coord = float2(1, coord.t/coord.s);"
-                              "return swap ? coord.ts : coord;",
-
-                    &normalizeHomogCoordFn);
-
     // Find the cubic's power basis coefficients.
     s->codeAppendf("float2x4 C = float4x4(-1,  3, -3,  1, "
                                          " 3, -6,  3,  0, "
@@ -48,6 +26,21 @@ void GrCCCubicShader::emitSetupCode(GrGLSLVertexGeoBuilder* s, const char* pts,
     s->codeAppend ("float D2 = -determinant(float2x2(C[0].xz, C[1].xz));");
     s->codeAppend ("float D1 = +determinant(float2x2(C));");
 
+    // Shift the exponents in D so the largest magnitude falls somewhere in 1..2. This protects us
+    // from overflow while solving for roots and KLM functionals.
+    s->codeAppend ("float Dmax = max(max(abs(D1), abs(D2)), abs(D3));");
+    s->codeAppend ("float norm;");
+    if (s->getProgramBuilder()->shaderCaps()->fpManipulationSupport()) {
+        s->codeAppend ("int exp;");
+        s->codeAppend ("frexp(Dmax, exp);");
+        s->codeAppend ("norm = ldexp(1, 1 - exp);");
+    } else {
+        s->codeAppend ("norm = 1/Dmax;"); // Dmax will not be 0 because we cull line cubics on CPU.
+    }
+    s->codeAppend ("D3 *= norm;");
+    s->codeAppend ("D2 *= norm;");
+    s->codeAppend ("D1 *= norm;");
+
     // Calculate the KLM matrix.
     s->declareGlobal(fKLMMatrix);
     s->codeAppend ("float discr = 3*D2*D2 - 4*D1*D3;");
@@ -56,10 +49,9 @@ void GrCCCubicShader::emitSetupCode(GrGLSLVertexGeoBuilder* s, const char* pts,
     s->codeAppend ("q = x*D2 + (D2 >= 0 ? q : -q);");
 
     s->codeAppend ("float2 l, m;");
-    s->codeAppendf("l.ts = %s(float2(q, 2*x * D1));", normalizeHomogCoordFn.c_str());
-    s->codeAppendf("m.ts = %s(float2(2, q) * (discr >= 0 ? float2(D3, 1) "
-                                                        ": float2(D2*D2 - D3*D1, D1)));",
-                                                        normalizeHomogCoordFn.c_str());
+    s->codeAppend ("l.ts = float2(q, 2*x * D1);");
+    s->codeAppend ("m.ts = float2(2, q) * (discr >= 0 ? float2(D3, 1) "
+                                                     ": float2(D2*D2 - D3*D1, D1));");
 
     s->codeAppend ("float4 K;");
     s->codeAppend ("float4 lm = l.sstt * m.stst;");
diff --git a/src/pathops/SkPathOpsCubic.cpp b/src/pathops/SkPathOpsCubic.cpp
index eb2ddb7c7d..8d3baa038b 100644
--- a/src/pathops/SkPathOpsCubic.cpp
+++ b/src/pathops/SkPathOpsCubic.cpp
@@ -253,9 +253,7 @@ int SkDCubic::ComplexBreak(const SkPoint pointsPtr[4], SkScalar* t) {
         case SkCubicType::kLoop: {
             const double &td = tt[0], &te = tt[1], &sd = ss[0], &se = ss[1];
             if (roughly_between(0, td, sd) && roughly_between(0, te, se)) {
-                SkASSERT(roughly_between(0, td/sd, 1) && roughly_between(0, te/se, 1));
                 t[0] = static_cast<SkScalar>((td * se + te * sd) / (2 * sd * se));
-                SkASSERT(roughly_between(0, *t, 1));
                 return (int) (t[0] > 0 && t[0] < 1);
             }
         }