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
14bool 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
100static bool is_unit_interval(SkScalar x)
101{
102    return x > 0 && x < SK_Scalar1;
103}
104
105static 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*/
135int 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*/
179static 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
195static 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
212static 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
220static 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
227void 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
239void 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
256static 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
268void 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
276void 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*/
295int 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
307static 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 */
315int 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 */
359int 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//
397int 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
449void 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
463static 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
471void 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
481static 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*/
508static 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
515static 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
524static 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
532void 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*/
553int 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
568static 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
586void 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
617void 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
662void 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
685static 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*/
698int 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
714int 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*/
741int 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
774int 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
789template <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
805static 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 */
841static 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
860static 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*/
907static 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*/
998static 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*/
1030int 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
1056int 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
1075bool 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
1162int 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*/
1195static 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*/
1213static 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
1258static 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
1282int 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