Home | History | Annotate | Download | only in core
      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 int SkChopQuadAtMaxCurvature(const SkPoint src[3], SkPoint dst[5])
    398 {
    399     SkScalar    Ax = src[1].fX - src[0].fX;
    400     SkScalar    Ay = src[1].fY - src[0].fY;
    401     SkScalar    Bx = src[0].fX - src[1].fX - src[1].fX + src[2].fX;
    402     SkScalar    By = src[0].fY - src[1].fY - src[1].fY + src[2].fY;
    403     SkScalar    t = 0;  // 0 means don't chop
    404 
    405 #ifdef SK_SCALAR_IS_FLOAT
    406     (void)valid_unit_divide(-(Ax * Bx + Ay * By), Bx * Bx + By * By, &t);
    407 #else
    408     // !!! should I use SkFloat here? seems like it
    409     Sk64    numer, denom, tmp;
    410 
    411     numer.setMul(Ax, -Bx);
    412     tmp.setMul(Ay, -By);
    413     numer.add(tmp);
    414 
    415     if (numer.isPos())  // do nothing if numer <= 0
    416     {
    417         denom.setMul(Bx, Bx);
    418         tmp.setMul(By, By);
    419         denom.add(tmp);
    420         SkASSERT(!denom.isNeg());
    421         if (numer < denom)
    422         {
    423             t = numer.getFixedDiv(denom);
    424             SkASSERT(t >= 0 && t <= SK_Fixed1);     // assert that we're numerically stable (ha!)
    425             if ((unsigned)t >= SK_Fixed1)           // runtime check for numerical stability
    426                 t = 0;  // ignore the chop
    427         }
    428     }
    429 #endif
    430 
    431     if (t == 0)
    432     {
    433         memcpy(dst, src, 3 * sizeof(SkPoint));
    434         return 1;
    435     }
    436     else
    437     {
    438         SkChopQuadAt(src, dst, t);
    439         return 2;
    440     }
    441 }
    442 
    443 #ifdef SK_SCALAR_IS_FLOAT
    444     #define SK_ScalarTwoThirds  (0.666666666f)
    445 #else
    446     #define SK_ScalarTwoThirds  ((SkFixed)(43691))
    447 #endif
    448 
    449 void SkConvertQuadToCubic(const SkPoint src[3], SkPoint dst[4]) {
    450     const SkScalar scale = SK_ScalarTwoThirds;
    451     dst[0] = src[0];
    452     dst[1].set(src[0].fX + SkScalarMul(src[1].fX - src[0].fX, scale),
    453                src[0].fY + SkScalarMul(src[1].fY - src[0].fY, scale));
    454     dst[2].set(src[2].fX + SkScalarMul(src[1].fX - src[2].fX, scale),
    455                src[2].fY + SkScalarMul(src[1].fY - src[2].fY, scale));
    456     dst[3] = src[2];
    457 }
    458 
    459 ////////////////////////////////////////////////////////////////////////////////////////
    460 ///// CUBICS // CUBICS // CUBICS // CUBICS // CUBICS // CUBICS // CUBICS // CUBICS /////
    461 ////////////////////////////////////////////////////////////////////////////////////////
    462 
    463 static void get_cubic_coeff(const SkScalar pt[], SkScalar coeff[4])
    464 {
    465     coeff[0] = pt[6] + 3*(pt[2] - pt[4]) - pt[0];
    466     coeff[1] = 3*(pt[4] - pt[2] - pt[2] + pt[0]);
    467     coeff[2] = 3*(pt[2] - pt[0]);
    468     coeff[3] = pt[0];
    469 }
    470 
    471 void SkGetCubicCoeff(const SkPoint pts[4], SkScalar cx[4], SkScalar cy[4])
    472 {
    473     SkASSERT(pts);
    474 
    475     if (cx)
    476         get_cubic_coeff(&pts[0].fX, cx);
    477     if (cy)
    478         get_cubic_coeff(&pts[0].fY, cy);
    479 }
    480 
    481 static SkScalar eval_cubic(const SkScalar src[], SkScalar t)
    482 {
    483     SkASSERT(src);
    484     SkASSERT(t >= 0 && t <= SK_Scalar1);
    485 
    486     if (t == 0)
    487         return src[0];
    488 
    489 #ifdef DIRECT_EVAL_OF_POLYNOMIALS
    490     SkScalar D = src[0];
    491     SkScalar A = src[6] + 3*(src[2] - src[4]) - D;
    492     SkScalar B = 3*(src[4] - src[2] - src[2] + D);
    493     SkScalar C = 3*(src[2] - D);
    494 
    495     return SkScalarMulAdd(SkScalarMulAdd(SkScalarMulAdd(A, t, B), t, C), t, D);
    496 #else
    497     SkScalar    ab = SkScalarInterp(src[0], src[2], t);
    498     SkScalar    bc = SkScalarInterp(src[2], src[4], t);
    499     SkScalar    cd = SkScalarInterp(src[4], src[6], t);
    500     SkScalar    abc = SkScalarInterp(ab, bc, t);
    501     SkScalar    bcd = SkScalarInterp(bc, cd, t);
    502     return SkScalarInterp(abc, bcd, t);
    503 #endif
    504 }
    505 
    506 /** return At^2 + Bt + C
    507 */
    508 static SkScalar eval_quadratic(SkScalar A, SkScalar B, SkScalar C, SkScalar t)
    509 {
    510     SkASSERT(t >= 0 && t <= SK_Scalar1);
    511 
    512     return SkScalarMulAdd(SkScalarMulAdd(A, t, B), t, C);
    513 }
    514 
    515 static SkScalar eval_cubic_derivative(const SkScalar src[], SkScalar t)
    516 {
    517     SkScalar A = src[6] + 3*(src[2] - src[4]) - src[0];
    518     SkScalar B = 2*(src[4] - 2 * src[2] + src[0]);
    519     SkScalar C = src[2] - src[0];
    520 
    521     return eval_quadratic(A, B, C, t);
    522 }
    523 
    524 static SkScalar eval_cubic_2ndDerivative(const SkScalar src[], SkScalar t)
    525 {
    526     SkScalar A = src[6] + 3*(src[2] - src[4]) - src[0];
    527     SkScalar B = src[4] - 2 * src[2] + src[0];
    528 
    529     return SkScalarMulAdd(A, t, B);
    530 }
    531 
    532 void SkEvalCubicAt(const SkPoint src[4], SkScalar t, SkPoint* loc, SkVector* tangent, SkVector* curvature)
    533 {
    534     SkASSERT(src);
    535     SkASSERT(t >= 0 && t <= SK_Scalar1);
    536 
    537     if (loc)
    538         loc->set(eval_cubic(&src[0].fX, t), eval_cubic(&src[0].fY, t));
    539     if (tangent)
    540         tangent->set(eval_cubic_derivative(&src[0].fX, t),
    541                      eval_cubic_derivative(&src[0].fY, t));
    542     if (curvature)
    543         curvature->set(eval_cubic_2ndDerivative(&src[0].fX, t),
    544                        eval_cubic_2ndDerivative(&src[0].fY, t));
    545 }
    546 
    547 /** Cubic'(t) = At^2 + Bt + C, where
    548     A = 3(-a + 3(b - c) + d)
    549     B = 6(a - 2b + c)
    550     C = 3(b - a)
    551     Solve for t, keeping only those that fit betwee 0 < t < 1
    552 */
    553 int SkFindCubicExtrema(SkScalar a, SkScalar b, SkScalar c, SkScalar d, SkScalar tValues[2])
    554 {
    555 #ifdef SK_SCALAR_IS_FIXED
    556     if (!is_not_monotonic(a, b, c, d))
    557         return 0;
    558 #endif
    559 
    560     // we divide A,B,C by 3 to simplify
    561     SkScalar A = d - a + 3*(b - c);
    562     SkScalar B = 2*(a - b - b + c);
    563     SkScalar C = b - a;
    564 
    565     return SkFindUnitQuadRoots(A, B, C, tValues);
    566 }
    567 
    568 static void interp_cubic_coords(const SkScalar* src, SkScalar* dst, SkScalar t)
    569 {
    570     SkScalar    ab = SkScalarInterp(src[0], src[2], t);
    571     SkScalar    bc = SkScalarInterp(src[2], src[4], t);
    572     SkScalar    cd = SkScalarInterp(src[4], src[6], t);
    573     SkScalar    abc = SkScalarInterp(ab, bc, t);
    574     SkScalar    bcd = SkScalarInterp(bc, cd, t);
    575     SkScalar    abcd = SkScalarInterp(abc, bcd, t);
    576 
    577     dst[0] = src[0];
    578     dst[2] = ab;
    579     dst[4] = abc;
    580     dst[6] = abcd;
    581     dst[8] = bcd;
    582     dst[10] = cd;
    583     dst[12] = src[6];
    584 }
    585 
    586 void SkChopCubicAt(const SkPoint src[4], SkPoint dst[7], SkScalar t)
    587 {
    588     SkASSERT(t > 0 && t < SK_Scalar1);
    589 
    590     interp_cubic_coords(&src[0].fX, &dst[0].fX, t);
    591     interp_cubic_coords(&src[0].fY, &dst[0].fY, t);
    592 }
    593 
    594 /*  http://code.google.com/p/skia/issues/detail?id=32
    595 
    596     This test code would fail when we didn't check the return result of
    597     valid_unit_divide in SkChopCubicAt(... tValues[], int roots). The reason is
    598     that after the first chop, the parameters to valid_unit_divide are equal
    599     (thanks to finite float precision and rounding in the subtracts). Thus
    600     even though the 2nd tValue looks < 1.0, after we renormalize it, we end
    601     up with 1.0, hence the need to check and just return the last cubic as
    602     a degenerate clump of 4 points in the sampe place.
    603 
    604     static void test_cubic() {
    605         SkPoint src[4] = {
    606             { 556.25000, 523.03003 },
    607             { 556.23999, 522.96002 },
    608             { 556.21997, 522.89001 },
    609             { 556.21997, 522.82001 }
    610         };
    611         SkPoint dst[10];
    612         SkScalar tval[] = { 0.33333334f, 0.99999994f };
    613         SkChopCubicAt(src, dst, tval, 2);
    614     }
    615  */
    616 
    617 void SkChopCubicAt(const SkPoint src[4], SkPoint dst[], const SkScalar tValues[], int roots)
    618 {
    619 #ifdef SK_DEBUG
    620     {
    621         for (int i = 0; i < roots - 1; i++)
    622         {
    623             SkASSERT(is_unit_interval(tValues[i]));
    624             SkASSERT(is_unit_interval(tValues[i+1]));
    625             SkASSERT(tValues[i] < tValues[i+1]);
    626         }
    627     }
    628 #endif
    629 
    630     if (dst)
    631     {
    632         if (roots == 0) // nothing to chop
    633             memcpy(dst, src, 4*sizeof(SkPoint));
    634         else
    635         {
    636             SkScalar    t = tValues[0];
    637             SkPoint     tmp[4];
    638 
    639             for (int i = 0; i < roots; i++)
    640             {
    641                 SkChopCubicAt(src, dst, t);
    642                 if (i == roots - 1)
    643                     break;
    644 
    645                 dst += 3;
    646                 // have src point to the remaining cubic (after the chop)
    647                 memcpy(tmp, dst, 4 * sizeof(SkPoint));
    648                 src = tmp;
    649 
    650                 // watch out in case the renormalized t isn't in range
    651                 if (!valid_unit_divide(tValues[i+1] - tValues[i],
    652                                        SK_Scalar1 - tValues[i], &t)) {
    653                     // if we can't, just create a degenerate cubic
    654                     dst[4] = dst[5] = dst[6] = src[3];
    655                     break;
    656                 }
    657             }
    658         }
    659     }
    660 }
    661 
    662 void SkChopCubicAtHalf(const SkPoint src[4], SkPoint dst[7])
    663 {
    664     SkScalar x01 = SkScalarAve(src[0].fX, src[1].fX);
    665     SkScalar y01 = SkScalarAve(src[0].fY, src[1].fY);
    666     SkScalar x12 = SkScalarAve(src[1].fX, src[2].fX);
    667     SkScalar y12 = SkScalarAve(src[1].fY, src[2].fY);
    668     SkScalar x23 = SkScalarAve(src[2].fX, src[3].fX);
    669     SkScalar y23 = SkScalarAve(src[2].fY, src[3].fY);
    670 
    671     SkScalar x012 = SkScalarAve(x01, x12);
    672     SkScalar y012 = SkScalarAve(y01, y12);
    673     SkScalar x123 = SkScalarAve(x12, x23);
    674     SkScalar y123 = SkScalarAve(y12, y23);
    675 
    676     dst[0] = src[0];
    677     dst[1].set(x01, y01);
    678     dst[2].set(x012, y012);
    679     dst[3].set(SkScalarAve(x012, x123), SkScalarAve(y012, y123));
    680     dst[4].set(x123, y123);
    681     dst[5].set(x23, y23);
    682     dst[6] = src[3];
    683 }
    684 
    685 static void flatten_double_cubic_extrema(SkScalar coords[14])
    686 {
    687     coords[4] = coords[8] = coords[6];
    688 }
    689 
    690 /** Given 4 points on a cubic bezier, chop it into 1, 2, 3 beziers such that
    691     the resulting beziers are monotonic in Y. This is called by the scan converter.
    692     Depending on what is returned, dst[] is treated as follows
    693     0   dst[0..3] is the original cubic
    694     1   dst[0..3] and dst[3..6] are the two new cubics
    695     2   dst[0..3], dst[3..6], dst[6..9] are the three new cubics
    696     If dst == null, it is ignored and only the count is returned.
    697 */
    698 int SkChopCubicAtYExtrema(const SkPoint src[4], SkPoint dst[10]) {
    699     SkScalar    tValues[2];
    700     int         roots = SkFindCubicExtrema(src[0].fY, src[1].fY, src[2].fY,
    701                                            src[3].fY, tValues);
    702 
    703     SkChopCubicAt(src, dst, tValues, roots);
    704     if (dst && roots > 0) {
    705         // we do some cleanup to ensure our Y extrema are flat
    706         flatten_double_cubic_extrema(&dst[0].fY);
    707         if (roots == 2) {
    708             flatten_double_cubic_extrema(&dst[3].fY);
    709         }
    710     }
    711     return roots;
    712 }
    713 
    714 int SkChopCubicAtXExtrema(const SkPoint src[4], SkPoint dst[10]) {
    715     SkScalar    tValues[2];
    716     int         roots = SkFindCubicExtrema(src[0].fX, src[1].fX, src[2].fX,
    717                                            src[3].fX, tValues);
    718 
    719     SkChopCubicAt(src, dst, tValues, roots);
    720     if (dst && roots > 0) {
    721         // we do some cleanup to ensure our Y extrema are flat
    722         flatten_double_cubic_extrema(&dst[0].fX);
    723         if (roots == 2) {
    724             flatten_double_cubic_extrema(&dst[3].fX);
    725         }
    726     }
    727     return roots;
    728 }
    729 
    730 /** http://www.faculty.idc.ac.il/arik/quality/appendixA.html
    731 
    732     Inflection means that curvature is zero.
    733     Curvature is [F' x F''] / [F'^3]
    734     So we solve F'x X F''y - F'y X F''y == 0
    735     After some canceling of the cubic term, we get
    736     A = b - a
    737     B = c - 2b + a
    738     C = d - 3c + 3b - a
    739     (BxCy - ByCx)t^2 + (AxCy - AyCx)t + AxBy - AyBx == 0
    740 */
    741 int SkFindCubicInflections(const SkPoint src[4], SkScalar tValues[])
    742 {
    743     SkScalar    Ax = src[1].fX - src[0].fX;
    744     SkScalar    Ay = src[1].fY - src[0].fY;
    745     SkScalar    Bx = src[2].fX - 2 * src[1].fX + src[0].fX;
    746     SkScalar    By = src[2].fY - 2 * src[1].fY + src[0].fY;
    747     SkScalar    Cx = src[3].fX + 3 * (src[1].fX - src[2].fX) - src[0].fX;
    748     SkScalar    Cy = src[3].fY + 3 * (src[1].fY - src[2].fY) - src[0].fY;
    749     int         count;
    750 
    751 #ifdef SK_SCALAR_IS_FLOAT
    752     count = SkFindUnitQuadRoots(Bx*Cy - By*Cx, Ax*Cy - Ay*Cx, Ax*By - Ay*Bx, tValues);
    753 #else
    754     Sk64    A, B, C, tmp;
    755 
    756     A.setMul(Bx, Cy);
    757     tmp.setMul(By, Cx);
    758     A.sub(tmp);
    759 
    760     B.setMul(Ax, Cy);
    761     tmp.setMul(Ay, Cx);
    762     B.sub(tmp);
    763 
    764     C.setMul(Ax, By);
    765     tmp.setMul(Ay, Bx);
    766     C.sub(tmp);
    767 
    768     count = Sk64FindFixedQuadRoots(A, B, C, tValues);
    769 #endif
    770 
    771     return count;
    772 }
    773 
    774 int SkChopCubicAtInflections(const SkPoint src[], SkPoint dst[10])
    775 {
    776     SkScalar    tValues[2];
    777     int         count = SkFindCubicInflections(src, tValues);
    778 
    779     if (dst)
    780     {
    781         if (count == 0)
    782             memcpy(dst, src, 4 * sizeof(SkPoint));
    783         else
    784             SkChopCubicAt(src, dst, tValues, count);
    785     }
    786     return count + 1;
    787 }
    788 
    789 template <typename T> void bubble_sort(T array[], int count)
    790 {
    791     for (int i = count - 1; i > 0; --i)
    792         for (int j = i; j > 0; --j)
    793             if (array[j] < array[j-1])
    794             {
    795                 T   tmp(array[j]);
    796                 array[j] = array[j-1];
    797                 array[j-1] = tmp;
    798             }
    799 }
    800 
    801 #include "SkFP.h"
    802 
    803 // newton refinement
    804 #if 0
    805 static SkScalar refine_cubic_root(const SkFP coeff[4], SkScalar root)
    806 {
    807     //  x1 = x0 - f(t) / f'(t)
    808 
    809     SkFP    T = SkScalarToFloat(root);
    810     SkFP    N, D;
    811 
    812     // f' = 3*coeff[0]*T^2 + 2*coeff[1]*T + coeff[2]
    813     D = SkFPMul(SkFPMul(coeff[0], SkFPMul(T,T)), 3);
    814     D = SkFPAdd(D, SkFPMulInt(SkFPMul(coeff[1], T), 2));
    815     D = SkFPAdd(D, coeff[2]);
    816 
    817     if (D == 0)
    818         return root;
    819 
    820     // f = coeff[0]*T^3 + coeff[1]*T^2 + coeff[2]*T + coeff[3]
    821     N = SkFPMul(SkFPMul(SkFPMul(T, T), T), coeff[0]);
    822     N = SkFPAdd(N, SkFPMul(SkFPMul(T, T), coeff[1]));
    823     N = SkFPAdd(N, SkFPMul(T, coeff[2]));
    824     N = SkFPAdd(N, coeff[3]);
    825 
    826     if (N)
    827     {
    828         SkScalar delta = SkFPToScalar(SkFPDiv(N, D));
    829 
    830         if (delta)
    831             root -= delta;
    832     }
    833     return root;
    834 }
    835 #endif
    836 
    837 /**
    838  *  Given an array and count, remove all pair-wise duplicates from the array,
    839  *  keeping the existing sorting, and return the new count
    840  */
    841 static int collaps_duplicates(float array[], int count) {
    842     int n = 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     {
   1067         if (count == 0)
   1068             memcpy(dst, src, 4 * sizeof(SkPoint));
   1069         else
   1070             SkChopCubicAt(src, dst, tValues, count);
   1071     }
   1072     return count + 1;
   1073 }
   1074 
   1075 bool SkXRayCrossesMonotonicCubic(const SkXRay& pt, const SkPoint cubic[4], bool* ambiguous) {
   1076     if (ambiguous) {
   1077         *ambiguous = false;
   1078     }
   1079 
   1080     // Find the minimum and maximum y of the extrema, which are the
   1081     // first and last points since this cubic is monotonic
   1082     SkScalar min_y = SkMinScalar(cubic[0].fY, cubic[3].fY);
   1083     SkScalar max_y = SkMaxScalar(cubic[0].fY, cubic[3].fY);
   1084 
   1085     if (pt.fY == cubic[0].fY
   1086         || pt.fY < min_y
   1087         || pt.fY > max_y) {
   1088         // The query line definitely does not cross the curve
   1089         if (ambiguous) {
   1090             *ambiguous = (pt.fY == cubic[0].fY);
   1091         }
   1092         return false;
   1093     }
   1094 
   1095     bool pt_at_extremum = (pt.fY == cubic[3].fY);
   1096 
   1097     SkScalar min_x =
   1098         SkMinScalar(
   1099             SkMinScalar(
   1100                 SkMinScalar(cubic[0].fX, cubic[1].fX),
   1101                 cubic[2].fX),
   1102             cubic[3].fX);
   1103     if (pt.fX < min_x) {
   1104         // The query line definitely crosses the curve
   1105         if (ambiguous) {
   1106             *ambiguous = pt_at_extremum;
   1107         }
   1108         return true;
   1109     }
   1110 
   1111     SkScalar max_x =
   1112         SkMaxScalar(
   1113             SkMaxScalar(
   1114                 SkMaxScalar(cubic[0].fX, cubic[1].fX),
   1115                 cubic[2].fX),
   1116             cubic[3].fX);
   1117     if (pt.fX > max_x) {
   1118         // The query line definitely does not cross the curve
   1119         return false;
   1120     }
   1121 
   1122     // Do a binary search to find the parameter value which makes y as
   1123     // close as possible to the query point. See whether the query
   1124     // line's origin is to the left of the associated x coordinate.
   1125 
   1126     // kMaxIter is chosen as the number of mantissa bits for a float,
   1127     // since there's no way we are going to get more precision by
   1128     // iterating more times than that.
   1129     const int kMaxIter = 23;
   1130     SkPoint eval;
   1131     int iter = 0;
   1132     SkScalar upper_t;
   1133     SkScalar lower_t;
   1134     // Need to invert direction of t parameter if cubic goes up
   1135     // instead of down
   1136     if (cubic[3].fY > cubic[0].fY) {
   1137         upper_t = SK_Scalar1;
   1138         lower_t = SkFloatToScalar(0);
   1139     } else {
   1140         upper_t = SkFloatToScalar(0);
   1141         lower_t = SK_Scalar1;
   1142     }
   1143     do {
   1144         SkScalar t = SkScalarAve(upper_t, lower_t);
   1145         SkEvalCubicAt(cubic, t, &eval, NULL, NULL);
   1146         if (pt.fY > eval.fY) {
   1147             lower_t = t;
   1148         } else {
   1149             upper_t = t;
   1150         }
   1151     } while (++iter < kMaxIter
   1152              && !SkScalarNearlyZero(eval.fY - pt.fY));
   1153     if (pt.fX <= eval.fX) {
   1154         if (ambiguous) {
   1155             *ambiguous = pt_at_extremum;
   1156         }
   1157         return true;
   1158     }
   1159     return false;
   1160 }
   1161 
   1162 int SkNumXRayCrossingsForCubic(const SkXRay& pt, const SkPoint cubic[4], bool* ambiguous) {
   1163     int num_crossings = 0;
   1164     SkPoint monotonic_cubics[10];
   1165     int num_monotonic_cubics = SkChopCubicAtYExtrema(cubic, monotonic_cubics);
   1166     if (ambiguous) {
   1167         *ambiguous = false;
   1168     }
   1169     bool locally_ambiguous;
   1170     if (SkXRayCrossesMonotonicCubic(pt, &monotonic_cubics[0], &locally_ambiguous))
   1171         ++num_crossings;
   1172     if (ambiguous) {
   1173         *ambiguous |= locally_ambiguous;
   1174     }
   1175     if (num_monotonic_cubics > 0)
   1176         if (SkXRayCrossesMonotonicCubic(pt, &monotonic_cubics[3], &locally_ambiguous))
   1177             ++num_crossings;
   1178     if (ambiguous) {
   1179         *ambiguous |= locally_ambiguous;
   1180     }
   1181     if (num_monotonic_cubics > 1)
   1182         if (SkXRayCrossesMonotonicCubic(pt, &monotonic_cubics[6], &locally_ambiguous))
   1183             ++num_crossings;
   1184     if (ambiguous) {
   1185         *ambiguous |= locally_ambiguous;
   1186     }
   1187     return num_crossings;
   1188 }
   1189 
   1190 ////////////////////////////////////////////////////////////////////////////////
   1191 
   1192 /*  Find t value for quadratic [a, b, c] = d.
   1193     Return 0 if there is no solution within [0, 1)
   1194 */
   1195 static SkScalar quad_solve(SkScalar a, SkScalar b, SkScalar c, SkScalar d)
   1196 {
   1197     // At^2 + Bt + C = d
   1198     SkScalar A = a - 2 * b + c;
   1199     SkScalar B = 2 * (b - a);
   1200     SkScalar C = a - d;
   1201 
   1202     SkScalar    roots[2];
   1203     int         count = SkFindUnitQuadRoots(A, B, C, roots);
   1204 
   1205     SkASSERT(count <= 1);
   1206     return count == 1 ? roots[0] : 0;
   1207 }
   1208 
   1209 /*  given a quad-curve and a point (x,y), chop the quad at that point and return
   1210     the new quad's offCurve point. Should only return false if the computed pos
   1211     is the start of the curve (i.e. root == 0)
   1212 */
   1213 static bool quad_pt2OffCurve(const SkPoint quad[3], SkScalar x, SkScalar y, SkPoint* offCurve)
   1214 {
   1215     const SkScalar* base;
   1216     SkScalar        value;
   1217 
   1218     if (SkScalarAbs(x) < SkScalarAbs(y)) {
   1219         base = &quad[0].fX;
   1220         value = x;
   1221     } else {
   1222         base = &quad[0].fY;
   1223         value = y;
   1224     }
   1225 
   1226     // note: this returns 0 if it thinks value is out of range, meaning the
   1227     // root might return something outside of [0, 1)
   1228     SkScalar t = quad_solve(base[0], base[2], base[4], value);
   1229 
   1230     if (t > 0)
   1231     {
   1232         SkPoint tmp[5];
   1233         SkChopQuadAt(quad, tmp, t);
   1234         *offCurve = tmp[1];
   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             *offCurve = quad[1];
   1252             return true;
   1253         }
   1254     }
   1255     return false;
   1256 }
   1257 
   1258 static const SkPoint gQuadCirclePts[kSkBuildQuadArcStorage] = {
   1259     { SK_Scalar1,           0               },
   1260     { SK_Scalar1,           SK_ScalarTanPIOver8 },
   1261     { SK_ScalarRoot2Over2,  SK_ScalarRoot2Over2 },
   1262     { SK_ScalarTanPIOver8,  SK_Scalar1          },
   1263 
   1264     { 0,                    SK_Scalar1      },
   1265     { -SK_ScalarTanPIOver8, SK_Scalar1  },
   1266     { -SK_ScalarRoot2Over2, SK_ScalarRoot2Over2 },
   1267     { -SK_Scalar1,          SK_ScalarTanPIOver8 },
   1268 
   1269     { -SK_Scalar1,          0               },
   1270     { -SK_Scalar1,          -SK_ScalarTanPIOver8    },
   1271     { -SK_ScalarRoot2Over2, -SK_ScalarRoot2Over2    },
   1272     { -SK_ScalarTanPIOver8, -SK_Scalar1     },
   1273 
   1274     { 0,                    -SK_Scalar1     },
   1275     { SK_ScalarTanPIOver8,  -SK_Scalar1     },
   1276     { SK_ScalarRoot2Over2,  -SK_ScalarRoot2Over2    },
   1277     { SK_Scalar1,           -SK_ScalarTanPIOver8    },
   1278 
   1279     { SK_Scalar1,           0   }
   1280 };
   1281 
   1282 int SkBuildQuadArc(const SkVector& uStart, const SkVector& uStop,
   1283                    SkRotationDirection dir, const SkMatrix* userMatrix,
   1284                    SkPoint quadPoints[])
   1285 {
   1286     // rotate by x,y so that uStart is (1.0)
   1287     SkScalar x = SkPoint::DotProduct(uStart, uStop);
   1288     SkScalar y = SkPoint::CrossProduct(uStart, uStop);
   1289 
   1290     SkScalar absX = SkScalarAbs(x);
   1291     SkScalar absY = SkScalarAbs(y);
   1292 
   1293     int pointCount;
   1294 
   1295     // check for (effectively) coincident vectors
   1296     // this can happen if our angle is nearly 0 or nearly 180 (y == 0)
   1297     // ... we use the dot-prod to distinguish between 0 and 180 (x > 0)
   1298     if (absY <= SK_ScalarNearlyZero && x > 0 &&
   1299         ((y >= 0 && kCW_SkRotationDirection == dir) ||
   1300          (y <= 0 && kCCW_SkRotationDirection == dir))) {
   1301 
   1302         // just return the start-point
   1303         quadPoints[0].set(SK_Scalar1, 0);
   1304         pointCount = 1;
   1305     } else {
   1306         if (dir == kCCW_SkRotationDirection)
   1307             y = -y;
   1308 
   1309         // what octant (quadratic curve) is [xy] in?
   1310         int oct = 0;
   1311         bool sameSign = true;
   1312 
   1313         if (0 == y)
   1314         {
   1315             oct = 4;        // 180
   1316             SkASSERT(SkScalarAbs(x + SK_Scalar1) <= SK_ScalarNearlyZero);
   1317         }
   1318         else if (0 == x)
   1319         {
   1320             SkASSERT(absY - SK_Scalar1 <= SK_ScalarNearlyZero);
   1321             if (y > 0)
   1322                 oct = 2;    // 90
   1323             else
   1324                 oct = 6;    // 270
   1325         }
   1326         else
   1327         {
   1328             if (y < 0)
   1329                 oct += 4;
   1330             if ((x < 0) != (y < 0))
   1331             {
   1332                 oct += 2;
   1333                 sameSign = false;
   1334             }
   1335             if ((absX < absY) == sameSign)
   1336                 oct += 1;
   1337         }
   1338 
   1339         int wholeCount = oct << 1;
   1340         memcpy(quadPoints, gQuadCirclePts, (wholeCount + 1) * sizeof(SkPoint));
   1341 
   1342         const SkPoint* arc = &gQuadCirclePts[wholeCount];
   1343         if (quad_pt2OffCurve(arc, x, y, &quadPoints[wholeCount + 1]))
   1344         {
   1345             quadPoints[wholeCount + 2].set(x, y);
   1346             wholeCount += 2;
   1347         }
   1348         pointCount = wholeCount + 1;
   1349     }
   1350 
   1351     // now handle counter-clockwise and the initial unitStart rotation
   1352     SkMatrix    matrix;
   1353     matrix.setSinCos(uStart.fY, uStart.fX);
   1354     if (dir == kCCW_SkRotationDirection) {
   1355         matrix.preScale(SK_Scalar1, -SK_Scalar1);
   1356     }
   1357     if (userMatrix) {
   1358         matrix.postConcat(*userMatrix);
   1359     }
   1360     matrix.mapPoints(quadPoints, pointCount);
   1361     return pointCount;
   1362 }
   1363 
   1364