Implement Conics in SkCurveMeasure

src/utils/SkCurveMeasure.cpp

 /*
  * 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"
 #include "SkGeometry.h"
 
 // for abs
 #include <cmath>
 
-#define UNIMPLEMENTED SkDEBUGF(("%s:%d unimplemented\n", __FILE__, __LINE__))
+#define UNIMPLEMENTED SkDEBUGFAILF("%s:%d unimplemented\n", __FILE__, __LINE__)
 
 /// Used inside SkCurveMeasure::getTime's Newton's iteration
 static inline SkPoint evaluate(const SkPoint pts[4], SkSegType segType,
                                SkScalar t) {
     SkPoint pos;
     switch (segType) {
         case kQuad_SegType:
             pos = SkEvalQuadAt(pts, t);
             break;
         case kLine_SegType:
@@ skipping 33 matching lines @@
             SkConic conic(pts, pts[3].x());
             conic.evalAt(t, nullptr, &tan);
         }
             break;
         default:
             UNIMPLEMENTED;
     }
 
     return tan;
 }
+
+/*
+ * Derivation of conic formulae based on parameterization from
+ * http://citeseerx.ist.psu.edu/viewdoc/
+ *   download?doi=10.1.1.44.5740&rep=rep1&type=ps
+ * referenced in SkGeometry.cpp
+ *
+ * Usage: below and in the constructor for ArcLengthIntegrator
+ *
+ * Mathematica:
+ *
+ * conic[t_] := (A*(1 - t)^2 + 2*w*B*(1 - t)*t + C*t^2)/
+ *      (1 - 2*(1 - w)*t + 2*(1 - w)*t^2)
+ *
+ * conic[0]
+ *
+ * >>> A
+ *
+ * conic[1]
+ *
+ * >>> C
+ *
+ * Simplify[Derivative[1][conic][t]]
+ *
+ * >>> (2*(C*t*(1 + t*(-1 + w)) - A*(-1 + t)*(t*(-1 + w) - w) + B*(1 - 2*t)*w))/
+ *     (1 + 2*t*(-1 + w) - 2*t^2*(-1 + w))^2
+ *
+ * HornerForm[Numerator[Simplify[Derivative[1][conic][t]]], t]
+ *
+ * >>> -2*A*w + 2*B*w + t*(-2*A + 2*C + 4*A*w - 4*B*w +
+ *               t*(2*A - 2*C - 2*A*w + 2*C*w))
+ *
+ * HornerForm[Denominator[Simplify[Derivative[1][conic][t]]], t]
+ *
+ * >>> 1 + t*(-4 + 4*w + t*(8 - 12*w + 4*w^2 +
+ *               t*(-8 + 16*w - 8*w^2 + t*(4 - 8*w + 4*w^2))))
+ *
+ */
+
 /// Used in ArcLengthIntegrator::computeLength
 static inline Sk8f evaluateDerivativeLength(const Sk8f& ts,
-                                            const Sk8f (&xCoeff)[3],
-                                            const Sk8f (&yCoeff)[3],
+                                            const Sk8f (&fXCoeff)[3],
+                                            const Sk8f (&fYCoeff)[3],
+                                            const SkScalar w,
                                             const SkSegType segType) {
     Sk8f x;
     Sk8f y;
     switch (segType) {
         case kQuad_SegType:
-            x = xCoeff[0]*ts + xCoeff[1];
-            y = yCoeff[0]*ts + yCoeff[1];
+            x = fXCoeff[0]*ts + fXCoeff[1];
+            y = fYCoeff[0]*ts + fYCoeff[1];
             break;
         case kCubic_SegType:
-            x = (xCoeff[0]*ts + xCoeff[1])*ts + xCoeff[2];
-            y = (yCoeff[0]*ts + yCoeff[1])*ts + yCoeff[2];
+            x = (fXCoeff[0]*ts + fXCoeff[1])*ts + fXCoeff[2];
+            y = (fYCoeff[0]*ts + fYCoeff[1])*ts + fYCoeff[2];
             break;
-        case kConic_SegType:
-            UNIMPLEMENTED;
+        case kConic_SegType: {
+            Sk8f xnum = (fXCoeff[0]*ts + fXCoeff[1])*ts + fXCoeff[2];
+            Sk8f ynum = (fYCoeff[0]*ts + fYCoeff[1])*ts + fYCoeff[2];
+
+            SkScalar w484 = 4 - 8*w + 4*w*w;
+            Sk8f denom =  ts*(-4 + 4*w + ts*(8 - 12*w + 4*w*w + ts*(-2*w484
+                              + ts*(w484)))) + 1;

[caryclark, 2016/08/19 11:48:29]

More factoring is possible here, e.g. -(-4 + 4*w) + w484

[Harry Stern, 2016/08/19 16:31:28]

Oh clever, I didn't think of that. Done.

+
+            x = xnum / denom;
+            y = ynum / denom;
+        }
             break;
         default:
             UNIMPLEMENTED;
     }
 
     x = x * x;
     y = y * y;
 
     return (x + y).sqrt();
 }
 
 ArcLengthIntegrator::ArcLengthIntegrator(const SkPoint* pts, SkSegType segType)
     : fSegType(segType) {
     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
-            xCoeff[0] = Sk8f(2*(Ax - 2*Bx + Cx));
-            xCoeff[1] = Sk8f(2*(Bx - Ax));
+            fXCoeff[0] = Sk8f(2*(Ax - 2*Bx + Cx));
+            fXCoeff[1] = Sk8f(2*(Bx - Ax));
 
-            yCoeff[0] = Sk8f(2*(Ay - 2*By + Cy));
-            yCoeff[1] = Sk8f(2*(By - Ay));
+            fYCoeff[0] = Sk8f(2*(Ay - 2*By + Cy));
+            fYCoeff[1] = Sk8f(2*(By - Ay));
         }
             break;
         case kCubic_SegType:
         {
             float Ax = pts[0].x();
             float Bx = pts[1].x();
             float Cx = pts[2].x();
             float Dx = pts[3].x();
             float Ay = pts[0].y();
             float By = pts[1].y();
             float Cy = pts[2].y();
             float Dy = pts[3].y();
 
             // precompute coefficients for derivative
-            xCoeff[0] = Sk8f(3*(-Ax + 3*(Bx - Cx) + Dx));
-            xCoeff[1] = Sk8f(6*(Ax - 2*Bx + Cx));
-            xCoeff[2] = Sk8f(3*(-Ax + Bx));
+            fXCoeff[0] = Sk8f(3*(-Ax + 3*(Bx - Cx) + Dx));
+            fXCoeff[1] = Sk8f(6*(Ax - 2*Bx + Cx));
+            fXCoeff[2] = Sk8f(3*(-Ax + Bx));
 
-            yCoeff[0] = Sk8f(3*(-Ay + 3*(By - Cy) + Dy));
-            yCoeff[1] = Sk8f(6*(Ay - 2*By + Cy));
-            yCoeff[2] = Sk8f(3*(-Ay + By));
+            fYCoeff[0] = Sk8f(3*(-Ay + 3*(By - Cy) + Dy));
+            fYCoeff[1] = Sk8f(6*(Ay - 2*By + Cy));
+            fYCoeff[2] = Sk8f(3*(-Ay + By));
         }
             break;
