SkGeometry.cpp revision 945a139553a9c9da03766213661d7f5fd6ed3042
1/* libs/graphics/sgl/SkGeometry.cpp 2** 3** Copyright 2006, The Android Open Source Project 4** 5** Licensed under the Apache License, Version 2.0 (the "License"); 6** you may not use this file except in compliance with the License. 7** You may obtain a copy of the License at 8** 9** http://www.apache.org/licenses/LICENSE-2.0 10** 11** Unless required by applicable law or agreed to in writing, software 12** distributed under the License is distributed on an "AS IS" BASIS, 13** WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14** See the License for the specific language governing permissions and 15** limitations under the License. 16*/ 17 18#include "SkGeometry.h" 19#include "Sk64.h" 20#include "SkMatrix.h" 21 22bool SkXRayCrossesLine(const SkXRay& pt, const SkPoint pts[2]) { 23 // Determine quick discards. 24 // Consider query line going exactly through point 0 to not 25 // intersect, for symmetry with SkXRayCrossesMonotonicCubic. 26 if (pt.fY == pts[0].fY) 27 return false; 28 if (pt.fY < pts[0].fY && pt.fY < pts[1].fY) 29 return false; 30 if (pt.fY > pts[0].fY && pt.fY > pts[1].fY) 31 return false; 32 if (pt.fX > pts[0].fX && pt.fX > pts[1].fX) 33 return false; 34 // Determine degenerate cases 35 if (SkScalarNearlyZero(pts[0].fY - pts[1].fY)) 36 return false; 37 if (SkScalarNearlyZero(pts[0].fX - pts[1].fX)) 38 // We've already determined the query point lies within the 39 // vertical range of the line segment. 40 return pt.fX <= pts[0].fX; 41 // Full line segment evaluation 42 SkScalar delta_y = pts[1].fY - pts[0].fY; 43 SkScalar delta_x = pts[1].fX - pts[0].fX; 44 SkScalar slope = SkScalarDiv(delta_y, delta_x); 45 SkScalar b = pts[0].fY - SkScalarMul(slope, pts[0].fX); 46 // Solve for x coordinate at y = pt.fY 47 SkScalar x = SkScalarDiv(pt.fY - b, slope); 48 return pt.fX <= x; 49} 50 51/** If defined, this makes eval_quad and eval_cubic do more setup (sometimes 52 involving integer multiplies by 2 or 3, but fewer calls to SkScalarMul. 53 May also introduce overflow of fixed when we compute our setup. 54*/ 55#ifdef SK_SCALAR_IS_FIXED 56 #define DIRECT_EVAL_OF_POLYNOMIALS 57#endif 58 59//////////////////////////////////////////////////////////////////////// 60 61#ifdef SK_SCALAR_IS_FIXED 62 static int is_not_monotonic(int a, int b, int c, int d) 63 { 64 return (((a - b) | (b - c) | (c - d)) & ((b - a) | (c - b) | (d - c))) >> 31; 65 } 66 67 static int is_not_monotonic(int a, int b, int c) 68 { 69 return (((a - b) | (b - c)) & ((b - a) | (c - b))) >> 31; 70 } 71#else 72 static int is_not_monotonic(float a, float b, float c) 73 { 74 float ab = a - b; 75 float bc = b - c; 76 if (ab < 0) 77 bc = -bc; 78 return ab == 0 || bc < 0; 79 } 80#endif 81 82//////////////////////////////////////////////////////////////////////// 83 84static bool is_unit_interval(SkScalar x) 85{ 86 return x > 0 && x < SK_Scalar1; 87} 88 89static int valid_unit_divide(SkScalar numer, SkScalar denom, SkScalar* ratio) 90{ 91 SkASSERT(ratio); 92 93 if (numer < 0) 94 { 95 numer = -numer; 96 denom = -denom; 97 } 98 99 if (denom == 0 || numer == 0 || numer >= denom) 100 return 0; 101 102 SkScalar r = SkScalarDiv(numer, denom); 103 SkASSERT(r >= 0 && r < SK_Scalar1); 104 if (r == 0) // catch underflow if numer <<<< denom 105 return 0; 106 *ratio = r; 107 return 1; 108} 109 110/** From Numerical Recipes in C. 111 112 Q = -1/2 (B + sign(B) sqrt[B*B - 4*A*C]) 113 x1 = Q / A 114 x2 = C / Q 115*/ 116int SkFindUnitQuadRoots(SkScalar A, SkScalar B, SkScalar C, SkScalar roots[2]) 117{ 118 SkASSERT(roots); 119 120 if (A == 0) 121 return valid_unit_divide(-C, B, roots); 122 123 SkScalar* r = roots; 124 125#ifdef SK_SCALAR_IS_FLOAT 126 float R = B*B - 4*A*C; 127 if (R < 0) // complex roots 128 return 0; 129 R = sk_float_sqrt(R); 130#else 131 Sk64 RR, tmp; 132 133 RR.setMul(B,B); 134 tmp.setMul(A,C); 135 tmp.shiftLeft(2); 136 RR.sub(tmp); 137 if (RR.isNeg()) 138 return 0; 139 SkFixed R = RR.getSqrt(); 140#endif 141 142 SkScalar Q = (B < 0) ? -(B-R)/2 : -(B+R)/2; 143 r += valid_unit_divide(Q, A, r); 144 r += valid_unit_divide(C, Q, r); 145 if (r - roots == 2) 146 { 147 if (roots[0] > roots[1]) 148 SkTSwap<SkScalar>(roots[0], roots[1]); 149 else if (roots[0] == roots[1]) // nearly-equal? 150 r -= 1; // skip the double root 151 } 152 return (int)(r - roots); 153} 154 155#ifdef SK_SCALAR_IS_FIXED 156/** Trim A/B/C down so that they are all <= 32bits 157 and then call SkFindUnitQuadRoots() 158*/ 159static int Sk64FindFixedQuadRoots(const Sk64& A, const Sk64& B, const Sk64& C, SkFixed roots[2]) 160{ 161 int na = A.shiftToMake32(); 162 int nb = B.shiftToMake32(); 163 int nc = C.shiftToMake32(); 164 165 int shift = SkMax32(na, SkMax32(nb, nc)); 166 SkASSERT(shift >= 0); 167 168 return SkFindUnitQuadRoots(A.getShiftRight(shift), B.getShiftRight(shift), C.getShiftRight(shift), roots); 169} 170#endif 171 172///////////////////////////////////////////////////////////////////////////////////// 173///////////////////////////////////////////////////////////////////////////////////// 174 175static SkScalar eval_quad(const SkScalar src[], SkScalar t) 176{ 177 SkASSERT(src); 178 SkASSERT(t >= 0 && t <= SK_Scalar1); 179 180#ifdef DIRECT_EVAL_OF_POLYNOMIALS 181 SkScalar C = src[0]; 182 SkScalar A = src[4] - 2 * src[2] + C; 183 SkScalar B = 2 * (src[2] - C); 184 return SkScalarMulAdd(SkScalarMulAdd(A, t, B), t, C); 185#else 186 SkScalar ab = SkScalarInterp(src[0], src[2], t); 187 SkScalar bc = SkScalarInterp(src[2], src[4], t); 188 return SkScalarInterp(ab, bc, t); 189#endif 190} 191 192static SkScalar eval_quad_derivative(const SkScalar src[], SkScalar t) 193{ 194 SkScalar A = src[4] - 2 * src[2] + src[0]; 195 SkScalar B = src[2] - src[0]; 196 197 return 2 * SkScalarMulAdd(A, t, B); 198} 199 200static SkScalar eval_quad_derivative_at_half(const SkScalar src[]) 201{ 202 SkScalar A = src[4] - 2 * src[2] + src[0]; 203 SkScalar B = src[2] - src[0]; 204 return A + 2 * B; 205} 206 207void SkEvalQuadAt(const SkPoint src[3], SkScalar t, SkPoint* pt, SkVector* tangent) 208{ 209 SkASSERT(src); 210 SkASSERT(t >= 0 && t <= SK_Scalar1); 211 212 if (pt) 213 pt->set(eval_quad(&src[0].fX, t), eval_quad(&src[0].fY, t)); 214 if (tangent) 215 tangent->set(eval_quad_derivative(&src[0].fX, t), 216 eval_quad_derivative(&src[0].fY, t)); 217} 218 219void SkEvalQuadAtHalf(const SkPoint src[3], SkPoint* pt, SkVector* tangent) 220{ 221 SkASSERT(src); 222 223 if (pt) 224 { 225 SkScalar x01 = SkScalarAve(src[0].fX, src[1].fX); 226 SkScalar y01 = SkScalarAve(src[0].fY, src[1].fY); 227 SkScalar x12 = SkScalarAve(src[1].fX, src[2].fX); 228 SkScalar y12 = SkScalarAve(src[1].fY, src[2].fY); 229 pt->set(SkScalarAve(x01, x12), SkScalarAve(y01, y12)); 230 } 231 if (tangent) 232 tangent->set(eval_quad_derivative_at_half(&src[0].fX), 233 eval_quad_derivative_at_half(&src[0].fY)); 234} 235 236static void interp_quad_coords(const SkScalar* src, SkScalar* dst, SkScalar t) 237{ 238 SkScalar ab = SkScalarInterp(src[0], src[2], t); 239 SkScalar bc = SkScalarInterp(src[2], src[4], t); 240 241 dst[0] = src[0]; 242 dst[2] = ab; 243 dst[4] = SkScalarInterp(ab, bc, t); 244 dst[6] = bc; 245 dst[8] = src[4]; 246} 247 248void SkChopQuadAt(const SkPoint src[3], SkPoint dst[5], SkScalar t) 249{ 250 SkASSERT(t > 0 && t < SK_Scalar1); 251 252 interp_quad_coords(&src[0].fX, &dst[0].fX, t); 253 interp_quad_coords(&src[0].fY, &dst[0].fY, t); 254} 255 256void SkChopQuadAtHalf(const SkPoint src[3], SkPoint dst[5]) 257{ 258 SkScalar x01 = SkScalarAve(src[0].fX, src[1].fX); 259 SkScalar y01 = SkScalarAve(src[0].fY, src[1].fY); 260 SkScalar x12 = SkScalarAve(src[1].fX, src[2].fX); 261 SkScalar y12 = SkScalarAve(src[1].fY, src[2].fY); 262 263 dst[0] = src[0]; 264 dst[1].set(x01, y01); 265 dst[2].set(SkScalarAve(x01, x12), SkScalarAve(y01, y12)); 266 dst[3].set(x12, y12); 267 dst[4] = src[2]; 268} 269 270/** Quad'(t) = At + B, where 271 A = 2(a - 2b + c) 272 B = 2(b - a) 273 Solve for t, only if it fits between 0 < t < 1 274*/ 275int SkFindQuadExtrema(SkScalar a, SkScalar b, SkScalar c, SkScalar tValue[1]) 276{ 277 /* At + B == 0 278 t = -B / A 279 */ 280#ifdef SK_SCALAR_IS_FIXED 281 return is_not_monotonic(a, b, c) && valid_unit_divide(a - b, a - b - b + c, tValue); 282#else 283 return valid_unit_divide(a - b, a - b - b + c, tValue); 284#endif 285} 286 287static inline void flatten_double_quad_extrema(SkScalar coords[14]) 288{ 289 coords[2] = coords[6] = coords[4]; 290} 291 292/* Returns 0 for 1 quad, and 1 for two quads, either way the answer is 293 stored in dst[]. Guarantees that the 1/2 quads will be monotonic. 294 */ 295int SkChopQuadAtYExtrema(const SkPoint src[3], SkPoint dst[5]) 296{ 297 SkASSERT(src); 298 SkASSERT(dst); 299 300#if 0 301 static bool once = true; 302 if (once) 303 { 304 once = false; 305 SkPoint s[3] = { 0, 26398, 0, 26331, 0, 20621428 }; 306 SkPoint d[6]; 307 308 int n = SkChopQuadAtYExtrema(s, d); 309 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); 310 } 311#endif 312 313 SkScalar a = src[0].fY; 314 SkScalar b = src[1].fY; 315 SkScalar c = src[2].fY; 316 317 if (is_not_monotonic(a, b, c)) 318 { 319 SkScalar tValue; 320 if (valid_unit_divide(a - b, a - b - b + c, &tValue)) 321 { 322 SkChopQuadAt(src, dst, tValue); 323 flatten_double_quad_extrema(&dst[0].fY); 324 return 1; 325 } 326 // if we get here, we need to force dst to be monotonic, even though 327 // we couldn't compute a unit_divide value (probably underflow). 328 b = SkScalarAbs(a - b) < SkScalarAbs(b - c) ? a : c; 329 } 330 dst[0].set(src[0].fX, a); 331 dst[1].set(src[1].fX, b); 332 dst[2].set(src[2].fX, c); 333 return 0; 334} 335 336/* Returns 0 for 1 quad, and 1 for two quads, either way the answer is 337 stored in dst[]. Guarantees that the 1/2 quads will be monotonic. 338 */ 339int SkChopQuadAtXExtrema(const SkPoint src[3], SkPoint dst[5]) 340{ 341 SkASSERT(src); 342 SkASSERT(dst); 343 344 SkScalar a = src[0].fX; 345 SkScalar b = src[1].fX; 346 SkScalar c = src[2].fX; 347 348 if (is_not_monotonic(a, b, c)) { 349 SkScalar tValue; 350 if (valid_unit_divide(a - b, a - b - b + c, &tValue)) { 351 SkChopQuadAt(src, dst, tValue); 352 flatten_double_quad_extrema(&dst[0].fX); 353 return 1; 354 } 355 // if we get here, we need to force dst to be monotonic, even though 356 // we couldn't compute a unit_divide value (probably underflow). 357 b = SkScalarAbs(a - b) < SkScalarAbs(b - c) ? a : c; 358 } 359 dst[0].set(a, src[0].fY); 360 dst[1].set(b, src[1].fY); 361 dst[2].set(c, src[2].fY); 362 return 0; 363} 364 365// F(t) = a (1 - t) ^ 2 + 2 b t (1 - t) + c t ^ 2 366// F'(t) = 2 (b - a) + 2 (a - 2b + c) t 367// F''(t) = 2 (a - 2b + c) 368// 369// A = 2 (b - a) 370// B = 2 (a - 2b + c) 371// 372// Maximum curvature for a quadratic means solving 373// Fx' Fx'' + Fy' Fy'' = 0 374// 375// t = - (Ax Bx + Ay By) / (Bx ^ 2 + By ^ 2) 376// 377int SkChopQuadAtMaxCurvature(const SkPoint src[3], SkPoint dst[5]) 378{ 379 SkScalar Ax = src[1].fX - src[0].fX; 380 SkScalar Ay = src[1].fY - src[0].fY; 381 SkScalar Bx = src[0].fX - src[1].fX - src[1].fX + src[2].fX; 382 SkScalar By = src[0].fY - src[1].fY - src[1].fY + src[2].fY; 383 SkScalar t = 0; // 0 means don't chop 384 385#ifdef SK_SCALAR_IS_FLOAT 386 (void)valid_unit_divide(-(Ax * Bx + Ay * By), Bx * Bx + By * By, &t); 387#else 388 // !!! should I use SkFloat here? seems like it 389 Sk64 numer, denom, tmp; 390 391 numer.setMul(Ax, -Bx); 392 tmp.setMul(Ay, -By); 393 numer.add(tmp); 394 395 if (numer.isPos()) // do nothing if numer <= 0 396 { 397 denom.setMul(Bx, Bx); 398 tmp.setMul(By, By); 399 denom.add(tmp); 400 SkASSERT(!denom.isNeg()); 401 if (numer < denom) 402 { 403 t = numer.getFixedDiv(denom); 404 SkASSERT(t >= 0 && t <= SK_Fixed1); // assert that we're numerically stable (ha!) 405 if ((unsigned)t >= SK_Fixed1) // runtime check for numerical stability 406 t = 0; // ignore the chop 407 } 408 } 409#endif 410 411 if (t == 0) 412 { 413 memcpy(dst, src, 3 * sizeof(SkPoint)); 414 return 1; 415 } 416 else 417 { 418 SkChopQuadAt(src, dst, t); 419 return 2; 420 } 421} 422 423void SkConvertQuadToCubic(const SkPoint src[3], SkPoint dst[4]) 424{ 425 SkScalar two = SkIntToScalar(2); 426 SkScalar one_third = SkScalarDiv(SK_Scalar1, SkIntToScalar(3)); 427 dst[0].set(src[0].fX, src[0].fY); 428 dst[1].set( 429 SkScalarMul(SkScalarMulAdd(src[1].fX, two, src[0].fX), one_third), 430 SkScalarMul(SkScalarMulAdd(src[1].fY, two, src[0].fY), one_third)); 431 dst[2].set( 432 SkScalarMul(SkScalarMulAdd(src[1].fX, two, src[2].fX), one_third), 433 SkScalarMul(SkScalarMulAdd(src[1].fY, two, src[2].fY), one_third)); 434 dst[3].set(src[2].fX, src[2].fY); 435} 436 437//////////////////////////////////////////////////////////////////////////////////////// 438///// CUBICS // CUBICS // CUBICS // CUBICS // CUBICS // CUBICS // CUBICS // CUBICS ///// 439//////////////////////////////////////////////////////////////////////////////////////// 440 441static void get_cubic_coeff(const SkScalar pt[], SkScalar coeff[4]) 442{ 443 coeff[0] = pt[6] + 3*(pt[2] - pt[4]) - pt[0]; 444 coeff[1] = 3*(pt[4] - pt[2] - pt[2] + pt[0]); 445 coeff[2] = 3*(pt[2] - pt[0]); 446 coeff[3] = pt[0]; 447} 448 449void SkGetCubicCoeff(const SkPoint pts[4], SkScalar cx[4], SkScalar cy[4]) 450{ 451 SkASSERT(pts); 452 453 if (cx) 454 get_cubic_coeff(&pts[0].fX, cx); 455 if (cy) 456 get_cubic_coeff(&pts[0].fY, cy); 457} 458 459static SkScalar eval_cubic(const SkScalar src[], SkScalar t) 460{ 461 SkASSERT(src); 462 SkASSERT(t >= 0 && t <= SK_Scalar1); 463 464 if (t == 0) 465 return src[0]; 466 467#ifdef DIRECT_EVAL_OF_POLYNOMIALS 468 SkScalar D = src[0]; 469 SkScalar A = src[6] + 3*(src[2] - src[4]) - D; 470 SkScalar B = 3*(src[4] - src[2] - src[2] + D); 471 SkScalar C = 3*(src[2] - D); 472 473 return SkScalarMulAdd(SkScalarMulAdd(SkScalarMulAdd(A, t, B), t, C), t, D); 474#else 475 SkScalar ab = SkScalarInterp(src[0], src[2], t); 476 SkScalar bc = SkScalarInterp(src[2], src[4], t); 477 SkScalar cd = SkScalarInterp(src[4], src[6], t); 478 SkScalar abc = SkScalarInterp(ab, bc, t); 479 SkScalar bcd = SkScalarInterp(bc, cd, t); 480 return SkScalarInterp(abc, bcd, t); 481#endif 482} 483 484/** return At^2 + Bt + C 485*/ 486static SkScalar eval_quadratic(SkScalar A, SkScalar B, SkScalar C, SkScalar t) 487{ 488 SkASSERT(t >= 0 && t <= SK_Scalar1); 489 490 return SkScalarMulAdd(SkScalarMulAdd(A, t, B), t, C); 491} 492 493static SkScalar eval_cubic_derivative(const SkScalar src[], SkScalar t) 494{ 495 SkScalar A = src[6] + 3*(src[2] - src[4]) - src[0]; 496 SkScalar B = 2*(src[4] - 2 * src[2] + src[0]); 497 SkScalar C = src[2] - src[0]; 498 499 return eval_quadratic(A, B, C, t); 500} 501 502static SkScalar eval_cubic_2ndDerivative(const SkScalar src[], SkScalar t) 503{ 504 SkScalar A = src[6] + 3*(src[2] - src[4]) - src[0]; 505 SkScalar B = src[4] - 2 * src[2] + src[0]; 506 507 return SkScalarMulAdd(A, t, B); 508} 509 510void SkEvalCubicAt(const SkPoint src[4], SkScalar t, SkPoint* loc, SkVector* tangent, SkVector* curvature) 511{ 512 SkASSERT(src); 513 SkASSERT(t >= 0 && t <= SK_Scalar1); 514 515 if (loc) 516 loc->set(eval_cubic(&src[0].fX, t), eval_cubic(&src[0].fY, t)); 517 if (tangent) 518 tangent->set(eval_cubic_derivative(&src[0].fX, t), 519 eval_cubic_derivative(&src[0].fY, t)); 520 if (curvature) 521 curvature->set(eval_cubic_2ndDerivative(&src[0].fX, t), 522 eval_cubic_2ndDerivative(&src[0].fY, t)); 523} 524 525/** Cubic'(t) = At^2 + Bt + C, where 526 A = 3(-a + 3(b - c) + d) 527 B = 6(a - 2b + c) 528 C = 3(b - a) 529 Solve for t, keeping only those that fit betwee 0 < t < 1 530*/ 531int SkFindCubicExtrema(SkScalar a, SkScalar b, SkScalar c, SkScalar d, SkScalar tValues[2]) 532{ 533#ifdef SK_SCALAR_IS_FIXED 534 if (!is_not_monotonic(a, b, c, d)) 535 return 0; 536#endif 537 538 // we divide A,B,C by 3 to simplify 539 SkScalar A = d - a + 3*(b - c); 540 SkScalar B = 2*(a - b - b + c); 541 SkScalar C = b - a; 542 543 return SkFindUnitQuadRoots(A, B, C, tValues); 544} 545 546static void interp_cubic_coords(const SkScalar* src, SkScalar* dst, SkScalar t) 547{ 548 SkScalar ab = SkScalarInterp(src[0], src[2], t); 549 SkScalar bc = SkScalarInterp(src[2], src[4], t); 550 SkScalar cd = SkScalarInterp(src[4], src[6], t); 551 SkScalar abc = SkScalarInterp(ab, bc, t); 552 SkScalar bcd = SkScalarInterp(bc, cd, t); 553 SkScalar abcd = SkScalarInterp(abc, bcd, t); 554 555 dst[0] = src[0]; 556 dst[2] = ab; 557 dst[4] = abc; 558 dst[6] = abcd; 559 dst[8] = bcd; 560 dst[10] = cd; 561 dst[12] = src[6]; 562} 563 564void SkChopCubicAt(const SkPoint src[4], SkPoint dst[7], SkScalar t) 565{ 566 SkASSERT(t > 0 && t < SK_Scalar1); 567 568 interp_cubic_coords(&src[0].fX, &dst[0].fX, t); 569 interp_cubic_coords(&src[0].fY, &dst[0].fY, t); 570} 571 572/* http://code.google.com/p/skia/issues/detail?id=32 573 574 This test code would fail when we didn't check the return result of 575 valid_unit_divide in SkChopCubicAt(... tValues[], int roots). The reason is 576 that after the first chop, the parameters to valid_unit_divide are equal 577 (thanks to finite float precision and rounding in the subtracts). Thus 578 even though the 2nd tValue looks < 1.0, after we renormalize it, we end 579 up with 1.0, hence the need to check and just return the last cubic as 580 a degenerate clump of 4 points in the sampe place. 581 582 static void test_cubic() { 583 SkPoint src[4] = { 584 { 556.25000, 523.03003 }, 585 { 556.23999, 522.96002 }, 586 { 556.21997, 522.89001 }, 587 { 556.21997, 522.82001 } 588 }; 589 SkPoint dst[10]; 590 SkScalar tval[] = { 0.33333334f, 0.99999994f }; 591 SkChopCubicAt(src, dst, tval, 2); 592 } 593 */ 594 595void SkChopCubicAt(const SkPoint src[4], SkPoint dst[], const SkScalar tValues[], int roots) 596{ 597#ifdef SK_DEBUG 598 { 599 for (int i = 0; i < roots - 1; i++) 600 { 601 SkASSERT(is_unit_interval(tValues[i])); 602 SkASSERT(is_unit_interval(tValues[i+1])); 603 SkASSERT(tValues[i] < tValues[i+1]); 604 } 605 } 606#endif 607 608 if (dst) 609 { 610 if (roots == 0) // nothing to chop 611 memcpy(dst, src, 4*sizeof(SkPoint)); 612 else 613 { 614 SkScalar t = tValues[0]; 615 SkPoint tmp[4]; 616 617 for (int i = 0; i < roots; i++) 618 { 619 SkChopCubicAt(src, dst, t); 620 if (i == roots - 1) 621 break; 622 623 dst += 3; 624 // have src point to the remaining cubic (after the chop) 625 memcpy(tmp, dst, 4 * sizeof(SkPoint)); 626 src = tmp; 627 628 // watch out in case the renormalized t isn't in range 629 if (!valid_unit_divide(tValues[i+1] - tValues[i], 630 SK_Scalar1 - tValues[i], &t)) { 631 // if we can't, just create a degenerate cubic 632 dst[4] = dst[5] = dst[6] = src[3]; 633 break; 634 } 635 } 636 } 637 } 638} 639 640void SkChopCubicAtHalf(const SkPoint src[4], SkPoint dst[7]) 641{ 642 SkScalar x01 = SkScalarAve(src[0].fX, src[1].fX); 643 SkScalar y01 = SkScalarAve(src[0].fY, src[1].fY); 644 SkScalar x12 = SkScalarAve(src[1].fX, src[2].fX); 645 SkScalar y12 = SkScalarAve(src[1].fY, src[2].fY); 646 SkScalar x23 = SkScalarAve(src[2].fX, src[3].fX); 647 SkScalar y23 = SkScalarAve(src[2].fY, src[3].fY); 648 649 SkScalar x012 = SkScalarAve(x01, x12); 650 SkScalar y012 = SkScalarAve(y01, y12); 651 SkScalar x123 = SkScalarAve(x12, x23); 652 SkScalar y123 = SkScalarAve(y12, y23); 653 654 dst[0] = src[0]; 655 dst[1].set(x01, y01); 656 dst[2].set(x012, y012); 657 dst[3].set(SkScalarAve(x012, x123), SkScalarAve(y012, y123)); 658 dst[4].set(x123, y123); 659 dst[5].set(x23, y23); 660 dst[6] = src[3]; 661} 662 663static void flatten_double_cubic_extrema(SkScalar coords[14]) 664{ 665 coords[4] = coords[8] = coords[6]; 666} 667 668/** Given 4 points on a cubic bezier, chop it into 1, 2, 3 beziers such that 669 the resulting beziers are monotonic in Y. This is called by the scan converter. 670 Depending on what is returned, dst[] is treated as follows 671 0 dst[0..3] is the original cubic 672 1 dst[0..3] and dst[3..6] are the two new cubics 673 2 dst[0..3], dst[3..6], dst[6..9] are the three new cubics 674 If dst == null, it is ignored and only the count is returned. 675*/ 676int SkChopCubicAtYExtrema(const SkPoint src[4], SkPoint dst[10]) { 677 SkScalar tValues[2]; 678 int roots = SkFindCubicExtrema(src[0].fY, src[1].fY, src[2].fY, 679 src[3].fY, tValues); 680 681 SkChopCubicAt(src, dst, tValues, roots); 682 if (dst && roots > 0) { 683 // we do some cleanup to ensure our Y extrema are flat 684 flatten_double_cubic_extrema(&dst[0].fY); 685 if (roots == 2) { 686 flatten_double_cubic_extrema(&dst[3].fY); 687 } 688 } 689 return roots; 690} 691 692int SkChopCubicAtXExtrema(const SkPoint src[4], SkPoint dst[10]) { 693 SkScalar tValues[2]; 694 int roots = SkFindCubicExtrema(src[0].fX, src[1].fX, src[2].fX, 695 src[3].fX, tValues); 696 697 SkChopCubicAt(src, dst, tValues, roots); 698 if (dst && roots > 0) { 699 // we do some cleanup to ensure our Y extrema are flat 700 flatten_double_cubic_extrema(&dst[0].fX); 701 if (roots == 2) { 702 flatten_double_cubic_extrema(&dst[3].fX); 703 } 704 } 705 return roots; 706} 707 708/** http://www.faculty.idc.ac.il/arik/quality/appendixA.html 709 710 Inflection means that curvature is zero. 711 Curvature is [F' x F''] / [F'^3] 712 So we solve F'x X F''y - F'y X F''y == 0 713 After some canceling of the cubic term, we get 714 A = b - a 715 B = c - 2b + a 716 C = d - 3c + 3b - a 717 (BxCy - ByCx)t^2 + (AxCy - AyCx)t + AxBy - AyBx == 0 718*/ 719int SkFindCubicInflections(const SkPoint src[4], SkScalar tValues[]) 720{ 721 SkScalar Ax = src[1].fX - src[0].fX; 722 SkScalar Ay = src[1].fY - src[0].fY; 723 SkScalar Bx = src[2].fX - 2 * src[1].fX + src[0].fX; 724 SkScalar By = src[2].fY - 2 * src[1].fY + src[0].fY; 725 SkScalar Cx = src[3].fX + 3 * (src[1].fX - src[2].fX) - src[0].fX; 726 SkScalar Cy = src[3].fY + 3 * (src[1].fY - src[2].fY) - src[0].fY; 727 int count; 728 729#ifdef SK_SCALAR_IS_FLOAT 730 count = SkFindUnitQuadRoots(Bx*Cy - By*Cx, Ax*Cy - Ay*Cx, Ax*By - Ay*Bx, tValues); 731#else 732 Sk64 A, B, C, tmp; 733 734 A.setMul(Bx, Cy); 735 tmp.setMul(By, Cx); 736 A.sub(tmp); 737 738 B.setMul(Ax, Cy); 739 tmp.setMul(Ay, Cx); 740 B.sub(tmp); 741 742 C.setMul(Ax, By); 743 tmp.setMul(Ay, Bx); 744 C.sub(tmp); 745 746 count = Sk64FindFixedQuadRoots(A, B, C, tValues); 747#endif 748 749 return count; 750} 751 752int SkChopCubicAtInflections(const SkPoint src[], SkPoint dst[10]) 753{ 754 SkScalar tValues[2]; 755 int count = SkFindCubicInflections(src, tValues); 756 757 if (dst) 758 { 759 if (count == 0) 760 memcpy(dst, src, 4 * sizeof(SkPoint)); 761 else 762 SkChopCubicAt(src, dst, tValues, count); 763 } 764 return count + 1; 765} 766 767template <typename T> void bubble_sort(T array[], int count) 768{ 769 for (int i = count - 1; i > 0; --i) 770 for (int j = i; j > 0; --j) 771 if (array[j] < array[j-1]) 772 { 773 T tmp(array[j]); 774 array[j] = array[j-1]; 775 array[j-1] = tmp; 776 } 777} 778 779#include "SkFP.h" 780 781// newton refinement 782#if 0 783static SkScalar refine_cubic_root(const SkFP coeff[4], SkScalar root) 784{ 785 // x1 = x0 - f(t) / f'(t) 786 787 SkFP T = SkScalarToFloat(root); 788 SkFP N, D; 789 790 // f' = 3*coeff[0]*T^2 + 2*coeff[1]*T + coeff[2] 791 D = SkFPMul(SkFPMul(coeff[0], SkFPMul(T,T)), 3); 792 D = SkFPAdd(D, SkFPMulInt(SkFPMul(coeff[1], T), 2)); 793 D = SkFPAdd(D, coeff[2]); 794 795 if (D == 0) 796 return root; 797 798 // f = coeff[0]*T^3 + coeff[1]*T^2 + coeff[2]*T + coeff[3] 799 N = SkFPMul(SkFPMul(SkFPMul(T, T), T), coeff[0]); 800 N = SkFPAdd(N, SkFPMul(SkFPMul(T, T), coeff[1])); 801 N = SkFPAdd(N, SkFPMul(T, coeff[2])); 802 N = SkFPAdd(N, coeff[3]); 803 804 if (N) 805 { 806 SkScalar delta = SkFPToScalar(SkFPDiv(N, D)); 807 808 if (delta) 809 root -= delta; 810 } 811 return root; 812} 813#endif 814 815#if defined _WIN32 && _MSC_VER >= 1300 && defined SK_SCALAR_IS_FIXED // disable warning : unreachable code if building fixed point for windows desktop 816#pragma warning ( disable : 4702 ) 817#endif 818 819/* Solve coeff(t) == 0, returning the number of roots that 820 lie withing 0 < t < 1. 821 coeff[0]t^3 + coeff[1]t^2 + coeff[2]t + coeff[3] 822*/ 823static int solve_cubic_polynomial(const SkFP coeff[4], SkScalar tValues[3]) 824{ 825#ifndef SK_SCALAR_IS_FLOAT 826 return 0; // this is not yet implemented for software float 827#endif 828 829 if (SkScalarNearlyZero(coeff[0])) // we're just a quadratic 830 { 831 return SkFindUnitQuadRoots(coeff[1], coeff[2], coeff[3], tValues); 832 } 833 834 SkFP a, b, c, Q, R; 835 836 { 837 SkASSERT(coeff[0] != 0); 838 839 SkFP inva = SkFPInvert(coeff[0]); 840 a = SkFPMul(coeff[1], inva); 841 b = SkFPMul(coeff[2], inva); 842 c = SkFPMul(coeff[3], inva); 843 } 844 Q = SkFPDivInt(SkFPSub(SkFPMul(a,a), SkFPMulInt(b, 3)), 9); 845// R = (2*a*a*a - 9*a*b + 27*c) / 54; 846 R = SkFPMulInt(SkFPMul(SkFPMul(a, a), a), 2); 847 R = SkFPSub(R, SkFPMulInt(SkFPMul(a, b), 9)); 848 R = SkFPAdd(R, SkFPMulInt(c, 27)); 849 R = SkFPDivInt(R, 54); 850 851 SkFP Q3 = SkFPMul(SkFPMul(Q, Q), Q); 852 SkFP R2MinusQ3 = SkFPSub(SkFPMul(R,R), Q3); 853 SkFP adiv3 = SkFPDivInt(a, 3); 854 855 SkScalar* roots = tValues; 856 SkScalar r; 857 858 if (SkFPLT(R2MinusQ3, 0)) // we have 3 real roots 859 { 860#ifdef SK_SCALAR_IS_FLOAT 861 float theta = sk_float_acos(R / sk_float_sqrt(Q3)); 862 float neg2RootQ = -2 * sk_float_sqrt(Q); 863 864 r = neg2RootQ * sk_float_cos(theta/3) - adiv3; 865 if (is_unit_interval(r)) 866 *roots++ = r; 867 868 r = neg2RootQ * sk_float_cos((theta + 2*SK_ScalarPI)/3) - adiv3; 869 if (is_unit_interval(r)) 870 *roots++ = r; 871 872 r = neg2RootQ * sk_float_cos((theta - 2*SK_ScalarPI)/3) - adiv3; 873 if (is_unit_interval(r)) 874 *roots++ = r; 875 876 // now sort the roots 877 bubble_sort(tValues, (int)(roots - tValues)); 878#endif 879 } 880 else // we have 1 real root 881 { 882 SkFP A = SkFPAdd(SkFPAbs(R), SkFPSqrt(R2MinusQ3)); 883 A = SkFPCubeRoot(A); 884 if (SkFPGT(R, 0)) 885 A = SkFPNeg(A); 886 887 if (A != 0) 888 A = SkFPAdd(A, SkFPDiv(Q, A)); 889 r = SkFPToScalar(SkFPSub(A, adiv3)); 890 if (is_unit_interval(r)) 891 *roots++ = r; 892 } 893 894 return (int)(roots - tValues); 895} 896 897/* Looking for F' dot F'' == 0 898 899 A = b - a 900 B = c - 2b + a 901 C = d - 3c + 3b - a 902 903 F' = 3Ct^2 + 6Bt + 3A 904 F'' = 6Ct + 6B 905 906 F' dot F'' -> CCt^3 + 3BCt^2 + (2BB + CA)t + AB 907*/ 908static void formulate_F1DotF2(const SkScalar src[], SkFP coeff[4]) 909{ 910 SkScalar a = src[2] - src[0]; 911 SkScalar b = src[4] - 2 * src[2] + src[0]; 912 SkScalar c = src[6] + 3 * (src[2] - src[4]) - src[0]; 913 914 SkFP A = SkScalarToFP(a); 915 SkFP B = SkScalarToFP(b); 916 SkFP C = SkScalarToFP(c); 917 918 coeff[0] = SkFPMul(C, C); 919 coeff[1] = SkFPMulInt(SkFPMul(B, C), 3); 920 coeff[2] = SkFPMulInt(SkFPMul(B, B), 2); 921 coeff[2] = SkFPAdd(coeff[2], SkFPMul(C, A)); 922 coeff[3] = SkFPMul(A, B); 923} 924 925// EXPERIMENTAL: can set this to zero to accept all t-values 0 < t < 1 926//#define kMinTValueForChopping (SK_Scalar1 / 256) 927#define kMinTValueForChopping 0 928 929/* Looking for F' dot F'' == 0 930 931 A = b - a 932 B = c - 2b + a 933 C = d - 3c + 3b - a 934 935 F' = 3Ct^2 + 6Bt + 3A 936 F'' = 6Ct + 6B 937 938 F' dot F'' -> CCt^3 + 3BCt^2 + (2BB + CA)t + AB 939*/ 940int SkFindCubicMaxCurvature(const SkPoint src[4], SkScalar tValues[3]) 941{ 942 SkFP coeffX[4], coeffY[4]; 943 int i; 944 945 formulate_F1DotF2(&src[0].fX, coeffX); 946 formulate_F1DotF2(&src[0].fY, coeffY); 947 948 for (i = 0; i < 4; i++) 949 coeffX[i] = SkFPAdd(coeffX[i],coeffY[i]); 950 951 SkScalar t[3]; 952 int count = solve_cubic_polynomial(coeffX, t); 953 int maxCount = 0; 954 955 // now remove extrema where the curvature is zero (mins) 956 // !!!! need a test for this !!!! 957 for (i = 0; i < count; i++) 958 { 959 // if (not_min_curvature()) 960 if (t[i] > kMinTValueForChopping && t[i] < SK_Scalar1 - kMinTValueForChopping) 961 tValues[maxCount++] = t[i]; 962 } 963 return maxCount; 964} 965 966int SkChopCubicAtMaxCurvature(const SkPoint src[4], SkPoint dst[13], SkScalar tValues[3]) 967{ 968 SkScalar t_storage[3]; 969 970 if (tValues == NULL) 971 tValues = t_storage; 972 973 int count = SkFindCubicMaxCurvature(src, tValues); 974 975 if (dst) 976 { 977 if (count == 0) 978 memcpy(dst, src, 4 * sizeof(SkPoint)); 979 else 980 SkChopCubicAt(src, dst, tValues, count); 981 } 982 return count + 1; 983} 984 985bool SkXRayCrossesMonotonicCubic(const SkXRay& pt, const SkPoint cubic[4]) { 986 // Find the minimum and maximum y of the extrema, which are the 987 // first and last points since this cubic is monotonic 988 SkScalar min_y = SkMinScalar(cubic[0].fY, cubic[3].fY); 989 SkScalar max_y = SkMaxScalar(cubic[0].fY, cubic[3].fY); 990 991 if (pt.fY == cubic[0].fY 992 || pt.fY < min_y 993 || pt.fY > max_y) { 994 // The query line definitely does not cross the curve 995 return false; 996 } 997 998 SkScalar min_x = 999 SkMinScalar( 1000 SkMinScalar( 1001 SkMinScalar(cubic[0].fX, cubic[1].fX), 1002 cubic[2].fX), 1003 cubic[3].fX); 1004 if (pt.fX < min_x) { 1005 // The query line definitely crosses the curve 1006 return true; 1007 } 1008 1009 SkScalar max_x = 1010 SkMaxScalar( 1011 SkMaxScalar( 1012 SkMaxScalar(cubic[0].fX, cubic[1].fX), 1013 cubic[2].fX), 1014 cubic[3].fX); 1015 if (pt.fX > max_x) { 1016 // The query line definitely does not cross the curve 1017 return false; 1018 } 1019 1020 // Do a binary search to find the parameter value which makes y as 1021 // close as possible to the query point. See whether the query 1022 // line's origin is to the left of the associated x coordinate. 1023 1024 // kMaxIter is chosen as the number of mantissa bits for a float, 1025 // since there's no way we are going to get more precision by 1026 // iterating more times than that. 1027 const int kMaxIter = 23; 1028 SkPoint eval; 1029 int iter = 0; 1030 SkScalar upper_t; 1031 SkScalar lower_t; 1032 // Need to invert direction of t parameter if cubic goes up 1033 // instead of down 1034 if (cubic[3].fY > cubic[0].fY) { 1035 upper_t = SK_Scalar1; 1036 lower_t = SkFloatToScalar(0); 1037 } else { 1038 upper_t = SkFloatToScalar(0); 1039 lower_t = SK_Scalar1; 1040 } 1041 do { 1042 SkScalar t = SkScalarAve(upper_t, lower_t); 1043 SkEvalCubicAt(cubic, t, &eval, NULL, NULL); 1044 if (pt.fY > eval.fY) { 1045 lower_t = t; 1046 } else { 1047 upper_t = t; 1048 } 1049 } while (++iter < kMaxIter 1050 && !SkScalarNearlyZero(eval.fY - pt.fY)); 1051 if (pt.fX <= eval.fX) { 1052 return true; 1053 } 1054 return false; 1055} 1056 1057int SkNumXRayCrossingsForCubic(const SkXRay& pt, const SkPoint cubic[4]) { 1058 int num_crossings = 0; 1059 SkPoint monotonic_cubics[10]; 1060 int num_monotonic_cubics = SkChopCubicAtYExtrema(cubic, monotonic_cubics); 1061 if (SkXRayCrossesMonotonicCubic(pt, &monotonic_cubics[0])) 1062 ++num_crossings; 1063 if (num_monotonic_cubics > 0) 1064 if (SkXRayCrossesMonotonicCubic(pt, &monotonic_cubics[3])) 1065 ++num_crossings; 1066 if (num_monotonic_cubics > 1) 1067 if (SkXRayCrossesMonotonicCubic(pt, &monotonic_cubics[6])) 1068 ++num_crossings; 1069 return num_crossings; 1070} 1071 1072//////////////////////////////////////////////////////////////////////////////// 1073 1074/* Find t value for quadratic [a, b, c] = d. 1075 Return 0 if there is no solution within [0, 1) 1076*/ 1077static SkScalar quad_solve(SkScalar a, SkScalar b, SkScalar c, SkScalar d) 1078{ 1079 // At^2 + Bt + C = d 1080 SkScalar A = a - 2 * b + c; 1081 SkScalar B = 2 * (b - a); 1082 SkScalar C = a - d; 1083 1084 SkScalar roots[2]; 1085 int count = SkFindUnitQuadRoots(A, B, C, roots); 1086 1087 SkASSERT(count <= 1); 1088 return count == 1 ? roots[0] : 0; 1089} 1090 1091/* given a quad-curve and a point (x,y), chop the quad at that point and return 1092 the new quad's offCurve point. Should only return false if the computed pos 1093 is the start of the curve (i.e. root == 0) 1094*/ 1095static bool quad_pt2OffCurve(const SkPoint quad[3], SkScalar x, SkScalar y, SkPoint* offCurve) 1096{ 1097 const SkScalar* base; 1098 SkScalar value; 1099 1100 if (SkScalarAbs(x) < SkScalarAbs(y)) { 1101 base = &quad[0].fX; 1102 value = x; 1103 } else { 1104 base = &quad[0].fY; 1105 value = y; 1106 } 1107 1108 // note: this returns 0 if it thinks value is out of range, meaning the 1109 // root might return something outside of [0, 1) 1110 SkScalar t = quad_solve(base[0], base[2], base[4], value); 1111 1112 if (t > 0) 1113 { 1114 SkPoint tmp[5]; 1115 SkChopQuadAt(quad, tmp, t); 1116 *offCurve = tmp[1]; 1117 return true; 1118 } else { 1119 /* t == 0 means either the value triggered a root outside of [0, 1) 1120 For our purposes, we can ignore the <= 0 roots, but we want to 1121 catch the >= 1 roots (which given our caller, will basically mean 1122 a root of 1, give-or-take numerical instability). If we are in the 1123 >= 1 case, return the existing offCurve point. 1124 1125 The test below checks to see if we are close to the "end" of the 1126 curve (near base[4]). Rather than specifying a tolerance, I just 1127 check to see if value is on to the right/left of the middle point 1128 (depending on the direction/sign of the end points). 1129 */ 1130 if ((base[0] < base[4] && value > base[2]) || 1131 (base[0] > base[4] && value < base[2])) // should root have been 1 1132 { 1133 *offCurve = quad[1]; 1134 return true; 1135 } 1136 } 1137 return false; 1138} 1139 1140static const SkPoint gQuadCirclePts[kSkBuildQuadArcStorage] = { 1141 { SK_Scalar1, 0 }, 1142 { SK_Scalar1, SK_ScalarTanPIOver8 }, 1143 { SK_ScalarRoot2Over2, SK_ScalarRoot2Over2 }, 1144 { SK_ScalarTanPIOver8, SK_Scalar1 }, 1145 1146 { 0, SK_Scalar1 }, 1147 { -SK_ScalarTanPIOver8, SK_Scalar1 }, 1148 { -SK_ScalarRoot2Over2, SK_ScalarRoot2Over2 }, 1149 { -SK_Scalar1, SK_ScalarTanPIOver8 }, 1150 1151 { -SK_Scalar1, 0 }, 1152 { -SK_Scalar1, -SK_ScalarTanPIOver8 }, 1153 { -SK_ScalarRoot2Over2, -SK_ScalarRoot2Over2 }, 1154 { -SK_ScalarTanPIOver8, -SK_Scalar1 }, 1155 1156 { 0, -SK_Scalar1 }, 1157 { SK_ScalarTanPIOver8, -SK_Scalar1 }, 1158 { SK_ScalarRoot2Over2, -SK_ScalarRoot2Over2 }, 1159 { SK_Scalar1, -SK_ScalarTanPIOver8 }, 1160 1161 { SK_Scalar1, 0 } 1162}; 1163 1164int SkBuildQuadArc(const SkVector& uStart, const SkVector& uStop, 1165 SkRotationDirection dir, const SkMatrix* userMatrix, 1166 SkPoint quadPoints[]) 1167{ 1168 // rotate by x,y so that uStart is (1.0) 1169 SkScalar x = SkPoint::DotProduct(uStart, uStop); 1170 SkScalar y = SkPoint::CrossProduct(uStart, uStop); 1171 1172 SkScalar absX = SkScalarAbs(x); 1173 SkScalar absY = SkScalarAbs(y); 1174 1175 int pointCount; 1176 1177 // check for (effectively) coincident vectors 1178 // this can happen if our angle is nearly 0 or nearly 180 (y == 0) 1179 // ... we use the dot-prod to distinguish between 0 and 180 (x > 0) 1180 if (absY <= SK_ScalarNearlyZero && x > 0 && 1181 ((y >= 0 && kCW_SkRotationDirection == dir) || 1182 (y <= 0 && kCCW_SkRotationDirection == dir))) { 1183 1184 // just return the start-point 1185 quadPoints[0].set(SK_Scalar1, 0); 1186 pointCount = 1; 1187 } else { 1188 if (dir == kCCW_SkRotationDirection) 1189 y = -y; 1190 1191 // what octant (quadratic curve) is [xy] in? 1192 int oct = 0; 1193 bool sameSign = true; 1194 1195 if (0 == y) 1196 { 1197 oct = 4; // 180 1198 SkASSERT(SkScalarAbs(x + SK_Scalar1) <= SK_ScalarNearlyZero); 1199 } 1200 else if (0 == x) 1201 { 1202 SkASSERT(absY - SK_Scalar1 <= SK_ScalarNearlyZero); 1203 if (y > 0) 1204 oct = 2; // 90 1205 else 1206 oct = 6; // 270 1207 } 1208 else 1209 { 1210 if (y < 0) 1211 oct += 4; 1212 if ((x < 0) != (y < 0)) 1213 { 1214 oct += 2; 1215 sameSign = false; 1216 } 1217 if ((absX < absY) == sameSign) 1218 oct += 1; 1219 } 1220 1221 int wholeCount = oct << 1; 1222 memcpy(quadPoints, gQuadCirclePts, (wholeCount + 1) * sizeof(SkPoint)); 1223 1224 const SkPoint* arc = &gQuadCirclePts[wholeCount]; 1225 if (quad_pt2OffCurve(arc, x, y, &quadPoints[wholeCount + 1])) 1226 { 1227 quadPoints[wholeCount + 2].set(x, y); 1228 wholeCount += 2; 1229 } 1230 pointCount = wholeCount + 1; 1231 } 1232 1233 // now handle counter-clockwise and the initial unitStart rotation 1234 SkMatrix matrix; 1235 matrix.setSinCos(uStart.fY, uStart.fX); 1236 if (dir == kCCW_SkRotationDirection) { 1237 matrix.preScale(SK_Scalar1, -SK_Scalar1); 1238 } 1239 if (userMatrix) { 1240 matrix.postConcat(*userMatrix); 1241 } 1242 matrix.mapPoints(quadPoints, pointCount); 1243 return pointCount; 1244} 1245 1246