1 2/* 3 * Copyright 2006 The Android Open Source Project 4 * 5 * Use of this source code is governed by a BSD-style license that can be 6 * found in the LICENSE file. 7 */ 8 9 10#include "SkGeometry.h" 11#include "Sk64.h" 12#include "SkMatrix.h" 13 14bool SkXRayCrossesLine(const SkXRay& pt, const SkPoint pts[2], bool* ambiguous) { 15 if (ambiguous) { 16 *ambiguous = false; 17 } 18 // Determine quick discards. 19 // Consider query line going exactly through point 0 to not 20 // intersect, for symmetry with SkXRayCrossesMonotonicCubic. 21 if (pt.fY == pts[0].fY) { 22 if (ambiguous) { 23 *ambiguous = true; 24 } 25 return false; 26 } 27 if (pt.fY < pts[0].fY && pt.fY < pts[1].fY) 28 return false; 29 if (pt.fY > pts[0].fY && pt.fY > pts[1].fY) 30 return false; 31 if (pt.fX > pts[0].fX && pt.fX > pts[1].fX) 32 return false; 33 // Determine degenerate cases 34 if (SkScalarNearlyZero(pts[0].fY - pts[1].fY)) 35 return false; 36 if (SkScalarNearlyZero(pts[0].fX - pts[1].fX)) { 37 // We've already determined the query point lies within the 38 // vertical range of the line segment. 39 if (pt.fX <= pts[0].fX) { 40 if (ambiguous) { 41 *ambiguous = (pt.fY == pts[1].fY); 42 } 43 return true; 44 } 45 return false; 46 } 47 // Ambiguity check 48 if (pt.fY == pts[1].fY) { 49 if (pt.fX <= pts[1].fX) { 50 if (ambiguous) { 51 *ambiguous = true; 52 } 53 return true; 54 } 55 return false; 56 } 57 // Full line segment evaluation 58 SkScalar delta_y = pts[1].fY - pts[0].fY; 59 SkScalar delta_x = pts[1].fX - pts[0].fX; 60 SkScalar slope = SkScalarDiv(delta_y, delta_x); 61 SkScalar b = pts[0].fY - SkScalarMul(slope, pts[0].fX); 62 // Solve for x coordinate at y = pt.fY 63 SkScalar x = SkScalarDiv(pt.fY - b, slope); 64 return pt.fX <= x; 65} 66 67/** If defined, this makes eval_quad and eval_cubic do more setup (sometimes 68 involving integer multiplies by 2 or 3, but fewer calls to SkScalarMul. 69 May also introduce overflow of fixed when we compute our setup. 70*/ 71#ifdef SK_SCALAR_IS_FIXED 72 #define DIRECT_EVAL_OF_POLYNOMIALS 73#endif 74 75//////////////////////////////////////////////////////////////////////// 76 77#ifdef SK_SCALAR_IS_FIXED 78 static int is_not_monotonic(int a, int b, int c, int d) 79 { 80 return (((a - b) | (b - c) | (c - d)) & ((b - a) | (c - b) | (d - c))) >> 31; 81 } 82 83 static int is_not_monotonic(int a, int b, int c) 84 { 85 return (((a - b) | (b - c)) & ((b - a) | (c - b))) >> 31; 86 } 87#else 88 static int is_not_monotonic(float a, float b, float c) 89 { 90 float ab = a - b; 91 float bc = b - c; 92 if (ab < 0) 93 bc = -bc; 94 return ab == 0 || bc < 0; 95 } 96#endif 97 98//////////////////////////////////////////////////////////////////////// 99 100static bool is_unit_interval(SkScalar x) 101{ 102 return x > 0 && x < SK_Scalar1; 103} 104 105static int valid_unit_divide(SkScalar numer, SkScalar denom, SkScalar* ratio) 106{ 107 SkASSERT(ratio); 108 109 if (numer < 0) 110 { 111 numer = -numer; 112 denom = -denom; 113 } 114 115 if (denom == 0 || numer == 0 || numer >= denom) 116 return 0; 117 118 SkScalar r = SkScalarDiv(numer, denom); 119 if (SkScalarIsNaN(r)) { 120 return 0; 121 } 122 SkASSERT(r >= 0 && r < SK_Scalar1); 123 if (r == 0) // catch underflow if numer <<<< denom 124 return 0; 125 *ratio = r; 126 return 1; 127} 128 129/** From Numerical Recipes in C. 130 131 Q = -1/2 (B + sign(B) sqrt[B*B - 4*A*C]) 132 x1 = Q / A 133 x2 = C / Q 134*/ 135int SkFindUnitQuadRoots(SkScalar A, SkScalar B, SkScalar C, SkScalar roots[2]) 136{ 137 SkASSERT(roots); 138 139 if (A == 0) 140 return valid_unit_divide(-C, B, roots); 141 142 SkScalar* r = roots; 143 144#ifdef SK_SCALAR_IS_FLOAT 145 float R = B*B - 4*A*C; 146 if (R < 0 || SkScalarIsNaN(R)) { // complex roots 147 return 0; 148 } 149 R = sk_float_sqrt(R); 150#else 151 Sk64 RR, tmp; 152 153 RR.setMul(B,B); 154 tmp.setMul(A,C); 155 tmp.shiftLeft(2); 156 RR.sub(tmp); 157 if (RR.isNeg()) 158 return 0; 159 SkFixed R = RR.getSqrt(); 160#endif 161 162 SkScalar Q = (B < 0) ? -(B-R)/2 : -(B+R)/2; 163 r += valid_unit_divide(Q, A, r); 164 r += valid_unit_divide(C, Q, r); 165 if (r - roots == 2) 166 { 167 if (roots[0] > roots[1]) 168 SkTSwap<SkScalar>(roots[0], roots[1]); 169 else if (roots[0] == roots[1]) // nearly-equal? 170 r -= 1; // skip the double root 171 } 172 return (int)(r - roots); 173} 174 175#ifdef SK_SCALAR_IS_FIXED 176/** Trim A/B/C down so that they are all <= 32bits 177 and then call SkFindUnitQuadRoots() 178*/ 179static int Sk64FindFixedQuadRoots(const Sk64& A, const Sk64& B, const Sk64& C, SkFixed roots[2]) 180{ 181 int na = A.shiftToMake32(); 182 int nb = B.shiftToMake32(); 183 int nc = C.shiftToMake32(); 184 185 int shift = SkMax32(na, SkMax32(nb, nc)); 186 SkASSERT(shift >= 0); 187 188 return SkFindUnitQuadRoots(A.getShiftRight(shift), B.getShiftRight(shift), C.getShiftRight(shift), roots); 189} 190#endif 191 192///////////////////////////////////////////////////////////////////////////////////// 193///////////////////////////////////////////////////////////////////////////////////// 194 195static SkScalar eval_quad(const SkScalar src[], SkScalar t) 196{ 197 SkASSERT(src); 198 SkASSERT(t >= 0 && t <= SK_Scalar1); 199 200#ifdef DIRECT_EVAL_OF_POLYNOMIALS 201 SkScalar C = src[0]; 202 SkScalar A = src[4] - 2 * src[2] + C; 203 SkScalar B = 2 * (src[2] - C); 204 return SkScalarMulAdd(SkScalarMulAdd(A, t, B), t, C); 205#else 206 SkScalar ab = SkScalarInterp(src[0], src[2], t); 207 SkScalar bc = SkScalarInterp(src[2], src[4], t); 208 return SkScalarInterp(ab, bc, t); 209#endif 210} 211 212static SkScalar eval_quad_derivative(const SkScalar src[], SkScalar t) 213{ 214 SkScalar A = src[4] - 2 * src[2] + src[0]; 215 SkScalar B = src[2] - src[0]; 216 217 return 2 * SkScalarMulAdd(A, t, B); 218} 219 220static SkScalar eval_quad_derivative_at_half(const SkScalar src[]) 221{ 222 SkScalar A = src[4] - 2 * src[2] + src[0]; 223 SkScalar B = src[2] - src[0]; 224 return A + 2 * B; 225} 226 227void SkEvalQuadAt(const SkPoint src[3], SkScalar t, SkPoint* pt, SkVector* tangent) 228{ 229 SkASSERT(src); 230 SkASSERT(t >= 0 && t <= SK_Scalar1); 231 232 if (pt) 233 pt->set(eval_quad(&src[0].fX, t), eval_quad(&src[0].fY, t)); 234 if (tangent) 235 tangent->set(eval_quad_derivative(&src[0].fX, t), 236 eval_quad_derivative(&src[0].fY, t)); 237} 238 239void SkEvalQuadAtHalf(const SkPoint src[3], SkPoint* pt, SkVector* tangent) 240{ 241 SkASSERT(src); 242 243 if (pt) 244 { 245 SkScalar x01 = SkScalarAve(src[0].fX, src[1].fX); 246 SkScalar y01 = SkScalarAve(src[0].fY, src[1].fY); 247 SkScalar x12 = SkScalarAve(src[1].fX, src[2].fX); 248 SkScalar y12 = SkScalarAve(src[1].fY, src[2].fY); 249 pt->set(SkScalarAve(x01, x12), SkScalarAve(y01, y12)); 250 } 251 if (tangent) 252 tangent->set(eval_quad_derivative_at_half(&src[0].fX), 253 eval_quad_derivative_at_half(&src[0].fY)); 254} 255 256static void interp_quad_coords(const SkScalar* src, SkScalar* dst, SkScalar t) 257{ 258 SkScalar ab = SkScalarInterp(src[0], src[2], t); 259 SkScalar bc = SkScalarInterp(src[2], src[4], t); 260 261 dst[0] = src[0]; 262 dst[2] = ab; 263 dst[4] = SkScalarInterp(ab, bc, t); 264 dst[6] = bc; 265 dst[8] = src[4]; 266} 267 268void SkChopQuadAt(const SkPoint src[3], SkPoint dst[5], SkScalar t) 269{ 270 SkASSERT(t > 0 && t < SK_Scalar1); 271 272 interp_quad_coords(&src[0].fX, &dst[0].fX, t); 273 interp_quad_coords(&src[0].fY, &dst[0].fY, t); 274} 275 276void SkChopQuadAtHalf(const SkPoint src[3], SkPoint dst[5]) 277{ 278 SkScalar x01 = SkScalarAve(src[0].fX, src[1].fX); 279 SkScalar y01 = SkScalarAve(src[0].fY, src[1].fY); 280 SkScalar x12 = SkScalarAve(src[1].fX, src[2].fX); 281 SkScalar y12 = SkScalarAve(src[1].fY, src[2].fY); 282 283 dst[0] = src[0]; 284 dst[1].set(x01, y01); 285 dst[2].set(SkScalarAve(x01, x12), SkScalarAve(y01, y12)); 286 dst[3].set(x12, y12); 287 dst[4] = src[2]; 288} 289 290/** Quad'(t) = At + B, where 291 A = 2(a - 2b + c) 292 B = 2(b - a) 293 Solve for t, only if it fits between 0 < t < 1 294*/ 295int SkFindQuadExtrema(SkScalar a, SkScalar b, SkScalar c, SkScalar tValue[1]) 296{ 297 /* At + B == 0 298 t = -B / A 299 */ 300#ifdef SK_SCALAR_IS_FIXED 301 return is_not_monotonic(a, b, c) && valid_unit_divide(a - b, a - b - b + c, tValue); 302#else 303 return valid_unit_divide(a - b, a - b - b + c, tValue); 304#endif 305} 306 307static inline void flatten_double_quad_extrema(SkScalar coords[14]) 308{ 309 coords[2] = coords[6] = coords[4]; 310} 311 312/* Returns 0 for 1 quad, and 1 for two quads, either way the answer is 313 stored in dst[]. Guarantees that the 1/2 quads will be monotonic. 314 */ 315int SkChopQuadAtYExtrema(const SkPoint src[3], SkPoint dst[5]) 316{ 317 SkASSERT(src); 318 SkASSERT(dst); 319 320#if 0 321 static bool once = true; 322 if (once) 323 { 324 once = false; 325 SkPoint s[3] = { 0, 26398, 0, 26331, 0, 20621428 }; 326 SkPoint d[6]; 327 328 int n = SkChopQuadAtYExtrema(s, d); 329 SkDebugf("chop=%d, Y=[%x %x %x %x %x %x]\n", n, d[0].fY, d[1].fY, d[2].fY, d[3].fY, d[4].fY, d[5].fY); 330 } 331#endif 332 333 SkScalar a = src[0].fY; 334 SkScalar b = src[1].fY; 335 SkScalar c = src[2].fY; 336 337 if (is_not_monotonic(a, b, c)) 338 { 339 SkScalar tValue; 340 if (valid_unit_divide(a - b, a - b - b + c, &tValue)) 341 { 342 SkChopQuadAt(src, dst, tValue); 343 flatten_double_quad_extrema(&dst[0].fY); 344 return 1; 345 } 346 // if we get here, we need to force dst to be monotonic, even though 347 // we couldn't compute a unit_divide value (probably underflow). 348 b = SkScalarAbs(a - b) < SkScalarAbs(b - c) ? a : c; 349 } 350 dst[0].set(src[0].fX, a); 351 dst[1].set(src[1].fX, b); 352 dst[2].set(src[2].fX, c); 353 return 0; 354} 355 356/* Returns 0 for 1 quad, and 1 for two quads, either way the answer is 357 stored in dst[]. Guarantees that the 1/2 quads will be monotonic. 358 */ 359int SkChopQuadAtXExtrema(const SkPoint src[3], SkPoint dst[5]) 360{ 361 SkASSERT(src); 362 SkASSERT(dst); 363 364 SkScalar a = src[0].fX; 365 SkScalar b = src[1].fX; 366 SkScalar c = src[2].fX; 367 368 if (is_not_monotonic(a, b, c)) { 369 SkScalar tValue; 370 if (valid_unit_divide(a - b, a - b - b + c, &tValue)) { 371 SkChopQuadAt(src, dst, tValue); 372 flatten_double_quad_extrema(&dst[0].fX); 373 return 1; 374 } 375 // if we get here, we need to force dst to be monotonic, even though 376 // we couldn't compute a unit_divide value (probably underflow). 377 b = SkScalarAbs(a - b) < SkScalarAbs(b - c) ? a : c; 378 } 379 dst[0].set(a, src[0].fY); 380 dst[1].set(b, src[1].fY); 381 dst[2].set(c, src[2].fY); 382 return 0; 383} 384 385// F(t) = a (1 - t) ^ 2 + 2 b t (1 - t) + c t ^ 2 386// F'(t) = 2 (b - a) + 2 (a - 2b + c) t 387// F''(t) = 2 (a - 2b + c) 388// 389// A = 2 (b - a) 390// B = 2 (a - 2b + c) 391// 392// Maximum curvature for a quadratic means solving 393// Fx' Fx'' + Fy' Fy'' = 0 394// 395// t = - (Ax Bx + Ay By) / (Bx ^ 2 + By ^ 2) 396// 397float SkFindQuadMaxCurvature(const SkPoint src[3]) { 398 SkScalar Ax = src[1].fX - src[0].fX; 399 SkScalar Ay = src[1].fY - src[0].fY; 400 SkScalar Bx = src[0].fX - src[1].fX - src[1].fX + src[2].fX; 401 SkScalar By = src[0].fY - src[1].fY - src[1].fY + src[2].fY; 402 SkScalar t = 0; // 0 means don't chop 403 404#ifdef SK_SCALAR_IS_FLOAT 405 (void)valid_unit_divide(-(Ax * Bx + Ay * By), Bx * Bx + By * By, &t); 406#else 407 // !!! should I use SkFloat here? seems like it 408 Sk64 numer, denom, tmp; 409 410 numer.setMul(Ax, -Bx); 411 tmp.setMul(Ay, -By); 412 numer.add(tmp); 413 414 if (numer.isPos()) // do nothing if numer <= 0 415 { 416 denom.setMul(Bx, Bx); 417 tmp.setMul(By, By); 418 denom.add(tmp); 419 SkASSERT(!denom.isNeg()); 420 if (numer < denom) 421 { 422 t = numer.getFixedDiv(denom); 423 SkASSERT(t >= 0 && t <= SK_Fixed1); // assert that we're numerically stable (ha!) 424 if ((unsigned)t >= SK_Fixed1) // runtime check for numerical stability 425 t = 0; // ignore the chop 426 } 427 } 428#endif 429 return t; 430} 431 432int SkChopQuadAtMaxCurvature(const SkPoint src[3], SkPoint dst[5]) 433{ 434 SkScalar t = SkFindQuadMaxCurvature(src); 435 if (t == 0) { 436 memcpy(dst, src, 3 * sizeof(SkPoint)); 437 return 1; 438 } else { 439 SkChopQuadAt(src, dst, t); 440 return 2; 441 } 442} 443 444#ifdef SK_SCALAR_IS_FLOAT 445 #define SK_ScalarTwoThirds (0.666666666f) 446#else 447 #define SK_ScalarTwoThirds ((SkFixed)(43691)) 448#endif 449 450void SkConvertQuadToCubic(const SkPoint src[3], SkPoint dst[4]) { 451 const SkScalar scale = SK_ScalarTwoThirds; 452 dst[0] = src[0]; 453 dst[1].set(src[0].fX + SkScalarMul(src[1].fX - src[0].fX, scale), 454 src[0].fY + SkScalarMul(src[1].fY - src[0].fY, scale)); 455 dst[2].set(src[2].fX + SkScalarMul(src[1].fX - src[2].fX, scale), 456 src[2].fY + SkScalarMul(src[1].fY - src[2].fY, scale)); 457 dst[3] = src[2]; 458} 459 460//////////////////////////////////////////////////////////////////////////////////////// 461///// CUBICS // CUBICS // CUBICS // CUBICS // CUBICS // CUBICS // CUBICS // CUBICS ///// 462//////////////////////////////////////////////////////////////////////////////////////// 463 464static void get_cubic_coeff(const SkScalar pt[], SkScalar coeff[4]) 465{ 466 coeff[0] = pt[6] + 3*(pt[2] - pt[4]) - pt[0]; 467 coeff[1] = 3*(pt[4] - pt[2] - pt[2] + pt[0]); 468 coeff[2] = 3*(pt[2] - pt[0]); 469 coeff[3] = pt[0]; 470} 471 472void SkGetCubicCoeff(const SkPoint pts[4], SkScalar cx[4], SkScalar cy[4]) 473{ 474 SkASSERT(pts); 475 476 if (cx) 477 get_cubic_coeff(&pts[0].fX, cx); 478 if (cy) 479 get_cubic_coeff(&pts[0].fY, cy); 480} 481 482static SkScalar eval_cubic(const SkScalar src[], SkScalar t) 483{ 484 SkASSERT(src); 485 SkASSERT(t >= 0 && t <= SK_Scalar1); 486 487 if (t == 0) 488 return src[0]; 489 490#ifdef DIRECT_EVAL_OF_POLYNOMIALS 491 SkScalar D = src[0]; 492 SkScalar A = src[6] + 3*(src[2] - src[4]) - D; 493 SkScalar B = 3*(src[4] - src[2] - src[2] + D); 494 SkScalar C = 3*(src[2] - D); 495 496 return SkScalarMulAdd(SkScalarMulAdd(SkScalarMulAdd(A, t, B), t, C), t, D); 497#else 498 SkScalar ab = SkScalarInterp(src[0], src[2], t); 499 SkScalar bc = SkScalarInterp(src[2], src[4], t); 500 SkScalar cd = SkScalarInterp(src[4], src[6], t); 501 SkScalar abc = SkScalarInterp(ab, bc, t); 502 SkScalar bcd = SkScalarInterp(bc, cd, t); 503 return SkScalarInterp(abc, bcd, t); 504#endif 505} 506 507/** return At^2 + Bt + C 508*/ 509static SkScalar eval_quadratic(SkScalar A, SkScalar B, SkScalar C, SkScalar t) 510{ 511 SkASSERT(t >= 0 && t <= SK_Scalar1); 512 513 return SkScalarMulAdd(SkScalarMulAdd(A, t, B), t, C); 514} 515 516static SkScalar eval_cubic_derivative(const SkScalar src[], SkScalar t) 517{ 518 SkScalar A = src[6] + 3*(src[2] - src[4]) - src[0]; 519 SkScalar B = 2*(src[4] - 2 * src[2] + src[0]); 520 SkScalar C = src[2] - src[0]; 521 522 return eval_quadratic(A, B, C, t); 523} 524 525static SkScalar eval_cubic_2ndDerivative(const SkScalar src[], SkScalar t) 526{ 527 SkScalar A = src[6] + 3*(src[2] - src[4]) - src[0]; 528 SkScalar B = src[4] - 2 * src[2] + src[0]; 529 530 return SkScalarMulAdd(A, t, B); 531} 532 533void SkEvalCubicAt(const SkPoint src[4], SkScalar t, SkPoint* loc, SkVector* tangent, SkVector* curvature) 534{ 535 SkASSERT(src); 536 SkASSERT(t >= 0 && t <= SK_Scalar1); 537 538 if (loc) 539 loc->set(eval_cubic(&src[0].fX, t), eval_cubic(&src[0].fY, t)); 540 if (tangent) 541 tangent->set(eval_cubic_derivative(&src[0].fX, t), 542 eval_cubic_derivative(&src[0].fY, t)); 543 if (curvature) 544 curvature->set(eval_cubic_2ndDerivative(&src[0].fX, t), 545 eval_cubic_2ndDerivative(&src[0].fY, t)); 546} 547 548/** Cubic'(t) = At^2 + Bt + C, where 549 A = 3(-a + 3(b - c) + d) 550 B = 6(a - 2b + c) 551 C = 3(b - a) 552 Solve for t, keeping only those that fit betwee 0 < t < 1 553*/ 554int SkFindCubicExtrema(SkScalar a, SkScalar b, SkScalar c, SkScalar d, SkScalar tValues[2]) 555{ 556#ifdef SK_SCALAR_IS_FIXED 557 if (!is_not_monotonic(a, b, c, d)) 558 return 0; 559#endif 560 561 // we divide A,B,C by 3 to simplify 562 SkScalar A = d - a + 3*(b - c); 563 SkScalar B = 2*(a - b - b + c); 564 SkScalar C = b - a; 565 566 return SkFindUnitQuadRoots(A, B, C, tValues); 567} 568 569static void interp_cubic_coords(const SkScalar* src, SkScalar* dst, SkScalar t) 570{ 571 SkScalar ab = SkScalarInterp(src[0], src[2], t); 572 SkScalar bc = SkScalarInterp(src[2], src[4], t); 573 SkScalar cd = SkScalarInterp(src[4], src[6], t); 574 SkScalar abc = SkScalarInterp(ab, bc, t); 575 SkScalar bcd = SkScalarInterp(bc, cd, t); 576 SkScalar abcd = SkScalarInterp(abc, bcd, t); 577 578 dst[0] = src[0]; 579 dst[2] = ab; 580 dst[4] = abc; 581 dst[6] = abcd; 582 dst[8] = bcd; 583 dst[10] = cd; 584 dst[12] = src[6]; 585} 586 587void SkChopCubicAt(const SkPoint src[4], SkPoint dst[7], SkScalar t) 588{ 589 SkASSERT(t > 0 && t < SK_Scalar1); 590 591 interp_cubic_coords(&src[0].fX, &dst[0].fX, t); 592 interp_cubic_coords(&src[0].fY, &dst[0].fY, t); 593} 594 595/* http://code.google.com/p/skia/issues/detail?id=32 596 597 This test code would fail when we didn't check the return result of 598 valid_unit_divide in SkChopCubicAt(... tValues[], int roots). The reason is 599 that after the first chop, the parameters to valid_unit_divide are equal 600 (thanks to finite float precision and rounding in the subtracts). Thus 601 even though the 2nd tValue looks < 1.0, after we renormalize it, we end 602 up with 1.0, hence the need to check and just return the last cubic as 603 a degenerate clump of 4 points in the sampe place. 604 605 static void test_cubic() { 606 SkPoint src[4] = { 607 { 556.25000, 523.03003 }, 608 { 556.23999, 522.96002 }, 609 { 556.21997, 522.89001 }, 610 { 556.21997, 522.82001 } 611 }; 612 SkPoint dst[10]; 613 SkScalar tval[] = { 0.33333334f, 0.99999994f }; 614 SkChopCubicAt(src, dst, tval, 2); 615 } 616 */ 617 618void SkChopCubicAt(const SkPoint src[4], SkPoint dst[], const SkScalar tValues[], int roots) 619{ 620#ifdef SK_DEBUG 621 { 622 for (int i = 0; i < roots - 1; i++) 623 { 624 SkASSERT(is_unit_interval(tValues[i])); 625 SkASSERT(is_unit_interval(tValues[i+1])); 626 SkASSERT(tValues[i] < tValues[i+1]); 627 } 628 } 629#endif 630 631 if (dst) 632 { 633 if (roots == 0) // nothing to chop 634 memcpy(dst, src, 4*sizeof(SkPoint)); 635 else 636 { 637 SkScalar t = tValues[0]; 638 SkPoint tmp[4]; 639 640 for (int i = 0; i < roots; i++) 641 { 642 SkChopCubicAt(src, dst, t); 643 if (i == roots - 1) 644 break; 645 646 dst += 3; 647 // have src point to the remaining cubic (after the chop) 648 memcpy(tmp, dst, 4 * sizeof(SkPoint)); 649 src = tmp; 650 651 // watch out in case the renormalized t isn't in range 652 if (!valid_unit_divide(tValues[i+1] - tValues[i], 653 SK_Scalar1 - tValues[i], &t)) { 654 // if we can't, just create a degenerate cubic 655 dst[4] = dst[5] = dst[6] = src[3]; 656 break; 657 } 658 } 659 } 660 } 661} 662 663void SkChopCubicAtHalf(const SkPoint src[4], SkPoint dst[7]) 664{ 665 SkScalar x01 = SkScalarAve(src[0].fX, src[1].fX); 666 SkScalar y01 = SkScalarAve(src[0].fY, src[1].fY); 667 SkScalar x12 = SkScalarAve(src[1].fX, src[2].fX); 668 SkScalar y12 = SkScalarAve(src[1].fY, src[2].fY); 669 SkScalar x23 = SkScalarAve(src[2].fX, src[3].fX); 670 SkScalar y23 = SkScalarAve(src[2].fY, src[3].fY); 671 672 SkScalar x012 = SkScalarAve(x01, x12); 673 SkScalar y012 = SkScalarAve(y01, y12); 674 SkScalar x123 = SkScalarAve(x12, x23); 675 SkScalar y123 = SkScalarAve(y12, y23); 676 677 dst[0] = src[0]; 678 dst[1].set(x01, y01); 679 dst[2].set(x012, y012); 680 dst[3].set(SkScalarAve(x012, x123), SkScalarAve(y012, y123)); 681 dst[4].set(x123, y123); 682 dst[5].set(x23, y23); 683 dst[6] = src[3]; 684} 685 686static void flatten_double_cubic_extrema(SkScalar coords[14]) 687{ 688 coords[4] = coords[8] = coords[6]; 689} 690 691/** Given 4 points on a cubic bezier, chop it into 1, 2, 3 beziers such that 692 the resulting beziers are monotonic in Y. This is called by the scan converter. 693 Depending on what is returned, dst[] is treated as follows 694 0 dst[0..3] is the original cubic 695 1 dst[0..3] and dst[3..6] are the two new cubics 696 2 dst[0..3], dst[3..6], dst[6..9] are the three new cubics 697 If dst == null, it is ignored and only the count is returned. 698*/ 699int SkChopCubicAtYExtrema(const SkPoint src[4], SkPoint dst[10]) { 700 SkScalar tValues[2]; 701 int roots = SkFindCubicExtrema(src[0].fY, src[1].fY, src[2].fY, 702 src[3].fY, tValues); 703 704 SkChopCubicAt(src, dst, tValues, roots); 705 if (dst && roots > 0) { 706 // we do some cleanup to ensure our Y extrema are flat 707 flatten_double_cubic_extrema(&dst[0].fY); 708 if (roots == 2) { 709 flatten_double_cubic_extrema(&dst[3].fY); 710 } 711 } 712 return roots; 713} 714 715int SkChopCubicAtXExtrema(const SkPoint src[4], SkPoint dst[10]) { 716 SkScalar tValues[2]; 717 int roots = SkFindCubicExtrema(src[0].fX, src[1].fX, src[2].fX, 718 src[3].fX, tValues); 719 720 SkChopCubicAt(src, dst, tValues, roots); 721 if (dst && roots > 0) { 722 // we do some cleanup to ensure our Y extrema are flat 723 flatten_double_cubic_extrema(&dst[0].fX); 724 if (roots == 2) { 725 flatten_double_cubic_extrema(&dst[3].fX); 726 } 727 } 728 return roots; 729} 730 731/** http://www.faculty.idc.ac.il/arik/quality/appendixA.html 732 733 Inflection means that curvature is zero. 734 Curvature is [F' x F''] / [F'^3] 735 So we solve F'x X F''y - F'y X F''y == 0 736 After some canceling of the cubic term, we get 737 A = b - a 738 B = c - 2b + a 739 C = d - 3c + 3b - a 740 (BxCy - ByCx)t^2 + (AxCy - AyCx)t + AxBy - AyBx == 0 741*/ 742int SkFindCubicInflections(const SkPoint src[4], SkScalar tValues[]) 743{ 744 SkScalar Ax = src[1].fX - src[0].fX; 745 SkScalar Ay = src[1].fY - src[0].fY; 746 SkScalar Bx = src[2].fX - 2 * src[1].fX + src[0].fX; 747 SkScalar By = src[2].fY - 2 * src[1].fY + src[0].fY; 748 SkScalar Cx = src[3].fX + 3 * (src[1].fX - src[2].fX) - src[0].fX; 749 SkScalar Cy = src[3].fY + 3 * (src[1].fY - src[2].fY) - src[0].fY; 750 int count; 751 752#ifdef SK_SCALAR_IS_FLOAT 753 count = SkFindUnitQuadRoots(Bx*Cy - By*Cx, Ax*Cy - Ay*Cx, Ax*By - Ay*Bx, tValues); 754#else 755 Sk64 A, B, C, tmp; 756 757 A.setMul(Bx, Cy); 758 tmp.setMul(By, Cx); 759 A.sub(tmp); 760 761 B.setMul(Ax, Cy); 762 tmp.setMul(Ay, Cx); 763 B.sub(tmp); 764 765 C.setMul(Ax, By); 766 tmp.setMul(Ay, Bx); 767 C.sub(tmp); 768 769 count = Sk64FindFixedQuadRoots(A, B, C, tValues); 770#endif 771 772 return count; 773} 774 775int SkChopCubicAtInflections(const SkPoint src[], SkPoint dst[10]) 776{ 777 SkScalar tValues[2]; 778 int count = SkFindCubicInflections(src, tValues); 779 780 if (dst) 781 { 782 if (count == 0) 783 memcpy(dst, src, 4 * sizeof(SkPoint)); 784 else 785 SkChopCubicAt(src, dst, tValues, count); 786 } 787 return count + 1; 788} 789 790template <typename T> void bubble_sort(T array[], int count) 791{ 792 for (int i = count - 1; i > 0; --i) 793 for (int j = i; j > 0; --j) 794 if (array[j] < array[j-1]) 795 { 796 T tmp(array[j]); 797 array[j] = array[j-1]; 798 array[j-1] = tmp; 799 } 800} 801 802#include "SkFP.h" 803 804// newton refinement 805#if 0 806static SkScalar refine_cubic_root(const SkFP coeff[4], SkScalar root) 807{ 808 // x1 = x0 - f(t) / f'(t) 809 810 SkFP T = SkScalarToFloat(root); 811 SkFP N, D; 812 813 // f' = 3*coeff[0]*T^2 + 2*coeff[1]*T + coeff[2] 814 D = SkFPMul(SkFPMul(coeff[0], SkFPMul(T,T)), 3); 815 D = SkFPAdd(D, SkFPMulInt(SkFPMul(coeff[1], T), 2)); 816 D = SkFPAdd(D, coeff[2]); 817 818 if (D == 0) 819 return root; 820 821 // f = coeff[0]*T^3 + coeff[1]*T^2 + coeff[2]*T + coeff[3] 822 N = SkFPMul(SkFPMul(SkFPMul(T, T), T), coeff[0]); 823 N = SkFPAdd(N, SkFPMul(SkFPMul(T, T), coeff[1])); 824 N = SkFPAdd(N, SkFPMul(T, coeff[2])); 825 N = SkFPAdd(N, coeff[3]); 826 827 if (N) 828 { 829 SkScalar delta = SkFPToScalar(SkFPDiv(N, D)); 830 831 if (delta) 832 root -= delta; 833 } 834 return root; 835} 836#endif 837 838/** 839 * Given an array and count, remove all pair-wise duplicates from the array, 840 * keeping the existing sorting, and return the new count 841 */ 842static int collaps_duplicates(float array[], int count) { 843 for (int n = count; n > 1; --n) { 844 if (array[0] == array[1]) { 845 for (int i = 1; i < n; ++i) { 846 array[i - 1] = array[i]; 847 } 848 count -= 1; 849 } else { 850 array += 1; 851 } 852 } 853 return count; 854} 855 856#ifdef SK_DEBUG 857 858#define TEST_COLLAPS_ENTRY(array) array, SK_ARRAY_COUNT(array) 859 860static void test_collaps_duplicates() { 861 static bool gOnce; 862 if (gOnce) { return; } 863 gOnce = true; 864 const float src0[] = { 0 }; 865 const float src1[] = { 0, 0 }; 866 const float src2[] = { 0, 1 }; 867 const float src3[] = { 0, 0, 0 }; 868 const float src4[] = { 0, 0, 1 }; 869 const float src5[] = { 0, 1, 1 }; 870 const float src6[] = { 0, 1, 2 }; 871 const struct { 872 const float* fData; 873 int fCount; 874 int fCollapsedCount; 875 } data[] = { 876 { TEST_COLLAPS_ENTRY(src0), 1 }, 877 { TEST_COLLAPS_ENTRY(src1), 1 }, 878 { TEST_COLLAPS_ENTRY(src2), 2 }, 879 { TEST_COLLAPS_ENTRY(src3), 1 }, 880 { TEST_COLLAPS_ENTRY(src4), 2 }, 881 { TEST_COLLAPS_ENTRY(src5), 2 }, 882 { TEST_COLLAPS_ENTRY(src6), 3 }, 883 }; 884 for (size_t i = 0; i < SK_ARRAY_COUNT(data); ++i) { 885 float dst[3]; 886 memcpy(dst, data[i].fData, data[i].fCount * sizeof(dst[0])); 887 int count = collaps_duplicates(dst, data[i].fCount); 888 SkASSERT(data[i].fCollapsedCount == count); 889 for (int j = 1; j < count; ++j) { 890 SkASSERT(dst[j-1] < dst[j]); 891 } 892 } 893} 894#endif 895 896#if defined _WIN32 && _MSC_VER >= 1300 && defined SK_SCALAR_IS_FIXED // disable warning : unreachable code if building fixed point for windows desktop 897#pragma warning ( disable : 4702 ) 898#endif 899 900/* Solve coeff(t) == 0, returning the number of roots that 901 lie withing 0 < t < 1. 902 coeff[0]t^3 + coeff[1]t^2 + coeff[2]t + coeff[3] 903 904 Eliminates repeated roots (so that all tValues are distinct, and are always 905 in increasing order. 906*/ 907static int solve_cubic_polynomial(const SkFP coeff[4], SkScalar tValues[3]) 908{ 909#ifndef SK_SCALAR_IS_FLOAT 910 return 0; // this is not yet implemented for software float 911#endif 912 913 if (SkScalarNearlyZero(coeff[0])) // we're just a quadratic 914 { 915 return SkFindUnitQuadRoots(coeff[1], coeff[2], coeff[3], tValues); 916 } 917 918 SkFP a, b, c, Q, R; 919 920 { 921 SkASSERT(coeff[0] != 0); 922 923 SkFP inva = SkFPInvert(coeff[0]); 924 a = SkFPMul(coeff[1], inva); 925 b = SkFPMul(coeff[2], inva); 926 c = SkFPMul(coeff[3], inva); 927 } 928 Q = SkFPDivInt(SkFPSub(SkFPMul(a,a), SkFPMulInt(b, 3)), 9); 929// R = (2*a*a*a - 9*a*b + 27*c) / 54; 930 R = SkFPMulInt(SkFPMul(SkFPMul(a, a), a), 2); 931 R = SkFPSub(R, SkFPMulInt(SkFPMul(a, b), 9)); 932 R = SkFPAdd(R, SkFPMulInt(c, 27)); 933 R = SkFPDivInt(R, 54); 934 935 SkFP Q3 = SkFPMul(SkFPMul(Q, Q), Q); 936 SkFP R2MinusQ3 = SkFPSub(SkFPMul(R,R), Q3); 937 SkFP adiv3 = SkFPDivInt(a, 3); 938 939 SkScalar* roots = tValues; 940 SkScalar r; 941 942 if (SkFPLT(R2MinusQ3, 0)) // we have 3 real roots 943 { 944#ifdef SK_SCALAR_IS_FLOAT 945 float theta = sk_float_acos(R / sk_float_sqrt(Q3)); 946 float neg2RootQ = -2 * sk_float_sqrt(Q); 947 948 r = neg2RootQ * sk_float_cos(theta/3) - adiv3; 949 if (is_unit_interval(r)) 950 *roots++ = r; 951 952 r = neg2RootQ * sk_float_cos((theta + 2*SK_ScalarPI)/3) - adiv3; 953 if (is_unit_interval(r)) 954 *roots++ = r; 955 956 r = neg2RootQ * sk_float_cos((theta - 2*SK_ScalarPI)/3) - adiv3; 957 if (is_unit_interval(r)) 958 *roots++ = r; 959 960 SkDEBUGCODE(test_collaps_duplicates();) 961 962 // now sort the roots 963 int count = (int)(roots - tValues); 964 SkASSERT((unsigned)count <= 3); 965 bubble_sort(tValues, count); 966 count = collaps_duplicates(tValues, count); 967 roots = tValues + count; // so we compute the proper count below 968#endif 969 } 970 else // we have 1 real root 971 { 972 SkFP A = SkFPAdd(SkFPAbs(R), SkFPSqrt(R2MinusQ3)); 973 A = SkFPCubeRoot(A); 974 if (SkFPGT(R, 0)) 975 A = SkFPNeg(A); 976 977 if (A != 0) 978 A = SkFPAdd(A, SkFPDiv(Q, A)); 979 r = SkFPToScalar(SkFPSub(A, adiv3)); 980 if (is_unit_interval(r)) 981 *roots++ = r; 982 } 983 984 return (int)(roots - tValues); 985} 986 987/* Looking for F' dot F'' == 0 988 989 A = b - a 990 B = c - 2b + a 991 C = d - 3c + 3b - a 992 993 F' = 3Ct^2 + 6Bt + 3A 994 F'' = 6Ct + 6B 995 996 F' dot F'' -> CCt^3 + 3BCt^2 + (2BB + CA)t + AB 997*/ 998static void formulate_F1DotF2(const SkScalar src[], SkFP coeff[4]) 999{ 1000 SkScalar a = src[2] - src[0]; 1001 SkScalar b = src[4] - 2 * src[2] + src[0]; 1002 SkScalar c = src[6] + 3 * (src[2] - src[4]) - src[0]; 1003 1004 SkFP A = SkScalarToFP(a); 1005 SkFP B = SkScalarToFP(b); 1006 SkFP C = SkScalarToFP(c); 1007 1008 coeff[0] = SkFPMul(C, C); 1009 coeff[1] = SkFPMulInt(SkFPMul(B, C), 3); 1010 coeff[2] = SkFPMulInt(SkFPMul(B, B), 2); 1011 coeff[2] = SkFPAdd(coeff[2], SkFPMul(C, A)); 1012 coeff[3] = SkFPMul(A, B); 1013} 1014 1015// EXPERIMENTAL: can set this to zero to accept all t-values 0 < t < 1 1016//#define kMinTValueForChopping (SK_Scalar1 / 256) 1017#define kMinTValueForChopping 0 1018 1019/* Looking for F' dot F'' == 0 1020 1021 A = b - a 1022 B = c - 2b + a 1023 C = d - 3c + 3b - a 1024 1025 F' = 3Ct^2 + 6Bt + 3A 1026 F'' = 6Ct + 6B 1027 1028 F' dot F'' -> CCt^3 + 3BCt^2 + (2BB + CA)t + AB 1029*/ 1030int SkFindCubicMaxCurvature(const SkPoint src[4], SkScalar tValues[3]) 1031{ 1032 SkFP coeffX[4], coeffY[4]; 1033 int i; 1034 1035 formulate_F1DotF2(&src[0].fX, coeffX); 1036 formulate_F1DotF2(&src[0].fY, coeffY); 1037 1038 for (i = 0; i < 4; i++) 1039 coeffX[i] = SkFPAdd(coeffX[i],coeffY[i]); 1040 1041 SkScalar t[3]; 1042 int count = solve_cubic_polynomial(coeffX, t); 1043 int maxCount = 0; 1044 1045 // now remove extrema where the curvature is zero (mins) 1046 // !!!! need a test for this !!!! 1047 for (i = 0; i < count; i++) 1048 { 1049 // if (not_min_curvature()) 1050 if (t[i] > kMinTValueForChopping && t[i] < SK_Scalar1 - kMinTValueForChopping) 1051 tValues[maxCount++] = t[i]; 1052 } 1053 return maxCount; 1054} 1055 1056int SkChopCubicAtMaxCurvature(const SkPoint src[4], SkPoint dst[13], SkScalar tValues[3]) 1057{ 1058 SkScalar t_storage[3]; 1059 1060 if (tValues == NULL) 1061 tValues = t_storage; 1062 1063 int count = SkFindCubicMaxCurvature(src, tValues); 1064 1065 if (dst) { 1066 if (count == 0) 1067 memcpy(dst, src, 4 * sizeof(SkPoint)); 1068 else 1069 SkChopCubicAt(src, dst, tValues, count); 1070 } 1071 return count + 1; 1072} 1073 1074bool SkXRayCrossesMonotonicCubic(const SkXRay& pt, const SkPoint cubic[4], bool* ambiguous) { 1075 if (ambiguous) { 1076 *ambiguous = false; 1077 } 1078 1079 // Find the minimum and maximum y of the extrema, which are the 1080 // first and last points since this cubic is monotonic 1081 SkScalar min_y = SkMinScalar(cubic[0].fY, cubic[3].fY); 1082 SkScalar max_y = SkMaxScalar(cubic[0].fY, cubic[3].fY); 1083 1084 if (pt.fY == cubic[0].fY 1085 || pt.fY < min_y 1086 || pt.fY > max_y) { 1087 // The query line definitely does not cross the curve 1088 if (ambiguous) { 1089 *ambiguous = (pt.fY == cubic[0].fY); 1090 } 1091 return false; 1092 } 1093 1094 bool pt_at_extremum = (pt.fY == cubic[3].fY); 1095 1096 SkScalar min_x = 1097 SkMinScalar( 1098 SkMinScalar( 1099 SkMinScalar(cubic[0].fX, cubic[1].fX), 1100 cubic[2].fX), 1101 cubic[3].fX); 1102 if (pt.fX < min_x) { 1103 // The query line definitely crosses the curve 1104 if (ambiguous) { 1105 *ambiguous = pt_at_extremum; 1106 } 1107 return true; 1108 } 1109 1110 SkScalar max_x = 1111 SkMaxScalar( 1112 SkMaxScalar( 1113 SkMaxScalar(cubic[0].fX, cubic[1].fX), 1114 cubic[2].fX), 1115 cubic[3].fX); 1116 if (pt.fX > max_x) { 1117 // The query line definitely does not cross the curve 1118 return false; 1119 } 1120 1121 // Do a binary search to find the parameter value which makes y as 1122 // close as possible to the query point. See whether the query 1123 // line's origin is to the left of the associated x coordinate. 1124 1125 // kMaxIter is chosen as the number of mantissa bits for a float, 1126 // since there's no way we are going to get more precision by 1127 // iterating more times than that. 1128 const int kMaxIter = 23; 1129 SkPoint eval; 1130 int iter = 0; 1131 SkScalar upper_t; 1132 SkScalar lower_t; 1133 // Need to invert direction of t parameter if cubic goes up 1134 // instead of down 1135 if (cubic[3].fY > cubic[0].fY) { 1136 upper_t = SK_Scalar1; 1137 lower_t = SkFloatToScalar(0); 1138 } else { 1139 upper_t = SkFloatToScalar(0); 1140 lower_t = SK_Scalar1; 1141 } 1142 do { 1143 SkScalar t = SkScalarAve(upper_t, lower_t); 1144 SkEvalCubicAt(cubic, t, &eval, NULL, NULL); 1145 if (pt.fY > eval.fY) { 1146 lower_t = t; 1147 } else { 1148 upper_t = t; 1149 } 1150 } while (++iter < kMaxIter 1151 && !SkScalarNearlyZero(eval.fY - pt.fY)); 1152 if (pt.fX <= eval.fX) { 1153 if (ambiguous) { 1154 *ambiguous = pt_at_extremum; 1155 } 1156 return true; 1157 } 1158 return false; 1159} 1160 1161int SkNumXRayCrossingsForCubic(const SkXRay& pt, const SkPoint cubic[4], bool* ambiguous) { 1162 int num_crossings = 0; 1163 SkPoint monotonic_cubics[10]; 1164 int num_monotonic_cubics = SkChopCubicAtYExtrema(cubic, monotonic_cubics); 1165 if (ambiguous) { 1166 *ambiguous = false; 1167 } 1168 bool locally_ambiguous; 1169 if (SkXRayCrossesMonotonicCubic(pt, &monotonic_cubics[0], &locally_ambiguous)) 1170 ++num_crossings; 1171 if (ambiguous) { 1172 *ambiguous |= locally_ambiguous; 1173 } 1174 if (num_monotonic_cubics > 0) 1175 if (SkXRayCrossesMonotonicCubic(pt, &monotonic_cubics[3], &locally_ambiguous)) 1176 ++num_crossings; 1177 if (ambiguous) { 1178 *ambiguous |= locally_ambiguous; 1179 } 1180 if (num_monotonic_cubics > 1) 1181 if (SkXRayCrossesMonotonicCubic(pt, &monotonic_cubics[6], &locally_ambiguous)) 1182 ++num_crossings; 1183 if (ambiguous) { 1184 *ambiguous |= locally_ambiguous; 1185 } 1186 return num_crossings; 1187} 1188//////////////////////////////////////////////////////////////////////////////// 1189 1190/* Find t value for quadratic [a, b, c] = d. 1191 Return 0 if there is no solution within [0, 1) 1192*/ 1193static SkScalar quad_solve(SkScalar a, SkScalar b, SkScalar c, SkScalar d) 1194{ 1195 // At^2 + Bt + C = d 1196 SkScalar A = a - 2 * b + c; 1197 SkScalar B = 2 * (b - a); 1198 SkScalar C = a - d; 1199 1200 SkScalar roots[2]; 1201 int count = SkFindUnitQuadRoots(A, B, C, roots); 1202 1203 SkASSERT(count <= 1); 1204 return count == 1 ? roots[0] : 0; 1205} 1206 1207/* given a quad-curve and a point (x,y), chop the quad at that point and place 1208 the new off-curve point and endpoint into 'dest'. The new end point is used 1209 (rather than (x,y)) to compensate for numerical inaccuracies. 1210 Should only return false if the computed pos is the start of the curve 1211 (i.e. root == 0) 1212*/ 1213static bool truncate_last_curve(const SkPoint quad[3], SkScalar x, SkScalar y, SkPoint* dest) 1214{ 1215 const SkScalar* base; 1216 SkScalar value; 1217 1218 if (SkScalarAbs(x) < SkScalarAbs(y)) { 1219 base = &quad[0].fX; 1220 value = x; 1221 } else { 1222 base = &quad[0].fY; 1223 value = y; 1224 } 1225 1226 // note: this returns 0 if it thinks value is out of range, meaning the 1227 // root might return something outside of [0, 1) 1228 SkScalar t = quad_solve(base[0], base[2], base[4], value); 1229 1230 if (t > 0) 1231 { 1232 SkPoint tmp[5]; 1233 SkChopQuadAt(quad, tmp, t); 1234 dest[0] = tmp[1]; 1235 dest[1] = tmp[2]; 1236 return true; 1237 } else { 1238 /* t == 0 means either the value triggered a root outside of [0, 1) 1239 For our purposes, we can ignore the <= 0 roots, but we want to 1240 catch the >= 1 roots (which given our caller, will basically mean 1241 a root of 1, give-or-take numerical instability). If we are in the 1242 >= 1 case, return the existing offCurve point. 1243 1244 The test below checks to see if we are close to the "end" of the 1245 curve (near base[4]). Rather than specifying a tolerance, I just 1246 check to see if value is on to the right/left of the middle point 1247 (depending on the direction/sign of the end points). 1248 */ 1249 if ((base[0] < base[4] && value > base[2]) || 1250 (base[0] > base[4] && value < base[2])) // should root have been 1 1251 { 1252 dest[0] = quad[1]; 1253 dest[1].set(x, y); 1254 return true; 1255 } 1256 } 1257 return false; 1258} 1259 1260#ifdef SK_SCALAR_IS_FLOAT 1261 1262// Due to floating point issues (i.e., 1.0f - SK_ScalarRoot2Over2 != 1263// SK_ScalarRoot2Over2 - SK_ScalarTanPIOver8) a cruder root2over2 1264// approximation is required to make the quad circle points convex. The 1265// root of the problem is that with the root2over2 value in SkScalar.h 1266// the arcs really are ever so slightly concave. Some alternative fixes 1267// to this problem (besides just arbitrarily pushing out the mid-point as 1268// is done here) are: 1269// Adjust all the points (not just the middle) to both better approximate 1270// the curve and remain convex 1271// Switch over to using cubics rather then quads 1272// Use a different method to create the mid-point (e.g., compute 1273// the two side points, average them, then move it out as needed 1274#define SK_ScalarRoot2Over2_QuadCircle 0.7072f 1275 1276#else 1277 #define SK_ScalarRoot2Over2_QuadCircle SK_FixedRoot2Over2 1278#endif 1279 1280 1281static const SkPoint gQuadCirclePts[kSkBuildQuadArcStorage] = { 1282 { SK_Scalar1, 0 }, 1283 { SK_Scalar1, SK_ScalarTanPIOver8 }, 1284 { SK_ScalarRoot2Over2_QuadCircle, SK_ScalarRoot2Over2_QuadCircle }, 1285 { SK_ScalarTanPIOver8, SK_Scalar1 }, 1286 1287 { 0, SK_Scalar1 }, 1288 { -SK_ScalarTanPIOver8, SK_Scalar1 }, 1289 { -SK_ScalarRoot2Over2_QuadCircle, SK_ScalarRoot2Over2_QuadCircle }, 1290 { -SK_Scalar1, SK_ScalarTanPIOver8 }, 1291 1292 { -SK_Scalar1, 0 }, 1293 { -SK_Scalar1, -SK_ScalarTanPIOver8 }, 1294 { -SK_ScalarRoot2Over2_QuadCircle, -SK_ScalarRoot2Over2_QuadCircle }, 1295 { -SK_ScalarTanPIOver8, -SK_Scalar1 }, 1296 1297 { 0, -SK_Scalar1 }, 1298 { SK_ScalarTanPIOver8, -SK_Scalar1 }, 1299 { SK_ScalarRoot2Over2_QuadCircle, -SK_ScalarRoot2Over2_QuadCircle }, 1300 { SK_Scalar1, -SK_ScalarTanPIOver8 }, 1301 1302 { SK_Scalar1, 0 } 1303}; 1304 1305int SkBuildQuadArc(const SkVector& uStart, const SkVector& uStop, 1306 SkRotationDirection dir, const SkMatrix* userMatrix, 1307 SkPoint quadPoints[]) 1308{ 1309 // rotate by x,y so that uStart is (1.0) 1310 SkScalar x = SkPoint::DotProduct(uStart, uStop); 1311 SkScalar y = SkPoint::CrossProduct(uStart, uStop); 1312 1313 SkScalar absX = SkScalarAbs(x); 1314 SkScalar absY = SkScalarAbs(y); 1315 1316 int pointCount; 1317 1318 // check for (effectively) coincident vectors 1319 // this can happen if our angle is nearly 0 or nearly 180 (y == 0) 1320 // ... we use the dot-prod to distinguish between 0 and 180 (x > 0) 1321 if (absY <= SK_ScalarNearlyZero && x > 0 && 1322 ((y >= 0 && kCW_SkRotationDirection == dir) || 1323 (y <= 0 && kCCW_SkRotationDirection == dir))) { 1324 1325 // just return the start-point 1326 quadPoints[0].set(SK_Scalar1, 0); 1327 pointCount = 1; 1328 } else { 1329 if (dir == kCCW_SkRotationDirection) 1330 y = -y; 1331 1332 // what octant (quadratic curve) is [xy] in? 1333 int oct = 0; 1334 bool sameSign = true; 1335 1336 if (0 == y) 1337 { 1338 oct = 4; // 180 1339 SkASSERT(SkScalarAbs(x + SK_Scalar1) <= SK_ScalarNearlyZero); 1340 } 1341 else if (0 == x) 1342 { 1343 SkASSERT(absY - SK_Scalar1 <= SK_ScalarNearlyZero); 1344 if (y > 0) 1345 oct = 2; // 90 1346 else 1347 oct = 6; // 270 1348 } 1349 else 1350 { 1351 if (y < 0) 1352 oct += 4; 1353 if ((x < 0) != (y < 0)) 1354 { 1355 oct += 2; 1356 sameSign = false; 1357 } 1358 if ((absX < absY) == sameSign) 1359 oct += 1; 1360 } 1361 1362 int wholeCount = oct << 1; 1363 memcpy(quadPoints, gQuadCirclePts, (wholeCount + 1) * sizeof(SkPoint)); 1364 1365 const SkPoint* arc = &gQuadCirclePts[wholeCount]; 1366 if (truncate_last_curve(arc, x, y, &quadPoints[wholeCount + 1])) 1367 { 1368 wholeCount += 2; 1369 } 1370 pointCount = wholeCount + 1; 1371 } 1372 1373 // now handle counter-clockwise and the initial unitStart rotation 1374 SkMatrix matrix; 1375 matrix.setSinCos(uStart.fY, uStart.fX); 1376 if (dir == kCCW_SkRotationDirection) { 1377 matrix.preScale(SK_Scalar1, -SK_Scalar1); 1378 } 1379 if (userMatrix) { 1380 matrix.postConcat(*userMatrix); 1381 } 1382 matrix.mapPoints(quadPoints, pointCount); 1383 return pointCount; 1384} 1385 1386/////////////////////////////////////////////////////////////////////////////// 1387 1388// F = (A (1 - t)^2 + C t^2 + 2 B (1 - t) t w) 1389// ------------------------------------------ 1390// ((1 - t)^2 + t^2 + 2 (1 - t) t w) 1391// 1392// = {t^2 (P0 + P2 - 2 P1 w), t (-2 P0 + 2 P1 w), P0} 1393// ------------------------------------------------ 1394// {t^2 (2 - 2 w), t (-2 + 2 w), 1} 1395// 1396 1397// Take the parametric specification for the conic (either X or Y) and return 1398// in coeff[] the coefficients for the simple quadratic polynomial 1399// coeff[0] for t^2 1400// coeff[1] for t 1401// coeff[2] for constant term 1402// 1403static SkScalar conic_eval_pos(const SkScalar src[], SkScalar w, SkScalar t) { 1404 SkASSERT(src); 1405 SkASSERT(t >= 0 && t <= SK_Scalar1); 1406 1407 SkScalar src2w = SkScalarMul(src[2], w); 1408 SkScalar C = src[0]; 1409 SkScalar A = src[4] - 2 * src2w + C; 1410 SkScalar B = 2 * (src2w - C); 1411 SkScalar numer = SkScalarMulAdd(SkScalarMulAdd(A, t, B), t, C); 1412 1413 B = 2 * (w - SK_Scalar1); 1414 C = SK_Scalar1; 1415 A = -B; 1416 SkScalar denom = SkScalarMulAdd(SkScalarMulAdd(A, t, B), t, C); 1417 1418 return SkScalarDiv(numer, denom); 1419} 1420 1421// F' = 2 (C t (1 + t (-1 + w)) - A (-1 + t) (t (-1 + w) - w) + B (1 - 2 t) w) 1422// 1423// t^2 : (2 P0 - 2 P2 - 2 P0 w + 2 P2 w) 1424// t^1 : (-2 P0 + 2 P2 + 4 P0 w - 4 P1 w) 1425// t^0 : -2 P0 w + 2 P1 w 1426// 1427// We disregard magnitude, so we can freely ignore the denominator of F', and 1428// divide the numerator by 2 1429// 1430// coeff[0] for t^2 1431// coeff[1] for t^1 1432// coeff[2] for t^0 1433// 1434static void conic_deriv_coeff(const SkScalar src[], SkScalar w, SkScalar coeff[3]) { 1435 const SkScalar P20 = src[4] - src[0]; 1436 const SkScalar P10 = src[2] - src[0]; 1437 const SkScalar wP10 = w * P10; 1438 coeff[0] = w * P20 - P20; 1439 coeff[1] = P20 - 2 * wP10; 1440 coeff[2] = wP10; 1441} 1442 1443static SkScalar conic_eval_tan(const SkScalar coord[], SkScalar w, SkScalar t) { 1444 SkScalar coeff[3]; 1445 conic_deriv_coeff(coord, w, coeff); 1446 return t * (t * coeff[0] + coeff[1]) + coeff[2]; 1447} 1448 1449static bool conic_find_extrema(const SkScalar src[], SkScalar w, SkScalar* t) { 1450 SkScalar coeff[3]; 1451 conic_deriv_coeff(src, w, coeff); 1452 1453 SkScalar tValues[2]; 1454 int roots = SkFindUnitQuadRoots(coeff[0], coeff[1], coeff[2], tValues); 1455 SkASSERT(0 == roots || 1 == roots); 1456 1457 if (1 == roots) { 1458 *t = tValues[0]; 1459 return true; 1460 } 1461 return false; 1462} 1463 1464struct SkP3D { 1465 SkScalar fX, fY, fZ; 1466 1467 void set(SkScalar x, SkScalar y, SkScalar z) { 1468 fX = x; fY = y; fZ = z; 1469 } 1470 1471 void projectDown(SkPoint* dst) const { 1472 dst->set(fX / fZ, fY / fZ); 1473 } 1474}; 1475 1476// we just return the middle 3 points, since the first and last are dups of src 1477// 1478static void p3d_interp(const SkScalar src[3], SkScalar dst[3], SkScalar t) { 1479 SkScalar ab = SkScalarInterp(src[0], src[3], t); 1480 SkScalar bc = SkScalarInterp(src[3], src[6], t); 1481 dst[0] = ab; 1482 dst[3] = SkScalarInterp(ab, bc, t); 1483 dst[6] = bc; 1484} 1485 1486static void ratquad_mapTo3D(const SkPoint src[3], SkScalar w, SkP3D dst[]) { 1487 dst[0].set(src[0].fX * 1, src[0].fY * 1, 1); 1488 dst[1].set(src[1].fX * w, src[1].fY * w, w); 1489 dst[2].set(src[2].fX * 1, src[2].fY * 1, 1); 1490} 1491 1492void SkConic::evalAt(SkScalar t, SkPoint* pt, SkVector* tangent) const { 1493 SkASSERT(t >= 0 && t <= SK_Scalar1); 1494 1495 if (pt) { 1496 pt->set(conic_eval_pos(&fPts[0].fX, fW, t), 1497 conic_eval_pos(&fPts[0].fY, fW, t)); 1498 } 1499 if (tangent) { 1500 tangent->set(conic_eval_tan(&fPts[0].fX, fW, t), 1501 conic_eval_tan(&fPts[0].fY, fW, t)); 1502 } 1503} 1504 1505void SkConic::chopAt(SkScalar t, SkConic dst[2]) const { 1506 SkP3D tmp[3], tmp2[3]; 1507 1508 ratquad_mapTo3D(fPts, fW, tmp); 1509 1510 p3d_interp(&tmp[0].fX, &tmp2[0].fX, t); 1511 p3d_interp(&tmp[0].fY, &tmp2[0].fY, t); 1512 p3d_interp(&tmp[0].fZ, &tmp2[0].fZ, t); 1513 1514 dst[0].fPts[0] = fPts[0]; 1515 tmp2[0].projectDown(&dst[0].fPts[1]); 1516 tmp2[1].projectDown(&dst[0].fPts[2]); dst[1].fPts[0] = dst[0].fPts[2]; 1517 tmp2[2].projectDown(&dst[1].fPts[1]); 1518 dst[1].fPts[2] = fPts[2]; 1519 1520 // to put in "standard form", where w0 and w2 are both 1, we compute the 1521 // new w1 as sqrt(w1*w1/w0*w2) 1522 // or 1523 // w1 /= sqrt(w0*w2) 1524 // 1525 // However, in our case, we know that for dst[0], w0 == 1, and for dst[1], w2 == 1 1526 // 1527 SkScalar root = SkScalarSqrt(tmp2[1].fZ); 1528 dst[0].fW = tmp2[0].fZ / root; 1529 dst[1].fW = tmp2[2].fZ / root; 1530} 1531 1532static SkScalar subdivide_w_value(SkScalar w) { 1533 return SkScalarSqrt(SK_ScalarHalf + w * SK_ScalarHalf); 1534} 1535 1536void SkConic::chop(SkConic dst[2]) const { 1537 SkScalar scale = SkScalarInvert(SK_Scalar1 + fW); 1538 SkScalar p1x = fW * fPts[1].fX; 1539 SkScalar p1y = fW * fPts[1].fY; 1540 SkScalar mx = (fPts[0].fX + 2 * p1x + fPts[2].fX) * scale * SK_ScalarHalf; 1541 SkScalar my = (fPts[0].fY + 2 * p1y + fPts[2].fY) * scale * SK_ScalarHalf; 1542 1543 dst[0].fPts[0] = fPts[0]; 1544 dst[0].fPts[1].set((fPts[0].fX + p1x) * scale, 1545 (fPts[0].fY + p1y) * scale); 1546 dst[0].fPts[2].set(mx, my); 1547 1548 dst[1].fPts[0].set(mx, my); 1549 dst[1].fPts[1].set((p1x + fPts[2].fX) * scale, 1550 (p1y + fPts[2].fY) * scale); 1551 dst[1].fPts[2] = fPts[2]; 1552 1553 dst[0].fW = dst[1].fW = subdivide_w_value(fW); 1554} 1555 1556/* 1557 * "High order approximation of conic sections by quadratic splines" 1558 * by Michael Floater, 1993 1559 */ 1560#define AS_QUAD_ERROR_SETUP \ 1561 SkScalar a = fW - 1; \ 1562 SkScalar k = a / (4 * (2 + a)); \ 1563 SkScalar x = k * (fPts[0].fX - 2 * fPts[1].fX + fPts[2].fX); \ 1564 SkScalar y = k * (fPts[0].fY - 2 * fPts[1].fY + fPts[2].fY); 1565 1566void SkConic::computeAsQuadError(SkVector* err) const { 1567 AS_QUAD_ERROR_SETUP 1568 err->set(x, y); 1569} 1570 1571bool SkConic::asQuadTol(SkScalar tol) const { 1572 AS_QUAD_ERROR_SETUP 1573 return (x * x + y * y) <= tol * tol; 1574} 1575 1576int SkConic::computeQuadPOW2(SkScalar tol) const { 1577 AS_QUAD_ERROR_SETUP 1578 SkScalar error = SkScalarSqrt(x * x + y * y) - tol; 1579 1580 if (error <= 0) { 1581 return 0; 1582 } 1583 uint32_t ierr = (uint32_t)error; 1584 return (34 - SkCLZ(ierr)) >> 1; 1585} 1586 1587static SkPoint* subdivide(const SkConic& src, SkPoint pts[], int level) { 1588 SkASSERT(level >= 0); 1589 1590 if (0 == level) { 1591 memcpy(pts, &src.fPts[1], 2 * sizeof(SkPoint)); 1592 return pts + 2; 1593 } else { 1594 SkConic dst[2]; 1595 src.chop(dst); 1596 --level; 1597 pts = subdivide(dst[0], pts, level); 1598 return subdivide(dst[1], pts, level); 1599 } 1600} 1601 1602int SkConic::chopIntoQuadsPOW2(SkPoint pts[], int pow2) const { 1603 SkASSERT(pow2 >= 0); 1604 *pts = fPts[0]; 1605 SkDEBUGCODE(SkPoint* endPts =) subdivide(*this, pts + 1, pow2); 1606 SkASSERT(endPts - pts == (2 * (1 << pow2) + 1)); 1607 return 1 << pow2; 1608} 1609 1610bool SkConic::findXExtrema(SkScalar* t) const { 1611 return conic_find_extrema(&fPts[0].fX, fW, t); 1612} 1613 1614bool SkConic::findYExtrema(SkScalar* t) const { 1615 return conic_find_extrema(&fPts[0].fY, fW, t); 1616} 1617 1618bool SkConic::chopAtXExtrema(SkConic dst[2]) const { 1619 SkScalar t; 1620 if (this->findXExtrema(&t)) { 1621 this->chopAt(t, dst); 1622 // now clean-up the middle, since we know t was meant to be at 1623 // an X-extrema 1624 SkScalar value = dst[0].fPts[2].fX; 1625 dst[0].fPts[1].fX = value; 1626 dst[1].fPts[0].fX = value; 1627 dst[1].fPts[1].fX = value; 1628 return true; 1629 } 1630 return false; 1631} 1632 1633bool SkConic::chopAtYExtrema(SkConic dst[2]) const { 1634 SkScalar t; 1635 if (this->findYExtrema(&t)) { 1636 this->chopAt(t, dst); 1637 // now clean-up the middle, since we know t was meant to be at 1638 // an Y-extrema 1639 SkScalar value = dst[0].fPts[2].fY; 1640 dst[0].fPts[1].fY = value; 1641 dst[1].fPts[0].fY = value; 1642 dst[1].fPts[1].fY = value; 1643 return true; 1644 } 1645 return false; 1646} 1647 1648void SkConic::computeTightBounds(SkRect* bounds) const { 1649 SkPoint pts[4]; 1650 pts[0] = fPts[0]; 1651 pts[1] = fPts[2]; 1652 int count = 2; 1653 1654 SkScalar t; 1655 if (this->findXExtrema(&t)) { 1656 this->evalAt(t, &pts[count++]); 1657 } 1658 if (this->findYExtrema(&t)) { 1659 this->evalAt(t, &pts[count++]); 1660 } 1661 bounds->set(pts, count); 1662} 1663 1664void SkConic::computeFastBounds(SkRect* bounds) const { 1665 bounds->set(fPts, 3); 1666} 1667