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