SkGeometry.cpp revision c9adb05b64fa0bfadf9d1a782afcda470da68c9e
1/* 2 * Copyright 2006 The Android Open Source Project 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#include "SkGeometry.h" 9#include "SkMatrix.h" 10#include "SkNx.h" 11 12static Sk2s from_point(const SkPoint& point) { 13 return Sk2s::Load(&point.fX); 14} 15 16static SkPoint to_point(const Sk2s& x) { 17 SkPoint point; 18 x.store(&point.fX); 19 return point; 20} 21 22static SkVector to_vector(const Sk2s& x) { 23 SkVector vector; 24 x.store(&vector.fX); 25 return vector; 26} 27 28/** If defined, this makes eval_quad and eval_cubic do more setup (sometimes 29 involving integer multiplies by 2 or 3, but fewer calls to SkScalarMul. 30 May also introduce overflow of fixed when we compute our setup. 31*/ 32// #define DIRECT_EVAL_OF_POLYNOMIALS 33 34//////////////////////////////////////////////////////////////////////// 35 36static int is_not_monotonic(SkScalar a, SkScalar b, SkScalar c) { 37 SkScalar ab = a - b; 38 SkScalar bc = b - c; 39 if (ab < 0) { 40 bc = -bc; 41 } 42 return ab == 0 || bc < 0; 43} 44 45//////////////////////////////////////////////////////////////////////// 46 47static bool is_unit_interval(SkScalar x) { 48 return x > 0 && x < SK_Scalar1; 49} 50 51static int valid_unit_divide(SkScalar numer, SkScalar denom, SkScalar* ratio) { 52 SkASSERT(ratio); 53 54 if (numer < 0) { 55 numer = -numer; 56 denom = -denom; 57 } 58 59 if (denom == 0 || numer == 0 || numer >= denom) { 60 return 0; 61 } 62 63 SkScalar r = SkScalarDiv(numer, denom); 64 if (SkScalarIsNaN(r)) { 65 return 0; 66 } 67 SkASSERTF(r >= 0 && r < SK_Scalar1, "numer %f, denom %f, r %f", numer, denom, r); 68 if (r == 0) { // catch underflow if numer <<<< denom 69 return 0; 70 } 71 *ratio = r; 72 return 1; 73} 74 75/** From Numerical Recipes in C. 76 77 Q = -1/2 (B + sign(B) sqrt[B*B - 4*A*C]) 78 x1 = Q / A 79 x2 = C / Q 80*/ 81int SkFindUnitQuadRoots(SkScalar A, SkScalar B, SkScalar C, SkScalar roots[2]) { 82 SkASSERT(roots); 83 84 if (A == 0) { 85 return valid_unit_divide(-C, B, roots); 86 } 87 88 SkScalar* r = roots; 89 90 SkScalar R = B*B - 4*A*C; 91 if (R < 0 || SkScalarIsNaN(R)) { // complex roots 92 return 0; 93 } 94 R = SkScalarSqrt(R); 95 96 SkScalar Q = (B < 0) ? -(B-R)/2 : -(B+R)/2; 97 r += valid_unit_divide(Q, A, r); 98 r += valid_unit_divide(C, Q, r); 99 if (r - roots == 2) { 100 if (roots[0] > roots[1]) 101 SkTSwap<SkScalar>(roots[0], roots[1]); 102 else if (roots[0] == roots[1]) // nearly-equal? 103 r -= 1; // skip the double root 104 } 105 return (int)(r - roots); 106} 107 108/////////////////////////////////////////////////////////////////////////////// 109/////////////////////////////////////////////////////////////////////////////// 110 111static Sk2s quad_poly_eval(const Sk2s& A, const Sk2s& B, const Sk2s& C, const Sk2s& t) { 112 return (A * t + B) * t + C; 113} 114 115static SkScalar eval_quad(const SkScalar src[], SkScalar t) { 116 SkASSERT(src); 117 SkASSERT(t >= 0 && t <= SK_Scalar1); 118 119#ifdef DIRECT_EVAL_OF_POLYNOMIALS 120 SkScalar C = src[0]; 121 SkScalar A = src[4] - 2 * src[2] + C; 122 SkScalar B = 2 * (src[2] - C); 123 return SkScalarMulAdd(SkScalarMulAdd(A, t, B), t, C); 124#else 125 SkScalar ab = SkScalarInterp(src[0], src[2], t); 126 SkScalar bc = SkScalarInterp(src[2], src[4], t); 127 return SkScalarInterp(ab, bc, t); 128#endif 129} 130 131static SkScalar eval_quad_derivative(const SkScalar src[], SkScalar t) { 132 SkScalar A = src[4] - 2 * src[2] + src[0]; 133 SkScalar B = src[2] - src[0]; 134 135 return 2 * SkScalarMulAdd(A, t, B); 136} 137 138void SkEvalQuadAt(const SkPoint src[3], SkScalar t, SkPoint* pt, SkVector* tangent) { 139 SkASSERT(src); 140 SkASSERT(t >= 0 && t <= SK_Scalar1); 141 142 if (pt) { 143 pt->set(eval_quad(&src[0].fX, t), eval_quad(&src[0].fY, t)); 144 } 145 if (tangent) { 146 tangent->set(eval_quad_derivative(&src[0].fX, t), 147 eval_quad_derivative(&src[0].fY, t)); 148 } 149} 150 151SkPoint SkEvalQuadAt(const SkPoint src[3], SkScalar t) { 152 SkASSERT(src); 153 SkASSERT(t >= 0 && t <= SK_Scalar1); 154 155 const Sk2s t2(t); 156 157 Sk2s P0 = from_point(src[0]); 158 Sk2s P1 = from_point(src[1]); 159 Sk2s P2 = from_point(src[2]); 160 161 Sk2s B = P1 - P0; 162 Sk2s A = P2 - P1 - B; 163 164 return to_point((A * t2 + B+B) * t2 + P0); 165} 166 167SkVector SkEvalQuadTangentAt(const SkPoint src[3], SkScalar t) { 168 SkASSERT(src); 169 SkASSERT(t >= 0 && t <= SK_Scalar1); 170 171 Sk2s P0 = from_point(src[0]); 172 Sk2s P1 = from_point(src[1]); 173 Sk2s P2 = from_point(src[2]); 174 175 Sk2s B = P1 - P0; 176 Sk2s A = P2 - P1 - B; 177 Sk2s T = A * Sk2s(t) + B; 178 179 return to_vector(T + T); 180} 181 182static inline Sk2s interp(const Sk2s& v0, const Sk2s& v1, const Sk2s& t) { 183 return v0 + (v1 - v0) * t; 184} 185 186void SkChopQuadAt(const SkPoint src[3], SkPoint dst[5], SkScalar t) { 187 SkASSERT(t > 0 && t < SK_Scalar1); 188 189 Sk2s p0 = from_point(src[0]); 190 Sk2s p1 = from_point(src[1]); 191 Sk2s p2 = from_point(src[2]); 192 Sk2s tt(t); 193 194 Sk2s p01 = interp(p0, p1, tt); 195 Sk2s p12 = interp(p1, p2, tt); 196 197 dst[0] = to_point(p0); 198 dst[1] = to_point(p01); 199 dst[2] = to_point(interp(p01, p12, tt)); 200 dst[3] = to_point(p12); 201 dst[4] = to_point(p2); 202} 203 204void SkChopQuadAtHalf(const SkPoint src[3], SkPoint dst[5]) { 205 SkChopQuadAt(src, dst, 0.5f); return; 206} 207 208/** Quad'(t) = At + B, where 209 A = 2(a - 2b + c) 210 B = 2(b - a) 211 Solve for t, only if it fits between 0 < t < 1 212*/ 213int SkFindQuadExtrema(SkScalar a, SkScalar b, SkScalar c, SkScalar tValue[1]) { 214 /* At + B == 0 215 t = -B / A 216 */ 217 return valid_unit_divide(a - b, a - b - b + c, tValue); 218} 219 220static inline void flatten_double_quad_extrema(SkScalar coords[14]) { 221 coords[2] = coords[6] = coords[4]; 222} 223 224/* Returns 0 for 1 quad, and 1 for two quads, either way the answer is 225 stored in dst[]. Guarantees that the 1/2 quads will be monotonic. 226 */ 227int SkChopQuadAtYExtrema(const SkPoint src[3], SkPoint dst[5]) { 228 SkASSERT(src); 229 SkASSERT(dst); 230 231 SkScalar a = src[0].fY; 232 SkScalar b = src[1].fY; 233 SkScalar c = src[2].fY; 234 235 if (is_not_monotonic(a, b, c)) { 236 SkScalar tValue; 237 if (valid_unit_divide(a - b, a - b - b + c, &tValue)) { 238 SkChopQuadAt(src, dst, tValue); 239 flatten_double_quad_extrema(&dst[0].fY); 240 return 1; 241 } 242 // if we get here, we need to force dst to be monotonic, even though 243 // we couldn't compute a unit_divide value (probably underflow). 244 b = SkScalarAbs(a - b) < SkScalarAbs(b - c) ? a : c; 245 } 246 dst[0].set(src[0].fX, a); 247 dst[1].set(src[1].fX, b); 248 dst[2].set(src[2].fX, c); 249 return 0; 250} 251 252/* Returns 0 for 1 quad, and 1 for two quads, either way the answer is 253 stored in dst[]. Guarantees that the 1/2 quads will be monotonic. 254 */ 255int SkChopQuadAtXExtrema(const SkPoint src[3], SkPoint dst[5]) { 256 SkASSERT(src); 257 SkASSERT(dst); 258 259 SkScalar a = src[0].fX; 260 SkScalar b = src[1].fX; 261 SkScalar c = src[2].fX; 262 263 if (is_not_monotonic(a, b, c)) { 264 SkScalar tValue; 265 if (valid_unit_divide(a - b, a - b - b + c, &tValue)) { 266 SkChopQuadAt(src, dst, tValue); 267 flatten_double_quad_extrema(&dst[0].fX); 268 return 1; 269 } 270 // if we get here, we need to force dst to be monotonic, even though 271 // we couldn't compute a unit_divide value (probably underflow). 272 b = SkScalarAbs(a - b) < SkScalarAbs(b - c) ? a : c; 273 } 274 dst[0].set(a, src[0].fY); 275 dst[1].set(b, src[1].fY); 276 dst[2].set(c, src[2].fY); 277 return 0; 278} 279 280// F(t) = a (1 - t) ^ 2 + 2 b t (1 - t) + c t ^ 2 281// F'(t) = 2 (b - a) + 2 (a - 2b + c) t 282// F''(t) = 2 (a - 2b + c) 283// 284// A = 2 (b - a) 285// B = 2 (a - 2b + c) 286// 287// Maximum curvature for a quadratic means solving 288// Fx' Fx'' + Fy' Fy'' = 0 289// 290// t = - (Ax Bx + Ay By) / (Bx ^ 2 + By ^ 2) 291// 292SkScalar SkFindQuadMaxCurvature(const SkPoint src[3]) { 293 SkScalar Ax = src[1].fX - src[0].fX; 294 SkScalar Ay = src[1].fY - src[0].fY; 295 SkScalar Bx = src[0].fX - src[1].fX - src[1].fX + src[2].fX; 296 SkScalar By = src[0].fY - src[1].fY - src[1].fY + src[2].fY; 297 SkScalar t = 0; // 0 means don't chop 298 299 (void)valid_unit_divide(-(Ax * Bx + Ay * By), Bx * Bx + By * By, &t); 300 return t; 301} 302 303int SkChopQuadAtMaxCurvature(const SkPoint src[3], SkPoint dst[5]) { 304 SkScalar t = SkFindQuadMaxCurvature(src); 305 if (t == 0) { 306 memcpy(dst, src, 3 * sizeof(SkPoint)); 307 return 1; 308 } else { 309 SkChopQuadAt(src, dst, t); 310 return 2; 311 } 312} 313 314void SkConvertQuadToCubic(const SkPoint src[3], SkPoint dst[4]) { 315 Sk2s scale(SkDoubleToScalar(2.0 / 3.0)); 316 Sk2s s0 = from_point(src[0]); 317 Sk2s s1 = from_point(src[1]); 318 Sk2s s2 = from_point(src[2]); 319 320 dst[0] = src[0]; 321 dst[1] = to_point(s0 + (s1 - s0) * scale); 322 dst[2] = to_point(s2 + (s1 - s2) * scale); 323 dst[3] = src[2]; 324} 325 326////////////////////////////////////////////////////////////////////////////// 327///// CUBICS // CUBICS // CUBICS // CUBICS // CUBICS // CUBICS // CUBICS ///// 328////////////////////////////////////////////////////////////////////////////// 329 330static void get_cubic_coeff(const SkScalar pt[], SkScalar coeff[4]) { 331 coeff[0] = pt[6] + 3*(pt[2] - pt[4]) - pt[0]; 332 coeff[1] = 3*(pt[4] - pt[2] - pt[2] + pt[0]); 333 coeff[2] = 3*(pt[2] - pt[0]); 334 coeff[3] = pt[0]; 335} 336 337void SkGetCubicCoeff(const SkPoint pts[4], SkScalar cx[4], SkScalar cy[4]) { 338 SkASSERT(pts); 339 340 if (cx) { 341 get_cubic_coeff(&pts[0].fX, cx); 342 } 343 if (cy) { 344 get_cubic_coeff(&pts[0].fY, cy); 345 } 346} 347 348static SkScalar eval_cubic(const SkScalar src[], SkScalar t) { 349 SkASSERT(src); 350 SkASSERT(t >= 0 && t <= SK_Scalar1); 351 352 if (t == 0) { 353 return src[0]; 354 } 355 356#ifdef DIRECT_EVAL_OF_POLYNOMIALS 357 SkScalar D = src[0]; 358 SkScalar A = src[6] + 3*(src[2] - src[4]) - D; 359 SkScalar B = 3*(src[4] - src[2] - src[2] + D); 360 SkScalar C = 3*(src[2] - D); 361 362 return SkScalarMulAdd(SkScalarMulAdd(SkScalarMulAdd(A, t, B), t, C), t, D); 363#else 364 SkScalar ab = SkScalarInterp(src[0], src[2], t); 365 SkScalar bc = SkScalarInterp(src[2], src[4], t); 366 SkScalar cd = SkScalarInterp(src[4], src[6], t); 367 SkScalar abc = SkScalarInterp(ab, bc, t); 368 SkScalar bcd = SkScalarInterp(bc, cd, t); 369 return SkScalarInterp(abc, bcd, t); 370#endif 371} 372 373/** return At^2 + Bt + C 374*/ 375static SkScalar eval_quadratic(SkScalar A, SkScalar B, SkScalar C, SkScalar t) { 376 SkASSERT(t >= 0 && t <= SK_Scalar1); 377 378 return SkScalarMulAdd(SkScalarMulAdd(A, t, B), t, C); 379} 380 381static SkScalar eval_cubic_derivative(const SkScalar src[], SkScalar t) { 382 SkScalar A = src[6] + 3*(src[2] - src[4]) - src[0]; 383 SkScalar B = 2*(src[4] - 2 * src[2] + src[0]); 384 SkScalar C = src[2] - src[0]; 385 386 return eval_quadratic(A, B, C, t); 387} 388 389static SkScalar eval_cubic_2ndDerivative(const SkScalar src[], SkScalar t) { 390 SkScalar A = src[6] + 3*(src[2] - src[4]) - src[0]; 391 SkScalar B = src[4] - 2 * src[2] + src[0]; 392 393 return SkScalarMulAdd(A, t, B); 394} 395 396void SkEvalCubicAt(const SkPoint src[4], SkScalar t, SkPoint* loc, 397 SkVector* tangent, SkVector* curvature) { 398 SkASSERT(src); 399 SkASSERT(t >= 0 && t <= SK_Scalar1); 400 401 if (loc) { 402 loc->set(eval_cubic(&src[0].fX, t), eval_cubic(&src[0].fY, t)); 403 } 404 if (tangent) { 405 tangent->set(eval_cubic_derivative(&src[0].fX, t), 406 eval_cubic_derivative(&src[0].fY, t)); 407 } 408 if (curvature) { 409 curvature->set(eval_cubic_2ndDerivative(&src[0].fX, t), 410 eval_cubic_2ndDerivative(&src[0].fY, t)); 411 } 412} 413 414/** Cubic'(t) = At^2 + Bt + C, where 415 A = 3(-a + 3(b - c) + d) 416 B = 6(a - 2b + c) 417 C = 3(b - a) 418 Solve for t, keeping only those that fit betwee 0 < t < 1 419*/ 420int SkFindCubicExtrema(SkScalar a, SkScalar b, SkScalar c, SkScalar d, 421 SkScalar tValues[2]) { 422 // we divide A,B,C by 3 to simplify 423 SkScalar A = d - a + 3*(b - c); 424 SkScalar B = 2*(a - b - b + c); 425 SkScalar C = b - a; 426 427 return SkFindUnitQuadRoots(A, B, C, tValues); 428} 429 430void SkChopCubicAt(const SkPoint src[4], SkPoint dst[7], SkScalar t) { 431 SkASSERT(t > 0 && t < SK_Scalar1); 432 433 Sk2s p0 = from_point(src[0]); 434 Sk2s p1 = from_point(src[1]); 435 Sk2s p2 = from_point(src[2]); 436 Sk2s p3 = from_point(src[3]); 437 Sk2s tt(t); 438 439 Sk2s ab = interp(p0, p1, tt); 440 Sk2s bc = interp(p1, p2, tt); 441 Sk2s cd = interp(p2, p3, tt); 442 Sk2s abc = interp(ab, bc, tt); 443 Sk2s bcd = interp(bc, cd, tt); 444 Sk2s abcd = interp(abc, bcd, tt); 445 446 dst[0] = src[0]; 447 dst[1] = to_point(ab); 448 dst[2] = to_point(abc); 449 dst[3] = to_point(abcd); 450 dst[4] = to_point(bcd); 451 dst[5] = to_point(cd); 452 dst[6] = src[3]; 453} 454 455/* http://code.google.com/p/skia/issues/detail?id=32 456 457 This test code would fail when we didn't check the return result of 458 valid_unit_divide in SkChopCubicAt(... tValues[], int roots). The reason is 459 that after the first chop, the parameters to valid_unit_divide are equal 460 (thanks to finite float precision and rounding in the subtracts). Thus 461 even though the 2nd tValue looks < 1.0, after we renormalize it, we end 462 up with 1.0, hence the need to check and just return the last cubic as 463 a degenerate clump of 4 points in the sampe place. 464 465 static void test_cubic() { 466 SkPoint src[4] = { 467 { 556.25000, 523.03003 }, 468 { 556.23999, 522.96002 }, 469 { 556.21997, 522.89001 }, 470 { 556.21997, 522.82001 } 471 }; 472 SkPoint dst[10]; 473 SkScalar tval[] = { 0.33333334f, 0.99999994f }; 474 SkChopCubicAt(src, dst, tval, 2); 475 } 476 */ 477 478void SkChopCubicAt(const SkPoint src[4], SkPoint dst[], 479 const SkScalar tValues[], int roots) { 480#ifdef SK_DEBUG 481 { 482 for (int i = 0; i < roots - 1; i++) 483 { 484 SkASSERT(is_unit_interval(tValues[i])); 485 SkASSERT(is_unit_interval(tValues[i+1])); 486 SkASSERT(tValues[i] < tValues[i+1]); 487 } 488 } 489#endif 490 491 if (dst) { 492 if (roots == 0) { // nothing to chop 493 memcpy(dst, src, 4*sizeof(SkPoint)); 494 } else { 495 SkScalar t = tValues[0]; 496 SkPoint tmp[4]; 497 498 for (int i = 0; i < roots; i++) { 499 SkChopCubicAt(src, dst, t); 500 if (i == roots - 1) { 501 break; 502 } 503 504 dst += 3; 505 // have src point to the remaining cubic (after the chop) 506 memcpy(tmp, dst, 4 * sizeof(SkPoint)); 507 src = tmp; 508 509 // watch out in case the renormalized t isn't in range 510 if (!valid_unit_divide(tValues[i+1] - tValues[i], 511 SK_Scalar1 - tValues[i], &t)) { 512 // if we can't, just create a degenerate cubic 513 dst[4] = dst[5] = dst[6] = src[3]; 514 break; 515 } 516 } 517 } 518 } 519} 520 521void SkChopCubicAtHalf(const SkPoint src[4], SkPoint dst[7]) { 522 SkChopCubicAt(src, dst, 0.5f); 523} 524 525static void flatten_double_cubic_extrema(SkScalar coords[14]) { 526 coords[4] = coords[8] = coords[6]; 527} 528 529/** Given 4 points on a cubic bezier, chop it into 1, 2, 3 beziers such that 530 the resulting beziers are monotonic in Y. This is called by the scan 531 converter. Depending on what is returned, dst[] is treated as follows: 532 0 dst[0..3] is the original cubic 533 1 dst[0..3] and dst[3..6] are the two new cubics 534 2 dst[0..3], dst[3..6], dst[6..9] are the three new cubics 535 If dst == null, it is ignored and only the count is returned. 536*/ 537int SkChopCubicAtYExtrema(const SkPoint src[4], SkPoint dst[10]) { 538 SkScalar tValues[2]; 539 int roots = SkFindCubicExtrema(src[0].fY, src[1].fY, src[2].fY, 540 src[3].fY, tValues); 541 542 SkChopCubicAt(src, dst, tValues, roots); 543 if (dst && roots > 0) { 544 // we do some cleanup to ensure our Y extrema are flat 545 flatten_double_cubic_extrema(&dst[0].fY); 546 if (roots == 2) { 547 flatten_double_cubic_extrema(&dst[3].fY); 548 } 549 } 550 return roots; 551} 552 553int SkChopCubicAtXExtrema(const SkPoint src[4], SkPoint dst[10]) { 554 SkScalar tValues[2]; 555 int roots = SkFindCubicExtrema(src[0].fX, src[1].fX, src[2].fX, 556 src[3].fX, tValues); 557 558 SkChopCubicAt(src, dst, tValues, roots); 559 if (dst && roots > 0) { 560 // we do some cleanup to ensure our Y extrema are flat 561 flatten_double_cubic_extrema(&dst[0].fX); 562 if (roots == 2) { 563 flatten_double_cubic_extrema(&dst[3].fX); 564 } 565 } 566 return roots; 567} 568 569/** http://www.faculty.idc.ac.il/arik/quality/appendixA.html 570 571 Inflection means that curvature is zero. 572 Curvature is [F' x F''] / [F'^3] 573 So we solve F'x X F''y - F'y X F''y == 0 574 After some canceling of the cubic term, we get 575 A = b - a 576 B = c - 2b + a 577 C = d - 3c + 3b - a 578 (BxCy - ByCx)t^2 + (AxCy - AyCx)t + AxBy - AyBx == 0 579*/ 580int SkFindCubicInflections(const SkPoint src[4], SkScalar tValues[]) { 581 SkScalar Ax = src[1].fX - src[0].fX; 582 SkScalar Ay = src[1].fY - src[0].fY; 583 SkScalar Bx = src[2].fX - 2 * src[1].fX + src[0].fX; 584 SkScalar By = src[2].fY - 2 * src[1].fY + src[0].fY; 585 SkScalar Cx = src[3].fX + 3 * (src[1].fX - src[2].fX) - src[0].fX; 586 SkScalar Cy = src[3].fY + 3 * (src[1].fY - src[2].fY) - src[0].fY; 587 588 return SkFindUnitQuadRoots(Bx*Cy - By*Cx, 589 Ax*Cy - Ay*Cx, 590 Ax*By - Ay*Bx, 591 tValues); 592} 593 594int SkChopCubicAtInflections(const SkPoint src[], SkPoint dst[10]) { 595 SkScalar tValues[2]; 596 int count = SkFindCubicInflections(src, tValues); 597 598 if (dst) { 599 if (count == 0) { 600 memcpy(dst, src, 4 * sizeof(SkPoint)); 601 } else { 602 SkChopCubicAt(src, dst, tValues, count); 603 } 604 } 605 return count + 1; 606} 607 608// See http://http.developer.nvidia.com/GPUGems3/gpugems3_ch25.html (from the book GPU Gems 3) 609// discr(I) = d0^2 * (3*d1^2 - 4*d0*d2) 610// Classification: 611// discr(I) > 0 Serpentine 612// discr(I) = 0 Cusp 613// discr(I) < 0 Loop 614// d0 = d1 = 0 Quadratic 615// d0 = d1 = d2 = 0 Line 616// p0 = p1 = p2 = p3 Point 617static SkCubicType classify_cubic(const SkPoint p[4], const SkScalar d[3]) { 618 if (p[0] == p[1] && p[0] == p[2] && p[0] == p[3]) { 619 return kPoint_SkCubicType; 620 } 621 const SkScalar discr = d[0] * d[0] * (3.f * d[1] * d[1] - 4.f * d[0] * d[2]); 622 if (discr > SK_ScalarNearlyZero) { 623 return kSerpentine_SkCubicType; 624 } else if (discr < -SK_ScalarNearlyZero) { 625 return kLoop_SkCubicType; 626 } else { 627 if (0.f == d[0] && 0.f == d[1]) { 628 return (0.f == d[2] ? kLine_SkCubicType : kQuadratic_SkCubicType); 629 } else { 630 return kCusp_SkCubicType; 631 } 632 } 633} 634 635// Assumes the third component of points is 1. 636// Calcs p0 . (p1 x p2) 637static SkScalar calc_dot_cross_cubic(const SkPoint& p0, const SkPoint& p1, const SkPoint& p2) { 638 const SkScalar xComp = p0.fX * (p1.fY - p2.fY); 639 const SkScalar yComp = p0.fY * (p2.fX - p1.fX); 640 const SkScalar wComp = p1.fX * p2.fY - p1.fY * p2.fX; 641 return (xComp + yComp + wComp); 642} 643 644// Calc coefficients of I(s,t) where roots of I are inflection points of curve 645// I(s,t) = t*(3*d0*s^2 - 3*d1*s*t + d2*t^2) 646// d0 = a1 - 2*a2+3*a3 647// d1 = -a2 + 3*a3 648// d2 = 3*a3 649// a1 = p0 . (p3 x p2) 650// a2 = p1 . (p0 x p3) 651// a3 = p2 . (p1 x p0) 652// Places the values of d1, d2, d3 in array d passed in 653static void calc_cubic_inflection_func(const SkPoint p[4], SkScalar d[3]) { 654 SkScalar a1 = calc_dot_cross_cubic(p[0], p[3], p[2]); 655 SkScalar a2 = calc_dot_cross_cubic(p[1], p[0], p[3]); 656 SkScalar a3 = calc_dot_cross_cubic(p[2], p[1], p[0]); 657 658 // need to scale a's or values in later calculations will grow to high 659 SkScalar max = SkScalarAbs(a1); 660 max = SkMaxScalar(max, SkScalarAbs(a2)); 661 max = SkMaxScalar(max, SkScalarAbs(a3)); 662 max = 1.f/max; 663 a1 = a1 * max; 664 a2 = a2 * max; 665 a3 = a3 * max; 666 667 d[2] = 3.f * a3; 668 d[1] = d[2] - a2; 669 d[0] = d[1] - a2 + a1; 670} 671 672SkCubicType SkClassifyCubic(const SkPoint src[4], SkScalar d[3]) { 673 calc_cubic_inflection_func(src, d); 674 return classify_cubic(src, d); 675} 676 677template <typename T> void bubble_sort(T array[], int count) { 678 for (int i = count - 1; i > 0; --i) 679 for (int j = i; j > 0; --j) 680 if (array[j] < array[j-1]) 681 { 682 T tmp(array[j]); 683 array[j] = array[j-1]; 684 array[j-1] = tmp; 685 } 686} 687 688/** 689 * Given an array and count, remove all pair-wise duplicates from the array, 690 * keeping the existing sorting, and return the new count 691 */ 692static int collaps_duplicates(SkScalar array[], int count) { 693 for (int n = count; n > 1; --n) { 694 if (array[0] == array[1]) { 695 for (int i = 1; i < n; ++i) { 696 array[i - 1] = array[i]; 697 } 698 count -= 1; 699 } else { 700 array += 1; 701 } 702 } 703 return count; 704} 705 706#ifdef SK_DEBUG 707 708#define TEST_COLLAPS_ENTRY(array) array, SK_ARRAY_COUNT(array) 709 710static void test_collaps_duplicates() { 711 static bool gOnce; 712 if (gOnce) { return; } 713 gOnce = true; 714 const SkScalar src0[] = { 0 }; 715 const SkScalar src1[] = { 0, 0 }; 716 const SkScalar src2[] = { 0, 1 }; 717 const SkScalar src3[] = { 0, 0, 0 }; 718 const SkScalar src4[] = { 0, 0, 1 }; 719 const SkScalar src5[] = { 0, 1, 1 }; 720 const SkScalar src6[] = { 0, 1, 2 }; 721 const struct { 722 const SkScalar* fData; 723 int fCount; 724 int fCollapsedCount; 725 } data[] = { 726 { TEST_COLLAPS_ENTRY(src0), 1 }, 727 { TEST_COLLAPS_ENTRY(src1), 1 }, 728 { TEST_COLLAPS_ENTRY(src2), 2 }, 729 { TEST_COLLAPS_ENTRY(src3), 1 }, 730 { TEST_COLLAPS_ENTRY(src4), 2 }, 731 { TEST_COLLAPS_ENTRY(src5), 2 }, 732 { TEST_COLLAPS_ENTRY(src6), 3 }, 733 }; 734 for (size_t i = 0; i < SK_ARRAY_COUNT(data); ++i) { 735 SkScalar dst[3]; 736 memcpy(dst, data[i].fData, data[i].fCount * sizeof(dst[0])); 737 int count = collaps_duplicates(dst, data[i].fCount); 738 SkASSERT(data[i].fCollapsedCount == count); 739 for (int j = 1; j < count; ++j) { 740 SkASSERT(dst[j-1] < dst[j]); 741 } 742 } 743} 744#endif 745 746static SkScalar SkScalarCubeRoot(SkScalar x) { 747 return SkScalarPow(x, 0.3333333f); 748} 749 750/* Solve coeff(t) == 0, returning the number of roots that 751 lie withing 0 < t < 1. 752 coeff[0]t^3 + coeff[1]t^2 + coeff[2]t + coeff[3] 753 754 Eliminates repeated roots (so that all tValues are distinct, and are always 755 in increasing order. 756*/ 757static int solve_cubic_poly(const SkScalar coeff[4], SkScalar tValues[3]) { 758 if (SkScalarNearlyZero(coeff[0])) { // we're just a quadratic 759 return SkFindUnitQuadRoots(coeff[1], coeff[2], coeff[3], tValues); 760 } 761 762 SkScalar a, b, c, Q, R; 763 764 { 765 SkASSERT(coeff[0] != 0); 766 767 SkScalar inva = SkScalarInvert(coeff[0]); 768 a = coeff[1] * inva; 769 b = coeff[2] * inva; 770 c = coeff[3] * inva; 771 } 772 Q = (a*a - b*3) / 9; 773 R = (2*a*a*a - 9*a*b + 27*c) / 54; 774 775 SkScalar Q3 = Q * Q * Q; 776 SkScalar R2MinusQ3 = R * R - Q3; 777 SkScalar adiv3 = a / 3; 778 779 SkScalar* roots = tValues; 780 SkScalar r; 781 782 if (R2MinusQ3 < 0) { // we have 3 real roots 783 SkScalar theta = SkScalarACos(R / SkScalarSqrt(Q3)); 784 SkScalar neg2RootQ = -2 * SkScalarSqrt(Q); 785 786 r = neg2RootQ * SkScalarCos(theta/3) - adiv3; 787 if (is_unit_interval(r)) { 788 *roots++ = r; 789 } 790 r = neg2RootQ * SkScalarCos((theta + 2*SK_ScalarPI)/3) - adiv3; 791 if (is_unit_interval(r)) { 792 *roots++ = r; 793 } 794 r = neg2RootQ * SkScalarCos((theta - 2*SK_ScalarPI)/3) - adiv3; 795 if (is_unit_interval(r)) { 796 *roots++ = r; 797 } 798 SkDEBUGCODE(test_collaps_duplicates();) 799 800 // now sort the roots 801 int count = (int)(roots - tValues); 802 SkASSERT((unsigned)count <= 3); 803 bubble_sort(tValues, count); 804 count = collaps_duplicates(tValues, count); 805 roots = tValues + count; // so we compute the proper count below 806 } else { // we have 1 real root 807 SkScalar A = SkScalarAbs(R) + SkScalarSqrt(R2MinusQ3); 808 A = SkScalarCubeRoot(A); 809 if (R > 0) { 810 A = -A; 811 } 812 if (A != 0) { 813 A += Q / A; 814 } 815 r = A - adiv3; 816 if (is_unit_interval(r)) { 817 *roots++ = r; 818 } 819 } 820 821 return (int)(roots - tValues); 822} 823 824/* Looking for F' dot F'' == 0 825 826 A = b - a 827 B = c - 2b + a 828 C = d - 3c + 3b - a 829 830 F' = 3Ct^2 + 6Bt + 3A 831 F'' = 6Ct + 6B 832 833 F' dot F'' -> CCt^3 + 3BCt^2 + (2BB + CA)t + AB 834*/ 835static void formulate_F1DotF2(const SkScalar src[], SkScalar coeff[4]) { 836 SkScalar a = src[2] - src[0]; 837 SkScalar b = src[4] - 2 * src[2] + src[0]; 838 SkScalar c = src[6] + 3 * (src[2] - src[4]) - src[0]; 839 840 coeff[0] = c * c; 841 coeff[1] = 3 * b * c; 842 coeff[2] = 2 * b * b + c * a; 843 coeff[3] = a * b; 844} 845 846/* Looking for F' dot F'' == 0 847 848 A = b - a 849 B = c - 2b + a 850 C = d - 3c + 3b - a 851 852 F' = 3Ct^2 + 6Bt + 3A 853 F'' = 6Ct + 6B 854 855 F' dot F'' -> CCt^3 + 3BCt^2 + (2BB + CA)t + AB 856*/ 857int SkFindCubicMaxCurvature(const SkPoint src[4], SkScalar tValues[3]) { 858 SkScalar coeffX[4], coeffY[4]; 859 int i; 860 861 formulate_F1DotF2(&src[0].fX, coeffX); 862 formulate_F1DotF2(&src[0].fY, coeffY); 863 864 for (i = 0; i < 4; i++) { 865 coeffX[i] += coeffY[i]; 866 } 867 868 SkScalar t[3]; 869 int count = solve_cubic_poly(coeffX, t); 870 int maxCount = 0; 871 872 // now remove extrema where the curvature is zero (mins) 873 // !!!! need a test for this !!!! 874 for (i = 0; i < count; i++) { 875 // if (not_min_curvature()) 876 if (t[i] > 0 && t[i] < SK_Scalar1) { 877 tValues[maxCount++] = t[i]; 878 } 879 } 880 return maxCount; 881} 882 883int SkChopCubicAtMaxCurvature(const SkPoint src[4], SkPoint dst[13], 884 SkScalar tValues[3]) { 885 SkScalar t_storage[3]; 886 887 if (tValues == NULL) { 888 tValues = t_storage; 889 } 890 891 int count = SkFindCubicMaxCurvature(src, tValues); 892 893 if (dst) { 894 if (count == 0) { 895 memcpy(dst, src, 4 * sizeof(SkPoint)); 896 } else { 897 SkChopCubicAt(src, dst, tValues, count); 898 } 899 } 900 return count + 1; 901} 902 903/////////////////////////////////////////////////////////////////////////////// 904 905/* Find t value for quadratic [a, b, c] = d. 906 Return 0 if there is no solution within [0, 1) 907*/ 908static SkScalar quad_solve(SkScalar a, SkScalar b, SkScalar c, SkScalar d) { 909 // At^2 + Bt + C = d 910 SkScalar A = a - 2 * b + c; 911 SkScalar B = 2 * (b - a); 912 SkScalar C = a - d; 913 914 SkScalar roots[2]; 915 int count = SkFindUnitQuadRoots(A, B, C, roots); 916 917 SkASSERT(count <= 1); 918 return count == 1 ? roots[0] : 0; 919} 920 921/* given a quad-curve and a point (x,y), chop the quad at that point and place 922 the new off-curve point and endpoint into 'dest'. 923 Should only return false if the computed pos is the start of the curve 924 (i.e. root == 0) 925*/ 926static bool truncate_last_curve(const SkPoint quad[3], SkScalar x, SkScalar y, 927 SkPoint* dest) { 928 const SkScalar* base; 929 SkScalar value; 930 931 if (SkScalarAbs(x) < SkScalarAbs(y)) { 932 base = &quad[0].fX; 933 value = x; 934 } else { 935 base = &quad[0].fY; 936 value = y; 937 } 938 939 // note: this returns 0 if it thinks value is out of range, meaning the 940 // root might return something outside of [0, 1) 941 SkScalar t = quad_solve(base[0], base[2], base[4], value); 942 943 if (t > 0) { 944 SkPoint tmp[5]; 945 SkChopQuadAt(quad, tmp, t); 946 dest[0] = tmp[1]; 947 dest[1].set(x, y); 948 return true; 949 } else { 950 /* t == 0 means either the value triggered a root outside of [0, 1) 951 For our purposes, we can ignore the <= 0 roots, but we want to 952 catch the >= 1 roots (which given our caller, will basically mean 953 a root of 1, give-or-take numerical instability). If we are in the 954 >= 1 case, return the existing offCurve point. 955 956 The test below checks to see if we are close to the "end" of the 957 curve (near base[4]). Rather than specifying a tolerance, I just 958 check to see if value is on to the right/left of the middle point 959 (depending on the direction/sign of the end points). 960 */ 961 if ((base[0] < base[4] && value > base[2]) || 962 (base[0] > base[4] && value < base[2])) // should root have been 1 963 { 964 dest[0] = quad[1]; 965 dest[1].set(x, y); 966 return true; 967 } 968 } 969 return false; 970} 971 972static const SkPoint gQuadCirclePts[kSkBuildQuadArcStorage] = { 973// The mid point of the quadratic arc approximation is half way between the two 974// control points. The float epsilon adjustment moves the on curve point out by 975// two bits, distributing the convex test error between the round rect 976// approximation and the convex cross product sign equality test. 977#define SK_MID_RRECT_OFFSET \ 978 (SK_Scalar1 + SK_ScalarTanPIOver8 + FLT_EPSILON * 4) / 2 979 { SK_Scalar1, 0 }, 980 { SK_Scalar1, SK_ScalarTanPIOver8 }, 981 { SK_MID_RRECT_OFFSET, SK_MID_RRECT_OFFSET }, 982 { SK_ScalarTanPIOver8, SK_Scalar1 }, 983 984 { 0, SK_Scalar1 }, 985 { -SK_ScalarTanPIOver8, SK_Scalar1 }, 986 { -SK_MID_RRECT_OFFSET, SK_MID_RRECT_OFFSET }, 987 { -SK_Scalar1, SK_ScalarTanPIOver8 }, 988 989 { -SK_Scalar1, 0 }, 990 { -SK_Scalar1, -SK_ScalarTanPIOver8 }, 991 { -SK_MID_RRECT_OFFSET, -SK_MID_RRECT_OFFSET }, 992 { -SK_ScalarTanPIOver8, -SK_Scalar1 }, 993 994 { 0, -SK_Scalar1 }, 995 { SK_ScalarTanPIOver8, -SK_Scalar1 }, 996 { SK_MID_RRECT_OFFSET, -SK_MID_RRECT_OFFSET }, 997 { SK_Scalar1, -SK_ScalarTanPIOver8 }, 998 999 { SK_Scalar1, 0 } 1000#undef SK_MID_RRECT_OFFSET 1001}; 1002 1003int SkBuildQuadArc(const SkVector& uStart, const SkVector& uStop, 1004 SkRotationDirection dir, const SkMatrix* userMatrix, 1005 SkPoint quadPoints[]) { 1006 // rotate by x,y so that uStart is (1.0) 1007 SkScalar x = SkPoint::DotProduct(uStart, uStop); 1008 SkScalar y = SkPoint::CrossProduct(uStart, uStop); 1009 1010 SkScalar absX = SkScalarAbs(x); 1011 SkScalar absY = SkScalarAbs(y); 1012 1013 int pointCount; 1014 1015 // check for (effectively) coincident vectors 1016 // this can happen if our angle is nearly 0 or nearly 180 (y == 0) 1017 // ... we use the dot-prod to distinguish between 0 and 180 (x > 0) 1018 if (absY <= SK_ScalarNearlyZero && x > 0 && 1019 ((y >= 0 && kCW_SkRotationDirection == dir) || 1020 (y <= 0 && kCCW_SkRotationDirection == dir))) { 1021 1022 // just return the start-point 1023 quadPoints[0].set(SK_Scalar1, 0); 1024 pointCount = 1; 1025 } else { 1026 if (dir == kCCW_SkRotationDirection) { 1027 y = -y; 1028 } 1029 // what octant (quadratic curve) is [xy] in? 1030 int oct = 0; 1031 bool sameSign = true; 1032 1033 if (0 == y) { 1034 oct = 4; // 180 1035 SkASSERT(SkScalarAbs(x + SK_Scalar1) <= SK_ScalarNearlyZero); 1036 } else if (0 == x) { 1037 SkASSERT(absY - SK_Scalar1 <= SK_ScalarNearlyZero); 1038 oct = y > 0 ? 2 : 6; // 90 : 270 1039 } else { 1040 if (y < 0) { 1041 oct += 4; 1042 } 1043 if ((x < 0) != (y < 0)) { 1044 oct += 2; 1045 sameSign = false; 1046 } 1047 if ((absX < absY) == sameSign) { 1048 oct += 1; 1049 } 1050 } 1051 1052 int wholeCount = oct << 1; 1053 memcpy(quadPoints, gQuadCirclePts, (wholeCount + 1) * sizeof(SkPoint)); 1054 1055 const SkPoint* arc = &gQuadCirclePts[wholeCount]; 1056 if (truncate_last_curve(arc, x, y, &quadPoints[wholeCount + 1])) { 1057 wholeCount += 2; 1058 } 1059 pointCount = wholeCount + 1; 1060 } 1061 1062 // now handle counter-clockwise and the initial unitStart rotation 1063 SkMatrix matrix; 1064 matrix.setSinCos(uStart.fY, uStart.fX); 1065 if (dir == kCCW_SkRotationDirection) { 1066 matrix.preScale(SK_Scalar1, -SK_Scalar1); 1067 } 1068 if (userMatrix) { 1069 matrix.postConcat(*userMatrix); 1070 } 1071 matrix.mapPoints(quadPoints, pointCount); 1072 return pointCount; 1073} 1074 1075 1076/////////////////////////////////////////////////////////////////////////////// 1077// 1078// NURB representation for conics. Helpful explanations at: 1079// 1080// http://citeseerx.ist.psu.edu/viewdoc/ 1081// download?doi=10.1.1.44.5740&rep=rep1&type=ps 1082// and 1083// http://www.cs.mtu.edu/~shene/COURSES/cs3621/NOTES/spline/NURBS/RB-conics.html 1084// 1085// F = (A (1 - t)^2 + C t^2 + 2 B (1 - t) t w) 1086// ------------------------------------------ 1087// ((1 - t)^2 + t^2 + 2 (1 - t) t w) 1088// 1089// = {t^2 (P0 + P2 - 2 P1 w), t (-2 P0 + 2 P1 w), P0} 1090// ------------------------------------------------ 1091// {t^2 (2 - 2 w), t (-2 + 2 w), 1} 1092// 1093 1094static SkScalar conic_eval_pos(const SkScalar src[], SkScalar w, SkScalar t) { 1095 SkASSERT(src); 1096 SkASSERT(t >= 0 && t <= SK_Scalar1); 1097 1098 SkScalar src2w = SkScalarMul(src[2], w); 1099 SkScalar C = src[0]; 1100 SkScalar A = src[4] - 2 * src2w + C; 1101 SkScalar B = 2 * (src2w - C); 1102 SkScalar numer = SkScalarMulAdd(SkScalarMulAdd(A, t, B), t, C); 1103 1104 B = 2 * (w - SK_Scalar1); 1105 C = SK_Scalar1; 1106 A = -B; 1107 SkScalar denom = SkScalarMulAdd(SkScalarMulAdd(A, t, B), t, C); 1108 1109 return SkScalarDiv(numer, denom); 1110} 1111 1112// F' = 2 (C t (1 + t (-1 + w)) - A (-1 + t) (t (-1 + w) - w) + B (1 - 2 t) w) 1113// 1114// t^2 : (2 P0 - 2 P2 - 2 P0 w + 2 P2 w) 1115// t^1 : (-2 P0 + 2 P2 + 4 P0 w - 4 P1 w) 1116// t^0 : -2 P0 w + 2 P1 w 1117// 1118// We disregard magnitude, so we can freely ignore the denominator of F', and 1119// divide the numerator by 2 1120// 1121// coeff[0] for t^2 1122// coeff[1] for t^1 1123// coeff[2] for t^0 1124// 1125static void conic_deriv_coeff(const SkScalar src[], 1126 SkScalar w, 1127 SkScalar coeff[3]) { 1128 const SkScalar P20 = src[4] - src[0]; 1129 const SkScalar P10 = src[2] - src[0]; 1130 const SkScalar wP10 = w * P10; 1131 coeff[0] = w * P20 - P20; 1132 coeff[1] = P20 - 2 * wP10; 1133 coeff[2] = wP10; 1134} 1135 1136static SkScalar conic_eval_tan(const SkScalar coord[], SkScalar w, SkScalar t) { 1137 SkScalar coeff[3]; 1138 conic_deriv_coeff(coord, w, coeff); 1139 return t * (t * coeff[0] + coeff[1]) + coeff[2]; 1140} 1141 1142static bool conic_find_extrema(const SkScalar src[], SkScalar w, SkScalar* t) { 1143 SkScalar coeff[3]; 1144 conic_deriv_coeff(src, w, coeff); 1145 1146 SkScalar tValues[2]; 1147 int roots = SkFindUnitQuadRoots(coeff[0], coeff[1], coeff[2], tValues); 1148 SkASSERT(0 == roots || 1 == roots); 1149 1150 if (1 == roots) { 1151 *t = tValues[0]; 1152 return true; 1153 } 1154 return false; 1155} 1156 1157struct SkP3D { 1158 SkScalar fX, fY, fZ; 1159 1160 void set(SkScalar x, SkScalar y, SkScalar z) { 1161 fX = x; fY = y; fZ = z; 1162 } 1163 1164 void projectDown(SkPoint* dst) const { 1165 dst->set(fX / fZ, fY / fZ); 1166 } 1167}; 1168 1169// We only interpolate one dimension at a time (the first, at +0, +3, +6). 1170static void p3d_interp(const SkScalar src[7], SkScalar dst[7], SkScalar t) { 1171 SkScalar ab = SkScalarInterp(src[0], src[3], t); 1172 SkScalar bc = SkScalarInterp(src[3], src[6], t); 1173 dst[0] = ab; 1174 dst[3] = SkScalarInterp(ab, bc, t); 1175 dst[6] = bc; 1176} 1177 1178static void ratquad_mapTo3D(const SkPoint src[3], SkScalar w, SkP3D dst[]) { 1179 dst[0].set(src[0].fX * 1, src[0].fY * 1, 1); 1180 dst[1].set(src[1].fX * w, src[1].fY * w, w); 1181 dst[2].set(src[2].fX * 1, src[2].fY * 1, 1); 1182} 1183 1184void SkConic::evalAt(SkScalar t, SkPoint* pt, SkVector* tangent) const { 1185 SkASSERT(t >= 0 && t <= SK_Scalar1); 1186 1187 if (pt) { 1188 pt->set(conic_eval_pos(&fPts[0].fX, fW, t), 1189 conic_eval_pos(&fPts[0].fY, fW, t)); 1190 } 1191 if (tangent) { 1192 tangent->set(conic_eval_tan(&fPts[0].fX, fW, t), 1193 conic_eval_tan(&fPts[0].fY, fW, t)); 1194 } 1195} 1196 1197void SkConic::chopAt(SkScalar t, SkConic dst[2]) const { 1198 SkP3D tmp[3], tmp2[3]; 1199 1200 ratquad_mapTo3D(fPts, fW, tmp); 1201 1202 p3d_interp(&tmp[0].fX, &tmp2[0].fX, t); 1203 p3d_interp(&tmp[0].fY, &tmp2[0].fY, t); 1204 p3d_interp(&tmp[0].fZ, &tmp2[0].fZ, t); 1205 1206 dst[0].fPts[0] = fPts[0]; 1207 tmp2[0].projectDown(&dst[0].fPts[1]); 1208 tmp2[1].projectDown(&dst[0].fPts[2]); dst[1].fPts[0] = dst[0].fPts[2]; 1209 tmp2[2].projectDown(&dst[1].fPts[1]); 1210 dst[1].fPts[2] = fPts[2]; 1211 1212 // to put in "standard form", where w0 and w2 are both 1, we compute the 1213 // new w1 as sqrt(w1*w1/w0*w2) 1214 // or 1215 // w1 /= sqrt(w0*w2) 1216 // 1217 // However, in our case, we know that for dst[0]: 1218 // w0 == 1, and for dst[1], w2 == 1 1219 // 1220 SkScalar root = SkScalarSqrt(tmp2[1].fZ); 1221 dst[0].fW = tmp2[0].fZ / root; 1222 dst[1].fW = tmp2[2].fZ / root; 1223} 1224 1225static Sk2s times_2(const Sk2s& value) { 1226 return value + value; 1227} 1228 1229SkPoint SkConic::evalAt(SkScalar t) const { 1230 Sk2s p0 = from_point(fPts[0]); 1231 Sk2s p1 = from_point(fPts[1]); 1232 Sk2s p2 = from_point(fPts[2]); 1233 Sk2s tt(t); 1234 Sk2s ww(fW); 1235 Sk2s one(1); 1236 1237 Sk2s p1w = p1 * ww; 1238 Sk2s C = p0; 1239 Sk2s A = p2 - times_2(p1w) + p0; 1240 Sk2s B = times_2(p1w - C); 1241 Sk2s numer = quad_poly_eval(A, B, C, tt); 1242 1243 B = times_2(ww - one); 1244 A = -B; 1245 Sk2s denom = quad_poly_eval(A, B, one, tt); 1246 1247 return to_point(numer / denom); 1248} 1249 1250SkVector SkConic::evalTangentAt(SkScalar t) const { 1251 Sk2s p0 = from_point(fPts[0]); 1252 Sk2s p1 = from_point(fPts[1]); 1253 Sk2s p2 = from_point(fPts[2]); 1254 Sk2s ww(fW); 1255 1256 Sk2s p20 = p2 - p0; 1257 Sk2s p10 = p1 - p0; 1258 1259 Sk2s C = ww * p10; 1260 Sk2s A = ww * p20 - p20; 1261 Sk2s B = p20 - C - C; 1262 1263 return to_vector(quad_poly_eval(A, B, C, Sk2s(t))); 1264} 1265 1266static SkScalar subdivide_w_value(SkScalar w) { 1267 return SkScalarSqrt(SK_ScalarHalf + w * SK_ScalarHalf); 1268} 1269 1270static Sk2s twice(const Sk2s& value) { 1271 return value + value; 1272} 1273 1274void SkConic::chop(SkConic * SK_RESTRICT dst) const { 1275 Sk2s scale = Sk2s(SkScalarInvert(SK_Scalar1 + fW)); 1276 SkScalar newW = subdivide_w_value(fW); 1277 1278 Sk2s p0 = from_point(fPts[0]); 1279 Sk2s p1 = from_point(fPts[1]); 1280 Sk2s p2 = from_point(fPts[2]); 1281 Sk2s ww(fW); 1282 1283 Sk2s wp1 = ww * p1; 1284 Sk2s m = (p0 + twice(wp1) + p2) * scale * Sk2s(0.5f); 1285 1286 dst[0].fPts[0] = fPts[0]; 1287 dst[0].fPts[1] = to_point((p0 + wp1) * scale); 1288 dst[0].fPts[2] = dst[1].fPts[0] = to_point(m); 1289 dst[1].fPts[1] = to_point((wp1 + p2) * scale); 1290 dst[1].fPts[2] = fPts[2]; 1291 1292 dst[0].fW = dst[1].fW = newW; 1293} 1294 1295/* 1296 * "High order approximation of conic sections by quadratic splines" 1297 * by Michael Floater, 1993 1298 */ 1299#define AS_QUAD_ERROR_SETUP \ 1300 SkScalar a = fW - 1; \ 1301 SkScalar k = a / (4 * (2 + a)); \ 1302 SkScalar x = k * (fPts[0].fX - 2 * fPts[1].fX + fPts[2].fX); \ 1303 SkScalar y = k * (fPts[0].fY - 2 * fPts[1].fY + fPts[2].fY); 1304 1305void SkConic::computeAsQuadError(SkVector* err) const { 1306 AS_QUAD_ERROR_SETUP 1307 err->set(x, y); 1308} 1309 1310bool SkConic::asQuadTol(SkScalar tol) const { 1311 AS_QUAD_ERROR_SETUP 1312 return (x * x + y * y) <= tol * tol; 1313} 1314 1315// Limit the number of suggested quads to approximate a conic 1316#define kMaxConicToQuadPOW2 5 1317 1318int SkConic::computeQuadPOW2(SkScalar tol) const { 1319 if (tol < 0 || !SkScalarIsFinite(tol)) { 1320 return 0; 1321 } 1322 1323 AS_QUAD_ERROR_SETUP 1324 1325 SkScalar error = SkScalarSqrt(x * x + y * y); 1326 int pow2; 1327 for (pow2 = 0; pow2 < kMaxConicToQuadPOW2; ++pow2) { 1328 if (error <= tol) { 1329 break; 1330 } 1331 error *= 0.25f; 1332 } 1333 // float version -- using ceil gives the same results as the above. 1334 if (false) { 1335 SkScalar err = SkScalarSqrt(x * x + y * y); 1336 if (err <= tol) { 1337 return 0; 1338 } 1339 SkScalar tol2 = tol * tol; 1340 if (tol2 == 0) { 1341 return kMaxConicToQuadPOW2; 1342 } 1343 SkScalar fpow2 = SkScalarLog2((x * x + y * y) / tol2) * 0.25f; 1344 int altPow2 = SkScalarCeilToInt(fpow2); 1345 if (altPow2 != pow2) { 1346 SkDebugf("pow2 %d altPow2 %d fbits %g err %g tol %g\n", pow2, altPow2, fpow2, err, tol); 1347 } 1348 pow2 = altPow2; 1349 } 1350 return pow2; 1351} 1352 1353static SkPoint* subdivide(const SkConic& src, SkPoint pts[], int level) { 1354 SkASSERT(level >= 0); 1355 1356 if (0 == level) { 1357 memcpy(pts, &src.fPts[1], 2 * sizeof(SkPoint)); 1358 return pts + 2; 1359 } else { 1360 SkConic dst[2]; 1361 src.chop(dst); 1362 --level; 1363 pts = subdivide(dst[0], pts, level); 1364 return subdivide(dst[1], pts, level); 1365 } 1366} 1367 1368int SkConic::chopIntoQuadsPOW2(SkPoint pts[], int pow2) const { 1369 SkASSERT(pow2 >= 0); 1370 *pts = fPts[0]; 1371 SkDEBUGCODE(SkPoint* endPts =) subdivide(*this, pts + 1, pow2); 1372 SkASSERT(endPts - pts == (2 * (1 << pow2) + 1)); 1373 return 1 << pow2; 1374} 1375 1376bool SkConic::findXExtrema(SkScalar* t) const { 1377 return conic_find_extrema(&fPts[0].fX, fW, t); 1378} 1379 1380bool SkConic::findYExtrema(SkScalar* t) const { 1381 return conic_find_extrema(&fPts[0].fY, fW, t); 1382} 1383 1384bool SkConic::chopAtXExtrema(SkConic dst[2]) const { 1385 SkScalar t; 1386 if (this->findXExtrema(&t)) { 1387 this->chopAt(t, dst); 1388 // now clean-up the middle, since we know t was meant to be at 1389 // an X-extrema 1390 SkScalar value = dst[0].fPts[2].fX; 1391 dst[0].fPts[1].fX = value; 1392 dst[1].fPts[0].fX = value; 1393 dst[1].fPts[1].fX = value; 1394 return true; 1395 } 1396 return false; 1397} 1398 1399bool SkConic::chopAtYExtrema(SkConic dst[2]) const { 1400 SkScalar t; 1401 if (this->findYExtrema(&t)) { 1402 this->chopAt(t, dst); 1403 // now clean-up the middle, since we know t was meant to be at 1404 // an Y-extrema 1405 SkScalar value = dst[0].fPts[2].fY; 1406 dst[0].fPts[1].fY = value; 1407 dst[1].fPts[0].fY = value; 1408 dst[1].fPts[1].fY = value; 1409 return true; 1410 } 1411 return false; 1412} 1413 1414void SkConic::computeTightBounds(SkRect* bounds) const { 1415 SkPoint pts[4]; 1416 pts[0] = fPts[0]; 1417 pts[1] = fPts[2]; 1418 int count = 2; 1419 1420 SkScalar t; 1421 if (this->findXExtrema(&t)) { 1422 this->evalAt(t, &pts[count++]); 1423 } 1424 if (this->findYExtrema(&t)) { 1425 this->evalAt(t, &pts[count++]); 1426 } 1427 bounds->set(pts, count); 1428} 1429 1430void SkConic::computeFastBounds(SkRect* bounds) const { 1431 bounds->set(fPts, 3); 1432} 1433 1434bool SkConic::findMaxCurvature(SkScalar* t) const { 1435 // TODO: Implement me 1436 return false; 1437} 1438 1439SkScalar SkConic::TransformW(const SkPoint pts[], SkScalar w, 1440 const SkMatrix& matrix) { 1441 if (!matrix.hasPerspective()) { 1442 return w; 1443 } 1444 1445 SkP3D src[3], dst[3]; 1446 1447 ratquad_mapTo3D(pts, w, src); 1448 1449 matrix.mapHomogeneousPoints(&dst[0].fX, &src[0].fX, 3); 1450 1451 // w' = sqrt(w1*w1/w0*w2) 1452 SkScalar w0 = dst[0].fZ; 1453 SkScalar w1 = dst[1].fZ; 1454 SkScalar w2 = dst[2].fZ; 1455 w = SkScalarSqrt((w1 * w1) / (w0 * w2)); 1456 return w; 1457} 1458 1459int SkConic::BuildUnitArc(const SkVector& uStart, const SkVector& uStop, SkRotationDirection dir, 1460 const SkMatrix* userMatrix, SkConic dst[kMaxConicsForArc]) { 1461 // rotate by x,y so that uStart is (1.0) 1462 SkScalar x = SkPoint::DotProduct(uStart, uStop); 1463 SkScalar y = SkPoint::CrossProduct(uStart, uStop); 1464 1465 SkScalar absY = SkScalarAbs(y); 1466 1467 // check for (effectively) coincident vectors 1468 // this can happen if our angle is nearly 0 or nearly 180 (y == 0) 1469 // ... we use the dot-prod to distinguish between 0 and 180 (x > 0) 1470 if (absY <= SK_ScalarNearlyZero && x > 0 && ((y >= 0 && kCW_SkRotationDirection == dir) || 1471 (y <= 0 && kCCW_SkRotationDirection == dir))) { 1472 return 0; 1473 } 1474 1475 if (dir == kCCW_SkRotationDirection) { 1476 y = -y; 1477 } 1478 1479 // We decide to use 1-conic per quadrant of a circle. What quadrant does [xy] lie in? 1480 // 0 == [0 .. 90) 1481 // 1 == [90 ..180) 1482 // 2 == [180..270) 1483 // 3 == [270..360) 1484 // 1485 int quadrant = 0; 1486 if (0 == y) { 1487 quadrant = 2; // 180 1488 SkASSERT(SkScalarAbs(x + SK_Scalar1) <= SK_ScalarNearlyZero); 1489 } else if (0 == x) { 1490 SkASSERT(absY - SK_Scalar1 <= SK_ScalarNearlyZero); 1491 quadrant = y > 0 ? 1 : 3; // 90 : 270 1492 } else { 1493 if (y < 0) { 1494 quadrant += 2; 1495 } 1496 if ((x < 0) != (y < 0)) { 1497 quadrant += 1; 1498 } 1499 } 1500 1501 const SkPoint quadrantPts[] = { 1502 { 1, 0 }, { 1, 1 }, { 0, 1 }, { -1, 1 }, { -1, 0 }, { -1, -1 }, { 0, -1 }, { 1, -1 } 1503 }; 1504 const SkScalar quadrantWeight = SK_ScalarRoot2Over2; 1505 1506 int conicCount = quadrant; 1507 for (int i = 0; i < conicCount; ++i) { 1508 dst[i].set(&quadrantPts[i * 2], quadrantWeight); 1509 } 1510 1511 // Now compute any remaing (sub-90-degree) arc for the last conic 1512 const SkPoint finalP = { x, y }; 1513 const SkPoint& lastQ = quadrantPts[quadrant * 2]; // will already be a unit-vector 1514 const SkScalar dot = SkVector::DotProduct(lastQ, finalP); 1515 SkASSERT(0 <= dot && dot <= SK_Scalar1 + SK_ScalarNearlyZero); 1516 1517 if (dot < 1 - SK_ScalarNearlyZero) { 1518 SkVector offCurve = { lastQ.x() + x, lastQ.y() + y }; 1519 // compute the bisector vector, and then rescale to be the off-curve point. 1520 // we compute its length from cos(theta/2) = length / 1, using half-angle identity we get 1521 // length = sqrt(2 / (1 + cos(theta)). We already have cos() when to computed the dot. 1522 // This is nice, since our computed weight is cos(theta/2) as well! 1523 // 1524 const SkScalar cosThetaOver2 = SkScalarSqrt((1 + dot) / 2); 1525 offCurve.setLength(SkScalarInvert(cosThetaOver2)); 1526 dst[conicCount].set(lastQ, offCurve, finalP, cosThetaOver2); 1527 conicCount += 1; 1528 } 1529 1530 // now handle counter-clockwise and the initial unitStart rotation 1531 SkMatrix matrix; 1532 matrix.setSinCos(uStart.fY, uStart.fX); 1533 if (dir == kCCW_SkRotationDirection) { 1534 matrix.preScale(SK_Scalar1, -SK_Scalar1); 1535 } 1536 if (userMatrix) { 1537 matrix.postConcat(*userMatrix); 1538 } 1539 for (int i = 0; i < conicCount; ++i) { 1540 matrix.mapPoints(dst[i].fPts, 3); 1541 } 1542 return conicCount; 1543} 1544