add tests for cubicmap

check for some special cases:
- degenerate or simple cubic
- 0/0

Bug: skia:
Change-Id: Ie978caf9d862755d9693768695bf84eff9ec80e4
Reviewed-on: https://skia-review.googlesource.com/147112
Commit-Queue: Mike Reed <reed@google.com>
Commit-Queue: Florin Malita <fmalita@chromium.org>
Reviewed-by: Florin Malita <fmalita@chromium.org>
Auto-Submit: Mike Reed <reed@google.com>

gn/tests.gni

@@ -46,6 +46,7 @@ tests_sources = [
   "$_tests/ColorTest.cpp",
   "$_tests/CopySurfaceTest.cpp",
   "$_tests/CTest.cpp",
+  "$_tests/CubicMapTest.cpp",
   "$_tests/DashPathEffectTest.cpp",
   "$_tests/DataRefTest.cpp",
   "$_tests/DefaultPathRendererTest.cpp",

src/core/SkCubicMap.cpp

@@ -7,7 +7,12 @@
 
 #include "SkCubicMap.h"
 #include "SkNx.h"
+
+//#define CUBICMAP_TRACK_MAX_ERROR
+
+#ifdef CUBICMAP_TRACK_MAX_ERROR
 #include "../../src/pathops/SkPathOpsCubic.h"
+#endif
 
 static float eval_poly3(float a, float b, float c, float d, float t) {
     return ((a * t + b) * t + c) * t + d;
@@ -31,6 +36,19 @@ static bool valid(float r) {
 }
 #endif
 
+static inline bool nearly_zero(SkScalar x) {
+    SkASSERT(x >= 0);
+    return x <= 0.0000000001f;
+}
+
+static inline bool delta_nearly_zero(float delta) {
+    return sk_float_abs(delta) <= 0.0001f;
+}
+
+#ifdef CUBICMAP_TRACK_MAX_ERROR
+    static int max_iters;
+#endif
+
 /*
  *  TODO: will this be faster if we algebraically compute the polynomials for the numer and denom
  *        rather than compute them in parts?
@@ -39,29 +57,48 @@ static bool valid(float r) {
  *        iterations (divides)
  */
 static float solve_nice_cubic_halley(float A, float B, float C, float D) {
-    const int MAX_ITERS = 3;    // 3 is accurate to 0.0000112 (anecdotally)
-                                // 2 is accurate to 0.0188965 (anecdotally)
+#if 0
+    A = 3.99999;
+    B = -5.99999;
+    C = 2.99999;
+    D = -0.5;
+#endif
+    const int MAX_ITERS = 8;
     const float A3 = 3 * A;
     const float B2 = B + B;
 
     float t = guess_nice_cubic_root(A, B, C, D);
-    for (int iters = 0; iters < MAX_ITERS; ++iters) {
+    int iters = 0;
+    for (; iters < MAX_ITERS; ++iters) {
         float f = eval_poly3(A, B, C, D, t);    // f   = At^3 + Bt^2 + Ct + D
         float fp = eval_poly2(A3, B2, C, t);    // f'  = 3At^2 + 2Bt + C
         float fpp = eval_poly1(A3 + A3, B2, t); // f'' = 6At + 2B
 
         float numer = 2 * fp * f;
+        if (numer == 0) {
+            break;
+        }
         float denom = 2 * fp * fp - f * fpp;
         float delta = numer / denom;
+//      SkDebugf("[%d] delta %g t %g\n", iters, delta, t);
+        if (delta_nearly_zero(delta)) {
+            break;
+        }
         float new_t = t - delta;
         SkASSERT(valid(new_t));
         t = new_t;
     }
     SkASSERT(valid(t));
+#ifdef CUBICMAP_TRACK_MAX_ERROR
+    if (iters > max_iters) {
+        max_iters = iters;
+        SkDebugf("max_iters %d\n", max_iters);
+    }
+#endif
     return t;
 }
 
-#ifdef SK_DEBUG
+#ifdef CUBICMAP_TRACK_MAX_ERROR
 static float compute_slow(float A, float B, float C, float x) {
     double roots[3];
     SkDEBUGCODE(int count =) SkDCubic::RootsValidT(A, B, C, -x, roots);
@@ -73,11 +110,11 @@ static float max_err;
 #endif
 
 static float compute_t_from_x(float A, float B, float C, float x) {
-#ifdef SK_DEBUG
+#ifdef CUBICMAP_TRACK_MAX_ERROR
     float answer = compute_slow(A, B, C, x);
 #endif
     float answer2 = solve_nice_cubic_halley(A, B, C, -x);
-#ifdef SK_DEBUG
+#ifdef CUBICMAP_TRACK_MAX_ERROR
     float err = sk_float_abs(answer - answer2);
     if (err > max_err) {
         max_err = err;
@@ -89,11 +126,29 @@ static float compute_t_from_x(float A, float B, float C, float x) {
 
 float SkCubicMap::computeYFromX(float x) const {
     SkASSERT(valid(x));
-    float t = compute_t_from_x(fCoeff[0].fX, fCoeff[1].fX, fCoeff[2].fX, x);
+
+    if (nearly_zero(x) || nearly_zero(1 - x)) {
+        return x;
+    }
+    if (fType == kLine_Type) {
+        return x;
+    }
+    float t;
+    if (fType == kCubeRoot_Type) {
+        t = sk_float_pow(x / fCoeff[0].fX, 1.0f / 3);
+    } else {
+        t = compute_t_from_x(fCoeff[0].fX, fCoeff[1].fX, fCoeff[2].fX, x);
+    }
     float a = fCoeff[0].fY;
     float b = fCoeff[1].fY;
     float c = fCoeff[2].fY;
-    return ((a * t + b) * t + c) * t;
+    float y = ((a * t + b) * t + c) * t;
+    SkASSERT(y >= 0);
+    return std::min(y, 1.0f);
+}
+
+static inline bool coeff_nearly_zero(float delta) {
+    return sk_float_abs(delta) <= 0.0000001f;
 }
 
 void SkCubicMap::setPts(SkPoint p1, SkPoint p2) {
@@ -106,6 +161,13 @@ void SkCubicMap::setPts(SkPoint p1, SkPoint p2) {
     (Sk2s(1) + s1 - s2).store(&fCoeff[0]);
     (s2 - s1 - s1).store(&fCoeff[1]);
     s1.store(&fCoeff[2]);
+
+    fType = kSolver_Type;
+    if (SkScalarNearlyEqual(p1.fX, p1.fY) && SkScalarNearlyEqual(p2.fX, p2.fY)) {
+        fType = kLine_Type;
+    } else if (coeff_nearly_zero(fCoeff[1].fX) && coeff_nearly_zero(fCoeff[2].fX)) {
+        fType = kCubeRoot_Type;
+    }
 }
 
 SkPoint SkCubicMap::computeFromT(float t) const {

src/core/SkCubicMap.h

@@ -25,15 +25,21 @@ class SkCubicMap {
     SkCubicMap(SkPoint p1, SkPoint p2) {
         this->setPts(p1, p2);
     }
+
     void setPts(SkPoint p1, SkPoint p2);
 
     float computeYFromX(float x) const;
 
-    // is this needed?
     SkPoint computeFromT(float t) const;
 
 private:
+    enum Type {
+        kLine_Type,     // x == y
+        kCubeRoot_Type, // At^3 == x
+        kSolver_Type,   // general monotonic cubic solver
+    };
     SkPoint fCoeff[3];
+    Type    fType;
 };
 
 #endif

tests/CubicMapTest.cpp

@@ -0,0 +1,42 @@
+/*
+ * Copyright 2018 Google Inc.
+ *
+ * Use of this source code is governed by a BSD-style license that can be
+ * found in the LICENSE file.
+ */
+
+#include "SkCubicMap.h"
+#include "SkRandom.h"
+#include "Test.h"
+
+static bool nearly_le(SkScalar a, SkScalar b) {
+    return a <= b || SkScalarNearlyZero(a - b);
+}
+
+static void exercise_cubicmap(const SkCubicMap& cmap, skiatest::Reporter* reporter) {
+    SkScalar prev_y = 0;
+    for (SkScalar x = 0; x <= 1; x += 1.0f / 512) {
+        SkScalar y = cmap.computeYFromX(x);
+        if (y < 0 || y > 1 || !nearly_le(prev_y, y)) {
+            cmap.computeYFromX(x);
+            REPORTER_ASSERT(reporter, false);
+        }
+        prev_y = y;
+    }
+}
+
+DEF_TEST(CubicMap, r) {
+    const SkScalar values[] = {
+        0, 1, 0.25f, 0.75f, 0.5f, 1.0f/3, 2.0f/3, 0.0000001f, 0.999999f,
+    };
+
+    for (SkScalar x0 : values) {
+        for (SkScalar y0 : values) {
+            for (SkScalar x1 : values) {
+                for (SkScalar y1 : values) {
+                    exercise_cubicmap(SkCubicMap({ x0, y0 }, { x1, y1 }), r);
+                }
+            }
+        }
+    }
+}