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