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