Add initial CurveMeasure code

src/utils/SkCurveMeasure.cpp

Index: src/utils/SkCurveMeasure.cpp
diff --git a/src/utils/SkCurveMeasure.cpp b/src/utils/SkCurveMeasure.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a7be484a97ccc6a5b37527e3200be3fba112a0b4
--- /dev/null
+++ b/src/utils/SkCurveMeasure.cpp
@@ -0,0 +1,250 @@
+/*
+ * Copyright 2016 Google Inc.
+ *
+ * Use of this source code is governed by a BSD-style license that can be
+ * found in the LICENSE file.
+ */
+
+#include "SkCurveMeasure.h"
+
+// for abs
+#include <cmath>
+
+ArcLengthIntegrator::ArcLengthIntegrator(const SkPoint* pts, SkSegType segType) : fSegType(segType) {

[reed1, 2016/08/02 21:46:23]

exceed 100cols

[Harry Stern, 2016/08/02 22:05:07]

Done, and various other formatting stuff.

+    switch (fSegType) {
+        case kQuad_SegType: {
+            float Ax = pts[0].x();
+            float Bx = pts[1].x();
+            float Cx = pts[2].x();
+            float Ay = pts[0].y();
+            float By = pts[1].y();
+            float Cy = pts[2].y();
+
+            // precompute coefficients for derivative
+            A2BC_x = Sk8f(2.0f*(Ax - 2*Bx + Cx));
+            A2BC_y = Sk8f(2.0f*(Ay - 2*By + Cy));
+
+            AB_x = Sk8f(2.0f*(Bx - Ax));
+            AB_y = Sk8f(2.0f*(By - Ay));
+        }
+            break;
+        case kLine_SegType:
+            SkDebugf("Unimplemented");
+            break;
+        case kCubic_SegType:
+            SkDebugf("Unimplemented");
+            break;
+        case kConic_SegType:
+            SkDebugf("Unimplemented");
+            break;
+        default:
+            SkDebugf("Unimplemented");
+    }
+}
+
+// We use Gaussian quadrature (https://en.wikipedia.org/wiki/Gaussian_quadrature)
+// to approximate the arc length integral here, because it is amenable to SIMD.
+SkScalar ArcLengthIntegrator::computeLength(SkScalar t) {
+    SkScalar length = 0.0f;
+
+    Sk8f lengths = evaluateDerivativeLength(absc*t);
+    lengths = weights*lengths;
+    // is it faster or more accurate to sum and then multiply or vice versa?
+    lengths = lengths*(t*0.5f);
+
+    // Why does SkNx index with ints? does negative index mean something?
+    for (int i = 0; i < 8; i++) {
+        length += lengths[i];
+    }
+    return length;
+}
+
+Sk8f ArcLengthIntegrator::evaluateDerivativeLength(Sk8f ts) {

[Harry Stern, 2016/08/02 21:33:21]

Windows buildbots are giving me
"error C2719: 'ts': formal parameter with requested alignment of 16 won't be
aligned" on this line and in the header. What does this mean?

[reed1, 2016/08/02 21:46:23]

Can you pass ts by const& ?

[Harry Stern, 2016/08/02 22:05:08]

Mike Klein told me I might be able to use SK_VECTORCALL, and according to
https://msdn.microsoft.com/en-us/library/dn375768.aspx that might additionally
make it faster. Apparently __vectorcall makes it pass via registers. We'll see
if that makes it stop complaining.

[mtklein, 2016/08/02 22:32:35]

Actually I think I had a rather more nuanced response and I would like to
reiterate that making functions that work with SkNx types static and inline is
almost always the best option.  If we _can_ make this a static function, it'd
just take your argument by const&.  Is this virtual or anything?  Does it use
ArcLengthIntegrator's state?  Can you just pass a pointer to 8 floats instead of
this?

+    Sk8f x;
+    Sk8f y;
+    switch (fSegType) {
+        case kQuad_SegType:
+            x = A2BC_x*ts + AB_x;
+            y = A2BC_y*ts + AB_y;
+            break;
+        case kLine_SegType:
+            SkDebugf("Unimplemented");
+            break;
+        case kCubic_SegType:
+            SkDebugf("Unimplemented");
+            break;
+        case kConic_SegType:
+            SkDebugf("Unimplemented");
+            break;
+        default:
+            SkDebugf("Unimplemented");
+    }
+
+    x = x * x;
+    y = y * y;
+
+    return (x + y).sqrt();
+}
+
+
+SkCurveMeasure::SkCurveMeasure(const SkPoint* pts, SkSegType segType) : fSegType(segType) {
+	switch (fSegType) {
+        case SkSegType::kQuad_SegType:
+            for (size_t i = 0; i < 3; i++) {
+                fPts[i] = pts[i];
+            }
+            fIntegrator = ArcLengthIntegrator(fPts, fSegType);
+            break;
+        default:
+            SkDebugf("Unimplemented");
+            break;
+    }
+}
+
+SkScalar SkCurveMeasure::getLength() {
+    if (-1.0f == fLength) {
+        fLength = fIntegrator.computeLength(1.0f);
+    }
+    return fLength;
+}
+
+// Given an arc length targetLength, we want to determine what t
+// gives us the corresponding arc length along the curve.
+// We do this by letting the arc length integral := f(t) and
+// solving for the root of the equation f(t) - targetLength = 0
+// using Newton's method and lerp-bisection.
+// The computationally expensive parts are the integral approximation
+// at each step, and computing the derivative of the arc length integral,
+// which is equal to the length of the tangent (so we have to do a sqrt).
+
+SkScalar SkCurveMeasure::getTime(SkScalar targetLength) {
+    if (targetLength == 0.0f) {
+        return 0.0f;
+    }
+
+    SkScalar currentLength = getLength();
+
+    if (SkScalarNearlyEqual(targetLength, currentLength)) {
+        return 1.0f;
+    }
+
+    // initial estimate of t is percentage of total length
+    SkScalar currentT = targetLength / currentLength;
+    SkScalar prevT = -1.0f;
+    SkScalar newT;
+
+    SkScalar minT = 0.0f;
+    SkScalar maxT = 1.0f;
+
+    int iterations = 0;
+    while (iterations < kNewtonIters + kBisectIters) {
+        currentLength = fIntegrator.computeLength(currentT);
+        SkScalar lengthDiff = currentLength - targetLength;
+
+        // Update root bounds.
+        // If lengthDiff is positive, we have overshot the target, so
+        // we know the current t is an upper bound, and similarly
+        // for the lower bound.
+        if (lengthDiff > 0.0f) {
+            if (currentT < maxT) {
+                maxT = currentT;
+            }
+        } else {
+            if (currentT > minT) {
+                minT = currentT;
+            }
+        }
+
+        // We have a tolerance on both the absolute value of the difference and
+        // on the t value
+        // because we may not have enough precision in the t to get close enough
+        // in the length.
+        if ((std::abs(lengthDiff) < kTolerance) ||
+            (std::abs(prevT - currentT) < kTolerance)) {
+            break;
+        }
+
+        prevT = currentT;
+        if (iterations < kNewtonIters) {
+            // TODO(hstern) switch here on curve type.
+            // This is just newton's formula.
+            SkScalar dt = evaluateQuadDerivative(currentT).length();
+            newT = currentT - (lengthDiff / dt);
+
+            // If newT is out of bounds, bisect inside newton.
+            if ((newT < 0.0f) || (newT > 1.0f)) {
+                newT = (minT + maxT) * 0.5f;
+            }
+        } else if (iterations < kNewtonIters + kBisectIters) {
+            if (lengthDiff > 0.0f) {
+                maxT = currentT;
+            } else {
+                minT = currentT;
+            }
+            // TODO(hstern) do a lerp here instead of a bisection
+            newT = (minT + maxT) * 0.5f;
+        } else {
+            SkDebugf("%.7f %.7f didn't get close enough after bisection.\n",
+                    currentT, currentLength);
+            break;
+        }
+        currentT = newT;
+
+        SkASSERT(minT <= maxT);
+
+        iterations++;
+    }
+
+    // debug. is there an SKDEBUG or something for ifdefs?
+    fIters = iterations;
+
+    return currentT;
+}
+
+void SkCurveMeasure::getPosTan(SkScalar targetLength, SkPoint* pos,
+                               SkVector* tan) {
+    SkScalar t = getTime(targetLength);
+
+    if (pos) {
+        // TODO(hstern) switch here on curve type.
+        *pos = evaluateQuad(t);
+    }
+    if (tan) {
+        // TODO(hstern) switch here on curve type.
+        *tan = evaluateQuadDerivative(t);
+    }
+}
+
+// this is why I feel that the ArcLengthIntegrator should be combined
+// with some sort of evaluator that caches the constants computed from the
+// control points. this is basically the same code in ArcLengthIntegrator
+SkPoint SkCurveMeasure::evaluateQuad(SkScalar t) {
+    SkScalar ti = 1.0f - t;
+
+    SkScalar Ax = fPts[0].x();
+    SkScalar Bx = fPts[1].x();
+    SkScalar Cx = fPts[2].x();
+    SkScalar Ay = fPts[0].y();
+    SkScalar By = fPts[1].y();
+    SkScalar Cy = fPts[2].y();
+
+    SkScalar x = Ax*ti*ti + 2.0f*Bx*t*ti + Cx*t*t;
+    SkScalar y = Ay*ti*ti + 2.0f*By*t*ti + Cy*t*t;
+    return SkPoint::Make(x, y);
+}
+
+SkVector SkCurveMeasure::evaluateQuadDerivative(SkScalar t) {
+    SkScalar Ax = fPts[0].x();
+    SkScalar Bx = fPts[1].x();
+    SkScalar Cx = fPts[2].x();
+    SkScalar Ay = fPts[0].y();
+    SkScalar By = fPts[1].y();
+    SkScalar Cy = fPts[2].y();
+
+    SkScalar A2BCx = 2.0f*(Ax - 2*Bx + Cx);
+    SkScalar A2BCy = 2.0f*(Ay - 2*By + Cy);
+    SkScalar ABx = 2.0f*(Bx - Ax);
+    SkScalar ABy = 2.0f*(By - Ay);
+
+    return SkPoint::Make(A2BCx*t + ABx, A2BCy*t + ABy);
+}