Change cubicmap to eval directly (no table)

- gives us zero setup cost
- perfect accuracy (at least better than our 16-polyline approx)
- slower at runtime, but for skottie seems still way under the radar

Bug: skia:
Change-Id: Ic01f76184c359c65d6af9a3bedeac34a2ad7d3a6
Reviewed-on: https://skia-review.googlesource.com/146961
Reviewed-by: Mike Klein <mtklein@google.com>
Reviewed-by: Florin Malita <fmalita@chromium.org>
Commit-Queue: Mike Reed <reed@google.com>

gm/aaclip.cpp

@@ -11,52 +11,6 @@
 #include "SkPath.h"
 #include "SkMakeUnique.h"
 
-
-#include "SkCubicMap.h"
-
-static void test_cubic(SkCanvas* canvas) {
-    const SkPoint pts[] = {
-        { 0.333333f, 0.333333f }, { 0.666666f, 0.666666f },
-        { 1, 0 }, { 0, 1 },
-        { 0, 1 }, { 1, 0 },
-        { 0, 0 }, { 1, 1 },
-        { 1, 1 }, { 0, 0 },
-        { 0, 1 }, { 0, 1 },
-        { 1, 0 }, { 1, 0 },
-    };
-
-    SkPaint paint0, paint1;
-    paint0.setAntiAlias(true);  paint0.setStrokeWidth(3/256.0f); paint0.setColor(SK_ColorRED);
-    paint1.setAntiAlias(true);
-
-    SkCubicMap cmap;
-
-    canvas->translate(10, 266);
-    canvas->scale(256, -256);
-    for (size_t i = 0; i < SK_ARRAY_COUNT(pts); i += 2) {
-        cmap.setPts(pts[i], pts[i+1]);
-
-        const int N = 128;
-        SkPoint tmp0[N+1], tmp1[N+1], tmp2[N+1];
-        for (int j = 0; j <= N; ++j) {
-            float p = j * 1.0f / N;
-            tmp0[j] = cmap.computeFromT(p);
-            tmp1[j].set(p, cmap.computeYFromX(p));
-            tmp2[j].set(p, cmap.hackYFromX(p));
-        }
-
-        canvas->save();
-        canvas->drawPoints(SkCanvas::kPolygon_PointMode, N+1, tmp0, paint0);
-        canvas->drawPoints(SkCanvas::kPolygon_PointMode, N+1, tmp1, paint1);
-        canvas->translate(0, -1.2f);
-        canvas->drawPoints(SkCanvas::kPolygon_PointMode, N+1, tmp0, paint0);
-        canvas->drawPoints(SkCanvas::kPolygon_PointMode, N+1, tmp2, paint1);
-        canvas->restore();
-
-        canvas->translate(1.1f, 0);
-    }
-}
-
 static void do_draw(SkCanvas* canvas, const SkRect& r) {
     SkPaint paint;
     paint.setBlendMode(SkBlendMode::kSrc);
@@ -187,8 +141,6 @@ class AAClipGM : public skiagm::GM {
     }
 
     void onDraw(SkCanvas* canvas) override {
-        if (0) { test_cubic(canvas); return; }
-
         // Initial pixel-boundary-aligned draw
         draw_rect_tests(canvas);
 

gm/stringart.cpp

@@ -83,3 +83,75 @@ class StringArtGM : public skiagm::GM {
 };
 
 DEF_GM( return new StringArtGM; )
+
+/////////////////////////////////////////////////////////////////////////////////////////////////
+
+#if 0
+#include "Skottie.h"
+
+class SkottieGM : public skiagm::GM {
+    enum {
+        kWidth = 800,
+        kHeight = 600,
+    };
+
+    enum {
+        N = 100,
+    };
+    skottie::Animation* fAnims[N];
+    SkRect              fRects[N];
+    SkScalar            fDur;
+
+public:
+    SkottieGM() {
+        sk_bzero(fAnims, sizeof(fAnims));
+    }
+    ~SkottieGM() override {
+        for (auto anim : fAnims) {
+            SkSafeUnref(anim);
+        }
+    }
+
+protected:
+
+    SkString onShortName() override { return SkString("skottie"); }
+
+    SkISize onISize() override { return SkISize::Make(kWidth, kHeight); }
+
+    void init() {
+        SkRandom rand;
+        auto data = SkData::MakeFromFileName("/Users/reed/Downloads/maps_pinlet.json");
+   //     for (;;) skottie::Animation::Make((const char*)data->data(), data->size());
+        for (int i = 0; i < N; ++i) {
+            fAnims[i] = skottie::Animation::Make((const char*)data->data(), data->size()).release();
+            SkScalar x = rand.nextF() * kWidth;
+            SkScalar y = rand.nextF() * kHeight;
+            fRects[i].setXYWH(x, y, 400, 400);
+        }
+        fDur = fAnims[0]->duration();
+    }
+
+    void onDraw(SkCanvas* canvas) override {
+        if (!fAnims[0]) {
+            this->init();
+        }
+        canvas->drawColor(0xFFBBBBBB);
+        for (int i = 0; i < N; ++i) {
+            fAnims[0]->render(canvas, &fRects[i]);
+        }
+    }
+
+    bool onAnimate(const SkAnimTimer& timer) override {
+        SkScalar time = (float)(fmod(timer.secs(), fDur) / fDur);
+        for (auto anim : fAnims) {
+            anim->seek(time);
+        }
+        return true;
+    }
+
+private:
+    typedef GM INHERITED;
+};
+DEF_GM( return new SkottieGM; )
+#endif
+

src/core/SkCubicMap.cpp

@@ -9,75 +9,111 @@
 #include "SkNx.h"
 #include "../../src/pathops/SkPathOpsCubic.h"
 
-void SkCubicMap::setPts(SkPoint p1, SkPoint p2) {
-    Sk2s s1 = Sk2s::Load(&p1) * 3;
-    Sk2s s2 = Sk2s::Load(&p2) * 3;
-
-    s1 = Sk2s::Min(Sk2s::Max(s1, 0), 3);
-    s2 = Sk2s::Min(Sk2s::Max(s2, 0), 3);
-
-    (Sk2s(1) + s1 - s2).store(&fCoeff[0]);
-    (s2 - s1 - s1).store(&fCoeff[1]);
-    s1.store(&fCoeff[2]);
+static float eval_poly3(float a, float b, float c, float d, float t) {
+    return ((a * t + b) * t + c) * t + d;
+}
 
-    this->buildXTable();
+static float eval_poly2(float a, float b, float c, float t) {
+    return (a * t + b) * t + c;
 }
 
-SkPoint SkCubicMap::computeFromT(float t) const {
-    Sk2s a = Sk2s::Load(&fCoeff[0]);
-    Sk2s b = Sk2s::Load(&fCoeff[1]);
-    Sk2s c = Sk2s::Load(&fCoeff[2]);
+static float eval_poly1(float a, float b, float t) {
+    return a * t + b;
+}
 
-    SkPoint result;
-    (((a * t + b) * t + c) * t).store(&result);
-    return result;
+static float guess_nice_cubic_root(float A, float B, float C, float D) {
+    return -D;
 }
 
-float SkCubicMap::computeYFromX(float x) const {
-    x = SkTPin<float>(x, 0, 0.99999f) * kTableCount;
-    float ix = sk_float_floor(x);
-    int index = (int)ix;
-    SkASSERT((unsigned)index < SK_ARRAY_COUNT(fXTable));
-    return this->computeFromT(fXTable[index].fT0 + fXTable[index].fDT * (x - ix)).fY;
+#ifdef SK_DEBUG
+static bool valid(float r) {
+    return r >= 0 && r <= 1;
 }
+#endif
 
-float SkCubicMap::hackYFromX(float x) const {
-    x = SkTPin<float>(x, 0, 0.99999f) * kTableCount;
-    float ix = sk_float_floor(x);
-    int index = (int)ix;
-    SkASSERT((unsigned)index < SK_ARRAY_COUNT(fXTable));
-    return fXTable[index].fY0 + fXTable[index].fDY * (x - ix);
+/*
+ *  TODO: will this be faster if we algebraically compute the polynomials for the numer and denom
+ *        rather than compute them in parts?
+ *
+ *  TODO: investigate Householder's method, to see if we can get away with even fewer
+ *        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)
+    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) {
+        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;
+        float denom = 2 * fp * fp - f * fpp;
+        float delta = numer / denom;
+        float new_t = t - delta;
+        SkASSERT(valid(new_t));
+        t = new_t;
+    }
+    SkASSERT(valid(t));
+    return t;
 }
 
-static float compute_t_from_x(float A, float B, float C, float x) {
+#ifdef SK_DEBUG
+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);
     SkASSERT(count == 1);
     return (float)roots[0];
 }
 
-void SkCubicMap::buildXTable() {
-    float prevT = 0;
+static float max_err;
+#endif
 
-    const float dx = 1.0f / kTableCount;
-    float x = dx;
+static float compute_t_from_x(float A, float B, float C, float x) {
+#ifdef SK_DEBUG
+    float answer = compute_slow(A, B, C, x);
+#endif
+    float answer2 = solve_nice_cubic_halley(A, B, C, -x);
+#ifdef SK_DEBUG
+    float err = sk_float_abs(answer - answer2);
+    if (err > max_err) {
+        max_err = err;
+        SkDebugf("max error %g\n", max_err);
+    }
+#endif
+    return answer2;
+}
 
-    fXTable[0].fT0 = 0;
-    fXTable[0].fY0 = 0;
-    for (int i = 1; i < kTableCount; ++i) {
-        float t = compute_t_from_x(fCoeff[0].fX, fCoeff[1].fX, fCoeff[2].fX, x);
-        SkASSERT(t > prevT);
+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);
+    float a = fCoeff[0].fY;
+    float b = fCoeff[1].fY;
+    float c = fCoeff[2].fY;
+    return ((a * t + b) * t + c) * t;
+}
 
-        fXTable[i - 1].fDT = t - prevT;
-        fXTable[i].fT0 = t;
+void SkCubicMap::setPts(SkPoint p1, SkPoint p2) {
+    Sk2s s1 = Sk2s::Load(&p1) * 3;
+    Sk2s s2 = Sk2s::Load(&p2) * 3;
 
-        SkPoint p = this->computeFromT(t);
-        fXTable[i - 1].fDY = p.fY - fXTable[i - 1].fY0;
-        fXTable[i].fY0 = p.fY;
+    s1 = Sk2s::Min(Sk2s::Max(s1, 0), 3);
+    s2 = Sk2s::Min(Sk2s::Max(s2, 0), 3);
 
-        prevT = t;
-        x += dx;
-    }
-    fXTable[kTableCount - 1].fDT = 1 - prevT;
-    fXTable[kTableCount - 1].fDY = 1 - fXTable[kTableCount - 1].fY0;
+    (Sk2s(1) + s1 - s2).store(&fCoeff[0]);
+    (s2 - s1 - s1).store(&fCoeff[1]);
+    s1.store(&fCoeff[2]);
+}
+
+SkPoint SkCubicMap::computeFromT(float t) const {
+    Sk2s a = Sk2s::Load(&fCoeff[0]);
+    Sk2s b = Sk2s::Load(&fCoeff[1]);
+    Sk2s c = Sk2s::Load(&fCoeff[2]);
+
+    SkPoint result;
+    (((a * t + b) * t + c) * t).store(&result);
+    return result;
 }

src/core/SkCubicMap.h

@@ -10,33 +10,31 @@
 
 #include "SkPoint.h"
 
+/**
+ *  Fast evaluation of a cubic ease-in / ease-out curve. This is defined as a parametric cubic
+ *  curve inside the unit square.
+ *
+ *  pt[0] is implicitly { 0, 0 }
+ *  pt[3] is implicitly { 1, 1 }
+ *  pts[1,2] are inside the unit square
+ */
 class SkCubicMap {
 public:
-    void setPts(SkPoint p1, SkPoint p2);
-    void setPts(float x1, float y1, float x2, float y2) {
-        this->setPts({x1, y1}, {x2, y2});
+    SkCubicMap() {} // must call setPts() before using
+
+    SkCubicMap(SkPoint p1, SkPoint p2) {
+        this->setPts(p1, p2);
     }
+    void setPts(SkPoint p1, SkPoint p2);
 
-    SkPoint computeFromT(float t) const;
     float computeYFromX(float x) const;
 
-    // experimental
-    float hackYFromX(float x) const;
+    // is this needed?
+    SkPoint computeFromT(float t) const;
 
 private:
-    SkPoint fCoeff[4];
-    // x->t lookup
-    enum { kTableCount = 16 };
-    struct Rec {
-        float   fT0;
-        float   fDT;
-
-        float fY0;
-        float fDY;
-    };
-    Rec fXTable[kTableCount];
-
-    void buildXTable();
+    SkPoint fCoeff[3];
 };
+
 #endif