OLD | NEW |
| (Empty) |
1 /* | |
2 * Copyright 2012 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 #include "CubicUtilities.h" | |
8 #include "Extrema.h" | |
9 #include "LineUtilities.h" | |
10 #include "QuadraticUtilities.h" | |
11 | |
12 const int gPrecisionUnit = 256; // FIXME: arbitrary -- should try different valu
es in test framework | |
13 | |
14 // FIXME: cache keep the bounds and/or precision with the caller? | |
15 double calcPrecision(const Cubic& cubic) { | |
16 _Rect dRect; | |
17 dRect.setBounds(cubic); // OPTIMIZATION: just use setRawBounds ? | |
18 double width = dRect.right - dRect.left; | |
19 double height = dRect.bottom - dRect.top; | |
20 return (width > height ? width : height) / gPrecisionUnit; | |
21 } | |
22 | |
23 #ifdef SK_DEBUG | |
24 double calcPrecision(const Cubic& cubic, double t, double scale) { | |
25 Cubic part; | |
26 sub_divide(cubic, SkTMax(0., t - scale), SkTMin(1., t + scale), part); | |
27 return calcPrecision(part); | |
28 } | |
29 #endif | |
30 | |
31 bool clockwise(const Cubic& c) { | |
32 double sum = (c[0].x - c[3].x) * (c[0].y + c[3].y); | |
33 for (int idx = 0; idx < 3; ++idx){ | |
34 sum += (c[idx + 1].x - c[idx].x) * (c[idx + 1].y + c[idx].y); | |
35 } | |
36 return sum <= 0; | |
37 } | |
38 | |
39 void coefficients(const double* cubic, double& A, double& B, double& C, double&
D) { | |
40 A = cubic[6]; // d | |
41 B = cubic[4] * 3; // 3*c | |
42 C = cubic[2] * 3; // 3*b | |
43 D = cubic[0]; // a | |
44 A -= D - C + B; // A = -a + 3*b - 3*c + d | |
45 B += 3 * D - 2 * C; // B = 3*a - 6*b + 3*c | |
46 C -= 3 * D; // C = -3*a + 3*b | |
47 } | |
48 | |
49 bool controls_contained_by_ends(const Cubic& c) { | |
50 _Vector startTan = c[1] - c[0]; | |
51 if (startTan.x == 0 && startTan.y == 0) { | |
52 startTan = c[2] - c[0]; | |
53 } | |
54 _Vector endTan = c[2] - c[3]; | |
55 if (endTan.x == 0 && endTan.y == 0) { | |
56 endTan = c[1] - c[3]; | |
57 } | |
58 if (startTan.dot(endTan) >= 0) { | |
59 return false; | |
60 } | |
61 _Line startEdge = {c[0], c[0]}; | |
62 startEdge[1].x -= startTan.y; | |
63 startEdge[1].y += startTan.x; | |
64 _Line endEdge = {c[3], c[3]}; | |
65 endEdge[1].x -= endTan.y; | |
66 endEdge[1].y += endTan.x; | |
67 double leftStart1 = is_left(startEdge, c[1]); | |
68 if (leftStart1 * is_left(startEdge, c[2]) < 0) { | |
69 return false; | |
70 } | |
71 double leftEnd1 = is_left(endEdge, c[1]); | |
72 if (leftEnd1 * is_left(endEdge, c[2]) < 0) { | |
73 return false; | |
74 } | |
75 return leftStart1 * leftEnd1 >= 0; | |
76 } | |
77 | |
78 bool ends_are_extrema_in_x_or_y(const Cubic& c) { | |
79 return (between(c[0].x, c[1].x, c[3].x) && between(c[0].x, c[2].x, c[3].x)) | |
80 || (between(c[0].y, c[1].y, c[3].y) && between(c[0].y, c[2].y, c[3].
y)); | |
81 } | |
82 | |
83 bool monotonic_in_y(const Cubic& c) { | |
84 return between(c[0].y, c[1].y, c[3].y) && between(c[0].y, c[2].y, c[3].y); | |
85 } | |
86 | |
87 bool serpentine(const Cubic& c) { | |
88 if (!controls_contained_by_ends(c)) { | |
89 return false; | |
90 } | |
91 double wiggle = (c[0].x - c[2].x) * (c[0].y + c[2].y); | |
92 for (int idx = 0; idx < 2; ++idx){ | |
93 wiggle += (c[idx + 1].x - c[idx].x) * (c[idx + 1].y + c[idx].y); | |
94 } | |
95 double waggle = (c[1].x - c[3].x) * (c[1].y + c[3].y); | |
96 for (int idx = 1; idx < 3; ++idx){ | |
97 waggle += (c[idx + 1].x - c[idx].x) * (c[idx + 1].y + c[idx].y); | |
98 } | |
99 return wiggle * waggle < 0; | |
100 } | |
101 | |
102 // cubic roots | |
103 | |
104 const double PI = 4 * atan(1); | |
105 | |
106 // from SkGeometry.cpp (and Numeric Solutions, 5.6) | |
107 int cubicRootsValidT(double A, double B, double C, double D, double t[3]) { | |
108 #if 0 | |
109 if (approximately_zero(A)) { // we're just a quadratic | |
110 return quadraticRootsValidT(B, C, D, t); | |
111 } | |
112 double a, b, c; | |
113 { | |
114 double invA = 1 / A; | |
115 a = B * invA; | |
116 b = C * invA; | |
117 c = D * invA; | |
118 } | |
119 double a2 = a * a; | |
120 double Q = (a2 - b * 3) / 9; | |
121 double R = (2 * a2 * a - 9 * a * b + 27 * c) / 54; | |
122 double Q3 = Q * Q * Q; | |
123 double R2MinusQ3 = R * R - Q3; | |
124 double adiv3 = a / 3; | |
125 double* roots = t; | |
126 double r; | |
127 | |
128 if (R2MinusQ3 < 0) // we have 3 real roots | |
129 { | |
130 double theta = acos(R / sqrt(Q3)); | |
131 double neg2RootQ = -2 * sqrt(Q); | |
132 | |
133 r = neg2RootQ * cos(theta / 3) - adiv3; | |
134 if (is_unit_interval(r)) | |
135 *roots++ = r; | |
136 | |
137 r = neg2RootQ * cos((theta + 2 * PI) / 3) - adiv3; | |
138 if (is_unit_interval(r)) | |
139 *roots++ = r; | |
140 | |
141 r = neg2RootQ * cos((theta - 2 * PI) / 3) - adiv3; | |
142 if (is_unit_interval(r)) | |
143 *roots++ = r; | |
144 } | |
145 else // we have 1 real root | |
146 { | |
147 double A = fabs(R) + sqrt(R2MinusQ3); | |
148 A = cube_root(A); | |
149 if (R > 0) { | |
150 A = -A; | |
151 } | |
152 if (A != 0) { | |
153 A += Q / A; | |
154 } | |
155 r = A - adiv3; | |
156 if (is_unit_interval(r)) | |
157 *roots++ = r; | |
158 } | |
159 return (int)(roots - t); | |
160 #else | |
161 double s[3]; | |
162 int realRoots = cubicRootsReal(A, B, C, D, s); | |
163 int foundRoots = add_valid_ts(s, realRoots, t); | |
164 return foundRoots; | |
165 #endif | |
166 } | |
167 | |
168 int cubicRootsReal(double A, double B, double C, double D, double s[3]) { | |
169 #ifdef SK_DEBUG | |
170 // create a string mathematica understands | |
171 // GDB set print repe 15 # if repeated digits is a bother | |
172 // set print elements 400 # if line doesn't fit | |
173 char str[1024]; | |
174 bzero(str, sizeof(str)); | |
175 sprintf(str, "Solve[%1.19g x^3 + %1.19g x^2 + %1.19g x + %1.19g == 0, x]", A
, B, C, D); | |
176 mathematica_ize(str, sizeof(str)); | |
177 #if ONE_OFF_DEBUG && ONE_OFF_DEBUG_MATHEMATICA | |
178 SkDebugf("%s\n", str); | |
179 #endif | |
180 #endif | |
181 if (approximately_zero(A) | |
182 && approximately_zero_when_compared_to(A, B) | |
183 && approximately_zero_when_compared_to(A, C) | |
184 && approximately_zero_when_compared_to(A, D)) { // we're just a qua
dratic | |
185 return quadraticRootsReal(B, C, D, s); | |
186 } | |
187 if (approximately_zero_when_compared_to(D, A) | |
188 && approximately_zero_when_compared_to(D, B) | |
189 && approximately_zero_when_compared_to(D, C)) { // 0 is one root | |
190 int num = quadraticRootsReal(A, B, C, s); | |
191 for (int i = 0; i < num; ++i) { | |
192 if (approximately_zero(s[i])) { | |
193 return num; | |
194 } | |
195 } | |
196 s[num++] = 0; | |
197 return num; | |
198 } | |
199 if (approximately_zero(A + B + C + D)) { // 1 is one root | |
200 int num = quadraticRootsReal(A, A + B, -D, s); | |
201 for (int i = 0; i < num; ++i) { | |
202 if (AlmostEqualUlps(s[i], 1)) { | |
203 return num; | |
204 } | |
205 } | |
206 s[num++] = 1; | |
207 return num; | |
208 } | |
209 double a, b, c; | |
210 { | |
211 double invA = 1 / A; | |
212 a = B * invA; | |
213 b = C * invA; | |
214 c = D * invA; | |
215 } | |
216 double a2 = a * a; | |
217 double Q = (a2 - b * 3) / 9; | |
218 double R = (2 * a2 * a - 9 * a * b + 27 * c) / 54; | |
219 double R2 = R * R; | |
220 double Q3 = Q * Q * Q; | |
221 double R2MinusQ3 = R2 - Q3; | |
222 double adiv3 = a / 3; | |
223 double r; | |
224 double* roots = s; | |
225 #if 0 | |
226 if (approximately_zero_squared(R2MinusQ3) && AlmostEqualUlps(R2, Q3)) { | |
227 if (approximately_zero_squared(R)) {/* one triple solution */ | |
228 *roots++ = -adiv3; | |
229 } else { /* one single and one double solution */ | |
230 | |
231 double u = cube_root(-R); | |
232 *roots++ = 2 * u - adiv3; | |
233 *roots++ = -u - adiv3; | |
234 } | |
235 } | |
236 else | |
237 #endif | |
238 if (R2MinusQ3 < 0) // we have 3 real roots | |
239 { | |
240 double theta = acos(R / sqrt(Q3)); | |
241 double neg2RootQ = -2 * sqrt(Q); | |
242 | |
243 r = neg2RootQ * cos(theta / 3) - adiv3; | |
244 *roots++ = r; | |
245 | |
246 r = neg2RootQ * cos((theta + 2 * PI) / 3) - adiv3; | |
247 if (!AlmostEqualUlps(s[0], r)) { | |
248 *roots++ = r; | |
249 } | |
250 r = neg2RootQ * cos((theta - 2 * PI) / 3) - adiv3; | |
251 if (!AlmostEqualUlps(s[0], r) && (roots - s == 1 || !AlmostEqualUlps(s[1
], r))) { | |
252 *roots++ = r; | |
253 } | |
254 } | |
255 else // we have 1 real root | |
256 { | |
257 double sqrtR2MinusQ3 = sqrt(R2MinusQ3); | |
258 double A = fabs(R) + sqrtR2MinusQ3; | |
259 A = cube_root(A); | |
260 if (R > 0) { | |
261 A = -A; | |
262 } | |
263 if (A != 0) { | |
264 A += Q / A; | |
265 } | |
266 r = A - adiv3; | |
267 *roots++ = r; | |
268 if (AlmostEqualUlps(R2, Q3)) { | |
269 r = -A / 2 - adiv3; | |
270 if (!AlmostEqualUlps(s[0], r)) { | |
271 *roots++ = r; | |
272 } | |
273 } | |
274 } | |
275 return (int)(roots - s); | |
276 } | |
277 | |
278 // from http://www.cs.sunysb.edu/~qin/courses/geometry/4.pdf | |
279 // c(t) = a(1-t)^3 + 3bt(1-t)^2 + 3c(1-t)t^2 + dt^3 | |
280 // c'(t) = -3a(1-t)^2 + 3b((1-t)^2 - 2t(1-t)) + 3c(2t(1-t) - t^2) + 3dt^2 | |
281 // = 3(b-a)(1-t)^2 + 6(c-b)t(1-t) + 3(d-c)t^2 | |
282 static double derivativeAtT(const double* cubic, double t) { | |
283 double one_t = 1 - t; | |
284 double a = cubic[0]; | |
285 double b = cubic[2]; | |
286 double c = cubic[4]; | |
287 double d = cubic[6]; | |
288 return 3 * ((b - a) * one_t * one_t + 2 * (c - b) * t * one_t + (d - c) * t
* t); | |
289 } | |
290 | |
291 double dx_at_t(const Cubic& cubic, double t) { | |
292 return derivativeAtT(&cubic[0].x, t); | |
293 } | |
294 | |
295 double dy_at_t(const Cubic& cubic, double t) { | |
296 return derivativeAtT(&cubic[0].y, t); | |
297 } | |
298 | |
299 // OPTIMIZE? compute t^2, t(1-t), and (1-t)^2 and pass them to another version o
f derivative at t? | |
300 _Vector dxdy_at_t(const Cubic& cubic, double t) { | |
301 _Vector result = { derivativeAtT(&cubic[0].x, t), derivativeAtT(&cubic[0].y,
t) }; | |
302 return result; | |
303 } | |
304 | |
305 // OPTIMIZE? share code with formulate_F1DotF2 | |
306 int find_cubic_inflections(const Cubic& src, double tValues[]) | |
307 { | |
308 double Ax = src[1].x - src[0].x; | |
309 double Ay = src[1].y - src[0].y; | |
310 double Bx = src[2].x - 2 * src[1].x + src[0].x; | |
311 double By = src[2].y - 2 * src[1].y + src[0].y; | |
312 double Cx = src[3].x + 3 * (src[1].x - src[2].x) - src[0].x; | |
313 double Cy = src[3].y + 3 * (src[1].y - src[2].y) - src[0].y; | |
314 return quadraticRootsValidT(Bx * Cy - By * Cx, Ax * Cy - Ay * Cx, Ax * By -
Ay * Bx, tValues); | |
315 } | |
316 | |
317 static void formulate_F1DotF2(const double src[], double coeff[4]) | |
318 { | |
319 double a = src[2] - src[0]; | |
320 double b = src[4] - 2 * src[2] + src[0]; | |
321 double c = src[6] + 3 * (src[2] - src[4]) - src[0]; | |
322 coeff[0] = c * c; | |
323 coeff[1] = 3 * b * c; | |
324 coeff[2] = 2 * b * b + c * a; | |
325 coeff[3] = a * b; | |
326 } | |
327 | |
328 /* from SkGeometry.cpp | |
329 Looking for F' dot F'' == 0 | |
330 | |
331 A = b - a | |
332 B = c - 2b + a | |
333 C = d - 3c + 3b - a | |
334 | |
335 F' = 3Ct^2 + 6Bt + 3A | |
336 F'' = 6Ct + 6B | |
337 | |
338 F' dot F'' -> CCt^3 + 3BCt^2 + (2BB + CA)t + AB | |
339 */ | |
340 int find_cubic_max_curvature(const Cubic& src, double tValues[]) | |
341 { | |
342 double coeffX[4], coeffY[4]; | |
343 int i; | |
344 formulate_F1DotF2(&src[0].x, coeffX); | |
345 formulate_F1DotF2(&src[0].y, coeffY); | |
346 for (i = 0; i < 4; i++) { | |
347 coeffX[i] = coeffX[i] + coeffY[i]; | |
348 } | |
349 return cubicRootsValidT(coeffX[0], coeffX[1], coeffX[2], coeffX[3], tValues)
; | |
350 } | |
351 | |
352 | |
353 bool rotate(const Cubic& cubic, int zero, int index, Cubic& rotPath) { | |
354 double dy = cubic[index].y - cubic[zero].y; | |
355 double dx = cubic[index].x - cubic[zero].x; | |
356 if (approximately_zero(dy)) { | |
357 if (approximately_zero(dx)) { | |
358 return false; | |
359 } | |
360 memcpy(rotPath, cubic, sizeof(Cubic)); | |
361 return true; | |
362 } | |
363 for (int index = 0; index < 4; ++index) { | |
364 rotPath[index].x = cubic[index].x * dx + cubic[index].y * dy; | |
365 rotPath[index].y = cubic[index].y * dx - cubic[index].x * dy; | |
366 } | |
367 return true; | |
368 } | |
369 | |
370 #if 0 // unused for now | |
371 double secondDerivativeAtT(const double* cubic, double t) { | |
372 double a = cubic[0]; | |
373 double b = cubic[2]; | |
374 double c = cubic[4]; | |
375 double d = cubic[6]; | |
376 return (c - 2 * b + a) * (1 - t) + (d - 2 * c + b) * t; | |
377 } | |
378 #endif | |
379 | |
380 _Point top(const Cubic& cubic, double startT, double endT) { | |
381 Cubic sub; | |
382 sub_divide(cubic, startT, endT, sub); | |
383 _Point topPt = sub[0]; | |
384 if (topPt.y > sub[3].y || (topPt.y == sub[3].y && topPt.x > sub[3].x)) { | |
385 topPt = sub[3]; | |
386 } | |
387 double extremeTs[2]; | |
388 if (!monotonic_in_y(sub)) { | |
389 int roots = findExtrema(sub[0].y, sub[1].y, sub[2].y, sub[3].y, extremeT
s); | |
390 for (int index = 0; index < roots; ++index) { | |
391 _Point mid; | |
392 double t = startT + (endT - startT) * extremeTs[index]; | |
393 xy_at_t(cubic, t, mid.x, mid.y); | |
394 if (topPt.y > mid.y || (topPt.y == mid.y && topPt.x > mid.x)) { | |
395 topPt = mid; | |
396 } | |
397 } | |
398 } | |
399 return topPt; | |
400 } | |
401 | |
402 // OPTIMIZE: avoid computing the unused half | |
403 void xy_at_t(const Cubic& cubic, double t, double& x, double& y) { | |
404 _Point xy = xy_at_t(cubic, t); | |
405 if (&x) { | |
406 x = xy.x; | |
407 } | |
408 if (&y) { | |
409 y = xy.y; | |
410 } | |
411 } | |
412 | |
413 _Point xy_at_t(const Cubic& cubic, double t) { | |
414 double one_t = 1 - t; | |
415 double one_t2 = one_t * one_t; | |
416 double a = one_t2 * one_t; | |
417 double b = 3 * one_t2 * t; | |
418 double t2 = t * t; | |
419 double c = 3 * one_t * t2; | |
420 double d = t2 * t; | |
421 _Point result = {a * cubic[0].x + b * cubic[1].x + c * cubic[2].x + d * cubi
c[3].x, | |
422 a * cubic[0].y + b * cubic[1].y + c * cubic[2].y + d * cubic[3].y}; | |
423 return result; | |
424 } | |
OLD | NEW |