Add GrWangsFormula implementation for conics

Also add a unit test that the vectorized version equals the reference
implementation.

Bug: skia:10419
Change-Id: I4d165fd45532e9ec468565d0637fb769b51f5fcd
Reviewed-on: https://skia-review.googlesource.com/c/skia/+/345122
Commit-Queue: Tyler Denniston <tdenniston@google.com>
Reviewed-by: Chris Dalton <csmartdalton@google.com>

bench/TessellateBench.cpp

@@ -50,6 +50,19 @@ static SkPath make_cubic_path() {
     return path;
 }
 
+static SkPath make_conic_path() {
+    SkRandom rand;
+    SkPath path;
+    for (int i = 0; i < kNumCubicsInChalkboard / 40; ++i) {
+        for (int j = -10; j <= 10; j++) {
+            const float x = std::ldexp(rand.nextF(), (i % 18)) / 1e3f;
+            const float w = std::ldexp(1 + rand.nextF(), j);
+            path.conicTo(111.625f * x, 308.188f * x, 764.62f * x, -435.688f * x, w);
+        }
+    }
+    return path;
+}
+
 // This serves as a base class for benchmarking individual methods on GrPathTessellateOp.
 class PathTessellateBenchmark : public Benchmark {
 public:
@@ -137,6 +150,46 @@ DEF_PATH_TESS_BENCH(wangs_formula_cubic_log2_affine, make_cubic_path(),
     benchmark_wangs_formula_cubic_log2(fMatrix, fPath);
 }
 
+static void benchmark_wangs_formula_conic(const SkMatrix& matrix, const SkPath& path) {
+    // Conic version expects tolerance, not intolerance
+    constexpr float kTolerance = 4;
+    int sum = 0;
+    GrVectorXform xform(matrix);
+    for (auto [verb, pts, w] : SkPathPriv::Iterate(path)) {
+        if (verb == SkPathVerb::kConic) {
+            sum += GrWangsFormula::conic(kTolerance, pts, *w, xform);
+        }
+    }
+    // Don't let the compiler optimize away GrWangsFormula::conic.
+    if (sum <= 0) {
+        SK_ABORT("sum should be > 0.");
+    }
+}
+
+static void benchmark_wangs_formula_conic_log2(const SkMatrix& matrix, const SkPath& path) {
+    // Conic version expects tolerance, not intolerance
+    constexpr float kTolerance = 4;
+    int sum = 0;
+    GrVectorXform xform(matrix);
+    for (auto [verb, pts, w] : SkPathPriv::Iterate(path)) {
+        if (verb == SkPathVerb::kConic) {
+            sum += GrWangsFormula::conic_log2(kTolerance, pts, *w, xform);
+        }
+    }
+    // Don't let the compiler optimize away GrWangsFormula::conic.
+    if (sum <= 0) {
+        SK_ABORT("sum should be > 0.");
+    }
+}
+
+DEF_PATH_TESS_BENCH(wangs_formula_conic, make_conic_path(), SkMatrix::I()) {
+    benchmark_wangs_formula_conic(fMatrix, fPath);
+}
+
+DEF_PATH_TESS_BENCH(wangs_formula_conic_log2, make_conic_path(), SkMatrix::I()) {
+    benchmark_wangs_formula_conic_log2(fMatrix, fPath);
+}
+
 DEF_PATH_TESS_BENCH(middle_out_triangulation,
                     ToolUtils::make_star(SkRect::MakeWH(500, 500), kNumCubicsInChalkboard),
                     SkMatrix::I()) {

src/gpu/tessellate/GrWangsFormula.h

@@ -37,6 +37,14 @@ SK_ALWAYS_INLINE static float root4(float x) {
     return sqrtf(sqrtf(x));
 }
 
+// Returns nextlog2(sqrt(x)):
+//
+//   log2(sqrt(x)) == log2(x^(1/2)) == log2(x)/2 == log2(x)/log2(4) == log4(x)
+//
+SK_ALWAYS_INLINE static int nextlog4(float x) {
+    return (sk_float_nextlog2(x) + 1) >> 1;
+}
+
 // Returns nextlog2(sqrt(sqrt(x))):
 //
 //   log2(sqrt(sqrt(x))) == log2(x^(1/4)) == log2(x)/4 == log2(x)/log2(16) == log16(x)
@@ -116,6 +124,61 @@ SK_ALWAYS_INLINE static int worst_case_cubic_log2(float intolerance, float devWi
     return nextlog16(4*kk * (devWidth * devWidth + devHeight * devHeight));
 }
 
+// Returns Wang's formula specialized for a conic curve, raised to the second power.
+// Input points should be in projected space, and note tolerance parameter is not intolerance.
+//
+// This is not actually due to Wang, but is an analogue from (Theorem 3, corollary 1):
+//   J. Zheng, T. Sederberg. "Estimating Tessellation Parameter Intervals for
+//   Rational Curves and Surfaces." ACM Transactions on Graphics 19(1). 2000.
+SK_ALWAYS_INLINE static float conic_pow2(float tolerance, const SkPoint pts[], float w,
+                                         const GrVectorXform& vectorXform = GrVectorXform()) {
+    using grvx::dot, grvx::float2, grvx::float4, skvx::bit_pun;
+    float2 p0 = vectorXform(bit_pun<float2>(pts[0]));
+    float2 p1 = vectorXform(bit_pun<float2>(pts[1]));
+    float2 p2 = vectorXform(bit_pun<float2>(pts[2]));
+
+    // Compute center of bounding box in projected space
+    const float2 C = 0.5f * (skvx::min(skvx::min(p0, p1), p2) + skvx::max(skvx::max(p0, p1), p2));
+
+    // Translate by -C. This improves translation-invariance of the formula,
+    // see Sec. 3.3 of cited paper
+    p0 -= C;
+    p1 -= C;
+    p2 -= C;
+
+    // Compute max length
+    const float max_len = sqrtf(std::max(dot(p0, p0), std::max(dot(p1, p1), dot(p2, p2))));
+
+    // Compute forward differences
+    const float2 dp = grvx::fast_madd<2>(-2 * w, p1, p0) + p2;
+    const float dw = fabsf(1 - 2 * w + 1);
+
+    // Compute numerator and denominator for parametric step size of linearization
+    const float r_minus_eps = std::max(0.f, max_len - tolerance);
+    const float min_w = std::min(w, 1.f);
+    const float numer = sqrtf(grvx::dot(dp, dp)) + r_minus_eps * dw;
+    const float denom = 4 * min_w * tolerance;
+
+    // Number of segments = sqrt(numer / denom).
+    // This assumes parametric interval of curve being linearized is [t0,t1] = [0, 1].
+    // If not, the number of segments is (tmax - tmin) / sqrt(denom / numer).
+    return numer / denom;
+}
+
+// Returns the value of Wang's formula specialized for a conic curve.
+SK_ALWAYS_INLINE static float conic(float tolerance, const SkPoint pts[], float w,
+                                    const GrVectorXform& vectorXform = GrVectorXform()) {
+    return sqrtf(conic_pow2(tolerance, pts, w, vectorXform));
+}
+
+// Returns the log2 value of Wang's formula specialized for a conic curve, rounded up to the next
+// int.
+SK_ALWAYS_INLINE static int conic_log2(float tolerance, const SkPoint pts[], float w,
+                                       const GrVectorXform& vectorXform = GrVectorXform()) {
+    // nextlog4(x) == ceil(log2(sqrt(x)))
+    return nextlog4(conic_pow2(tolerance, pts, w, vectorXform));
+}
+
 }  // namespace GrWangsFormula
 
 #endif

tests/WangsFormulaTest.cpp

@@ -446,3 +446,55 @@ DEF_TEST(WangsFormula_conic_within_tol, r) {
                 maxExponent);
     }
 }
