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? 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 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 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 of 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, 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 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 * 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
