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 values in test framework
13
14 // FIXME: cache keep the bounds and/or precision with the caller?
calcPrecision(const Cubic & cubic)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 #if SK_DEBUG
calcPrecision(const Cubic & cubic,double t,double scale)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
clockwise(const Cubic & c)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
coefficients(const double * cubic,double & A,double & B,double & C,double & D)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
controls_contained_by_ends(const Cubic & c)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
ends_are_extrema_in_x_or_y(const Cubic & c)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
monotonic_in_y(const Cubic & c)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
serpentine(const Cubic & c)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)
cubicRootsValidT(double A,double B,double C,double D,double t[3])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
cubicRootsReal(double A,double B,double C,double D,double s[3])168 int cubicRootsReal(double A, double B, double C, double D, double s[3]) {
169 #if 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 quadratic
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
derivativeAtT(const double * cubic,double t)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
dx_at_t(const Cubic & cubic,double t)291 double dx_at_t(const Cubic& cubic, double t) {
292 return derivativeAtT(&cubic[0].x, t);
293 }
294
dy_at_t(const Cubic & cubic,double t)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 of derivative at t?
dxdy_at_t(const Cubic & cubic,double 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
find_cubic_inflections(const Cubic & src,double tValues[])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
formulate_F1DotF2(const double src[],double coeff[4])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 */
find_cubic_max_curvature(const Cubic & src,double tValues[])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
rotate(const Cubic & cubic,int zero,int index,Cubic & rotPath)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
top(const Cubic & cubic,double startT,double endT)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, extremeTs);
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
xy_at_t(const Cubic & cubic,double t,double & x,double & y)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
xy_at_t(const Cubic & cubic,double t)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 * cubic[3].x,
422 a * cubic[0].y + b * cubic[1].y + c * cubic[2].y + d * cubic[3].y};
423 return result;
424 }
425