xref: /aosp_15_r20/external/skia/tests/WangsFormulaTest.cpp (revision c8dee2aa9b3f27cf6c858bd81872bdeb2c07ed17)
1 /*
2  * Copyright 2020 Google Inc.
3  *
4  * Use of this source code is governed by a BSD-style license that can be
5  * found in the LICENSE file.
6  */
7 
8 #include "include/core/SkMatrix.h"
9 #include "include/core/SkPoint.h"
10 #include "include/core/SkRect.h"
11 #include "include/core/SkScalar.h"
12 #include "include/core/SkTypes.h"
13 #include "src/base/SkRandom.h"
14 #include "src/base/SkVx.h"
15 #include "src/core/SkGeometry.h"
16 #include "src/gpu/tessellate/Tessellation.h"
17 #include "src/gpu/tessellate/WangsFormula.h"
18 #include "tests/Test.h"
19 
20 #include <algorithm>
21 #include <cmath>
22 #include <cstring>
23 #include <functional>
24 #include <limits>
25 
26 namespace skgpu::tess {
27 
28 const SkPoint kSerp[4] = {
29         {285.625f, 499.687f}, {411.625f, 808.188f}, {1064.62f, 135.688f}, {1042.63f, 585.187f}};
30 
31 const SkPoint kLoop[4] = {
32         {635.625f, 614.687f}, {171.625f, 236.188f}, {1064.62f, 135.688f}, {516.625f, 570.187f}};
33 
34 const SkPoint kQuad[4] = {
35         {460.625f, 557.187f}, {707.121f, 209.688f}, {779.628f, 577.687f}};
36 
wangs_formula_quadratic_reference_impl(float precision,const SkPoint p[3])37 static float wangs_formula_quadratic_reference_impl(float precision, const SkPoint p[3]) {
38     float k = (2 * 1) / 8.f * precision;
39     return sqrtf(k * (p[0] - p[1]*2 + p[2]).length());
40 }
41 
wangs_formula_cubic_reference_impl(float precision,const SkPoint p[4])42 static float wangs_formula_cubic_reference_impl(float precision, const SkPoint p[4]) {
43     float k = (3 * 2) / 8.f * precision;
44     return sqrtf(k * std::max((p[0] - p[1]*2 + p[2]).length(),
45                               (p[1] - p[2]*2 + p[3]).length()));
46 }
47 
48 // Returns number of segments for linearized quadratic rational. This is an analogue
49 // to Wang's formula, taken from:
50 //
51 //   J. Zheng, T. Sederberg. "Estimating Tessellation Parameter Intervals for
52 //   Rational Curves and Surfaces." ACM Transactions on Graphics 19(1). 2000.
53 // See Thm 3, Corollary 1.
54 //
55 // Input points should be in projected space.
wangs_formula_conic_reference_impl(float precision,const SkPoint P[3],const float w)56 static float wangs_formula_conic_reference_impl(float precision,
57                                                 const SkPoint P[3],
58                                                 const float w) {
59     // Compute center of bounding box in projected space
60     float min_x = P[0].fX, max_x = min_x,
61           min_y = P[0].fY, max_y = min_y;
62     for (int i = 1; i < 3; i++) {
63         min_x = std::min(min_x, P[i].fX);
64         max_x = std::max(max_x, P[i].fX);
65         min_y = std::min(min_y, P[i].fY);
66         max_y = std::max(max_y, P[i].fY);
67     }
68     const SkPoint C = SkPoint::Make(0.5f * (min_x + max_x), 0.5f * (min_y + max_y));
69 
70     // Translate control points and compute max length
71     SkPoint tP[3] = {P[0] - C, P[1] - C, P[2] - C};
72     float max_len = 0;
73     for (int i = 0; i < 3; i++) {
74         max_len = std::max(max_len, tP[i].length());
75     }
76     SkASSERT(max_len > 0);
77 
78     // Compute delta = parametric step size of linearization
79     const float eps = 1 / precision;
80     const float r_minus_eps = std::max(0.f, max_len - eps);
81     const float min_w = std::min(w, 1.f);
82     const float numer = 4 * min_w * eps;
83     const float denom =
84             (tP[2] - tP[1] * 2 * w + tP[0]).length() + r_minus_eps * std::abs(1 - 2 * w + 1);
85     const float delta = sqrtf(numer / denom);
86 
87     // Return corresponding num segments in the interval [tmin,tmax]
88     constexpr float tmin = 0, tmax = 1;
89     SkASSERT(delta > 0);
90     return (tmax - tmin) / delta;
91 }
92 
for_random_matrices(SkRandom * rand,const std::function<void (const SkMatrix &)> & f)93 static void for_random_matrices(SkRandom* rand, const std::function<void(const SkMatrix&)>& f) {
94     SkMatrix m;
95     m.setIdentity();
96     f(m);
97 
98     for (int i = -10; i <= 30; ++i) {
99         for (int j = -10; j <= 30; ++j) {
100             m.setScaleX(std::ldexp(1 + rand->nextF(), i));
101             m.setSkewX(0);
102             m.setSkewY(0);
103             m.setScaleY(std::ldexp(1 + rand->nextF(), j));
104             f(m);
105 
106             m.setScaleX(std::ldexp(1 + rand->nextF(), i));
107             m.setSkewX(std::ldexp(1 + rand->nextF(), (j + i) / 2));
108             m.setSkewY(std::ldexp(1 + rand->nextF(), (j + i) / 2));
109             m.setScaleY(std::ldexp(1 + rand->nextF(), j));
110             f(m);
111         }
112     }
113 }
114 
for_random_beziers(int numPoints,SkRandom * rand,const std::function<void (const SkPoint[])> & f,int maxExponent=30)115 static void for_random_beziers(int numPoints, SkRandom* rand,
116                                const std::function<void(const SkPoint[])>& f,
117                                int maxExponent = 30) {
118     SkASSERT(numPoints <= 4);
119     SkPoint pts[4];
120     for (int i = -10; i <= maxExponent; ++i) {
121         for (int j = 0; j < numPoints; ++j) {
122             pts[j].set(std::ldexp(1 + rand->nextF(), i), std::ldexp(1 + rand->nextF(), i));
123         }
124         f(pts);
125     }
126 }
127 
128 // Ensure the optimized "*_log2" versions return the same value as ceil(std::log2(f)).
DEF_TEST(wangs_formula_log2,r)129 DEF_TEST(wangs_formula_log2, r) {
130     // Constructs a cubic such that the 'length' term in wang's formula == term.
131     //
132     //     f = sqrt(k * length(max(abs(p0 - p1*2 + p2),
133     //                             abs(p1 - p2*2 + p3))));
134     auto setupCubicLengthTerm = [](int seed, SkPoint pts[], float term) {
135         memset(pts, 0, sizeof(SkPoint) * 4);
136 
137         SkPoint term2d = (seed & 1) ?
138                 SkPoint::Make(term, 0) : SkPoint::Make(.5f, std::sqrt(3)/2) * term;
139         seed >>= 1;
140 
141         if (seed & 1) {
142             term2d.fX = -term2d.fX;
143         }
144         seed >>= 1;
145 
146         if (seed & 1) {
147             std::swap(term2d.fX, term2d.fY);
148         }
149         seed >>= 1;
150 
151         switch (seed % 4) {
152             case 0:
153                 pts[0] = term2d;
154                 pts[3] = term2d * .75f;
155                 return;
156             case 1:
157                 pts[1] = term2d * -.5f;
158                 return;
159             case 2:
160                 pts[1] = term2d * -.5f;
161                 return;
162             case 3:
163                 pts[3] = term2d;
164                 pts[0] = term2d * .75f;
165                 return;
166         }
167     };
168 
169     // Constructs a quadratic such that the 'length' term in wang's formula == term.
170     //
171     //     f = sqrt(k * length(p0 - p1*2 + p2));
172     auto setupQuadraticLengthTerm = [](int seed, SkPoint pts[], float term) {
173         memset(pts, 0, sizeof(SkPoint) * 3);
174 
175         SkPoint term2d = (seed & 1) ?
176                 SkPoint::Make(term, 0) : SkPoint::Make(.5f, std::sqrt(3)/2) * term;
177         seed >>= 1;
178 
179         if (seed & 1) {
180             term2d.fX = -term2d.fX;
181         }
182         seed >>= 1;
183 
184         if (seed & 1) {
185             std::swap(term2d.fX, term2d.fY);
186         }
187         seed >>= 1;
188 
189         switch (seed % 3) {
190             case 0:
191                 pts[0] = term2d;
192                 return;
193             case 1:
194                 pts[1] = term2d * -.5f;
195                 return;
196             case 2:
197                 pts[2] = term2d;
198                 return;
199         }
200     };
201 
202     // wangs_formula_cubic and wangs_formula_quadratic both use rsqrt instead of sqrt for speed.
203     // Linearization is all approximate anyway, so as long as we are within ~1/2 tessellation
204     // segment of the reference value we are good enough.
205     constexpr static float kTessellationTolerance = 1/128.f;
206 
207     for (int level = 0; level < 30; ++level) {
208         float epsilon = std::ldexp(SK_ScalarNearlyZero, level * 2);
209         SkPoint pts[4];
210 
211         {
212             // Test cubic boundaries.
213             //     f = sqrt(k * length(max(abs(p0 - p1*2 + p2),
214             //                             abs(p1 - p2*2 + p3))));
215             constexpr static float k = (3 * 2) / (8 * (1.f/kPrecision));
216             float x = std::ldexp(1, level * 2) / k;
217             setupCubicLengthTerm(level << 1, pts, x - epsilon);
218             float referenceValue = wangs_formula_cubic_reference_impl(kPrecision, pts);
219             REPORTER_ASSERT(r, std::ceil(std::log2(referenceValue)) == level);
220             float c = wangs_formula::cubic(kPrecision, pts);
221             REPORTER_ASSERT(r, SkScalarNearlyEqual(c/referenceValue, 1, kTessellationTolerance));
222             REPORTER_ASSERT(r, wangs_formula::cubic_log2(kPrecision, pts) == level);
223             setupCubicLengthTerm(level << 1, pts, x + epsilon);
224             referenceValue = wangs_formula_cubic_reference_impl(kPrecision, pts);
225             REPORTER_ASSERT(r, std::ceil(std::log2(referenceValue)) == level + 1);
226             c = wangs_formula::cubic(kPrecision, pts);
227             REPORTER_ASSERT(r, SkScalarNearlyEqual(c/referenceValue, 1, kTessellationTolerance));
228             REPORTER_ASSERT(r, wangs_formula::cubic_log2(kPrecision, pts) == level + 1);
229         }
230 
231         {
232             // Test quadratic boundaries.
233             //     f = std::sqrt(k * Length(p0 - p1*2 + p2));
234             constexpr static float k = 2 / (8 * (1.f/kPrecision));
235             float x = std::ldexp(1, level * 2) / k;
236             setupQuadraticLengthTerm(level << 1, pts, x - epsilon);
237             float referenceValue = wangs_formula_quadratic_reference_impl(kPrecision, pts);
238             REPORTER_ASSERT(r, std::ceil(std::log2(referenceValue)) == level);
239             float q = wangs_formula::quadratic(kPrecision, pts);
240             REPORTER_ASSERT(r, SkScalarNearlyEqual(q/referenceValue, 1, kTessellationTolerance));
241             REPORTER_ASSERT(r, wangs_formula::quadratic_log2(kPrecision, pts) == level);
242             setupQuadraticLengthTerm(level << 1, pts, x + epsilon);
243             referenceValue = wangs_formula_quadratic_reference_impl(kPrecision, pts);
244             REPORTER_ASSERT(r, std::ceil(std::log2(referenceValue)) == level+1);
245             q = wangs_formula::quadratic(kPrecision, pts);
246             REPORTER_ASSERT(r, SkScalarNearlyEqual(q/referenceValue, 1, kTessellationTolerance));
247             REPORTER_ASSERT(r, wangs_formula::quadratic_log2(kPrecision, pts) == level + 1);
248         }
249     }
250 
251     auto check_cubic_log2 = [&](const SkPoint* pts) {
252         float f = std::max(1.f, wangs_formula_cubic_reference_impl(kPrecision, pts));
253         int f_log2 = wangs_formula::cubic_log2(kPrecision, pts);
254         REPORTER_ASSERT(r, SkScalarCeilToInt(std::log2(f)) == f_log2);
255         float c = std::max(1.f, wangs_formula::cubic(kPrecision, pts));
256         REPORTER_ASSERT(r, SkScalarNearlyEqual(c/f, 1, kTessellationTolerance));
257     };
258 
259     auto check_quadratic_log2 = [&](const SkPoint* pts) {
260         float f = std::max(1.f, wangs_formula_quadratic_reference_impl(kPrecision, pts));
261         int f_log2 = wangs_formula::quadratic_log2(kPrecision, pts);
262         REPORTER_ASSERT(r, SkScalarCeilToInt(std::log2(f)) == f_log2);
263         float q = std::max(1.f, wangs_formula::quadratic(kPrecision, pts));
264         REPORTER_ASSERT(r, SkScalarNearlyEqual(q/f, 1, kTessellationTolerance));
265     };
266 
267     SkRandom rand;
268 
269     for_random_matrices(&rand, [&](const SkMatrix& m) {
270         SkPoint pts[4];
271         m.mapPoints(pts, kSerp, 4);
272         check_cubic_log2(pts);
273 
274         m.mapPoints(pts, kLoop, 4);
275         check_cubic_log2(pts);
276 
277         m.mapPoints(pts, kQuad, 3);
278         check_quadratic_log2(pts);
279     });
280 
281     for_random_beziers(4, &rand, [&](const SkPoint pts[]) {
282         check_cubic_log2(pts);
283     });
284 
285     for_random_beziers(3, &rand, [&](const SkPoint pts[]) {
286         check_quadratic_log2(pts);
287     });
288 }
289 
290 // Ensure using transformations gives the same result as pre-transforming all points.
DEF_TEST(wangs_formula_vectorXforms,r)291 DEF_TEST(wangs_formula_vectorXforms, r) {
292     auto check_cubic_log2_with_transform = [&](const SkPoint* pts, const SkMatrix& m){
293         SkPoint ptsXformed[4];
294         m.mapPoints(ptsXformed, pts, 4);
295         int expected = wangs_formula::cubic_log2(kPrecision, ptsXformed);
296         int actual = wangs_formula::cubic_log2(kPrecision, pts, wangs_formula::VectorXform(m));
297         REPORTER_ASSERT(r, actual == expected);
298     };
299 
300     auto check_quadratic_log2_with_transform = [&](const SkPoint* pts, const SkMatrix& m) {
301         SkPoint ptsXformed[3];
302         m.mapPoints(ptsXformed, pts, 3);
303         int expected = wangs_formula::quadratic_log2(kPrecision, ptsXformed);
304         int actual = wangs_formula::quadratic_log2(kPrecision, pts, wangs_formula::VectorXform(m));
305         REPORTER_ASSERT(r, actual == expected);
306     };
307 
308     SkRandom rand;
309 
310     for_random_matrices(&rand, [&](const SkMatrix& m) {
311         check_cubic_log2_with_transform(kSerp, m);
312         check_cubic_log2_with_transform(kLoop, m);
313         check_quadratic_log2_with_transform(kQuad, m);
314 
315         for_random_beziers(4, &rand, [&](const SkPoint pts[]) {
316             check_cubic_log2_with_transform(pts, m);
317         });
318 
319         for_random_beziers(3, &rand, [&](const SkPoint pts[]) {
320             check_quadratic_log2_with_transform(pts, m);
321         });
322     });
323 }
324 
DEF_TEST(wangs_formula_worst_case_cubic,r)325 DEF_TEST(wangs_formula_worst_case_cubic, r) {
326     {
327         SkPoint worstP[] = {{0,0}, {100,100}, {0,0}, {0,0}};
328         REPORTER_ASSERT(r, wangs_formula::worst_case_cubic(kPrecision, 100, 100) ==
329                            wangs_formula_cubic_reference_impl(kPrecision, worstP));
330         REPORTER_ASSERT(r, wangs_formula::worst_case_cubic_log2(kPrecision, 100, 100) ==
331                            wangs_formula::cubic_log2(kPrecision, worstP));
332     }
333     {
334         SkPoint worstP[] = {{100,100}, {100,100}, {200,200}, {100,100}};
335         REPORTER_ASSERT(r, wangs_formula::worst_case_cubic(kPrecision, 100, 100) ==
336                            wangs_formula_cubic_reference_impl(kPrecision, worstP));
337         REPORTER_ASSERT(r, wangs_formula::worst_case_cubic_log2(kPrecision, 100, 100) ==
338                            wangs_formula::cubic_log2(kPrecision, worstP));
339     }
340     auto check_worst_case_cubic = [&](const SkPoint* pts) {
341         SkRect bbox;
342         bbox.setBoundsNoCheck(pts, 4);
343         float worst = wangs_formula::worst_case_cubic(kPrecision, bbox.width(), bbox.height());
344         int worst_log2 = wangs_formula::worst_case_cubic_log2(kPrecision, bbox.width(),
345                                                                bbox.height());
346         float actual = wangs_formula_cubic_reference_impl(kPrecision, pts);
347         REPORTER_ASSERT(r, worst >= actual);
348         REPORTER_ASSERT(r, std::ceil(std::log2(std::max(1.f, worst))) == worst_log2);
349     };
350     SkRandom rand;
351     for (int i = 0; i < 100; ++i) {
352         for_random_beziers(4, &rand, [&](const SkPoint pts[]) {
353             check_worst_case_cubic(pts);
354         });
355     }
356     // Make sure overflow saturates at infinity (not NaN).
357     constexpr static float inf = std::numeric_limits<float>::infinity();
358     REPORTER_ASSERT(r, wangs_formula::worst_case_cubic_p4(kPrecision, inf, inf) == inf);
359     REPORTER_ASSERT(r, wangs_formula::worst_case_cubic(kPrecision, inf, inf) == inf);
360 }
361 
362 // Ensure Wang's formula for quads produces max error within tolerance.
DEF_TEST(wangs_formula_quad_within_tol,r)363 DEF_TEST(wangs_formula_quad_within_tol, r) {
364     // Wang's formula and the quad math starts to lose precision with very large
365     // coordinate values, so limit the magnitude a bit to prevent test failures
366     // due to loss of precision.
367     constexpr int maxExponent = 15;
368     SkRandom rand;
369     for_random_beziers(3, &rand, [&r](const SkPoint pts[]) {
370         const int nsegs = static_cast<int>(
371                 std::ceil(wangs_formula_quadratic_reference_impl(kPrecision, pts)));
372 
373         const float tdelta = 1.f / nsegs;
374         for (int j = 0; j < nsegs; ++j) {
375             const float tmin = j * tdelta, tmax = (j + 1) * tdelta;
376 
377             // Get section of quad in [tmin,tmax]
378             const SkPoint* sectionPts;
379             SkPoint tmp0[5];
380             SkPoint tmp1[5];
381             if (tmin == 0) {
382                 if (tmax == 1) {
383                     sectionPts = pts;
384                 } else {
385                     SkChopQuadAt(pts, tmp0, tmax);
386                     sectionPts = tmp0;
387                 }
388             } else {
389                 SkChopQuadAt(pts, tmp0, tmin);
390                 if (tmax == 1) {
391                     sectionPts = tmp0 + 2;
392                 } else {
393                     SkChopQuadAt(tmp0 + 2, tmp1, (tmax - tmin) / (1 - tmin));
394                     sectionPts = tmp1;
395                 }
396             }
397 
398             // For quads, max distance from baseline is always at t=0.5.
399             SkPoint p;
400             p = SkEvalQuadAt(sectionPts, 0.5f);
401 
402             // Get distance of p to baseline
403             const SkPoint n = {sectionPts[2].fY - sectionPts[0].fY,
404                                sectionPts[0].fX - sectionPts[2].fX};
405             const float d = std::abs((p - sectionPts[0]).dot(n)) / n.length();
406 
407             // Check distance is within specified tolerance
408             REPORTER_ASSERT(r, d <= (1.f / kPrecision) + SK_ScalarNearlyZero);
409         }
410     }, maxExponent);
411 }
412 
413 // Ensure the specialized version for rational quads reduces to regular Wang's
414 // formula when all weights are equal to one
DEF_TEST(wangs_formula_rational_quad_reduces,r)415 DEF_TEST(wangs_formula_rational_quad_reduces, r) {
416     constexpr static float kTessellationTolerance = 1 / 128.f;
417 
418     SkRandom rand;
419     for (int i = 0; i < 100; ++i) {
420         for_random_beziers(3, &rand, [&r](const SkPoint pts[]) {
421             const float rational_nsegs = wangs_formula::conic(kPrecision, pts, 1.f);
422             const float integral_nsegs = wangs_formula_quadratic_reference_impl(kPrecision, pts);
423             REPORTER_ASSERT(
424                     r, SkScalarNearlyEqual(rational_nsegs, integral_nsegs, kTessellationTolerance));
425         });
426     }
427 }
428 
429 // Ensure the rational quad version (used for conics) produces max error within tolerance.
DEF_TEST(wangs_formula_conic_within_tol,r)430 DEF_TEST(wangs_formula_conic_within_tol, r) {
431     constexpr int maxExponent = 24;
432 
433     // Single-precision functions in SkConic/SkGeometry lose too much accuracy with
434     // large-magnitude curves and large weights for this test to pass.
435     using Sk2d = skvx::Vec<2, double>;
436     const auto eval_conic = [](const SkPoint pts[3], float w, float t) -> Sk2d {
437         const auto eval = [](Sk2d A, Sk2d B, Sk2d C, float t) -> Sk2d {
438             return (A * t + B) * t + C;
439         };
440 
441         const Sk2d p0 = {pts[0].fX, pts[0].fY};
442         const Sk2d p1 = {pts[1].fX, pts[1].fY};
443         const Sk2d p1w = p1 * w;
444         const Sk2d p2 = {pts[2].fX, pts[2].fY};
445         Sk2d numer = eval(p2 - p1w * 2 + p0, (p1w - p0) * 2, p0, t);
446 
447         Sk2d denomC = {1, 1};
448         Sk2d denomB = {2 * (w - 1), 2 * (w - 1)};
449         Sk2d denomA = {-2 * (w - 1), -2 * (w - 1)};
450         Sk2d denom = eval(denomA, denomB, denomC, t);
451         return numer / denom;
452     };
453 
454     const auto dot = [](const Sk2d& a, const Sk2d& b) -> double {
455         return a[0] * b[0] + a[1] * b[1];
456     };
457 
458     const auto length = [](const Sk2d& p) -> double { return sqrt(p[0] * p[0] + p[1] * p[1]); };
459 
460     SkRandom rand;
461     for (int i = -10; i <= 10; ++i) {
462         const float w = std::ldexp(1 + rand.nextF(), i);
463         for_random_beziers(
464                 3, &rand,
465                 [&](const SkPoint pts[]) {
466                     const int nsegs = SkScalarCeilToInt(wangs_formula::conic(kPrecision, pts, w));
467 
468                     const float tdelta = 1.f / nsegs;
469                     for (int j = 0; j < nsegs; ++j) {
470                         const float tmin = j * tdelta, tmax = (j + 1) * tdelta,
471                                     tmid = 0.5f * (tmin + tmax);
472 
473                         Sk2d p0, p1, p2;
474                         p0 = eval_conic(pts, w, tmin);
475                         p1 = eval_conic(pts, w, tmid);
476                         p2 = eval_conic(pts, w, tmax);
477 
478                         // Get distance of p1 to baseline (p0, p2).
479                         const Sk2d n = {p2[1] - p0[1], p0[0] - p2[0]};
480                         SkASSERT(length(n) != 0);
481                         const double d = std::abs(dot(p1 - p0, n)) / length(n);
482 
483                         // Check distance is within tolerance
484                         REPORTER_ASSERT(r, d <= (1.0 / kPrecision) + SK_ScalarNearlyZero);
485                     }
486                 },
487                 maxExponent);
488     }
489 }
490 
491 // Ensure the vectorized conic version equals the reference implementation
DEF_TEST(wangs_formula_conic_matches_reference,r)492 DEF_TEST(wangs_formula_conic_matches_reference, r) {
493     SkRandom rand;
494     for (int i = -10; i <= 10; ++i) {
495         const float w = std::ldexp(1 + rand.nextF(), i);
496         for_random_beziers(3, &rand, [&r, w](const SkPoint pts[]) {
497             const float ref_nsegs = wangs_formula_conic_reference_impl(kPrecision, pts, w);
498             const float nsegs = wangs_formula::conic(kPrecision, pts, w);
499 
500             // Because the Gr version may implement the math differently for performance,
501             // allow different slack in the comparison based on the rough scale of the answer.
502             const float cmpThresh = ref_nsegs * (1.f / (1 << 20));
503             REPORTER_ASSERT(r, SkScalarNearlyEqual(ref_nsegs, nsegs, cmpThresh));
504         });
505     }
506 }
507 
508 // Ensure using transformations gives the same result as pre-transforming all points.
DEF_TEST(wangs_formula_conic_vectorXforms,r)509 DEF_TEST(wangs_formula_conic_vectorXforms, r) {
510     auto check_conic_with_transform = [&](const SkPoint* pts, float w, const SkMatrix& m) {
511         SkPoint ptsXformed[3];
512         m.mapPoints(ptsXformed, pts, 3);
513         float expected = wangs_formula::conic(kPrecision, ptsXformed, w);
514         float actual = wangs_formula::conic(kPrecision, pts, w, wangs_formula::VectorXform(m));
515         REPORTER_ASSERT(r, SkScalarNearlyEqual(actual, expected));
516     };
517 
518     SkRandom rand;
519     for (int i = -10; i <= 10; ++i) {
520         const float w = std::ldexp(1 + rand.nextF(), i);
521         for_random_beziers(3, &rand, [&](const SkPoint pts[]) {
522             check_conic_with_transform(pts, w, SkMatrix::I());
523             check_conic_with_transform(
524                     pts, w, SkMatrix::Scale(rand.nextRangeF(-10, 10), rand.nextRangeF(-10, 10)));
525 
526             // Random 2x2 matrix
527             SkMatrix m;
528             m.setScaleX(rand.nextRangeF(-10, 10));
529             m.setSkewX(rand.nextRangeF(-10, 10));
530             m.setSkewY(rand.nextRangeF(-10, 10));
531             m.setScaleY(rand.nextRangeF(-10, 10));
532             check_conic_with_transform(pts, w, m);
533         });
534     }
535 }
536 
DEF_TEST(wangs_formula_nextlog2,r)537 DEF_TEST(wangs_formula_nextlog2, r) {
538     REPORTER_ASSERT(r, 0b0'00000000'111'1111111111'1111111111 == (1u << 23) - 1u);
539     REPORTER_ASSERT(r, wangs_formula::nextlog2(-std::numeric_limits<float>::infinity()) == 0);
540     REPORTER_ASSERT(r, wangs_formula::nextlog2(-std::numeric_limits<float>::max()) == 0);
541     REPORTER_ASSERT(r, wangs_formula::nextlog2(-1000.0f) == 0);
542     REPORTER_ASSERT(r, wangs_formula::nextlog2(-0.1f) == 0);
543     REPORTER_ASSERT(r, wangs_formula::nextlog2(-std::numeric_limits<float>::min()) == 0);
544     REPORTER_ASSERT(r, wangs_formula::nextlog2(-std::numeric_limits<float>::denorm_min()) == 0);
545     REPORTER_ASSERT(r, wangs_formula::nextlog2(0.0f) == 0);
546     REPORTER_ASSERT(r, wangs_formula::nextlog2(std::numeric_limits<float>::denorm_min()) == 0);
547     REPORTER_ASSERT(r, wangs_formula::nextlog2(std::numeric_limits<float>::min()) == 0);
548     REPORTER_ASSERT(r, wangs_formula::nextlog2(0.1f) == 0);
549     REPORTER_ASSERT(r, wangs_formula::nextlog2(1.0f) == 0);
550     REPORTER_ASSERT(r, wangs_formula::nextlog2(1.1f) == 1);
551     REPORTER_ASSERT(r, wangs_formula::nextlog2(2.0f) == 1);
552     REPORTER_ASSERT(r, wangs_formula::nextlog2(2.1f) == 2);
553     REPORTER_ASSERT(r, wangs_formula::nextlog2(3.0f) == 2);
554     REPORTER_ASSERT(r, wangs_formula::nextlog2(3.1f) == 2);
555     REPORTER_ASSERT(r, wangs_formula::nextlog2(4.0f) == 2);
556     REPORTER_ASSERT(r, wangs_formula::nextlog2(4.1f) == 3);
557     REPORTER_ASSERT(r, wangs_formula::nextlog2(5.0f) == 3);
558     REPORTER_ASSERT(r, wangs_formula::nextlog2(5.1f) == 3);
559     REPORTER_ASSERT(r, wangs_formula::nextlog2(6.0f) == 3);
560     REPORTER_ASSERT(r, wangs_formula::nextlog2(6.1f) == 3);
561     REPORTER_ASSERT(r, wangs_formula::nextlog2(7.0f) == 3);
562     REPORTER_ASSERT(r, wangs_formula::nextlog2(7.1f) == 3);
563     REPORTER_ASSERT(r, wangs_formula::nextlog2(8.0f) == 3);
564     REPORTER_ASSERT(r, wangs_formula::nextlog2(8.1f) == 4);
565     REPORTER_ASSERT(r, wangs_formula::nextlog2(9.0f) == 4);
566     REPORTER_ASSERT(r, wangs_formula::nextlog2(9.1f) == 4);
567     REPORTER_ASSERT(r, wangs_formula::nextlog2(std::numeric_limits<float>::max()) == 128);
568     REPORTER_ASSERT(r, wangs_formula::nextlog2(std::numeric_limits<float>::infinity()) == 128);
569     REPORTER_ASSERT(r, wangs_formula::nextlog2(std::numeric_limits<float>::quiet_NaN()) == 0);
570     REPORTER_ASSERT(r, wangs_formula::nextlog2(-std::numeric_limits<float>::quiet_NaN()) == 0);
571 
572     for (int i = 0; i < 100; ++i) {
573         float pow2 = std::ldexp(1, i);
574         float epsilon = std::ldexp(SK_ScalarNearlyZero, i);
575         REPORTER_ASSERT(r, wangs_formula::nextlog2(pow2) == i);
576         REPORTER_ASSERT(r, wangs_formula::nextlog2(pow2 + epsilon) == i + 1);
577         REPORTER_ASSERT(r, wangs_formula::nextlog2(pow2 - epsilon) == i);
578     }
579 }
580 
581 }  // namespace skgpu::tess
582