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