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 skgpu_tessellate_WangsFormula_DEFINED
9 #define skgpu_tessellate_WangsFormula_DEFINED
10
11 #include "include/core/SkM44.h"
12 #include "include/core/SkMatrix.h"
13 #include "include/core/SkPoint.h"
14 #include "include/core/SkTypes.h"
15 #include "src/base/SkFloatBits.h"
16 #include "src/base/SkUtils.h"
17 #include "src/base/SkVx.h"
18
19 #include <math.h>
20 #include <algorithm>
21 #include <cstdint>
22 #include <limits>
23
24 #define AI [[maybe_unused]] SK_ALWAYS_INLINE
25
26 // Wang's formula gives the minimum number of evenly spaced (in the parametric sense) line segments
27 // that a bezier curve must be chopped into in order to guarantee all lines stay within a distance
28 // of "1/precision" pixels from the true curve. Its definition for a bezier curve of degree "n" is
29 // as follows:
30 //
31 // maxLength = max([length(p[i+2] - 2p[i+1] + p[i]) for (0 <= i <= n-2)])
32 // numParametricSegments = sqrt(maxLength * precision * n*(n - 1)/8)
33 //
34 // (Goldman, Ron. (2003). 5.6.3 Wang's Formula. "Pyramid Algorithms: A Dynamic Programming Approach
35 // to Curves and Surfaces for Geometric Modeling". Morgan Kaufmann Publishers.)
36 namespace skgpu::wangs_formula {
37
38 // Returns the value by which to multiply length in Wang's formula. (See above.)
length_term(float precision)39 template<int Degree> constexpr float length_term(float precision) {
40 return (Degree * (Degree - 1) / 8.f) * precision;
41 }
length_term_p2(float precision)42 template<int Degree> constexpr float length_term_p2(float precision) {
43 return ((Degree * Degree) * ((Degree - 1) * (Degree - 1)) / 64.f) * (precision * precision);
44 }
45
root4(float x)46 AI float root4(float x) {
47 return sqrtf(sqrtf(x));
48 }
49
50 // For finite positive x > 1, return ceil(log2(x)) otherwise, return 0.
51 // For +/- NaN return 0.
52 // For +infinity return 128.
53 // For -infinity return 0.
54 //
55 // nextlog2((-inf..1]) -> 0
56 // nextlog2((1..2]) -> 1
57 // nextlog2((2..4]) -> 2
58 // nextlog2((4..8]) -> 3
59 // ...
nextlog2(float x)60 AI int nextlog2(float x) {
61 if (x <= 1) {
62 return 0;
63 }
64
65 uint32_t bits = SkFloat2Bits(x);
66 static constexpr uint32_t kDigitsAfterBinaryPoint = std::numeric_limits<float>::digits - 1;
67
68 // The constant is a significand of all 1s -- 0b0'00000000'111'1111111111'111111111. So, if
69 // the significand of x is all 0s (and therefore an integer power of two) this will not
70 // increment the exponent, but if it is just one ULP above the power of two the carry will
71 // ripple into the exponent incrementing the exponent by 1.
72 bits += (1u << kDigitsAfterBinaryPoint) - 1u;
73
74 // Shift the exponent down, and adjust it by the exponent offset so that 2^0 is really 0 instead
75 // of 127. Remember that 1 was added to the exponent, if x is NaN, then the exponent will
76 // carry a 1 into the sign bit during the addition to bits. Be sure to strip off the sign bit.
77 // In addition, infinity is an exponent of all 1's, and a significand of all 0, so
78 // the exponent is not affected during the addition to bits, and the exponent remains all 1's.
79 const int exp = ((bits >> kDigitsAfterBinaryPoint) & 0b1111'1111) - 127;
80
81 // Return 0 for x <= 1.
82 return exp > 0 ? exp : 0;
83 }
84
85 // Returns nextlog2(sqrt(x)):
86 //
87 // log2(sqrt(x)) == log2(x^(1/2)) == log2(x)/2 == log2(x)/log2(4) == log4(x)
88 //
89 AI int nextlog4(float x) {
90 return (nextlog2(x) + 1) >> 1;
91 }
92
93 // Returns nextlog2(sqrt(sqrt(x))):
94 //
95 // log2(sqrt(sqrt(x))) == log2(x^(1/4)) == log2(x)/4 == log2(x)/log2(16) == log16(x)
96 //
97 AI int nextlog16(float x) {
98 return (nextlog2(x) + 3) >> 2;
99 }
100
101 // Represents the upper-left 2x2 matrix of an affine transform for applying to vectors:
102 //
103 // VectorXform(p1 - p0) == M * float3(p1, 1) - M * float3(p0, 1)
104 //
105 class VectorXform {
106 public:
107 AI VectorXform() : fC0{1.0f, 0.f}, fC1{0.f, 1.f} {}
108 AI explicit VectorXform(const SkMatrix& m) { *this = m; }
109 AI explicit VectorXform(const SkM44& m) { *this = m; }
110
111 AI VectorXform& operator=(const SkMatrix& m) {
112 SkASSERT(!m.hasPerspective());
113 fC0 = {m.rc(0,0), m.rc(1,0)};
114 fC1 = {m.rc(0,1), m.rc(1,1)};
115 return *this;
116 }
117 AI VectorXform& operator=(const SkM44& m) {
118 SkASSERT(m.rc(3,0) == 0.f && m.rc(3,1) == 0.f && m.rc(3,2) == 0.f && m.rc(3,3) == 1.f);
119 fC0 = {m.rc(0,0), m.rc(1,0)};
120 fC1 = {m.rc(0,1), m.rc(1,1)};
121 return *this;
122 }
123 AI skvx::float2 operator()(skvx::float2 vector) const {
124 return fC0 * vector.x() + fC1 * vector.y();
125 }
126 AI skvx::float4 operator()(skvx::float4 vectors) const {
127 return join(fC0 * vectors.x() + fC1 * vectors.y(),
128 fC0 * vectors.z() + fC1 * vectors.w());
129 }
130 private:
131 // First and second columns of 2x2 matrix
132 skvx::float2 fC0;
133 skvx::float2 fC1;
134 };
135
136 // Returns Wang's formula, raised to the 4th power, specialized for a quadratic curve.
137 AI float quadratic_p4(float precision,
138 skvx::float2 p0, skvx::float2 p1, skvx::float2 p2,
139 const VectorXform& vectorXform = VectorXform()) {
140 skvx::float2 v = -2*p1 + p0 + p2;
141 v = vectorXform(v);
142 skvx::float2 vv = v*v;
143 return (vv[0] + vv[1]) * length_term_p2<2>(precision);
144 }
145 AI float quadratic_p4(float precision,
146 const SkPoint pts[],
147 const VectorXform& vectorXform = VectorXform()) {
148 return quadratic_p4(precision,
149 sk_bit_cast<skvx::float2>(pts[0]),
150 sk_bit_cast<skvx::float2>(pts[1]),
151 sk_bit_cast<skvx::float2>(pts[2]),
152 vectorXform);
153 }
154
155 // Returns Wang's formula specialized for a quadratic curve.
156 AI float quadratic(float precision,
157 const SkPoint pts[],
158 const VectorXform& vectorXform = VectorXform()) {
159 return root4(quadratic_p4(precision, pts, vectorXform));
160 }
161
162 // Returns the log2 value of Wang's formula specialized for a quadratic curve, rounded up to the
163 // next int.
164 AI int quadratic_log2(float precision,
165 const SkPoint pts[],
166 const VectorXform& vectorXform = VectorXform()) {
167 // nextlog16(x) == ceil(log2(sqrt(sqrt(x))))
168 return nextlog16(quadratic_p4(precision, pts, vectorXform));
169 }
170
171 // Returns Wang's formula, raised to the 4th power, specialized for a cubic curve.
172 AI float cubic_p4(float precision,
173 skvx::float2 p0, skvx::float2 p1, skvx::float2 p2, skvx::float2 p3,
174 const VectorXform& vectorXform = VectorXform()) {
175 skvx::float4 p01{p0, p1};
176 skvx::float4 p12{p1, p2};
177 skvx::float4 p23{p2, p3};
178 skvx::float4 v = -2*p12 + p01 + p23;
179 v = vectorXform(v);
180 skvx::float4 vv = v*v;
181 return std::max(vv[0] + vv[1], vv[2] + vv[3]) * length_term_p2<3>(precision);
182 }
183 AI float cubic_p4(float precision,
184 const SkPoint pts[],
185 const VectorXform& vectorXform = VectorXform()) {
186 return cubic_p4(precision,
187 sk_bit_cast<skvx::float2>(pts[0]),
188 sk_bit_cast<skvx::float2>(pts[1]),
189 sk_bit_cast<skvx::float2>(pts[2]),
190 sk_bit_cast<skvx::float2>(pts[3]),
191 vectorXform);
192 }
193
194 // Returns Wang's formula specialized for a cubic curve.
195 AI float cubic(float precision,
196 const SkPoint pts[],
197 const VectorXform& vectorXform = VectorXform()) {
198 return root4(cubic_p4(precision, pts, vectorXform));
199 }
200
201 // Returns the log2 value of Wang's formula specialized for a cubic curve, rounded up to the next
202 // int.
203 AI int cubic_log2(float precision,
204 const SkPoint pts[],
205 const VectorXform& vectorXform = VectorXform()) {
206 // nextlog16(x) == ceil(log2(sqrt(sqrt(x))))
207 return nextlog16(cubic_p4(precision, pts, vectorXform));
208 }
209
210 // Returns the maximum number of line segments a cubic with the given device-space bounding box size
211 // would ever need to be divided into, raised to the 4th power. This is simply a special case of the
212 // cubic formula where we maximize its value by placing control points on specific corners of the
213 // bounding box.
214 AI float worst_case_cubic_p4(float precision, float devWidth, float devHeight) {
215 float kk = length_term_p2<3>(precision);
216 return 4*kk * (devWidth * devWidth + devHeight * devHeight);
217 }
218
219 // Returns the maximum number of line segments a cubic with the given device-space bounding box size
220 // would ever need to be divided into.
221 AI float worst_case_cubic(float precision, float devWidth, float devHeight) {
222 return root4(worst_case_cubic_p4(precision, devWidth, devHeight));
223 }
224
225 // Returns the maximum log2 number of line segments a cubic with the given device-space bounding box
226 // size would ever need to be divided into.
227 AI int worst_case_cubic_log2(float precision, float devWidth, float devHeight) {
228 // nextlog16(x) == ceil(log2(sqrt(sqrt(x))))
229 return nextlog16(worst_case_cubic_p4(precision, devWidth, devHeight));
230 }
231
232 // Returns Wang's formula specialized for a conic curve, raised to the second power.
233 // Input points should be in projected space.
234 //
235 // This is not actually due to Wang, but is an analogue from (Theorem 3, corollary 1):
236 // J. Zheng, T. Sederberg. "Estimating Tessellation Parameter Intervals for
237 // Rational Curves and Surfaces." ACM Transactions on Graphics 19(1). 2000.
238 AI float conic_p2(float precision,
239 skvx::float2 p0, skvx::float2 p1, skvx::float2 p2,
240 float w,
241 const VectorXform& vectorXform = VectorXform()) {
242 p0 = vectorXform(p0);
243 p1 = vectorXform(p1);
244 p2 = vectorXform(p2);
245
246 // Compute center of bounding box in projected space
247 const skvx::float2 C = 0.5f * (min(min(p0, p1), p2) + max(max(p0, p1), p2));
248
249 // Translate by -C. This improves translation-invariance of the formula,
250 // see Sec. 3.3 of cited paper
251 p0 -= C;
252 p1 -= C;
253 p2 -= C;
254
255 // Compute max length
256 const float max_len = sqrtf(std::max(dot(p0, p0), std::max(dot(p1, p1), dot(p2, p2))));
257
258
259 // Compute forward differences
260 const skvx::float2 dp = -2*w*p1 + p0 + p2;
261 const float dw = fabsf(-2 * w + 2);
262
263 // Compute numerator and denominator for parametric step size of linearization. Here, the
264 // epsilon referenced from the cited paper is 1/precision.
265 const float rp_minus_1 = std::max(0.f, max_len * precision - 1);
266 const float numer = sqrtf(dot(dp, dp)) * precision + rp_minus_1 * dw;
267 const float denom = 4 * std::min(w, 1.f);
268
269 // Number of segments = sqrt(numer / denom).
270 // This assumes parametric interval of curve being linearized is [t0,t1] = [0, 1].
271 // If not, the number of segments is (tmax - tmin) / sqrt(denom / numer).
272 return numer / denom;
273 }
274 AI float conic_p2(float precision,
275 const SkPoint pts[],
276 float w,
277 const VectorXform& vectorXform = VectorXform()) {
278 return conic_p2(precision,
279 sk_bit_cast<skvx::float2>(pts[0]),
280 sk_bit_cast<skvx::float2>(pts[1]),
281 sk_bit_cast<skvx::float2>(pts[2]),
282 w,
283 vectorXform);
284 }
285
286 // Returns the value of Wang's formula specialized for a conic curve.
287 AI float conic(float tolerance,
288 const SkPoint pts[],
289 float w,
290 const VectorXform& vectorXform = VectorXform()) {
291 return sqrtf(conic_p2(tolerance, pts, w, vectorXform));
292 }
293
294 // Returns the log2 value of Wang's formula specialized for a conic curve, rounded up to the next
295 // int.
296 AI int conic_log2(float tolerance,
297 const SkPoint pts[],
298 float w,
299 const VectorXform& vectorXform = VectorXform()) {
300 // nextlog4(x) == ceil(log2(sqrt(x)))
301 return nextlog4(conic_p2(tolerance, pts, w, vectorXform));
302 }
303
304 } // namespace skgpu::wangs_formula
305
306 #undef AI
307
308 #endif // skgpu_tessellate_WangsFormula_DEFINED
309