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     for (int n = count; n > 1; --n) {
    843         if (array[0] == array[1]) {
    844             for (int i = 1; i < n; ++i) {
    845                 array[i - 1] = array[i];
    846             }
    847             count -= 1;
    848         } else {
    849             array += 1;
    850         }
    851     }
    852     return count;
    853 }
    854 
    855 #ifdef SK_DEBUG
    856 
    857 #define TEST_COLLAPS_ENTRY(array)   array, SK_ARRAY_COUNT(array)
    858 
    859 static void test_collaps_duplicates() {
    860     static bool gOnce;
    861     if (gOnce) { return; }
    862     gOnce = true;
    863     const float src0[] = { 0 };
    864     const float src1[] = { 0, 0 };
    865     const float src2[] = { 0, 1 };
    866     const float src3[] = { 0, 0, 0 };
    867     const float src4[] = { 0, 0, 1 };
    868     const float src5[] = { 0, 1, 1 };
    869     const float src6[] = { 0, 1, 2 };
    870     const struct {
    871         const float* fData;
    872         int fCount;
    873         int fCollapsedCount;
    874     } data[] = {
    875         { TEST_COLLAPS_ENTRY(src0), 1 },
    876         { TEST_COLLAPS_ENTRY(src1), 1 },
    877         { TEST_COLLAPS_ENTRY(src2), 2 },
    878         { TEST_COLLAPS_ENTRY(src3), 1 },
    879         { TEST_COLLAPS_ENTRY(src4), 2 },
    880         { TEST_COLLAPS_ENTRY(src5), 2 },
    881         { TEST_COLLAPS_ENTRY(src6), 3 },
    882     };
    883     for (size_t i = 0; i < SK_ARRAY_COUNT(data); ++i) {
    884         float dst[3];
    885         memcpy(dst, data[i].fData, data[i].fCount * sizeof(dst[0]));
    886         int count = collaps_duplicates(dst, data[i].fCount);
    887         SkASSERT(data[i].fCollapsedCount == count);
    888         for (int j = 1; j < count; ++j) {
    889             SkASSERT(dst[j-1] < dst[j]);
    890         }
    891     }
    892 }
    893 #endif
    894 
    895 #if defined _WIN32 && _MSC_VER >= 1300  && defined SK_SCALAR_IS_FIXED // disable warning : unreachable code if building fixed point for windows desktop
    896 #pragma warning ( disable : 4702 )
    897 #endif
    898 
    899 /*  Solve coeff(t) == 0, returning the number of roots that
    900     lie withing 0 < t < 1.
    901     coeff[0]t^3 + coeff[1]t^2 + coeff[2]t + coeff[3]
    902 
    903     Eliminates repeated roots (so that all tValues are distinct, and are always
    904     in increasing order.
    905 */
    906 static int solve_cubic_polynomial(const SkFP coeff[4], SkScalar tValues[3])
    907 {
    908 #ifndef SK_SCALAR_IS_FLOAT
    909     return 0;   // this is not yet implemented for software float
    910 #endif
    911 
    912     if (SkScalarNearlyZero(coeff[0]))   // we're just a quadratic
    913     {
    914         return SkFindUnitQuadRoots(coeff[1], coeff[2], coeff[3], tValues);
    915     }
    916 
    917     SkFP    a, b, c, Q, R;
    918 
    919     {
    920         SkASSERT(coeff[0] != 0);
    921 
    922         SkFP inva = SkFPInvert(coeff[0]);
    923         a = SkFPMul(coeff[1], inva);
    924         b = SkFPMul(coeff[2], inva);
    925         c = SkFPMul(coeff[3], inva);
    926     }
    927     Q = SkFPDivInt(SkFPSub(SkFPMul(a,a), SkFPMulInt(b, 3)), 9);
    928 //  R = (2*a*a*a - 9*a*b + 27*c) / 54;
    929     R = SkFPMulInt(SkFPMul(SkFPMul(a, a), a), 2);
    930     R = SkFPSub(R, SkFPMulInt(SkFPMul(a, b), 9));
    931     R = SkFPAdd(R, SkFPMulInt(c, 27));
    932     R = SkFPDivInt(R, 54);
    933 
    934     SkFP Q3 = SkFPMul(SkFPMul(Q, Q), Q);
    935     SkFP R2MinusQ3 = SkFPSub(SkFPMul(R,R), Q3);
    936     SkFP adiv3 = SkFPDivInt(a, 3);
    937 
    938     SkScalar*   roots = tValues;
    939     SkScalar    r;
    940 
    941     if (SkFPLT(R2MinusQ3, 0))   // we have 3 real roots
    942     {
    943 #ifdef SK_SCALAR_IS_FLOAT
    944         float theta = sk_float_acos(R / sk_float_sqrt(Q3));
    945         float neg2RootQ = -2 * sk_float_sqrt(Q);
    946 
    947         r = neg2RootQ * sk_float_cos(theta/3) - adiv3;
    948         if (is_unit_interval(r))
    949             *roots++ = r;
    950 
    951         r = neg2RootQ * sk_float_cos((theta + 2*SK_ScalarPI)/3) - adiv3;
    952         if (is_unit_interval(r))
    953             *roots++ = r;
    954 
    955         r = neg2RootQ * sk_float_cos((theta - 2*SK_ScalarPI)/3) - adiv3;
    956         if (is_unit_interval(r))
    957             *roots++ = r;
    958 
    959         SkDEBUGCODE(test_collaps_duplicates();)
    960 
    961         // now sort the roots
    962         int count = (int)(roots - tValues);
    963         SkASSERT((unsigned)count <= 3);
    964         bubble_sort(tValues, count);
    965         count = collaps_duplicates(tValues, count);
    966         roots = tValues + count;    // so we compute the proper count below
    967 #endif
    968     }
    969     else                // we have 1 real root
    970     {
    971         SkFP A = SkFPAdd(SkFPAbs(R), SkFPSqrt(R2MinusQ3));
    972         A = SkFPCubeRoot(A);
    973         if (SkFPGT(R, 0))
    974             A = SkFPNeg(A);
    975 
    976         if (A != 0)
    977             A = SkFPAdd(A, SkFPDiv(Q, A));
    978         r = SkFPToScalar(SkFPSub(A, adiv3));
    979         if (is_unit_interval(r))
    980             *roots++ = r;
    981     }
    982 
    983     return (int)(roots - tValues);
    984 }
    985 
    986 /*  Looking for F' dot F'' == 0
    987 
    988     A = b - a
    989     B = c - 2b + a
    990     C = d - 3c + 3b - a
    991 
    992     F' = 3Ct^2 + 6Bt + 3A
    993     F'' = 6Ct + 6B
    994 
    995     F' dot F'' -> CCt^3 + 3BCt^2 + (2BB + CA)t + AB
    996 */
    997 static void formulate_F1DotF2(const SkScalar src[], SkFP coeff[4])
    998 {
    999     SkScalar    a = src[2] - src[0];
   1000     SkScalar    b = src[4] - 2 * src[2] + src[0];
   1001     SkScalar    c = src[6] + 3 * (src[2] - src[4]) - src[0];
   1002 
   1003     SkFP    A = SkScalarToFP(a);
   1004     SkFP    B = SkScalarToFP(b);
   1005     SkFP    C = SkScalarToFP(c);
   1006 
   1007     coeff[0] = SkFPMul(C, C);
   1008     coeff[1] = SkFPMulInt(SkFPMul(B, C), 3);
   1009     coeff[2] = SkFPMulInt(SkFPMul(B, B), 2);
   1010     coeff[2] = SkFPAdd(coeff[2], SkFPMul(C, A));
   1011     coeff[3] = SkFPMul(A, B);
   1012 }
   1013 
   1014 // EXPERIMENTAL: can set this to zero to accept all t-values 0 < t < 1
   1015 //#define kMinTValueForChopping (SK_Scalar1 / 256)
   1016 #define kMinTValueForChopping   0
   1017 
   1018 /*  Looking for F' dot F'' == 0
   1019 
   1020     A = b - a
   1021     B = c - 2b + a
   1022     C = d - 3c + 3b - a
   1023 
   1024     F' = 3Ct^2 + 6Bt + 3A
   1025     F'' = 6Ct + 6B
   1026 
   1027     F' dot F'' -> CCt^3 + 3BCt^2 + (2BB + CA)t + AB
   1028 */
   1029 int SkFindCubicMaxCurvature(const SkPoint src[4], SkScalar tValues[3])
   1030 {
   1031     SkFP    coeffX[4], coeffY[4];
   1032     int     i;
   1033 
   1034     formulate_F1DotF2(&src[0].fX, coeffX);
   1035     formulate_F1DotF2(&src[0].fY, coeffY);
   1036 
   1037     for (i = 0; i < 4; i++)
   1038         coeffX[i] = SkFPAdd(coeffX[i],coeffY[i]);
   1039 
   1040     SkScalar    t[3];
   1041     int         count = solve_cubic_polynomial(coeffX, t);
   1042     int         maxCount = 0;
   1043 
   1044     // now remove extrema where the curvature is zero (mins)
   1045     // !!!! need a test for this !!!!
   1046     for (i = 0; i < count; i++)
   1047     {
   1048         // if (not_min_curvature())
   1049         if (t[i] > kMinTValueForChopping && t[i] < SK_Scalar1 - kMinTValueForChopping)
   1050             tValues[maxCount++] = t[i];
   1051     }
   1052     return maxCount;
   1053 }
   1054 
   1055 int SkChopCubicAtMaxCurvature(const SkPoint src[4], SkPoint dst[13], SkScalar tValues[3])
   1056 {
   1057     SkScalar    t_storage[3];
   1058 
   1059     if (tValues == NULL)
   1060         tValues = t_storage;
   1061 
   1062     int count = SkFindCubicMaxCurvature(src, tValues);
   1063 
   1064     if (dst)
   1065     {
   1066         if (count == 0)
   1067             memcpy(dst, src, 4 * sizeof(SkPoint));
   1068         else
   1069             SkChopCubicAt(src, dst, tValues, count);
   1070     }
   1071     return count + 1;
   1072 }
   1073 
   1074 bool SkXRayCrossesMonotonicCubic(const SkXRay& pt, const SkPoint cubic[4], bool* ambiguous) {
   1075     if (ambiguous) {
   1076         *ambiguous = false;
   1077     }
   1078 
   1079     // Find the minimum and maximum y of the extrema, which are the
   1080     // first and last points since this cubic is monotonic
   1081     SkScalar min_y = SkMinScalar(cubic[0].fY, cubic[3].fY);
   1082     SkScalar max_y = SkMaxScalar(cubic[0].fY, cubic[3].fY);
   1083 
   1084     if (pt.fY == cubic[0].fY
   1085         || pt.fY < min_y
   1086         || pt.fY > max_y) {
   1087         // The query line definitely does not cross the curve
   1088         if (ambiguous) {
   1089             *ambiguous = (pt.fY == cubic[0].fY);
   1090         }
   1091         return false;
   1092     }
   1093 
   1094     bool pt_at_extremum = (pt.fY == cubic[3].fY);
   1095 
   1096     SkScalar min_x =
   1097         SkMinScalar(
   1098             SkMinScalar(
   1099                 SkMinScalar(cubic[0].fX, cubic[1].fX),
   1100                 cubic[2].fX),
   1101             cubic[3].fX);
   1102     if (pt.fX < min_x) {
   1103         // The query line definitely crosses the curve
   1104         if (ambiguous) {
   1105             *ambiguous = pt_at_extremum;
   1106         }
   1107         return true;
   1108     }
   1109 
   1110     SkScalar max_x =
   1111         SkMaxScalar(
   1112             SkMaxScalar(
   1113                 SkMaxScalar(cubic[0].fX, cubic[1].fX),
   1114                 cubic[2].fX),
   1115             cubic[3].fX);
   1116     if (pt.fX > max_x) {
   1117         // The query line definitely does not cross the curve
   1118         return false;
   1119     }
   1120 
   1121     // Do a binary search to find the parameter value which makes y as
   1122     // close as possible to the query point. See whether the query
   1123     // line's origin is to the left of the associated x coordinate.
   1124 
   1125     // kMaxIter is chosen as the number of mantissa bits for a float,
   1126     // since there's no way we are going to get more precision by
   1127     // iterating more times than that.
   1128     const int kMaxIter = 23;
   1129     SkPoint eval;
   1130     int iter = 0;
   1131     SkScalar upper_t;
   1132     SkScalar lower_t;
   1133     // Need to invert direction of t parameter if cubic goes up
   1134     // instead of down
   1135     if (cubic[3].fY > cubic[0].fY) {
   1136         upper_t = SK_Scalar1;
   1137         lower_t = SkFloatToScalar(0);
   1138     } else {
   1139         upper_t = SkFloatToScalar(0);
   1140         lower_t = SK_Scalar1;
   1141     }
   1142     do {
   1143         SkScalar t = SkScalarAve(upper_t, lower_t);
   1144         SkEvalCubicAt(cubic, t, &eval, NULL, NULL);
   1145         if (pt.fY > eval.fY) {
   1146             lower_t = t;
   1147         } else {
   1148             upper_t = t;
   1149         }
   1150     } while (++iter < kMaxIter
   1151              && !SkScalarNearlyZero(eval.fY - pt.fY));
   1152     if (pt.fX <= eval.fX) {
   1153         if (ambiguous) {
   1154             *ambiguous = pt_at_extremum;
   1155         }
   1156         return true;
   1157     }
   1158     return false;
   1159 }
   1160 
   1161 int SkNumXRayCrossingsForCubic(const SkXRay& pt, const SkPoint cubic[4], bool* ambiguous) {
   1162     int num_crossings = 0;
   1163     SkPoint monotonic_cubics[10];
   1164     int num_monotonic_cubics = SkChopCubicAtYExtrema(cubic, monotonic_cubics);
   1165     if (ambiguous) {
   1166         *ambiguous = false;
   1167     }
   1168     bool locally_ambiguous;
   1169     if (SkXRayCrossesMonotonicCubic(pt, &monotonic_cubics[0], &locally_ambiguous))
   1170         ++num_crossings;
   1171     if (ambiguous) {
   1172         *ambiguous |= locally_ambiguous;
   1173     }
   1174     if (num_monotonic_cubics > 0)
   1175         if (SkXRayCrossesMonotonicCubic(pt, &monotonic_cubics[3], &locally_ambiguous))
   1176             ++num_crossings;
   1177     if (ambiguous) {
   1178         *ambiguous |= locally_ambiguous;
   1179     }
   1180     if (num_monotonic_cubics > 1)
   1181         if (SkXRayCrossesMonotonicCubic(pt, &monotonic_cubics[6], &locally_ambiguous))
   1182             ++num_crossings;
   1183     if (ambiguous) {
   1184         *ambiguous |= locally_ambiguous;
   1185     }
   1186     return num_crossings;
   1187 }
   1188 
   1189 ////////////////////////////////////////////////////////////////////////////////
   1190 
   1191 /*  Find t value for quadratic [a, b, c] = d.
   1192     Return 0 if there is no solution within [0, 1)
   1193 */
   1194 static SkScalar quad_solve(SkScalar a, SkScalar b, SkScalar c, SkScalar d)
   1195 {
   1196     // At^2 + Bt + C = d
   1197     SkScalar A = a - 2 * b + c;
   1198     SkScalar B = 2 * (b - a);
   1199     SkScalar C = a - d;
   1200 
   1201     SkScalar    roots[2];
   1202     int         count = SkFindUnitQuadRoots(A, B, C, roots);
   1203 
   1204     SkASSERT(count <= 1);
   1205     return count == 1 ? roots[0] : 0;
   1206 }
   1207 
   1208 /*  given a quad-curve and a point (x,y), chop the quad at that point and return
   1209     the new quad's offCurve point. Should only return false if the computed pos
   1210     is the start of the curve (i.e. root == 0)
   1211 */
   1212 static bool quad_pt2OffCurve(const SkPoint quad[3], SkScalar x, SkScalar y, SkPoint* offCurve)
   1213 {
   1214     const SkScalar* base;
   1215     SkScalar        value;
   1216 
   1217     if (SkScalarAbs(x) < SkScalarAbs(y)) {
   1218         base = &quad[0].fX;
   1219         value = x;
   1220     } else {
   1221         base = &quad[0].fY;
   1222         value = y;
   1223     }
   1224 
   1225     // note: this returns 0 if it thinks value is out of range, meaning the
   1226     // root might return something outside of [0, 1)
   1227     SkScalar t = quad_solve(base[0], base[2], base[4], value);
   1228 
   1229     if (t > 0)
   1230     {
   1231         SkPoint tmp[5];
   1232         SkChopQuadAt(quad, tmp, t);
   1233         *offCurve = tmp[1];
   1234         return true;
   1235     } else {
   1236         /*  t == 0 means either the value triggered a root outside of [0, 1)
   1237             For our purposes, we can ignore the <= 0 roots, but we want to
   1238             catch the >= 1 roots (which given our caller, will basically mean
   1239             a root of 1, give-or-take numerical instability). If we are in the
   1240             >= 1 case, return the existing offCurve point.
   1241 
   1242             The test below checks to see if we are close to the "end" of the
   1243             curve (near base[4]). Rather than specifying a tolerance, I just
   1244             check to see if value is on to the right/left of the middle point
   1245             (depending on the direction/sign of the end points).
   1246         */
   1247         if ((base[0] < base[4] && value > base[2]) ||
   1248             (base[0] > base[4] && value < base[2]))   // should root have been 1
   1249         {
   1250             *offCurve = quad[1];
   1251             return true;
   1252         }
   1253     }
   1254     return false;
   1255 }
   1256 
   1257 #ifdef SK_SCALAR_IS_FLOAT
   1258 // Due to floating point issues (i.e., 1.0f - SK_ScalarRoot2Over2 !=
   1259 // SK_ScalarRoot2Over2 - SK_ScalarTanPIOver8) a cruder root2over2
   1260 // approximation is required to make the quad circle points convex. The
   1261 // root of the problem is that with the root2over2 value in SkScalar.h
   1262 // the arcs really are ever so slightly concave. Some alternative fixes
   1263 // to this problem (besides just arbitrarily pushing out the mid-point as
   1264 // is done here) are:
   1265 //    Adjust all the points (not just the middle) to both better approximate
   1266 //             the curve and remain convex
   1267 //    Switch over to using cubics rather then quads
   1268 //    Use a different method to create the mid-point (e.g., compute
   1269 //             the two side points, average them, then move it out as needed
   1270 #ifndef SK_IGNORE_CONVEX_QUAD_OPT
   1271     #define SK_ScalarRoot2Over2_QuadCircle    0.7072f
   1272 #else
   1273     #define SK_ScalarRoot2Over2_QuadCircle    SK_ScalarRoot2Over2
   1274 #endif
   1275 
   1276 #else
   1277     #define SK_ScalarRoot2Over2_QuadCircle    SK_FixedRoot2Over2
   1278 #endif
   1279 
   1280 
   1281 static const SkPoint gQuadCirclePts[kSkBuildQuadArcStorage] = {
   1282     { SK_Scalar1,                      0                                  },
   1283     { SK_Scalar1,                      SK_ScalarTanPIOver8                },
   1284     { SK_ScalarRoot2Over2_QuadCircle,  SK_ScalarRoot2Over2_QuadCircle     },
   1285     { SK_ScalarTanPIOver8,             SK_Scalar1                         },
   1286 
   1287     { 0,                               SK_Scalar1                         },
   1288     { -SK_ScalarTanPIOver8,            SK_Scalar1                         },
   1289     { -SK_ScalarRoot2Over2_QuadCircle, SK_ScalarRoot2Over2_QuadCircle     },
   1290     { -SK_Scalar1,                     SK_ScalarTanPIOver8                },
   1291 
   1292     { -SK_Scalar1,                     0                                  },
   1293     { -SK_Scalar1,                     -SK_ScalarTanPIOver8               },
   1294     { -SK_ScalarRoot2Over2_QuadCircle, -SK_ScalarRoot2Over2_QuadCircle    },
   1295     { -SK_ScalarTanPIOver8,            -SK_Scalar1                        },
   1296 
   1297     { 0,                               -SK_Scalar1                        },
   1298     { SK_ScalarTanPIOver8,             -SK_Scalar1                        },
   1299     { SK_ScalarRoot2Over2_QuadCircle,  -SK_ScalarRoot2Over2_QuadCircle    },
   1300     { SK_Scalar1,                      -SK_ScalarTanPIOver8               },
   1301 
   1302     { SK_Scalar1,                      0                                  }
   1303 };
   1304 
   1305 int SkBuildQuadArc(const SkVector& uStart, const SkVector& uStop,
   1306                    SkRotationDirection dir, const SkMatrix* userMatrix,
   1307                    SkPoint quadPoints[])
   1308 {
   1309     // rotate by x,y so that uStart is (1.0)
   1310     SkScalar x = SkPoint::DotProduct(uStart, uStop);
   1311     SkScalar y = SkPoint::CrossProduct(uStart, uStop);
   1312 
   1313     SkScalar absX = SkScalarAbs(x);
   1314     SkScalar absY = SkScalarAbs(y);
   1315 
   1316     int pointCount;
   1317 
   1318     // check for (effectively) coincident vectors
   1319     // this can happen if our angle is nearly 0 or nearly 180 (y == 0)
   1320     // ... we use the dot-prod to distinguish between 0 and 180 (x > 0)
   1321     if (absY <= SK_ScalarNearlyZero && x > 0 &&
   1322         ((y >= 0 && kCW_SkRotationDirection == dir) ||
   1323          (y <= 0 && kCCW_SkRotationDirection == dir))) {
   1324 
   1325         // just return the start-point
   1326         quadPoints[0].set(SK_Scalar1, 0);
   1327         pointCount = 1;
   1328     } else {
   1329         if (dir == kCCW_SkRotationDirection)
   1330             y = -y;
   1331 
   1332         // what octant (quadratic curve) is [xy] in?
   1333         int oct = 0;
   1334         bool sameSign = true;
   1335 
   1336         if (0 == y)
   1337         {
   1338             oct = 4;        // 180
   1339             SkASSERT(SkScalarAbs(x + SK_Scalar1) <= SK_ScalarNearlyZero);
   1340         }
   1341         else if (0 == x)
   1342         {
   1343             SkASSERT(absY - SK_Scalar1 <= SK_ScalarNearlyZero);
   1344             if (y > 0)
   1345                 oct = 2;    // 90
   1346             else
   1347                 oct = 6;    // 270
   1348         }
   1349         else
   1350         {
   1351             if (y < 0)
   1352                 oct += 4;
   1353             if ((x < 0) != (y < 0))
   1354             {
   1355                 oct += 2;
   1356                 sameSign = false;
   1357             }
   1358             if ((absX < absY) == sameSign)
   1359                 oct += 1;
   1360         }
   1361 
   1362         int wholeCount = oct << 1;
   1363         memcpy(quadPoints, gQuadCirclePts, (wholeCount + 1) * sizeof(SkPoint));
   1364 
   1365         const SkPoint* arc = &gQuadCirclePts[wholeCount];
   1366         if (quad_pt2OffCurve(arc, x, y, &quadPoints[wholeCount + 1]))
   1367         {
   1368             quadPoints[wholeCount + 2].set(x, y);
   1369             wholeCount += 2;
   1370         }
   1371         pointCount = wholeCount + 1;
   1372     }
   1373 
   1374     // now handle counter-clockwise and the initial unitStart rotation
   1375     SkMatrix    matrix;
   1376     matrix.setSinCos(uStart.fY, uStart.fX);
   1377     if (dir == kCCW_SkRotationDirection) {
   1378         matrix.preScale(SK_Scalar1, -SK_Scalar1);
   1379     }
   1380     if (userMatrix) {
   1381         matrix.postConcat(*userMatrix);
   1382     }
   1383     matrix.mapPoints(quadPoints, pointCount);
   1384     return pointCount;
   1385 }
   1386