SkGeometry.cpp revision 4bb50b22fce5c4031549976cf7b04d63cc22e624
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// 397int SkChopQuadAtMaxCurvature(const SkPoint src[3], SkPoint dst[5]) 398{ 399 SkScalar Ax = src[1].fX - src[0].fX; 400 SkScalar Ay = src[1].fY - src[0].fY; 401 SkScalar Bx = src[0].fX - src[1].fX - src[1].fX + src[2].fX; 402 SkScalar By = src[0].fY - src[1].fY - src[1].fY + src[2].fY; 403 SkScalar t = 0; // 0 means don't chop 404 405#ifdef SK_SCALAR_IS_FLOAT 406 (void)valid_unit_divide(-(Ax * Bx + Ay * By), Bx * Bx + By * By, &t); 407#else 408 // !!! should I use SkFloat here? seems like it 409 Sk64 numer, denom, tmp; 410 411 numer.setMul(Ax, -Bx); 412 tmp.setMul(Ay, -By); 413 numer.add(tmp); 414 415 if (numer.isPos()) // do nothing if numer <= 0 416 { 417 denom.setMul(Bx, Bx); 418 tmp.setMul(By, By); 419 denom.add(tmp); 420 SkASSERT(!denom.isNeg()); 421 if (numer < denom) 422 { 423 t = numer.getFixedDiv(denom); 424 SkASSERT(t >= 0 && t <= SK_Fixed1); // assert that we're numerically stable (ha!) 425 if ((unsigned)t >= SK_Fixed1) // runtime check for numerical stability 426 t = 0; // ignore the chop 427 } 428 } 429#endif 430 431 if (t == 0) 432 { 433 memcpy(dst, src, 3 * sizeof(SkPoint)); 434 return 1; 435 } 436 else 437 { 438 SkChopQuadAt(src, dst, t); 439 return 2; 440 } 441} 442 443#ifdef SK_SCALAR_IS_FLOAT 444 #define SK_ScalarTwoThirds (0.666666666f) 445#else 446 #define SK_ScalarTwoThirds ((SkFixed)(43691)) 447#endif 448 449void SkConvertQuadToCubic(const SkPoint src[3], SkPoint dst[4]) { 450 const SkScalar scale = SK_ScalarTwoThirds; 451 dst[0] = src[0]; 452 dst[1].set(src[0].fX + SkScalarMul(src[1].fX - src[0].fX, scale), 453 src[0].fY + SkScalarMul(src[1].fY - src[0].fY, scale)); 454 dst[2].set(src[2].fX + SkScalarMul(src[1].fX - src[2].fX, scale), 455 src[2].fY + SkScalarMul(src[1].fY - src[2].fY, scale)); 456 dst[3] = src[2]; 457} 458 459//////////////////////////////////////////////////////////////////////////////////////// 460///// CUBICS // CUBICS // CUBICS // CUBICS // CUBICS // CUBICS // CUBICS // CUBICS ///// 461//////////////////////////////////////////////////////////////////////////////////////// 462 463static void get_cubic_coeff(const SkScalar pt[], SkScalar coeff[4]) 464{ 465 coeff[0] = pt[6] + 3*(pt[2] - pt[4]) - pt[0]; 466 coeff[1] = 3*(pt[4] - pt[2] - pt[2] + pt[0]); 467 coeff[2] = 3*(pt[2] - pt[0]); 468 coeff[3] = pt[0]; 469} 470 471void SkGetCubicCoeff(const SkPoint pts[4], SkScalar cx[4], SkScalar cy[4]) 472{ 473 SkASSERT(pts); 474 475 if (cx) 476 get_cubic_coeff(&pts[0].fX, cx); 477 if (cy) 478 get_cubic_coeff(&pts[0].fY, cy); 479} 480 481static SkScalar eval_cubic(const SkScalar src[], SkScalar t) 482{ 483 SkASSERT(src); 484 SkASSERT(t >= 0 && t <= SK_Scalar1); 485 486 if (t == 0) 487 return src[0]; 488 489#ifdef DIRECT_EVAL_OF_POLYNOMIALS 490 SkScalar D = src[0]; 491 SkScalar A = src[6] + 3*(src[2] - src[4]) - D; 492 SkScalar B = 3*(src[4] - src[2] - src[2] + D); 493 SkScalar C = 3*(src[2] - D); 494 495 return SkScalarMulAdd(SkScalarMulAdd(SkScalarMulAdd(A, t, B), t, C), t, D); 496#else 497 SkScalar ab = SkScalarInterp(src[0], src[2], t); 498 SkScalar bc = SkScalarInterp(src[2], src[4], t); 499 SkScalar cd = SkScalarInterp(src[4], src[6], t); 500 SkScalar abc = SkScalarInterp(ab, bc, t); 501 SkScalar bcd = SkScalarInterp(bc, cd, t); 502 return SkScalarInterp(abc, bcd, t); 503#endif 504} 505 506/** return At^2 + Bt + C 507*/ 508static SkScalar eval_quadratic(SkScalar A, SkScalar B, SkScalar C, SkScalar t) 509{ 510 SkASSERT(t >= 0 && t <= SK_Scalar1); 511 512 return SkScalarMulAdd(SkScalarMulAdd(A, t, B), t, C); 513} 514 515static SkScalar eval_cubic_derivative(const SkScalar src[], SkScalar t) 516{ 517 SkScalar A = src[6] + 3*(src[2] - src[4]) - src[0]; 518 SkScalar B = 2*(src[4] - 2 * src[2] + src[0]); 519 SkScalar C = src[2] - src[0]; 520 521 return eval_quadratic(A, B, C, t); 522} 523 524static SkScalar eval_cubic_2ndDerivative(const SkScalar src[], SkScalar t) 525{ 526 SkScalar A = src[6] + 3*(src[2] - src[4]) - src[0]; 527 SkScalar B = src[4] - 2 * src[2] + src[0]; 528 529 return SkScalarMulAdd(A, t, B); 530} 531 532void SkEvalCubicAt(const SkPoint src[4], SkScalar t, SkPoint* loc, SkVector* tangent, SkVector* curvature) 533{ 534 SkASSERT(src); 535 SkASSERT(t >= 0 && t <= SK_Scalar1); 536 537 if (loc) 538 loc->set(eval_cubic(&src[0].fX, t), eval_cubic(&src[0].fY, t)); 539 if (tangent) 540 tangent->set(eval_cubic_derivative(&src[0].fX, t), 541 eval_cubic_derivative(&src[0].fY, t)); 542 if (curvature) 543 curvature->set(eval_cubic_2ndDerivative(&src[0].fX, t), 544 eval_cubic_2ndDerivative(&src[0].fY, t)); 545} 546 547/** Cubic'(t) = At^2 + Bt + C, where 548 A = 3(-a + 3(b - c) + d) 549 B = 6(a - 2b + c) 550 C = 3(b - a) 551 Solve for t, keeping only those that fit betwee 0 < t < 1 552*/ 553int SkFindCubicExtrema(SkScalar a, SkScalar b, SkScalar c, SkScalar d, SkScalar tValues[2]) 554{ 555#ifdef SK_SCALAR_IS_FIXED 556 if (!is_not_monotonic(a, b, c, d)) 557 return 0; 558#endif 559 560 // we divide A,B,C by 3 to simplify 561 SkScalar A = d - a + 3*(b - c); 562 SkScalar B = 2*(a - b - b + c); 563 SkScalar C = b - a; 564 565 return SkFindUnitQuadRoots(A, B, C, tValues); 566} 567 568static void interp_cubic_coords(const SkScalar* src, SkScalar* dst, SkScalar t) 569{ 570 SkScalar ab = SkScalarInterp(src[0], src[2], t); 571 SkScalar bc = SkScalarInterp(src[2], src[4], t); 572 SkScalar cd = SkScalarInterp(src[4], src[6], t); 573 SkScalar abc = SkScalarInterp(ab, bc, t); 574 SkScalar bcd = SkScalarInterp(bc, cd, t); 575 SkScalar abcd = SkScalarInterp(abc, bcd, t); 576 577 dst[0] = src[0]; 578 dst[2] = ab; 579 dst[4] = abc; 580 dst[6] = abcd; 581 dst[8] = bcd; 582 dst[10] = cd; 583 dst[12] = src[6]; 584} 585 586void SkChopCubicAt(const SkPoint src[4], SkPoint dst[7], SkScalar t) 587{ 588 SkASSERT(t > 0 && t < SK_Scalar1); 589 590 interp_cubic_coords(&src[0].fX, &dst[0].fX, t); 591 interp_cubic_coords(&src[0].fY, &dst[0].fY, t); 592} 593 594/* http://code.google.com/p/skia/issues/detail?id=32 595 596 This test code would fail when we didn't check the return result of 597 valid_unit_divide in SkChopCubicAt(... tValues[], int roots). The reason is 598 that after the first chop, the parameters to valid_unit_divide are equal 599 (thanks to finite float precision and rounding in the subtracts). Thus 600 even though the 2nd tValue looks < 1.0, after we renormalize it, we end 601 up with 1.0, hence the need to check and just return the last cubic as 602 a degenerate clump of 4 points in the sampe place. 603 604 static void test_cubic() { 605 SkPoint src[4] = { 606 { 556.25000, 523.03003 }, 607 { 556.23999, 522.96002 }, 608 { 556.21997, 522.89001 }, 609 { 556.21997, 522.82001 } 610 }; 611 SkPoint dst[10]; 612 SkScalar tval[] = { 0.33333334f, 0.99999994f }; 613 SkChopCubicAt(src, dst, tval, 2); 614 } 615 */ 616 617void SkChopCubicAt(const SkPoint src[4], SkPoint dst[], const SkScalar tValues[], int roots) 618{ 619#ifdef SK_DEBUG 620 { 621 for (int i = 0; i < roots - 1; i++) 622 { 623 SkASSERT(is_unit_interval(tValues[i])); 624 SkASSERT(is_unit_interval(tValues[i+1])); 625 SkASSERT(tValues[i] < tValues[i+1]); 626 } 627 } 628#endif 629 630 if (dst) 631 { 632 if (roots == 0) // nothing to chop 633 memcpy(dst, src, 4*sizeof(SkPoint)); 634 else 635 { 636 SkScalar t = tValues[0]; 637 SkPoint tmp[4]; 638 639 for (int i = 0; i < roots; i++) 640 { 641 SkChopCubicAt(src, dst, t); 642 if (i == roots - 1) 643 break; 644 645 dst += 3; 646 // have src point to the remaining cubic (after the chop) 647 memcpy(tmp, dst, 4 * sizeof(SkPoint)); 648 src = tmp; 649 650 // watch out in case the renormalized t isn't in range 651 if (!valid_unit_divide(tValues[i+1] - tValues[i], 652 SK_Scalar1 - tValues[i], &t)) { 653 // if we can't, just create a degenerate cubic 654 dst[4] = dst[5] = dst[6] = src[3]; 655 break; 656 } 657 } 658 } 659 } 660} 661 662void SkChopCubicAtHalf(const SkPoint src[4], SkPoint dst[7]) 663{ 664 SkScalar x01 = SkScalarAve(src[0].fX, src[1].fX); 665 SkScalar y01 = SkScalarAve(src[0].fY, src[1].fY); 666 SkScalar x12 = SkScalarAve(src[1].fX, src[2].fX); 667 SkScalar y12 = SkScalarAve(src[1].fY, src[2].fY); 668 SkScalar x23 = SkScalarAve(src[2].fX, src[3].fX); 669 SkScalar y23 = SkScalarAve(src[2].fY, src[3].fY); 670 671 SkScalar x012 = SkScalarAve(x01, x12); 672 SkScalar y012 = SkScalarAve(y01, y12); 673 SkScalar x123 = SkScalarAve(x12, x23); 674 SkScalar y123 = SkScalarAve(y12, y23); 675 676 dst[0] = src[0]; 677 dst[1].set(x01, y01); 678 dst[2].set(x012, y012); 679 dst[3].set(SkScalarAve(x012, x123), SkScalarAve(y012, y123)); 680 dst[4].set(x123, y123); 681 dst[5].set(x23, y23); 682 dst[6] = src[3]; 683} 684 685static void flatten_double_cubic_extrema(SkScalar coords[14]) 686{ 687 coords[4] = coords[8] = coords[6]; 688} 689 690/** Given 4 points on a cubic bezier, chop it into 1, 2, 3 beziers such that 691 the resulting beziers are monotonic in Y. This is called by the scan converter. 692 Depending on what is returned, dst[] is treated as follows 693 0 dst[0..3] is the original cubic 694 1 dst[0..3] and dst[3..6] are the two new cubics 695 2 dst[0..3], dst[3..6], dst[6..9] are the three new cubics 696 If dst == null, it is ignored and only the count is returned. 697*/ 698int SkChopCubicAtYExtrema(const SkPoint src[4], SkPoint dst[10]) { 699 SkScalar tValues[2]; 700 int roots = SkFindCubicExtrema(src[0].fY, src[1].fY, src[2].fY, 701 src[3].fY, tValues); 702 703 SkChopCubicAt(src, dst, tValues, roots); 704 if (dst && roots > 0) { 705 // we do some cleanup to ensure our Y extrema are flat 706 flatten_double_cubic_extrema(&dst[0].fY); 707 if (roots == 2) { 708 flatten_double_cubic_extrema(&dst[3].fY); 709 } 710 } 711 return roots; 712} 713 714int SkChopCubicAtXExtrema(const SkPoint src[4], SkPoint dst[10]) { 715 SkScalar tValues[2]; 716 int roots = SkFindCubicExtrema(src[0].fX, src[1].fX, src[2].fX, 717 src[3].fX, tValues); 718 719 SkChopCubicAt(src, dst, tValues, roots); 720 if (dst && roots > 0) { 721 // we do some cleanup to ensure our Y extrema are flat 722 flatten_double_cubic_extrema(&dst[0].fX); 723 if (roots == 2) { 724 flatten_double_cubic_extrema(&dst[3].fX); 725 } 726 } 727 return roots; 728} 729 730/** http://www.faculty.idc.ac.il/arik/quality/appendixA.html 731 732 Inflection means that curvature is zero. 733 Curvature is [F' x F''] / [F'^3] 734 So we solve F'x X F''y - F'y X F''y == 0 735 After some canceling of the cubic term, we get 736 A = b - a 737 B = c - 2b + a 738 C = d - 3c + 3b - a 739 (BxCy - ByCx)t^2 + (AxCy - AyCx)t + AxBy - AyBx == 0 740*/ 741int SkFindCubicInflections(const SkPoint src[4], SkScalar tValues[]) 742{ 743 SkScalar Ax = src[1].fX - src[0].fX; 744 SkScalar Ay = src[1].fY - src[0].fY; 745 SkScalar Bx = src[2].fX - 2 * src[1].fX + src[0].fX; 746 SkScalar By = src[2].fY - 2 * src[1].fY + src[0].fY; 747 SkScalar Cx = src[3].fX + 3 * (src[1].fX - src[2].fX) - src[0].fX; 748 SkScalar Cy = src[3].fY + 3 * (src[1].fY - src[2].fY) - src[0].fY; 749 int count; 750 751#ifdef SK_SCALAR_IS_FLOAT 752 count = SkFindUnitQuadRoots(Bx*Cy - By*Cx, Ax*Cy - Ay*Cx, Ax*By - Ay*Bx, tValues); 753#else 754 Sk64 A, B, C, tmp; 755 756 A.setMul(Bx, Cy); 757 tmp.setMul(By, Cx); 758 A.sub(tmp); 759 760 B.setMul(Ax, Cy); 761 tmp.setMul(Ay, Cx); 762 B.sub(tmp); 763 764 C.setMul(Ax, By); 765 tmp.setMul(Ay, Bx); 766 C.sub(tmp); 767 768 count = Sk64FindFixedQuadRoots(A, B, C, tValues); 769#endif 770 771 return count; 772} 773 774int SkChopCubicAtInflections(const SkPoint src[], SkPoint dst[10]) 775{ 776 SkScalar tValues[2]; 777 int count = SkFindCubicInflections(src, tValues); 778 779 if (dst) 780 { 781 if (count == 0) 782 memcpy(dst, src, 4 * sizeof(SkPoint)); 783 else 784 SkChopCubicAt(src, dst, tValues, count); 785 } 786 return count + 1; 787} 788 789template <typename T> void bubble_sort(T array[], int count) 790{ 791 for (int i = count - 1; i > 0; --i) 792 for (int j = i; j > 0; --j) 793 if (array[j] < array[j-1]) 794 { 795 T tmp(array[j]); 796 array[j] = array[j-1]; 797 array[j-1] = tmp; 798 } 799} 800 801#include "SkFP.h" 802 803// newton refinement 804#if 0 805static SkScalar refine_cubic_root(const SkFP coeff[4], SkScalar root) 806{ 807 // x1 = x0 - f(t) / f'(t) 808 809 SkFP T = SkScalarToFloat(root); 810 SkFP N, D; 811 812 // f' = 3*coeff[0]*T^2 + 2*coeff[1]*T + coeff[2] 813 D = SkFPMul(SkFPMul(coeff[0], SkFPMul(T,T)), 3); 814 D = SkFPAdd(D, SkFPMulInt(SkFPMul(coeff[1], T), 2)); 815 D = SkFPAdd(D, coeff[2]); 816 817 if (D == 0) 818 return root; 819 820 // f = coeff[0]*T^3 + coeff[1]*T^2 + coeff[2]*T + coeff[3] 821 N = SkFPMul(SkFPMul(SkFPMul(T, T), T), coeff[0]); 822 N = SkFPAdd(N, SkFPMul(SkFPMul(T, T), coeff[1])); 823 N = SkFPAdd(N, SkFPMul(T, coeff[2])); 824 N = SkFPAdd(N, coeff[3]); 825 826 if (N) 827 { 828 SkScalar delta = SkFPToScalar(SkFPDiv(N, D)); 829 830 if (delta) 831 root -= delta; 832 } 833 return root; 834} 835#endif 836 837/** 838 * Given an array and count, remove all pair-wise duplicates from the array, 839 * keeping the existing sorting, and return the new count 840 */ 841static int collaps_duplicates(float array[], int count) { 842 for (int n = count; n > 1; --n) { 843 if (array[0] == array[1]) { 844 for (int i = 1; i < n; ++i) { 845 array[i - 1] = array[i]; 846 } 847 count -= 1; 848 } else { 849 array += 1; 850 } 851 } 852 return count; 853} 854 855#ifdef SK_DEBUG 856 857#define TEST_COLLAPS_ENTRY(array) array, SK_ARRAY_COUNT(array) 858 859static void test_collaps_duplicates() { 860 static bool gOnce; 861 if (gOnce) { return; } 862 gOnce = true; 863 const float src0[] = { 0 }; 864 const float src1[] = { 0, 0 }; 865 const float src2[] = { 0, 1 }; 866 const float src3[] = { 0, 0, 0 }; 867 const float src4[] = { 0, 0, 1 }; 868 const float src5[] = { 0, 1, 1 }; 869 const float src6[] = { 0, 1, 2 }; 870 const struct { 871 const float* fData; 872 int fCount; 873 int fCollapsedCount; 874 } data[] = { 875 { TEST_COLLAPS_ENTRY(src0), 1 }, 876 { TEST_COLLAPS_ENTRY(src1), 1 }, 877 { TEST_COLLAPS_ENTRY(src2), 2 }, 878 { TEST_COLLAPS_ENTRY(src3), 1 }, 879 { TEST_COLLAPS_ENTRY(src4), 2 }, 880 { TEST_COLLAPS_ENTRY(src5), 2 }, 881 { TEST_COLLAPS_ENTRY(src6), 3 }, 882 }; 883 for (size_t i = 0; i < SK_ARRAY_COUNT(data); ++i) { 884 float dst[3]; 885 memcpy(dst, data[i].fData, data[i].fCount * sizeof(dst[0])); 886 int count = collaps_duplicates(dst, data[i].fCount); 887 SkASSERT(data[i].fCollapsedCount == count); 888 for (int j = 1; j < count; ++j) { 889 SkASSERT(dst[j-1] < dst[j]); 890 } 891 } 892} 893#endif 894 895#if defined _WIN32 && _MSC_VER >= 1300 && defined SK_SCALAR_IS_FIXED // disable warning : unreachable code if building fixed point for windows desktop 896#pragma warning ( disable : 4702 ) 897#endif 898 899/* Solve coeff(t) == 0, returning the number of roots that 900 lie withing 0 < t < 1. 901 coeff[0]t^3 + coeff[1]t^2 + coeff[2]t + coeff[3] 902 903 Eliminates repeated roots (so that all tValues are distinct, and are always 904 in increasing order. 905*/ 906static int solve_cubic_polynomial(const SkFP coeff[4], SkScalar tValues[3]) 907{ 908#ifndef SK_SCALAR_IS_FLOAT 909 return 0; // this is not yet implemented for software float 910#endif 911 912 if (SkScalarNearlyZero(coeff[0])) // we're just a quadratic 913 { 914 return SkFindUnitQuadRoots(coeff[1], coeff[2], coeff[3], tValues); 915 } 916 917 SkFP a, b, c, Q, R; 918 919 { 920 SkASSERT(coeff[0] != 0); 921 922 SkFP inva = SkFPInvert(coeff[0]); 923 a = SkFPMul(coeff[1], inva); 924 b = SkFPMul(coeff[2], inva); 925 c = SkFPMul(coeff[3], inva); 926 } 927 Q = SkFPDivInt(SkFPSub(SkFPMul(a,a), SkFPMulInt(b, 3)), 9); 928// R = (2*a*a*a - 9*a*b + 27*c) / 54; 929 R = SkFPMulInt(SkFPMul(SkFPMul(a, a), a), 2); 930 R = SkFPSub(R, SkFPMulInt(SkFPMul(a, b), 9)); 931 R = SkFPAdd(R, SkFPMulInt(c, 27)); 932 R = SkFPDivInt(R, 54); 933 934 SkFP Q3 = SkFPMul(SkFPMul(Q, Q), Q); 935 SkFP R2MinusQ3 = SkFPSub(SkFPMul(R,R), Q3); 936 SkFP adiv3 = SkFPDivInt(a, 3); 937 938 SkScalar* roots = tValues; 939 SkScalar r; 940 941 if (SkFPLT(R2MinusQ3, 0)) // we have 3 real roots 942 { 943#ifdef SK_SCALAR_IS_FLOAT 944 float theta = sk_float_acos(R / sk_float_sqrt(Q3)); 945 float neg2RootQ = -2 * sk_float_sqrt(Q); 946 947 r = neg2RootQ * sk_float_cos(theta/3) - adiv3; 948 if (is_unit_interval(r)) 949 *roots++ = r; 950 951 r = neg2RootQ * sk_float_cos((theta + 2*SK_ScalarPI)/3) - adiv3; 952 if (is_unit_interval(r)) 953 *roots++ = r; 954 955 r = neg2RootQ * sk_float_cos((theta - 2*SK_ScalarPI)/3) - adiv3; 956 if (is_unit_interval(r)) 957 *roots++ = r; 958 959 SkDEBUGCODE(test_collaps_duplicates();) 960 961 // now sort the roots 962 int count = (int)(roots - tValues); 963 SkASSERT((unsigned)count <= 3); 964 bubble_sort(tValues, count); 965 count = collaps_duplicates(tValues, count); 966 roots = tValues + count; // so we compute the proper count below 967#endif 968 } 969 else // we have 1 real root 970 { 971 SkFP A = SkFPAdd(SkFPAbs(R), SkFPSqrt(R2MinusQ3)); 972 A = SkFPCubeRoot(A); 973 if (SkFPGT(R, 0)) 974 A = SkFPNeg(A); 975 976 if (A != 0) 977 A = SkFPAdd(A, SkFPDiv(Q, A)); 978 r = SkFPToScalar(SkFPSub(A, adiv3)); 979 if (is_unit_interval(r)) 980 *roots++ = r; 981 } 982 983 return (int)(roots - tValues); 984} 985 986/* Looking for F' dot F'' == 0 987 988 A = b - a 989 B = c - 2b + a 990 C = d - 3c + 3b - a 991 992 F' = 3Ct^2 + 6Bt + 3A 993 F'' = 6Ct + 6B 994 995 F' dot F'' -> CCt^3 + 3BCt^2 + (2BB + CA)t + AB 996*/ 997static void formulate_F1DotF2(const SkScalar src[], SkFP coeff[4]) 998{ 999 SkScalar a = src[2] - src[0]; 1000 SkScalar b = src[4] - 2 * src[2] + src[0]; 1001 SkScalar c = src[6] + 3 * (src[2] - src[4]) - src[0]; 1002 1003 SkFP A = SkScalarToFP(a); 1004 SkFP B = SkScalarToFP(b); 1005 SkFP C = SkScalarToFP(c); 1006 1007 coeff[0] = SkFPMul(C, C); 1008 coeff[1] = SkFPMulInt(SkFPMul(B, C), 3); 1009 coeff[2] = SkFPMulInt(SkFPMul(B, B), 2); 1010 coeff[2] = SkFPAdd(coeff[2], SkFPMul(C, A)); 1011 coeff[3] = SkFPMul(A, B); 1012} 1013 1014// EXPERIMENTAL: can set this to zero to accept all t-values 0 < t < 1 1015//#define kMinTValueForChopping (SK_Scalar1 / 256) 1016#define kMinTValueForChopping 0 1017 1018/* Looking for F' dot F'' == 0 1019 1020 A = b - a 1021 B = c - 2b + a 1022 C = d - 3c + 3b - a 1023 1024 F' = 3Ct^2 + 6Bt + 3A 1025 F'' = 6Ct + 6B 1026 1027 F' dot F'' -> CCt^3 + 3BCt^2 + (2BB + CA)t + AB 1028*/ 1029int SkFindCubicMaxCurvature(const SkPoint src[4], SkScalar tValues[3]) 1030{ 1031 SkFP coeffX[4], coeffY[4]; 1032 int i; 1033 1034 formulate_F1DotF2(&src[0].fX, coeffX); 1035 formulate_F1DotF2(&src[0].fY, coeffY); 1036 1037 for (i = 0; i < 4; i++) 1038 coeffX[i] = SkFPAdd(coeffX[i],coeffY[i]); 1039 1040 SkScalar t[3]; 1041 int count = solve_cubic_polynomial(coeffX, t); 1042 int maxCount = 0; 1043 1044 // now remove extrema where the curvature is zero (mins) 1045 // !!!! need a test for this !!!! 1046 for (i = 0; i < count; i++) 1047 { 1048 // if (not_min_curvature()) 1049 if (t[i] > kMinTValueForChopping && t[i] < SK_Scalar1 - kMinTValueForChopping) 1050 tValues[maxCount++] = t[i]; 1051 } 1052 return maxCount; 1053} 1054 1055int SkChopCubicAtMaxCurvature(const SkPoint src[4], SkPoint dst[13], SkScalar tValues[3]) 1056{ 1057 SkScalar t_storage[3]; 1058 1059 if (tValues == NULL) 1060 tValues = t_storage; 1061 1062 int count = SkFindCubicMaxCurvature(src, tValues); 1063 1064 if (dst) 1065 { 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 1191/* Find t value for quadratic [a, b, c] = d. 1192 Return 0 if there is no solution within [0, 1) 1193*/ 1194static SkScalar quad_solve(SkScalar a, SkScalar b, SkScalar c, SkScalar d) 1195{ 1196 // At^2 + Bt + C = d 1197 SkScalar A = a - 2 * b + c; 1198 SkScalar B = 2 * (b - a); 1199 SkScalar C = a - d; 1200 1201 SkScalar roots[2]; 1202 int count = SkFindUnitQuadRoots(A, B, C, roots); 1203 1204 SkASSERT(count <= 1); 1205 return count == 1 ? roots[0] : 0; 1206} 1207 1208/* given a quad-curve and a point (x,y), chop the quad at that point and return 1209 the new quad's offCurve point. Should only return false if the computed pos 1210 is the start of the curve (i.e. root == 0) 1211*/ 1212static bool quad_pt2OffCurve(const SkPoint quad[3], SkScalar x, SkScalar y, SkPoint* offCurve) 1213{ 1214 const SkScalar* base; 1215 SkScalar value; 1216 1217 if (SkScalarAbs(x) < SkScalarAbs(y)) { 1218 base = &quad[0].fX; 1219 value = x; 1220 } else { 1221 base = &quad[0].fY; 1222 value = y; 1223 } 1224 1225 // note: this returns 0 if it thinks value is out of range, meaning the 1226 // root might return something outside of [0, 1) 1227 SkScalar t = quad_solve(base[0], base[2], base[4], value); 1228 1229 if (t > 0) 1230 { 1231 SkPoint tmp[5]; 1232 SkChopQuadAt(quad, tmp, t); 1233 *offCurve = tmp[1]; 1234 return true; 1235 } else { 1236 /* t == 0 means either the value triggered a root outside of [0, 1) 1237 For our purposes, we can ignore the <= 0 roots, but we want to 1238 catch the >= 1 roots (which given our caller, will basically mean 1239 a root of 1, give-or-take numerical instability). If we are in the 1240 >= 1 case, return the existing offCurve point. 1241 1242 The test below checks to see if we are close to the "end" of the 1243 curve (near base[4]). Rather than specifying a tolerance, I just 1244 check to see if value is on to the right/left of the middle point 1245 (depending on the direction/sign of the end points). 1246 */ 1247 if ((base[0] < base[4] && value > base[2]) || 1248 (base[0] > base[4] && value < base[2])) // should root have been 1 1249 { 1250 *offCurve = quad[1]; 1251 return true; 1252 } 1253 } 1254 return false; 1255} 1256 1257#ifdef SK_SCALAR_IS_FLOAT 1258 1259// Due to floating point issues (i.e., 1.0f - SK_ScalarRoot2Over2 != 1260// SK_ScalarRoot2Over2 - SK_ScalarTanPIOver8) a cruder root2over2 1261// approximation is required to make the quad circle points convex. The 1262// root of the problem is that with the root2over2 value in SkScalar.h 1263// the arcs really are ever so slightly concave. Some alternative fixes 1264// to this problem (besides just arbitrarily pushing out the mid-point as 1265// is done here) are: 1266// Adjust all the points (not just the middle) to both better approximate 1267// the curve and remain convex 1268// Switch over to using cubics rather then quads 1269// Use a different method to create the mid-point (e.g., compute 1270// the two side points, average them, then move it out as needed 1271#define SK_ScalarRoot2Over2_QuadCircle 0.7072f 1272 1273#else 1274 #define SK_ScalarRoot2Over2_QuadCircle SK_FixedRoot2Over2 1275#endif 1276 1277 1278static const SkPoint gQuadCirclePts[kSkBuildQuadArcStorage] = { 1279 { SK_Scalar1, 0 }, 1280 { SK_Scalar1, SK_ScalarTanPIOver8 }, 1281 { SK_ScalarRoot2Over2_QuadCircle, SK_ScalarRoot2Over2_QuadCircle }, 1282 { SK_ScalarTanPIOver8, SK_Scalar1 }, 1283 1284 { 0, SK_Scalar1 }, 1285 { -SK_ScalarTanPIOver8, SK_Scalar1 }, 1286 { -SK_ScalarRoot2Over2_QuadCircle, SK_ScalarRoot2Over2_QuadCircle }, 1287 { -SK_Scalar1, SK_ScalarTanPIOver8 }, 1288 1289 { -SK_Scalar1, 0 }, 1290 { -SK_Scalar1, -SK_ScalarTanPIOver8 }, 1291 { -SK_ScalarRoot2Over2_QuadCircle, -SK_ScalarRoot2Over2_QuadCircle }, 1292 { -SK_ScalarTanPIOver8, -SK_Scalar1 }, 1293 1294 { 0, -SK_Scalar1 }, 1295 { SK_ScalarTanPIOver8, -SK_Scalar1 }, 1296 { SK_ScalarRoot2Over2_QuadCircle, -SK_ScalarRoot2Over2_QuadCircle }, 1297 { SK_Scalar1, -SK_ScalarTanPIOver8 }, 1298 1299 { SK_Scalar1, 0 } 1300}; 1301 1302int SkBuildQuadArc(const SkVector& uStart, const SkVector& uStop, 1303 SkRotationDirection dir, const SkMatrix* userMatrix, 1304 SkPoint quadPoints[]) 1305{ 1306 // rotate by x,y so that uStart is (1.0) 1307 SkScalar x = SkPoint::DotProduct(uStart, uStop); 1308 SkScalar y = SkPoint::CrossProduct(uStart, uStop); 1309 1310 SkScalar absX = SkScalarAbs(x); 1311 SkScalar absY = SkScalarAbs(y); 1312 1313 int pointCount; 1314 1315 // check for (effectively) coincident vectors 1316 // this can happen if our angle is nearly 0 or nearly 180 (y == 0) 1317 // ... we use the dot-prod to distinguish between 0 and 180 (x > 0) 1318 if (absY <= SK_ScalarNearlyZero && x > 0 && 1319 ((y >= 0 && kCW_SkRotationDirection == dir) || 1320 (y <= 0 && kCCW_SkRotationDirection == dir))) { 1321 1322 // just return the start-point 1323 quadPoints[0].set(SK_Scalar1, 0); 1324 pointCount = 1; 1325 } else { 1326 if (dir == kCCW_SkRotationDirection) 1327 y = -y; 1328 1329 // what octant (quadratic curve) is [xy] in? 1330 int oct = 0; 1331 bool sameSign = true; 1332 1333 if (0 == y) 1334 { 1335 oct = 4; // 180 1336 SkASSERT(SkScalarAbs(x + SK_Scalar1) <= SK_ScalarNearlyZero); 1337 } 1338 else if (0 == x) 1339 { 1340 SkASSERT(absY - SK_Scalar1 <= SK_ScalarNearlyZero); 1341 if (y > 0) 1342 oct = 2; // 90 1343 else 1344 oct = 6; // 270 1345 } 1346 else 1347 { 1348 if (y < 0) 1349 oct += 4; 1350 if ((x < 0) != (y < 0)) 1351 { 1352 oct += 2; 1353 sameSign = false; 1354 } 1355 if ((absX < absY) == sameSign) 1356 oct += 1; 1357 } 1358 1359 int wholeCount = oct << 1; 1360 memcpy(quadPoints, gQuadCirclePts, (wholeCount + 1) * sizeof(SkPoint)); 1361 1362 const SkPoint* arc = &gQuadCirclePts[wholeCount]; 1363 if (quad_pt2OffCurve(arc, x, y, &quadPoints[wholeCount + 1])) 1364 { 1365 quadPoints[wholeCount + 2].set(x, y); 1366 wholeCount += 2; 1367 } 1368 pointCount = wholeCount + 1; 1369 } 1370 1371 // now handle counter-clockwise and the initial unitStart rotation 1372 SkMatrix matrix; 1373 matrix.setSinCos(uStart.fY, uStart.fX); 1374 if (dir == kCCW_SkRotationDirection) { 1375 matrix.preScale(SK_Scalar1, -SK_Scalar1); 1376 } 1377 if (userMatrix) { 1378 matrix.postConcat(*userMatrix); 1379 } 1380 matrix.mapPoints(quadPoints, pointCount); 1381 return pointCount; 1382} 1383 1384/////////////////////////////////////////////////////////////////////////////// 1385 1386static SkScalar eval_ratquad(const SkScalar src[], SkScalar w, SkScalar t) { 1387 SkASSERT(src); 1388 SkASSERT(t >= 0 && t <= SK_Scalar1); 1389 1390 SkScalar src2w = SkScalarMul(src[2], w); 1391 SkScalar C = src[0]; 1392 SkScalar A = src[4] - 2 * src2w + C; 1393 SkScalar B = 2 * (src2w - C); 1394 SkScalar numer = SkScalarMulAdd(SkScalarMulAdd(A, t, B), t, C); 1395 1396 B = 2 * (w - SK_Scalar1); 1397 C = SK_Scalar1; 1398 A = -B; 1399 SkScalar denom = SkScalarMulAdd(SkScalarMulAdd(A, t, B), t, C); 1400 1401 return SkScalarDiv(numer, denom); 1402} 1403 1404struct SkP3D { 1405 SkScalar fX, fY, fZ; 1406 1407 void set(SkScalar x, SkScalar y, SkScalar z) { 1408 fX = x; fY = y; fZ = z; 1409 } 1410 1411 void projectDown(SkPoint* dst) const { 1412 dst->set(fX / fZ, fY / fZ); 1413 } 1414}; 1415 1416// we just return the middle 3 points, since the first and last are dups of src 1417// 1418static void p3d_interp(const SkScalar src[3], SkScalar dst[3], SkScalar t) { 1419 SkScalar ab = SkScalarInterp(src[0], src[3], t); 1420 SkScalar bc = SkScalarInterp(src[3], src[6], t); 1421 dst[0] = ab; 1422 dst[3] = SkScalarInterp(ab, bc, t); 1423 dst[6] = bc; 1424} 1425 1426static void ratquad_mapTo3D(const SkPoint src[3], SkScalar w, SkP3D dst[]) { 1427 dst[0].set(src[0].fX * 1, src[0].fY * 1, 1); 1428 dst[1].set(src[1].fX * w, src[1].fY * w, w); 1429 dst[2].set(src[2].fX * 1, src[2].fY * 1, 1); 1430} 1431 1432void SkRationalQuad::evalAt(SkScalar t, SkPoint* pt) const { 1433 SkASSERT(t >= 0 && t <= SK_Scalar1); 1434 1435 if (pt) { 1436 pt->set(eval_ratquad(&fPts[0].fX, fW, t), 1437 eval_ratquad(&fPts[0].fY, fW, t)); 1438 } 1439} 1440 1441void SkRationalQuad::chopAt(SkScalar t, SkRationalQuad dst[2]) const { 1442 SkP3D tmp[3], tmp2[3]; 1443 1444 ratquad_mapTo3D(fPts, fW, tmp); 1445 1446 p3d_interp(&tmp[0].fX, &tmp2[0].fX, t); 1447 p3d_interp(&tmp[0].fY, &tmp2[0].fY, t); 1448 p3d_interp(&tmp[0].fZ, &tmp2[0].fZ, t); 1449 1450 dst[0].fPts[0] = fPts[0]; 1451 tmp2[0].projectDown(&dst[0].fPts[1]); 1452 tmp2[1].projectDown(&dst[0].fPts[2]); dst[1].fPts[0] = dst[0].fPts[2]; 1453 tmp2[2].projectDown(&dst[1].fPts[1]); 1454 dst[1].fPts[2] = fPts[2]; 1455 1456 dst[0].fW = tmp2[0].fZ; // ????? 1457 dst[1].fW = tmp2[2].fZ; // ????? 1458} 1459