+
+// Ensure the vectorized conic version equals the reference implementation
+DEF_TEST(WangsFormula_conic_matches_reference, r) {
+    constexpr static float kTolerance = 1.f / kIntolerance;
+
+    SkRandom rand;
+    for (int i = -10; i <= 10; ++i) {
+        const float w = std::ldexp(1 + rand.nextF(), i);
+        for_random_beziers(3, &rand, [&r, w](const SkPoint pts[]) {
+            const SkPoint projPts[3] = {pts[0], pts[1] * (1.f / w), pts[2]};
+            const float ref_nsegs = wangs_formula_conic_reference_impl(kIntolerance, projPts, w);
+            const float nsegs = GrWangsFormula::conic(kTolerance, projPts, w);
+
+            // Because the Gr version may implement the math differently for performance,
+            // allow different slack in the comparison based on the rough scale of the answer.
+            const float cmpThresh = ref_nsegs * (1.f / (1 << 20));
+            REPORTER_ASSERT(r, SkScalarNearlyEqual(ref_nsegs, nsegs, cmpThresh));
+        });
+    }
+}
+
+// Ensure using transformations gives the same result as pre-transforming all points.
+DEF_TEST(WangsFormula_conic_vectorXforms, r) {
+    constexpr static float kTolerance = 1.f / kIntolerance;
+
+    auto check_conic_with_transform = [&](const SkPoint* pts, float w, const SkMatrix& m) {
+        const SkPoint projPts[3] = {pts[0], pts[1] * (1.f / w), pts[2]};
+        SkPoint ptsXformed[3];
+        m.mapPoints(ptsXformed, projPts, 3);
+        float expected = GrWangsFormula::conic(kTolerance, ptsXformed, w);
+        float actual = GrWangsFormula::conic(kTolerance, projPts, w, GrVectorXform(m));
+        REPORTER_ASSERT(r, SkScalarNearlyEqual(actual, expected));
+    };
+
+    SkRandom rand;
+    for (int i = -10; i <= 10; ++i) {
+        const float w = std::ldexp(1 + rand.nextF(), i);
+        for_random_beziers(3, &rand, [&](const SkPoint pts[]) {
+            check_conic_with_transform(pts, w, SkMatrix::I());
+            check_conic_with_transform(
+                    pts, w, SkMatrix::Scale(rand.nextRangeF(-10, 10), rand.nextRangeF(-10, 10)));
+
+            // Random 2x2 matrix
+            SkMatrix m;
+            m.setScaleX(rand.nextRangeF(-10, 10));
+            m.setSkewX(rand.nextRangeF(-10, 10));
+            m.setSkewY(rand.nextRangeF(-10, 10));
+            m.setScaleY(rand.nextRangeF(-10, 10));
+            check_conic_with_transform(pts, w, m);
+        });
+    }
+}