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

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