use Sk2s for conics

src/core/SkGeometry.cpp

 /*
  * Copyright 2006 The Android Open Source Project
  *
  * Use of this source code is governed by a BSD-style license that can be
  * found in the LICENSE file.
  */
 
 #include "SkGeometry.h"
 #include "SkMatrix.h"
 #include "Sk2x.h"
 
+static Sk2s from_point(const SkPoint& point) {
+    return Sk2s::Load(&point.fX);
+}
+
+static SkPoint to_point(const Sk2s& x) {
+    SkPoint point;
+    x.store(&point.fX);
+    return point;
+}
+
+static SkVector to_vector(const Sk2s& x) {
+    SkVector vector;
+    x.store(&vector.fX);
+    return vector;
+}
+
+#if 0
+static Sk2s divide(const Sk2s& numer, const Sk2s& denom) {
+    SkScalar numerStorage[2], denomStorage[2];
+    numer.store(numerStorage);
+    denom.store(denomStorage);
+    numerStorage[0] /= denomStorage[0];
+    numerStorage[1] /= denomStorage[1];
+    return Sk2s::Load(numerStorage);
+}
+#endif
+
 /** If defined, this makes eval_quad and eval_cubic do more setup (sometimes
     involving integer multiplies by 2 or 3, but fewer calls to SkScalarMul.
     May also introduce overflow of fixed when we compute our setup.
 */
 //    #define DIRECT_EVAL_OF_POLYNOMIALS
 
 ////////////////////////////////////////////////////////////////////////
 
 static int is_not_monotonic(SkScalar a, SkScalar b, SkScalar c) {
     SkScalar ab = a - b;
@@ skipping 63 matching lines @@
             SkTSwap<SkScalar>(roots[0], roots[1]);
         else if (roots[0] == roots[1])  // nearly-equal?
             r -= 1; // skip the double root
     }
     return (int)(r - roots);
 }
 
 ///////////////////////////////////////////////////////////////////////////////
 ///////////////////////////////////////////////////////////////////////////////
 
