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