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 #ifndef GrWangsFormula_DEFINED
9 #define GrWangsFormula_DEFINED
10
11 #include "include/core/SkPoint.h"
12 #include "include/private/SkFloatingPoint.h"
13 #include "src/gpu/GrVx.h"
14 #include "src/gpu/tessellate/GrVectorXform.h"
15
16 // Wang's formula gives the minimum number of evenly spaced (in the parametric sense) line segments
17 // that a bezier curve must be chopped into in order to guarantee all lines stay within a distance
18 // of "1/precision" pixels from the true curve. Its definition for a bezier curve of degree "n" is
19 // as follows:
20 //
21 // maxLength = max([length(p[i+2] - 2p[i+1] + p[i]) for (0 <= i <= n-2)])
22 // numParametricSegments = sqrt(maxLength * precision * n*(n - 1)/8)
23 //
24 // (Goldman, Ron. (2003). 5.6.3 Wang's Formula. "Pyramid Algorithms: A Dynamic Programming Approach
25 // to Curves and Surfaces for Geometric Modeling". Morgan Kaufmann Publishers.)
26 namespace GrWangsFormula {
27
28 // Returns the value by which to multiply length in Wang's formula. (See above.)
length_term(float precision)29 template<int Degree> constexpr float length_term(float precision) {
30 return (Degree * (Degree - 1) / 8.f) * precision;
31 }
length_term_pow2(float precision)32 template<int Degree> constexpr float length_term_pow2(float precision) {
33 return ((Degree * Degree) * ((Degree - 1) * (Degree - 1)) / 64.f) * (precision * precision);
34 }
35
root4(float x)36 SK_ALWAYS_INLINE static float root4(float x) {
37 return sqrtf(sqrtf(x));
38 }
39
40 // Returns nextlog2(sqrt(x)):
41 //
42 // log2(sqrt(x)) == log2(x^(1/2)) == log2(x)/2 == log2(x)/log2(4) == log4(x)
43 //
nextlog4(float x)44 SK_ALWAYS_INLINE static int nextlog4(float x) {
45 return (sk_float_nextlog2(x) + 1) >> 1;
46 }
47
48 // Returns nextlog2(sqrt(sqrt(x))):
49 //
50 // log2(sqrt(sqrt(x))) == log2(x^(1/4)) == log2(x)/4 == log2(x)/log2(16) == log16(x)
51 //
nextlog16(float x)52 SK_ALWAYS_INLINE static int nextlog16(float x) {
53 return (sk_float_nextlog2(x) + 3) >> 2;
54 }
55
56 // Returns Wang's formula, raised to the 4th power, specialized for a quadratic curve.
57 SK_ALWAYS_INLINE static float quadratic_pow4(float precision, const SkPoint pts[],
58 const GrVectorXform& vectorXform = GrVectorXform()) {
59 using grvx::float2, skvx::bit_pun;
60 float2 p0 = bit_pun<float2>(pts[0]);
61 float2 p1 = bit_pun<float2>(pts[1]);
62 float2 p2 = bit_pun<float2>(pts[2]);
63 float2 v = grvx::fast_madd<2>(-2, p1, p0) + p2;
64 v = vectorXform(v);
65 float2 vv = v*v;
66 return (vv[0] + vv[1]) * length_term_pow2<2>(precision);
67 }
68
69 // Returns Wang's formula specialized for a quadratic curve.
70 SK_ALWAYS_INLINE static float quadratic(float precision, const SkPoint pts[],
71 const GrVectorXform& vectorXform = GrVectorXform()) {
72 return root4(quadratic_pow4(precision, pts, vectorXform));
73 }
74
75 // Returns the log2 value of Wang's formula specialized for a quadratic curve, rounded up to the
76 // next int.
77 SK_ALWAYS_INLINE static int quadratic_log2(float precision, const SkPoint pts[],
78 const GrVectorXform& vectorXform = GrVectorXform()) {
79 // nextlog16(x) == ceil(log2(sqrt(sqrt(x))))
80 return nextlog16(quadratic_pow4(precision, pts, vectorXform));
81 }
82
83 // Returns Wang's formula, raised to the 4th power, specialized for a cubic curve.
84 SK_ALWAYS_INLINE static float cubic_pow4(float precision, const SkPoint pts[],
85 const GrVectorXform& vectorXform = GrVectorXform()) {
86 using grvx::float4;
87 float4 p01 = float4::Load(pts);
88 float4 p12 = float4::Load(pts + 1);
89 float4 p23 = float4::Load(pts + 2);
90 float4 v = grvx::fast_madd<4>(-2, p12, p01) + p23;
91 v = vectorXform(v);
92 float4 vv = v*v;
93 return std::max(vv[0] + vv[1], vv[2] + vv[3]) * length_term_pow2<3>(precision);
94 }
95
96 // Returns Wang's formula specialized for a cubic curve.
97 SK_ALWAYS_INLINE static float cubic(float precision, const SkPoint pts[],
98 const GrVectorXform& vectorXform = GrVectorXform()) {
99 return root4(cubic_pow4(precision, pts, vectorXform));
100 }
101
102 // Returns the log2 value of Wang's formula specialized for a cubic curve, rounded up to the next
103 // int.
104 SK_ALWAYS_INLINE static int cubic_log2(float precision, const SkPoint pts[],
105 const GrVectorXform& vectorXform = GrVectorXform()) {
106 // nextlog16(x) == ceil(log2(sqrt(sqrt(x))))
107 return nextlog16(cubic_pow4(precision, pts, vectorXform));
108 }
109
110 // Returns the maximum number of line segments a cubic with the given device-space bounding box size
111 // would ever need to be divided into. This is simply a special case of the cubic formula where we
112 // maximize its value by placing control points on specific corners of the bounding box.
worst_case_cubic(float precision,float devWidth,float devHeight)113 SK_ALWAYS_INLINE static float worst_case_cubic(float precision, float devWidth, float devHeight) {
114 float k = length_term<3>(precision);
115 return sqrtf(2*k * SkVector::Length(devWidth, devHeight));
116 }
117
118 // Returns the maximum log2 number of line segments a cubic with the given device-space bounding box
119 // size would ever need to be divided into.
worst_case_cubic_log2(float precision,float devWidth,float devHeight)120 SK_ALWAYS_INLINE static int worst_case_cubic_log2(float precision, float devWidth,
121 float devHeight) {
122 float kk = length_term_pow2<3>(precision);
123 // nextlog16(x) == ceil(log2(sqrt(sqrt(x))))
124 return nextlog16(4*kk * (devWidth * devWidth + devHeight * devHeight));
125 }
126
127 // Returns Wang's formula specialized for a conic curve, raised to the second power.
128 // Input points should be in projected space, and note tolerance parameter is not "precision".
129 //
130 // This is not actually due to Wang, but is an analogue from (Theorem 3, corollary 1):
131 // J. Zheng, T. Sederberg. "Estimating Tessellation Parameter Intervals for
132 // Rational Curves and Surfaces." ACM Transactions on Graphics 19(1). 2000.
133 SK_ALWAYS_INLINE static float conic_pow2(float tolerance, const SkPoint pts[], float w,
134 const GrVectorXform& vectorXform = GrVectorXform()) {
135 using grvx::dot, grvx::float2, grvx::float4, skvx::bit_pun;
136 float2 p0 = vectorXform(bit_pun<float2>(pts[0]));
137 float2 p1 = vectorXform(bit_pun<float2>(pts[1]));
138 float2 p2 = vectorXform(bit_pun<float2>(pts[2]));
139
140 // Compute center of bounding box in projected space
141 const float2 C = 0.5f * (skvx::min(skvx::min(p0, p1), p2) + skvx::max(skvx::max(p0, p1), p2));
142
143 // Translate by -C. This improves translation-invariance of the formula,
144 // see Sec. 3.3 of cited paper
145 p0 -= C;
146 p1 -= C;
147 p2 -= C;
148
149 // Compute max length
150 const float max_len = sqrtf(std::max(dot(p0, p0), std::max(dot(p1, p1), dot(p2, p2))));
151
152 // Compute forward differences
153 const float2 dp = grvx::fast_madd<2>(-2 * w, p1, p0) + p2;
154 const float dw = fabsf(1 - 2 * w + 1);
155
156 // Compute numerator and denominator for parametric step size of linearization
157 const float r_minus_eps = std::max(0.f, max_len - tolerance);
158 const float min_w = std::min(w, 1.f);
159 const float numer = sqrtf(grvx::dot(dp, dp)) + r_minus_eps * dw;
160 const float denom = 4 * min_w * tolerance;
161
162 // Number of segments = sqrt(numer / denom).
163 // This assumes parametric interval of curve being linearized is [t0,t1] = [0, 1].
164 // If not, the number of segments is (tmax - tmin) / sqrt(denom / numer).
165 return numer / denom;
166 }
167
168 // Returns the value of Wang's formula specialized for a conic curve.
169 SK_ALWAYS_INLINE static float conic(float tolerance, const SkPoint pts[], float w,
170 const GrVectorXform& vectorXform = GrVectorXform()) {
171 return sqrtf(conic_pow2(tolerance, pts, w, vectorXform));
172 }
173
174 // Returns the log2 value of Wang's formula specialized for a conic curve, rounded up to the next
175 // int.
176 SK_ALWAYS_INLINE static int conic_log2(float tolerance, const SkPoint pts[], float w,
177 const GrVectorXform& vectorXform = GrVectorXform()) {
178 // nextlog4(x) == ceil(log2(sqrt(x)))
179 return nextlog4(conic_pow2(tolerance, pts, w, vectorXform));
180 }
181
182 } // namespace GrWangsFormula
183
184 #endif
185