+static Sk2s quad_poly_eval(const Sk2s& A, const Sk2s& B, const Sk2s& C, const Sk2s& t) {
+    return (A * t + B) * t + C;
+}
+
 static SkScalar eval_quad(const SkScalar src[], SkScalar t) {
     SkASSERT(src);
     SkASSERT(t >= 0 && t <= SK_Scalar1);
 
 #ifdef DIRECT_EVAL_OF_POLYNOMIALS
     SkScalar    C = src[0];
     SkScalar    A = src[4] - 2 * src[2] + C;
     SkScalar    B = 2 * (src[2] - C);
     return SkScalarMulAdd(SkScalarMulAdd(A, t, B), t, C);
 #else
@@ skipping 20 matching lines @@
     if (tangent) {
         tangent->set(eval_quad_derivative(&src[0].fX, t),
                      eval_quad_derivative(&src[0].fY, t));
     }
 }
 
 SkPoint SkEvalQuadAt(const SkPoint src[3], SkScalar t) {
     SkASSERT(src);
     SkASSERT(t >= 0 && t <= SK_Scalar1);
 
-    const Sk2f t2(t);
+    const Sk2s t2(t);
 
-    Sk2f P0 = Sk2f::Load(&src[0].fX);
-    Sk2f P1 = Sk2f::Load(&src[1].fX);
-    Sk2f P2 = Sk2f::Load(&src[2].fX);
+    Sk2s P0 = from_point(src[0]);
+    Sk2s P1 = from_point(src[1]);
+    Sk2s P2 = from_point(src[2]);
 
-    Sk2f B = P1 - P0;
-    Sk2f A = P2 - P1 - B;
+    Sk2s B = P1 - P0;
+    Sk2s A = P2 - P1 - B;
 
-    SkPoint result;
-    ((A * t2 + B+B) * t2 + P0).store(&result.fX);
-    return result;
+    return to_point((A * t2 + B+B) * t2 + P0);
 }
 
 SkVector SkEvalQuadTangentAt(const SkPoint src[3], SkScalar t) {
     SkASSERT(src);
     SkASSERT(t >= 0 && t <= SK_Scalar1);
 
-    Sk2f P0 = Sk2f::Load(&src[0].fX);
-    Sk2f P1 = Sk2f::Load(&src[1].fX);
-    Sk2f P2 = Sk2f::Load(&src[2].fX);
+    Sk2s P0 = from_point(src[0]);
+    Sk2s P1 = from_point(src[1]);
+    Sk2s P2 = from_point(src[2]);
 
-    Sk2f B = P1 - P0;
-    Sk2f A = P2 - P1 - B;
-    Sk2f T = A * Sk2f(t) + B;
+    Sk2s B = P1 - P0;
+    Sk2s A = P2 - P1 - B;
+    Sk2s T = A * Sk2s(t) + B;
 
-    SkVector result;
-    (T + T).store(&result.fX);
-    return result;
+    return to_vector(T + T);
 }
 
 static void interp_quad_coords(const SkScalar* src, SkScalar* dst, SkScalar t) {
     SkScalar    ab = SkScalarInterp(src[0], src[2], t);
     SkScalar    bc = SkScalarInterp(src[2], src[4], t);
 
     dst[0] = src[0];
     dst[2] = ab;
     dst[4] = SkScalarInterp(ab, bc, t);
     dst[6] = bc;
     dst[8] = src[4];
 }
 
 void SkChopQuadAt(const SkPoint src[3], SkPoint dst[5], SkScalar t) {
     SkASSERT(t > 0 && t < SK_Scalar1);
 
     interp_quad_coords(&src[0].fX, &dst[0].fX, t);
     interp_quad_coords(&src[0].fY, &dst[0].fY, t);
 }
 
 static inline Sk2s interp(const Sk2s& v0, const Sk2s& v1, const Sk2s& t) {
     return v0 + (v1 - v0) * t;
 }
 
 void SkChopQuadAt2(const SkPoint src[3], SkPoint dst[5], SkScalar t) {
     SkASSERT(t > 0 && t < SK_Scalar1);
 
-    Sk2s p0 = Sk2f::Load(&src[0].fX);
-    Sk2s p1 = Sk2f::Load(&src[1].fX);
-    Sk2s p2 = Sk2f::Load(&src[2].fX);
+    Sk2s p0 = from_point(src[0]);
+    Sk2s p1 = from_point(src[1]);
+    Sk2s p2 = from_point(src[2]);
     Sk2s tt = Sk2s(t);
     
     Sk2s p01 = interp(p0, p1, tt);
     Sk2s p12 = interp(p1, p2, tt);
 
-    p0.store(&dst[0].fX);
-    p01.store(&dst[1].fX);
-    interp(p01, p12, tt).store(&dst[2].fX);
-    p12.store(&dst[3].fX);
-    p2.store(&dst[4].fX);
+    dst[0] = to_point(p0);
+    dst[1] = to_point(p01);
+    dst[2] = to_point(interp(p01, p12, tt));
+    dst[3] = to_point(p12);
+    dst[4] = to_point(p2);
 }
 
 void SkChopQuadAtHalf(const SkPoint src[3], SkPoint dst[5]) {
     SkScalar x01 = SkScalarAve(src[0].fX, src[1].fX);
     SkScalar y01 = SkScalarAve(src[0].fY, src[1].fY);
     SkScalar x12 = SkScalarAve(src[1].fX, src[2].fX);
     SkScalar y12 = SkScalarAve(src[1].fY, src[2].fY);
 
     dst[0] = src[0];
     dst[1].set(x01, y01);
@@ skipping 1030 matching lines @@
     // w1 /= sqrt(w0*w2)
     //
     // However, in our case, we know that for dst[0]:
     //     w0 == 1, and for dst[1], w2 == 1
     //
     SkScalar root = SkScalarSqrt(tmp2[1].fZ);
     dst[0].fW = tmp2[0].fZ / root;
     dst[1].fW = tmp2[2].fZ / root;
 }
 
