SkGeometry.cpp revision fc25abdabff76f913fb9d4f373418c10a1eca92b
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 263static inline void force_quad_monotonic_in_y(SkPoint pts[3]) 264{ 265 // zap pts[1].fY to the nearest value 266 SkScalar ab = SkScalarAbs(pts[0].fY - pts[1].fY); 267 SkScalar bc = SkScalarAbs(pts[1].fY - pts[2].fY); 268 pts[1].fY = ab < bc ? pts[0].fY : pts[2].fY; 269} 270 271/* Returns 0 for 1 quad, and 1 for two quads, either way the answer is 272 stored in dst[]. Guarantees that the 1/2 quads will be monotonic. 273*/ 274int SkChopQuadAtYExtrema(const SkPoint src[3], SkPoint dst[5]) 275{ 276 SkASSERT(src); 277 SkASSERT(dst); 278 279#if 0 280 static bool once = true; 281 if (once) 282 { 283 once = false; 284 SkPoint s[3] = { 0, 26398, 0, 26331, 0, 20621428 }; 285 SkPoint d[6]; 286 287 int n = SkChopQuadAtYExtrema(s, d); 288 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); 289 } 290#endif 291 292 SkScalar a = src[0].fY; 293 SkScalar b = src[1].fY; 294 SkScalar c = src[2].fY; 295 296 if (is_not_monotonic(a, b, c)) 297 { 298 SkScalar tValue; 299 if (valid_unit_divide(a - b, a - b - b + c, &tValue)) 300 { 301 SkChopQuadAt(src, dst, tValue); 302 flatten_double_quad_extrema(&dst[0].fY); 303 return 1; 304 } 305 // if we get here, we need to force dst to be monotonic, even though 306 // we couldn't compute a unit_divide value (probably underflow). 307 b = SkScalarAbs(a - b) < SkScalarAbs(b - c) ? a : c; 308 } 309 dst[0].set(src[0].fX, a); 310 dst[1].set(src[1].fX, b); 311 dst[2].set(src[2].fX, c); 312 return 0; 313} 314 315// F(t) = a (1 - t) ^ 2 + 2 b t (1 - t) + c t ^ 2 316// F'(t) = 2 (b - a) + 2 (a - 2b + c) t 317// F''(t) = 2 (a - 2b + c) 318// 319// A = 2 (b - a) 320// B = 2 (a - 2b + c) 321// 322// Maximum curvature for a quadratic means solving 323// Fx' Fx'' + Fy' Fy'' = 0 324// 325// t = - (Ax Bx + Ay By) / (Bx ^ 2 + By ^ 2) 326// 327int SkChopQuadAtMaxCurvature(const SkPoint src[3], SkPoint dst[5]) 328{ 329 SkScalar Ax = src[1].fX - src[0].fX; 330 SkScalar Ay = src[1].fY - src[0].fY; 331 SkScalar Bx = src[0].fX - src[1].fX - src[1].fX + src[2].fX; 332 SkScalar By = src[0].fY - src[1].fY - src[1].fY + src[2].fY; 333 SkScalar t = 0; // 0 means don't chop 334 335#ifdef SK_SCALAR_IS_FLOAT 336 (void)valid_unit_divide(-(Ax * Bx + Ay * By), Bx * Bx + By * By, &t); 337#else 338 // !!! should I use SkFloat here? seems like it 339 Sk64 numer, denom, tmp; 340 341 numer.setMul(Ax, -Bx); 342 tmp.setMul(Ay, -By); 343 numer.add(tmp); 344 345 if (numer.isPos()) // do nothing if numer <= 0 346 { 347 denom.setMul(Bx, Bx); 348 tmp.setMul(By, By); 349 denom.add(tmp); 350 SkASSERT(!denom.isNeg()); 351 if (numer < denom) 352 { 353 t = numer.getFixedDiv(denom); 354 SkASSERT(t >= 0 && t <= SK_Fixed1); // assert that we're numerically stable (ha!) 355 if ((unsigned)t >= SK_Fixed1) // runtime check for numerical stability 356 t = 0; // ignore the chop 357 } 358 } 359#endif 360 361 if (t == 0) 362 { 363 memcpy(dst, src, 3 * sizeof(SkPoint)); 364 return 1; 365 } 366 else 367 { 368 SkChopQuadAt(src, dst, t); 369 return 2; 370 } 371} 372 373//////////////////////////////////////////////////////////////////////////////////////// 374///// CUBICS // CUBICS // CUBICS // CUBICS // CUBICS // CUBICS // CUBICS // CUBICS ///// 375//////////////////////////////////////////////////////////////////////////////////////// 376 377static void get_cubic_coeff(const SkScalar pt[], SkScalar coeff[4]) 378{ 379 coeff[0] = pt[6] + 3*(pt[2] - pt[4]) - pt[0]; 380 coeff[1] = 3*(pt[4] - pt[2] - pt[2] + pt[0]); 381 coeff[2] = 3*(pt[2] - pt[0]); 382 coeff[3] = pt[0]; 383} 384 385void SkGetCubicCoeff(const SkPoint pts[4], SkScalar cx[4], SkScalar cy[4]) 386{ 387 SkASSERT(pts); 388 389 if (cx) 390 get_cubic_coeff(&pts[0].fX, cx); 391 if (cy) 392 get_cubic_coeff(&pts[0].fY, cy); 393} 394 395static SkScalar eval_cubic(const SkScalar src[], SkScalar t) 396{ 397 SkASSERT(src); 398 SkASSERT(t >= 0 && t <= SK_Scalar1); 399 400 if (t == 0) 401 return src[0]; 402 403#ifdef DIRECT_EVAL_OF_POLYNOMIALS 404 SkScalar D = src[0]; 405 SkScalar A = src[6] + 3*(src[2] - src[4]) - D; 406 SkScalar B = 3*(src[4] - src[2] - src[2] + D); 407 SkScalar C = 3*(src[2] - D); 408 409 return SkScalarMulAdd(SkScalarMulAdd(SkScalarMulAdd(A, t, B), t, C), t, D); 410#else 411 SkScalar ab = SkScalarInterp(src[0], src[2], t); 412 SkScalar bc = SkScalarInterp(src[2], src[4], t); 413 SkScalar cd = SkScalarInterp(src[4], src[6], t); 414 SkScalar abc = SkScalarInterp(ab, bc, t); 415 SkScalar bcd = SkScalarInterp(bc, cd, t); 416 return SkScalarInterp(abc, bcd, t); 417#endif 418} 419 420/** return At^2 + Bt + C 421*/ 422static SkScalar eval_quadratic(SkScalar A, SkScalar B, SkScalar C, SkScalar t) 423{ 424 SkASSERT(t >= 0 && t <= SK_Scalar1); 425 426 return SkScalarMulAdd(SkScalarMulAdd(A, t, B), t, C); 427} 428 429static SkScalar eval_cubic_derivative(const SkScalar src[], SkScalar t) 430{ 431 SkScalar A = src[6] + 3*(src[2] - src[4]) - src[0]; 432 SkScalar B = 2*(src[4] - 2 * src[2] + src[0]); 433 SkScalar C = src[2] - src[0]; 434 435 return eval_quadratic(A, B, C, t); 436} 437 438static SkScalar eval_cubic_2ndDerivative(const SkScalar src[], SkScalar t) 439{ 440 SkScalar A = src[6] + 3*(src[2] - src[4]) - src[0]; 441 SkScalar B = src[4] - 2 * src[2] + src[0]; 442 443 return SkScalarMulAdd(A, t, B); 444} 445 446void SkEvalCubicAt(const SkPoint src[4], SkScalar t, SkPoint* loc, SkVector* tangent, SkVector* curvature) 447{ 448 SkASSERT(src); 449 SkASSERT(t >= 0 && t <= SK_Scalar1); 450 451 if (loc) 452 loc->set(eval_cubic(&src[0].fX, t), eval_cubic(&src[0].fY, t)); 453 if (tangent) 454 tangent->set(eval_cubic_derivative(&src[0].fX, t), 455 eval_cubic_derivative(&src[0].fY, t)); 456 if (curvature) 457 curvature->set(eval_cubic_2ndDerivative(&src[0].fX, t), 458 eval_cubic_2ndDerivative(&src[0].fY, t)); 459} 460 461/** Cubic'(t) = At^2 + Bt + C, where 462 A = 3(-a + 3(b - c) + d) 463 B = 6(a - 2b + c) 464 C = 3(b - a) 465 Solve for t, keeping only those that fit betwee 0 < t < 1 466*/ 467int SkFindCubicExtrema(SkScalar a, SkScalar b, SkScalar c, SkScalar d, SkScalar tValues[2]) 468{ 469#ifdef SK_SCALAR_IS_FIXED 470 if (!is_not_monotonic(a, b, c, d)) 471 return 0; 472#endif 473 474 // we divide A,B,C by 3 to simplify 475 SkScalar A = d - a + 3*(b - c); 476 SkScalar B = 2*(a - b - b + c); 477 SkScalar C = b - a; 478 479 return SkFindUnitQuadRoots(A, B, C, tValues); 480} 481 482static void interp_cubic_coords(const SkScalar* src, SkScalar* dst, SkScalar t) 483{ 484 SkScalar ab = SkScalarInterp(src[0], src[2], t); 485 SkScalar bc = SkScalarInterp(src[2], src[4], t); 486 SkScalar cd = SkScalarInterp(src[4], src[6], t); 487 SkScalar abc = SkScalarInterp(ab, bc, t); 488 SkScalar bcd = SkScalarInterp(bc, cd, t); 489 SkScalar abcd = SkScalarInterp(abc, bcd, t); 490 491 dst[0] = src[0]; 492 dst[2] = ab; 493 dst[4] = abc; 494 dst[6] = abcd; 495 dst[8] = bcd; 496 dst[10] = cd; 497 dst[12] = src[6]; 498} 499 500void SkChopCubicAt(const SkPoint src[4], SkPoint dst[7], SkScalar t) 501{ 502 SkASSERT(t > 0 && t < SK_Scalar1); 503 504 interp_cubic_coords(&src[0].fX, &dst[0].fX, t); 505 interp_cubic_coords(&src[0].fY, &dst[0].fY, t); 506} 507 508void SkChopCubicAt(const SkPoint src[4], SkPoint dst[], const SkScalar tValues[], int roots) 509{ 510#ifdef SK_DEBUG 511 { 512 for (int i = 0; i < roots - 1; i++) 513 { 514 SkASSERT(is_unit_interval(tValues[i])); 515 SkASSERT(is_unit_interval(tValues[i+1])); 516 SkASSERT(tValues[i] < tValues[i+1]); 517 } 518 } 519#endif 520 521 if (dst) 522 { 523 if (roots == 0) // nothing to chop 524 memcpy(dst, src, 4*sizeof(SkPoint)); 525 else 526 { 527 SkScalar t = tValues[0]; 528 SkPoint tmp[4]; 529 530 for (int i = 0; i < roots; i++) 531 { 532 SkChopCubicAt(src, dst, t); 533 if (i == roots - 1) 534 break; 535 536 SkDEBUGCODE(int valid =) valid_unit_divide(tValues[i+1] - tValues[i], SK_Scalar1 - tValues[i], &t); 537 SkASSERT(valid); 538 539 dst += 3; 540 memcpy(tmp, dst, 4 * sizeof(SkPoint)); 541 src = tmp; 542 } 543 } 544 } 545} 546 547void SkChopCubicAtHalf(const SkPoint src[4], SkPoint dst[7]) 548{ 549 SkScalar x01 = SkScalarAve(src[0].fX, src[1].fX); 550 SkScalar y01 = SkScalarAve(src[0].fY, src[1].fY); 551 SkScalar x12 = SkScalarAve(src[1].fX, src[2].fX); 552 SkScalar y12 = SkScalarAve(src[1].fY, src[2].fY); 553 SkScalar x23 = SkScalarAve(src[2].fX, src[3].fX); 554 SkScalar y23 = SkScalarAve(src[2].fY, src[3].fY); 555 556 SkScalar x012 = SkScalarAve(x01, x12); 557 SkScalar y012 = SkScalarAve(y01, y12); 558 SkScalar x123 = SkScalarAve(x12, x23); 559 SkScalar y123 = SkScalarAve(y12, y23); 560 561 dst[0] = src[0]; 562 dst[1].set(x01, y01); 563 dst[2].set(x012, y012); 564 dst[3].set(SkScalarAve(x012, x123), SkScalarAve(y012, y123)); 565 dst[4].set(x123, y123); 566 dst[5].set(x23, y23); 567 dst[6] = src[3]; 568} 569 570static void flatten_double_cubic_extrema(SkScalar coords[14]) 571{ 572 coords[4] = coords[8] = coords[6]; 573} 574 575/** Given 4 points on a cubic bezier, chop it into 1, 2, 3 beziers such that 576 the resulting beziers are monotonic in Y. This is called by the scan converter. 577 Depending on what is returned, dst[] is treated as follows 578 0 dst[0..3] is the original cubic 579 1 dst[0..3] and dst[3..6] are the two new cubics 580 2 dst[0..3], dst[3..6], dst[6..9] are the three new cubics 581 If dst == null, it is ignored and only the count is returned. 582*/ 583int SkChopCubicAtYExtrema(const SkPoint src[4], SkPoint dst[10]) 584{ 585 SkScalar tValues[2]; 586 int roots = SkFindCubicExtrema(src[0].fY, src[1].fY, src[2].fY, src[3].fY, tValues); 587 588 SkChopCubicAt(src, dst, tValues, roots); 589 if (dst && roots > 0) 590 { 591 // we do some cleanup to ensure our Y extrema are flat 592 flatten_double_cubic_extrema(&dst[0].fY); 593 if (roots == 2) 594 flatten_double_cubic_extrema(&dst[3].fY); 595 } 596 return roots; 597} 598 599/** http://www.faculty.idc.ac.il/arik/quality/appendixA.html 600 601 Inflection means that curvature is zero. 602 Curvature is [F' x F''] / [F'^3] 603 So we solve F'x X F''y - F'y X F''y == 0 604 After some canceling of the cubic term, we get 605 A = b - a 606 B = c - 2b + a 607 C = d - 3c + 3b - a 608 (BxCy - ByCx)t^2 + (AxCy - AyCx)t + AxBy - AyBx == 0 609*/ 610int SkFindCubicInflections(const SkPoint src[4], SkScalar tValues[]) 611{ 612 SkScalar Ax = src[1].fX - src[0].fX; 613 SkScalar Ay = src[1].fY - src[0].fY; 614 SkScalar Bx = src[2].fX - 2 * src[1].fX + src[0].fX; 615 SkScalar By = src[2].fY - 2 * src[1].fY + src[0].fY; 616 SkScalar Cx = src[3].fX + 3 * (src[1].fX - src[2].fX) - src[0].fX; 617 SkScalar Cy = src[3].fY + 3 * (src[1].fY - src[2].fY) - src[0].fY; 618 int count; 619 620#ifdef SK_SCALAR_IS_FLOAT 621 count = SkFindUnitQuadRoots(Bx*Cy - By*Cx, Ax*Cy - Ay*Cx, Ax*By - Ay*Bx, tValues); 622#else 623 Sk64 A, B, C, tmp; 624 625 A.setMul(Bx, Cy); 626 tmp.setMul(By, Cx); 627 A.sub(tmp); 628 629 B.setMul(Ax, Cy); 630 tmp.setMul(Ay, Cx); 631 B.sub(tmp); 632 633 C.setMul(Ax, By); 634 tmp.setMul(Ay, Bx); 635 C.sub(tmp); 636 637 count = Sk64FindFixedQuadRoots(A, B, C, tValues); 638#endif 639 640 return count; 641} 642 643int SkChopCubicAtInflections(const SkPoint src[], SkPoint dst[10]) 644{ 645 SkScalar tValues[2]; 646 int count = SkFindCubicInflections(src, tValues); 647 648 if (dst) 649 { 650 if (count == 0) 651 memcpy(dst, src, 4 * sizeof(SkPoint)); 652 else 653 SkChopCubicAt(src, dst, tValues, count); 654 } 655 return count + 1; 656} 657 658template <typename T> void bubble_sort(T array[], int count) 659{ 660 for (int i = count - 1; i > 0; --i) 661 for (int j = i; j > 0; --j) 662 if (array[j] < array[j-1]) 663 { 664 T tmp(array[j]); 665 array[j] = array[j-1]; 666 array[j-1] = tmp; 667 } 668} 669 670#include "SkFP.h" 671 672// newton refinement 673#if 0 674static SkScalar refine_cubic_root(const SkFP coeff[4], SkScalar root) 675{ 676 // x1 = x0 - f(t) / f'(t) 677 678 SkFP T = SkScalarToFloat(root); 679 SkFP N, D; 680 681 // f' = 3*coeff[0]*T^2 + 2*coeff[1]*T + coeff[2] 682 D = SkFPMul(SkFPMul(coeff[0], SkFPMul(T,T)), 3); 683 D = SkFPAdd(D, SkFPMulInt(SkFPMul(coeff[1], T), 2)); 684 D = SkFPAdd(D, coeff[2]); 685 686 if (D == 0) 687 return root; 688 689 // f = coeff[0]*T^3 + coeff[1]*T^2 + coeff[2]*T + coeff[3] 690 N = SkFPMul(SkFPMul(SkFPMul(T, T), T), coeff[0]); 691 N = SkFPAdd(N, SkFPMul(SkFPMul(T, T), coeff[1])); 692 N = SkFPAdd(N, SkFPMul(T, coeff[2])); 693 N = SkFPAdd(N, coeff[3]); 694 695 if (N) 696 { 697 SkScalar delta = SkFPToScalar(SkFPDiv(N, D)); 698 699 if (delta) 700 root -= delta; 701 } 702 return root; 703} 704#endif 705 706#if defined _WIN32 && _MSC_VER >= 1300 && defined SK_SCALAR_IS_FIXED // disable warning : unreachable code if building fixed point for windows desktop 707#pragma warning ( disable : 4702 ) 708#endif 709 710/* Solve coeff(t) == 0, returning the number of roots that 711 lie withing 0 < t < 1. 712 coeff[0]t^3 + coeff[1]t^2 + coeff[2]t + coeff[3] 713*/ 714static int solve_cubic_polynomial(const SkFP coeff[4], SkScalar tValues[3]) 715{ 716#ifndef SK_SCALAR_IS_FLOAT 717 return 0; // this is not yet implemented for software float 718#endif 719 720 if (SkScalarNearlyZero(coeff[0])) // we're just a quadratic 721 { 722 return SkFindUnitQuadRoots(coeff[1], coeff[2], coeff[3], tValues); 723 } 724 725 SkFP a, b, c, Q, R; 726 727 { 728 SkASSERT(coeff[0] != 0); 729 730 SkFP inva = SkFPInvert(coeff[0]); 731 a = SkFPMul(coeff[1], inva); 732 b = SkFPMul(coeff[2], inva); 733 c = SkFPMul(coeff[3], inva); 734 } 735 Q = SkFPDivInt(SkFPSub(SkFPMul(a,a), SkFPMulInt(b, 3)), 9); 736// R = (2*a*a*a - 9*a*b + 27*c) / 54; 737 R = SkFPMulInt(SkFPMul(SkFPMul(a, a), a), 2); 738 R = SkFPSub(R, SkFPMulInt(SkFPMul(a, b), 9)); 739 R = SkFPAdd(R, SkFPMulInt(c, 27)); 740 R = SkFPDivInt(R, 54); 741 742 SkFP Q3 = SkFPMul(SkFPMul(Q, Q), Q); 743 SkFP R2MinusQ3 = SkFPSub(SkFPMul(R,R), Q3); 744 SkFP adiv3 = SkFPDivInt(a, 3); 745 746 SkScalar* roots = tValues; 747 SkScalar r; 748 749 if (SkFPLT(R2MinusQ3, 0)) // we have 3 real roots 750 { 751#ifdef SK_SCALAR_IS_FLOAT 752 float theta = sk_float_acos(R / sk_float_sqrt(Q3)); 753 float neg2RootQ = -2 * sk_float_sqrt(Q); 754 755 r = neg2RootQ * sk_float_cos(theta/3) - adiv3; 756 if (is_unit_interval(r)) 757 *roots++ = r; 758 759 r = neg2RootQ * sk_float_cos((theta + 2*SK_ScalarPI)/3) - adiv3; 760 if (is_unit_interval(r)) 761 *roots++ = r; 762 763 r = neg2RootQ * sk_float_cos((theta - 2*SK_ScalarPI)/3) - adiv3; 764 if (is_unit_interval(r)) 765 *roots++ = r; 766 767 // now sort the roots 768 bubble_sort(tValues, (int)(roots - tValues)); 769#endif 770 } 771 else // we have 1 real root 772 { 773 SkFP A = SkFPAdd(SkFPAbs(R), SkFPSqrt(R2MinusQ3)); 774 A = SkFPCubeRoot(A); 775 if (SkFPGT(R, 0)) 776 A = SkFPNeg(A); 777 778 if (A != 0) 779 A = SkFPAdd(A, SkFPDiv(Q, A)); 780 r = SkFPToScalar(SkFPSub(A, adiv3)); 781 if (is_unit_interval(r)) 782 *roots++ = r; 783 } 784 785 return (int)(roots - tValues); 786} 787 788/* Looking for F' dot F'' == 0 789 790 A = b - a 791 B = c - 2b + a 792 C = d - 3c + 3b - a 793 794 F' = 3Ct^2 + 6Bt + 3A 795 F'' = 6Ct + 6B 796 797 F' dot F'' -> CCt^3 + 3BCt^2 + (2BB + CA)t + AB 798*/ 799static void formulate_F1DotF2(const SkScalar src[], SkFP coeff[4]) 800{ 801 SkScalar a = src[2] - src[0]; 802 SkScalar b = src[4] - 2 * src[2] + src[0]; 803 SkScalar c = src[6] + 3 * (src[2] - src[4]) - src[0]; 804 805 SkFP A = SkScalarToFP(a); 806 SkFP B = SkScalarToFP(b); 807 SkFP C = SkScalarToFP(c); 808 809 coeff[0] = SkFPMul(C, C); 810 coeff[1] = SkFPMulInt(SkFPMul(B, C), 3); 811 coeff[2] = SkFPMulInt(SkFPMul(B, B), 2); 812 coeff[2] = SkFPAdd(coeff[2], SkFPMul(C, A)); 813 coeff[3] = SkFPMul(A, B); 814} 815 816// EXPERIMENTAL: can set this to zero to accept all t-values 0 < t < 1 817//#define kMinTValueForChopping (SK_Scalar1 / 256) 818#define kMinTValueForChopping 0 819 820/* Looking for F' dot F'' == 0 821 822 A = b - a 823 B = c - 2b + a 824 C = d - 3c + 3b - a 825 826 F' = 3Ct^2 + 6Bt + 3A 827 F'' = 6Ct + 6B 828 829 F' dot F'' -> CCt^3 + 3BCt^2 + (2BB + CA)t + AB 830*/ 831int SkFindCubicMaxCurvature(const SkPoint src[4], SkScalar tValues[3]) 832{ 833 SkFP coeffX[4], coeffY[4]; 834 int i; 835 836 formulate_F1DotF2(&src[0].fX, coeffX); 837 formulate_F1DotF2(&src[0].fY, coeffY); 838 839 for (i = 0; i < 4; i++) 840 coeffX[i] = SkFPAdd(coeffX[i],coeffY[i]); 841 842 SkScalar t[3]; 843 int count = solve_cubic_polynomial(coeffX, t); 844 int maxCount = 0; 845 846 // now remove extrema where the curvature is zero (mins) 847 // !!!! need a test for this !!!! 848 for (i = 0; i < count; i++) 849 { 850 // if (not_min_curvature()) 851 if (t[i] > kMinTValueForChopping && t[i] < SK_Scalar1 - kMinTValueForChopping) 852 tValues[maxCount++] = t[i]; 853 } 854 return maxCount; 855} 856 857int SkChopCubicAtMaxCurvature(const SkPoint src[4], SkPoint dst[13], SkScalar tValues[3]) 858{ 859 SkScalar t_storage[3]; 860 861 if (tValues == NULL) 862 tValues = t_storage; 863 864 int count = SkFindCubicMaxCurvature(src, tValues); 865 866 if (dst) 867 { 868 if (count == 0) 869 memcpy(dst, src, 4 * sizeof(SkPoint)); 870 else 871 SkChopCubicAt(src, dst, tValues, count); 872 } 873 return count + 1; 874} 875 876//////////////////////////////////////////////////////////////////////////////// 877 878/* Find t value for quadratic [a, b, c] = d. 879 Return 0 if there is no solution within [0, 1) 880*/ 881static SkScalar quad_solve(SkScalar a, SkScalar b, SkScalar c, SkScalar d) 882{ 883 // At^2 + Bt + C = d 884 SkScalar A = a - 2 * b + c; 885 SkScalar B = 2 * (b - a); 886 SkScalar C = a - d; 887 888 SkScalar roots[2]; 889 int count = SkFindUnitQuadRoots(A, B, C, roots); 890 891 SkASSERT(count <= 1); 892 return count == 1 ? roots[0] : 0; 893} 894 895/* given a quad-curve and a point (x,y), chop the quad at that point and return 896 the new quad's offCurve point. Should only return false if the computed pos 897 is the start of the curve (i.e. root == 0) 898*/ 899static bool quad_pt2OffCurve(const SkPoint quad[3], SkScalar x, SkScalar y, SkPoint* offCurve) 900{ 901 const SkScalar* base; 902 SkScalar value; 903 904 if (SkScalarAbs(x) < SkScalarAbs(y)) { 905 base = &quad[0].fX; 906 value = x; 907 } else { 908 base = &quad[0].fY; 909 value = y; 910 } 911 912 // note: this returns 0 if it thinks value is out of range, meaning the 913 // root might return something outside of [0, 1) 914 SkScalar t = quad_solve(base[0], base[2], base[4], value); 915 916 if (t > 0) 917 { 918 SkPoint tmp[5]; 919 SkChopQuadAt(quad, tmp, t); 920 *offCurve = tmp[1]; 921 return true; 922 } else { 923 /* t == 0 means either the value triggered a root outside of [0, 1) 924 For our purposes, we can ignore the <= 0 roots, but we want to 925 catch the >= 1 roots (which given our caller, will basically mean 926 a root of 1, give-or-take numerical instability). If we are in the 927 >= 1 case, return the existing offCurve point. 928 929 The test below checks to see if we are close to the "end" of the 930 curve (near base[4]). Rather than specifying a tolerance, I just 931 check to see if value is on to the right/left of the middle point 932 (depending on the direction/sign of the end points). 933 */ 934 if ((base[0] < base[4] && value > base[2]) || 935 (base[0] > base[4] && value < base[2])) // should root have been 1 936 { 937 *offCurve = quad[1]; 938 return true; 939 } 940 } 941 return false; 942} 943 944static const SkPoint gQuadCirclePts[kSkBuildQuadArcStorage] = { 945 { SK_Scalar1, 0 }, 946 { SK_Scalar1, SK_ScalarTanPIOver8 }, 947 { SK_ScalarRoot2Over2, SK_ScalarRoot2Over2 }, 948 { SK_ScalarTanPIOver8, SK_Scalar1 }, 949 950 { 0, SK_Scalar1 }, 951 { -SK_ScalarTanPIOver8, SK_Scalar1 }, 952 { -SK_ScalarRoot2Over2, SK_ScalarRoot2Over2 }, 953 { -SK_Scalar1, SK_ScalarTanPIOver8 }, 954 955 { -SK_Scalar1, 0 }, 956 { -SK_Scalar1, -SK_ScalarTanPIOver8 }, 957 { -SK_ScalarRoot2Over2, -SK_ScalarRoot2Over2 }, 958 { -SK_ScalarTanPIOver8, -SK_Scalar1 }, 959 960 { 0, -SK_Scalar1 }, 961 { SK_ScalarTanPIOver8, -SK_Scalar1 }, 962 { SK_ScalarRoot2Over2, -SK_ScalarRoot2Over2 }, 963 { SK_Scalar1, -SK_ScalarTanPIOver8 }, 964 965 { SK_Scalar1, 0 } 966}; 967 968int SkBuildQuadArc(const SkVector& uStart, const SkVector& uStop, 969 SkRotationDirection dir, const SkMatrix* userMatrix, 970 SkPoint quadPoints[]) 971{ 972 // rotate by x,y so that uStart is (1.0) 973 SkScalar x = SkPoint::DotProduct(uStart, uStop); 974 SkScalar y = SkPoint::CrossProduct(uStart, uStop); 975 976 SkScalar absX = SkScalarAbs(x); 977 SkScalar absY = SkScalarAbs(y); 978 979 int pointCount; 980 981 // check for (effectively) coincident vectors 982 // this can happen if our angle is nearly 0 or nearly 180 (y == 0) 983 // ... we use the dot-prod to distinguish between 0 and 180 (x > 0) 984 if (absY <= SK_ScalarNearlyZero && x > 0 && 985 ((y >= 0 && kCW_SkRotationDirection == dir) || 986 (y <= 0 && kCCW_SkRotationDirection == dir))) { 987 988 // just return the start-point 989 quadPoints[0].set(SK_Scalar1, 0); 990 pointCount = 1; 991 } else { 992 if (dir == kCCW_SkRotationDirection) 993 y = -y; 994 995 // what octant (quadratic curve) is [xy] in? 996 int oct = 0; 997 bool sameSign = true; 998 999 if (0 == y) 1000 { 1001 oct = 4; // 180 1002 SkASSERT(SkScalarAbs(x + SK_Scalar1) <= SK_ScalarNearlyZero); 1003 } 1004 else if (0 == x) 1005 { 1006 SkASSERT(absY - SK_Scalar1 <= SK_ScalarNearlyZero); 1007 if (y > 0) 1008 oct = 2; // 90 1009 else 1010 oct = 6; // 270 1011 } 1012 else 1013 { 1014 if (y < 0) 1015 oct += 4; 1016 if ((x < 0) != (y < 0)) 1017 { 1018 oct += 2; 1019 sameSign = false; 1020 } 1021 if ((absX < absY) == sameSign) 1022 oct += 1; 1023 } 1024 1025 int wholeCount = oct << 1; 1026 memcpy(quadPoints, gQuadCirclePts, (wholeCount + 1) * sizeof(SkPoint)); 1027 1028 const SkPoint* arc = &gQuadCirclePts[wholeCount]; 1029 if (quad_pt2OffCurve(arc, x, y, &quadPoints[wholeCount + 1])) 1030 { 1031 quadPoints[wholeCount + 2].set(x, y); 1032 wholeCount += 2; 1033 } 1034 pointCount = wholeCount + 1; 1035 } 1036 1037 // now handle counter-clockwise and the initial unitStart rotation 1038 SkMatrix matrix; 1039 matrix.setSinCos(uStart.fY, uStart.fX); 1040 if (dir == kCCW_SkRotationDirection) { 1041 matrix.preScale(SK_Scalar1, -SK_Scalar1); 1042 } 1043 if (userMatrix) { 1044 matrix.postConcat(*userMatrix); 1045 } 1046 matrix.mapPoints(quadPoints, pointCount); 1047 return pointCount; 1048} 1049 1050 1051///////////////////////////////////////////////////////////////////////////////////////// 1052///////////////////////////////////////////////////////////////////////////////////////// 1053 1054#ifdef SK_DEBUG 1055 1056void SkGeometry::UnitTest() 1057{ 1058#ifdef SK_SUPPORT_UNITTEST 1059 SkPoint pts[3], dst[5]; 1060 1061 pts[0].set(0, 0); 1062 pts[1].set(100, 50); 1063 pts[2].set(0, 100); 1064 1065 int count = SkChopQuadAtMaxCurvature(pts, dst); 1066 SkASSERT(count == 1 || count == 2); 1067#endif 1068} 1069 1070#endif 1071 1072 1073