src/third_party/skia/src/core/SkGeometry.cpp

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