+static Sk2s times_2(const Sk2s& value) {
+    return value + value;
+}
+
+SkPoint SkConic::evalAt(SkScalar t) const {
+    Sk2s p0 = from_point(fPts[0]);
+    Sk2s p1 = from_point(fPts[1]);
+    Sk2s p2 = from_point(fPts[2]);
+    Sk2s tt = Sk2s(t);
+    Sk2s ww = Sk2s(fW);
+    Sk2s one = Sk2s(1);
+
+    Sk2s p1w = p1 * ww;
+    Sk2s C = p0;
+    Sk2s A = p2 - times_2(p1w) + p0;
+    Sk2s B = times_2(p1w - C);
+    Sk2s numer = quad_poly_eval(A, B, C, tt);
+
+    B = times_2(ww - one);
+    A = -B;
+    Sk2s denom = quad_poly_eval(A, B, one, tt);
+
+    return to_point(numer / denom);
+}
+
+SkVector SkConic::evalTangentAt(SkScalar t) const {
+    Sk2s p0 = from_point(fPts[0]);
+    Sk2s p1 = from_point(fPts[1]);
+    Sk2s p2 = from_point(fPts[2]);
+    Sk2s ww = Sk2s(fW);
+
+    Sk2s p20 = p2 - p0;
+    Sk2s p10 = p1 - p0;
+
+    Sk2s C = ww * p10;
+    Sk2s A = ww * p20 - p20;
+    Sk2s B = p20 - C - C;
+
+    return to_vector(quad_poly_eval(A, B, C, Sk2s(t)));
+#if 0
+    static void conic_deriv_coeff(const SkScalar src[],
+                                  SkScalar w,
+                                  SkScalar coeff[3]) {
+        const SkScalar P20 = src[4] - src[0];
+        const SkScalar P10 = src[2] - src[0];
+        const SkScalar wP10 = w * P10;
+        coeff[0] = w * P20 - P20;
+        coeff[1] = P20 - 2 * wP10;
+        coeff[2] = wP10;
+    }
+
+    static SkScalar conic_eval_tan(const SkScalar coord[], SkScalar w, SkScalar t) {
+        SkScalar coeff[3];
+        conic_deriv_coeff(coord, w, coeff);
+        return t * (t * coeff[0] + coeff[1]) + coeff[2];
+    }
+#endif
+}
+
 static SkScalar subdivide_w_value(SkScalar w) {
     return SkScalarSqrt(SK_ScalarHalf + w * SK_ScalarHalf);
 }
 
 void SkConic::chop(SkConic dst[2]) const {
     SkScalar scale = SkScalarInvert(SK_Scalar1 + fW);
     SkScalar p1x = fW * fPts[1].fX;
     SkScalar p1y = fW * fPts[1].fY;
     SkScalar mx = (fPts[0].fX + 2 * p1x + fPts[2].fX) * scale * SK_ScalarHalf;
     SkScalar my = (fPts[0].fY + 2 * p1y + fPts[2].fY) * scale * SK_ScalarHalf;
 
     dst[0].fPts[0] = fPts[0];
     dst[0].fPts[1].set((fPts[0].fX + p1x) * scale,
                        (fPts[0].fY + p1y) * scale);
     dst[0].fPts[2].set(mx, my);
 
     dst[1].fPts[0].set(mx, my);
     dst[1].fPts[1].set((p1x + fPts[2].fX) * scale,
                        (p1y + fPts[2].fY) * scale);
     dst[1].fPts[2] = fPts[2];
 
     dst[0].fW = dst[1].fW = subdivide_w_value(fW);
 }
 
+void SkConic::chop2(SkConic * SK_RESTRICT dst) const {
+    Sk2s scale(SkScalarInvert(SK_Scalar1 + fW));
+//    Sk2s scale = Sk2s(SK_Scalar1 + fW).invert();
+    SkScalar newW = subdivide_w_value(fW);
+
+    Sk2s p0 = from_point(fPts[0]);
+    Sk2s p1 = from_point(fPts[1]);
+    Sk2s p2 = from_point(fPts[2]);
+    Sk2s ww = Sk2s(fW);
+    Sk2s half = Sk2s(0.5f);
+
+    Sk2s wp1 = ww * p1;
+    Sk2s m = ((p0 + wp1 + wp1 + p2) * half) * scale;
+
+    dst[0].fPts[0] = fPts[0];
+    dst[0].fPts[1] = to_point((p0 + wp1) * scale);
+    dst[0].fPts[2] = dst[1].fPts[0] = to_point(m);
+    dst[1].fPts[1] = to_point((wp1 + p2) * scale);
+    dst[1].fPts[2] = fPts[2];
+
+    dst[0].fW = dst[1].fW = newW;
+}
+
 /*
  *  "High order approximation of conic sections by quadratic splines"
  *      by Michael Floater, 1993
  */
 #define AS_QUAD_ERROR_SETUP                                         \
     SkScalar a = fW - 1;                                            \
     SkScalar k = a / (4 * (2 + a));                                 \
     SkScalar x = k * (fPts[0].fX - 2 * fPts[1].fX + fPts[2].fX);    \
     SkScalar y = k * (fPts[0].fY - 2 * fPts[1].fY + fPts[2].fY);
 
@@ skipping 229 matching lines @@
         matrix.preScale(SK_Scalar1, -SK_Scalar1);
     }
     if (userMatrix) {
         matrix.postConcat(*userMatrix);
     }
     for (int i = 0; i < conicCount; ++i) {
         matrix.mapPoints(dst[i].fPts, 3);
     }
     return conicCount;
 }