-        case kConic_SegType:
-            UNIMPLEMENTED;
+        case kConic_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();
+
+            fConicW = pts[3].x();
+
+            SkScalar w = fConicW;
+
+            // precompute coefficients for derivative
+            fXCoeff[0] = Sk8f(2*(Ax - Cx - Ax*w + Cx*w));
+            fXCoeff[1] = Sk8f(-2*(Ax - Cx - 2*(Ax*w - Bx*w)));
+            fXCoeff[2] = Sk8f(-2*(Ax*w - Bx*w));
+
+            fYCoeff[0] = Sk8f(2*(Ay - Cy - Ay*w + Cy*w));
+            fYCoeff[1] = Sk8f(-2*(Ay - Cy - 2*(Ay*w - By*w)));
+            fYCoeff[2] = Sk8f(-2*(Ay*w - By*w));

[caryclark, 2016/08/19 11:48:29]

Can come elements, e.g. -2*(Ax*w - Bx*w) be factored out here?

[Harry Stern, 2016/08/19 16:31:28]

Sure. The only reason I didn't do it initially is because this is called at
construction, so it only happens once per segment of a path.

Done.

+        }
             break;
         default:
             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, xCoeff, yCoeff, fSegType);
-    lengths = weights*lengths;
+    Sk8f lengths =
+        evaluateDerivativeLength(absc * t, fXCoeff, fYCoeff, fConicW, fSegType);
+    lengths = weights * lengths;
     // is it faster or more accurate to sum and then multiply or vice versa?
     lengths = lengths*(t*0.5f);

[caryclark, 2016/08/19 11:48:30]

Thanks for adding spaces around the 'weight * lengths' above. Can you add spaces
everywhere on both sides of multiplies?

[Harry Stern, 2016/08/19 16:31:28]

In a previous review mtklein basically said not to do so everywhere, though I
thought it was fine for the above two.

He gave the reasoning that multiplication binds tighter than addition,
essentially.

I really don't have a preference either way.

 
     // Why does SkNx index with ints? does negative index mean something?
     for (int i = 0; i < 8; i++) {
         length += lengths[i];
     }
     return length;
 }
 
 SkCurveMeasure::SkCurveMeasure(const SkPoint* pts, SkSegType segType)
     : fSegType(segType) {
@@ skipping 43 matching lines @@
 // 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 (targetLength > currentLength || (SkScalarNearlyEqual(targetLength, currentLength))) {
+    if (targetLength > currentLength ||
+        (SkScalarNearlyEqual(targetLength, currentLength))) {

[caryclark, 2016/08/19 11:48:29]

parens unneeded

targetLength - some_epsilon >= currentLength

[Harry Stern, 2016/08/19 16:31:28]

Per discussion we're going to leave it like it is.

         return 1.0f;
     }
     if (kLine_SegType == fSegType) {
         return targetLength / currentLength;
     }
 
     // initial estimate of t is percentage of total length
     SkScalar currentT = targetLength / currentLength;
     SkScalar prevT = -1.0f;
     SkScalar newT;
@@ skipping 72 matching lines @@
     if (time) {
         *time = t;
     }
     if (pos) {
         *pos = evaluate(fPts, fSegType, t);
     }
     if (tan) {
         *tan = evaluateDerivative(fPts, fSegType, t);
     }
 }