1/*
2 * Copyright 2012 Google Inc.
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#include "SkGeometry.h"
8#include "SkLineParameters.h"
9#include "SkPathOpsConic.h"
10#include "SkPathOpsCubic.h"
11#include "SkPathOpsCurve.h"
12#include "SkPathOpsLine.h"
13#include "SkPathOpsQuad.h"
14#include "SkPathOpsRect.h"
15#include "SkTSort.h"
16
17const int SkDCubic::gPrecisionUnit = 256;  // FIXME: test different values in test framework
18
19void SkDCubic::align(int endIndex, int ctrlIndex, SkDPoint* dstPt) const {
20    if (fPts[endIndex].fX == fPts[ctrlIndex].fX) {
21        dstPt->fX = fPts[endIndex].fX;
22    }
23    if (fPts[endIndex].fY == fPts[ctrlIndex].fY) {
24        dstPt->fY = fPts[endIndex].fY;
25    }
26}
27
28// give up when changing t no longer moves point
29// also, copy point rather than recompute it when it does change
30double SkDCubic::binarySearch(double min, double max, double axisIntercept,
31        SearchAxis xAxis) const {
32    double t = (min + max) / 2;
33    double step = (t - min) / 2;
34    SkDPoint cubicAtT = ptAtT(t);
35    double calcPos = (&cubicAtT.fX)[xAxis];
36    double calcDist = calcPos - axisIntercept;
37    do {
38        double priorT = t - step;
39        SkASSERT(priorT >= min);
40        SkDPoint lessPt = ptAtT(priorT);
41        if (approximately_equal_half(lessPt.fX, cubicAtT.fX)
42                && approximately_equal_half(lessPt.fY, cubicAtT.fY)) {
43            return -1;  // binary search found no point at this axis intercept
44        }
45        double lessDist = (&lessPt.fX)[xAxis] - axisIntercept;
46#if DEBUG_CUBIC_BINARY_SEARCH
47        SkDebugf("t=%1.9g calc=%1.9g dist=%1.9g step=%1.9g less=%1.9g\n", t, calcPos, calcDist,
48                step, lessDist);
49#endif
50        double lastStep = step;
51        step /= 2;
52        if (calcDist > 0 ? calcDist > lessDist : calcDist < lessDist) {
53            t = priorT;
54        } else {
55            double nextT = t + lastStep;
56            if (nextT > max) {
57                return -1;
58            }
59            SkDPoint morePt = ptAtT(nextT);
60            if (approximately_equal_half(morePt.fX, cubicAtT.fX)
61                    && approximately_equal_half(morePt.fY, cubicAtT.fY)) {
62                return -1;  // binary search found no point at this axis intercept
63            }
64            double moreDist = (&morePt.fX)[xAxis] - axisIntercept;
65            if (calcDist > 0 ? calcDist <= moreDist : calcDist >= moreDist) {
66                continue;
67            }
68            t = nextT;
69        }
70        SkDPoint testAtT = ptAtT(t);
71        cubicAtT = testAtT;
72        calcPos = (&cubicAtT.fX)[xAxis];
73        calcDist = calcPos - axisIntercept;
74    } while (!approximately_equal(calcPos, axisIntercept));
75    return t;
76}
77
78// FIXME: cache keep the bounds and/or precision with the caller?
79double SkDCubic::calcPrecision() const {
80    SkDRect dRect;
81    dRect.setBounds(*this);  // OPTIMIZATION: just use setRawBounds ?
82    double width = dRect.fRight - dRect.fLeft;
83    double height = dRect.fBottom - dRect.fTop;
84    return (width > height ? width : height) / gPrecisionUnit;
85}
86
87
88/* classic one t subdivision */
89static void interp_cubic_coords(const double* src, double* dst, double t) {
90    double ab = SkDInterp(src[0], src[2], t);
91    double bc = SkDInterp(src[2], src[4], t);
92    double cd = SkDInterp(src[4], src[6], t);
93    double abc = SkDInterp(ab, bc, t);
94    double bcd = SkDInterp(bc, cd, t);
95    double abcd = SkDInterp(abc, bcd, t);
96
97    dst[0] = src[0];
98    dst[2] = ab;
99    dst[4] = abc;
100    dst[6] = abcd;
101    dst[8] = bcd;
102    dst[10] = cd;
103    dst[12] = src[6];
104}
105
106SkDCubicPair SkDCubic::chopAt(double t) const {
107    SkDCubicPair dst;
108    if (t == 0.5) {
109        dst.pts[0] = fPts[0];
110        dst.pts[1].fX = (fPts[0].fX + fPts[1].fX) / 2;
111        dst.pts[1].fY = (fPts[0].fY + fPts[1].fY) / 2;
112        dst.pts[2].fX = (fPts[0].fX + 2 * fPts[1].fX + fPts[2].fX) / 4;
113        dst.pts[2].fY = (fPts[0].fY + 2 * fPts[1].fY + fPts[2].fY) / 4;
114        dst.pts[3].fX = (fPts[0].fX + 3 * (fPts[1].fX + fPts[2].fX) + fPts[3].fX) / 8;
115        dst.pts[3].fY = (fPts[0].fY + 3 * (fPts[1].fY + fPts[2].fY) + fPts[3].fY) / 8;
116        dst.pts[4].fX = (fPts[1].fX + 2 * fPts[2].fX + fPts[3].fX) / 4;
117        dst.pts[4].fY = (fPts[1].fY + 2 * fPts[2].fY + fPts[3].fY) / 4;
118        dst.pts[5].fX = (fPts[2].fX + fPts[3].fX) / 2;
119        dst.pts[5].fY = (fPts[2].fY + fPts[3].fY) / 2;
120        dst.pts[6] = fPts[3];
121        return dst;
122    }
123    interp_cubic_coords(&fPts[0].fX, &dst.pts[0].fX, t);
124    interp_cubic_coords(&fPts[0].fY, &dst.pts[0].fY, t);
125    return dst;
126}
127
128void SkDCubic::Coefficients(const double* src, double* A, double* B, double* C, double* D) {
129    *A = src[6];  // d
130    *B = src[4] * 3;  // 3*c
131    *C = src[2] * 3;  // 3*b
132    *D = src[0];  // a
133    *A -= *D - *C + *B;     // A =   -a + 3*b - 3*c + d
134    *B += 3 * *D - 2 * *C;  // B =  3*a - 6*b + 3*c
135    *C -= 3 * *D;           // C = -3*a + 3*b
136}
137
138bool SkDCubic::endsAreExtremaInXOrY() const {
139    return (between(fPts[0].fX, fPts[1].fX, fPts[3].fX)
140            && between(fPts[0].fX, fPts[2].fX, fPts[3].fX))
141            || (between(fPts[0].fY, fPts[1].fY, fPts[3].fY)
142            && between(fPts[0].fY, fPts[2].fY, fPts[3].fY));
143}
144
145// Do a quick reject by rotating all points relative to a line formed by
146// a pair of one cubic's points. If the 2nd cubic's points
147// are on the line or on the opposite side from the 1st cubic's 'odd man', the
148// curves at most intersect at the endpoints.
149/* if returning true, check contains true if cubic's hull collapsed, making the cubic linear
150   if returning false, check contains true if the the cubic pair have only the end point in common
151*/
152bool SkDCubic::hullIntersects(const SkDPoint* pts, int ptCount, bool* isLinear) const {
153    bool linear = true;
154    char hullOrder[4];
155    int hullCount = convexHull(hullOrder);
156    int end1 = hullOrder[0];
157    int hullIndex = 0;
158    const SkDPoint* endPt[2];
159    endPt[0] = &fPts[end1];
160    do {
161        hullIndex = (hullIndex + 1) % hullCount;
162        int end2 = hullOrder[hullIndex];
163        endPt[1] = &fPts[end2];
164        double origX = endPt[0]->fX;
165        double origY = endPt[0]->fY;
166        double adj = endPt[1]->fX - origX;
167        double opp = endPt[1]->fY - origY;
168        int oddManMask = other_two(end1, end2);
169        int oddMan = end1 ^ oddManMask;
170        double sign = (fPts[oddMan].fY - origY) * adj - (fPts[oddMan].fX - origX) * opp;
171        int oddMan2 = end2 ^ oddManMask;
172        double sign2 = (fPts[oddMan2].fY - origY) * adj - (fPts[oddMan2].fX - origX) * opp;
173        if (sign * sign2 < 0) {
174            continue;
175        }
176        if (approximately_zero(sign)) {
177            sign = sign2;
178            if (approximately_zero(sign)) {
179                continue;
180            }
181        }
182        linear = false;
183        bool foundOutlier = false;
184        for (int n = 0; n < ptCount; ++n) {
185            double test = (pts[n].fY - origY) * adj - (pts[n].fX - origX) * opp;
186            if (test * sign > 0 && !precisely_zero(test)) {
187                foundOutlier = true;
188                break;
189            }
190        }
191        if (!foundOutlier) {
192            return false;
193        }
194        endPt[0] = endPt[1];
195        end1 = end2;
196    } while (hullIndex);
197    *isLinear = linear;
198    return true;
199}
200
201bool SkDCubic::hullIntersects(const SkDCubic& c2, bool* isLinear) const {
202    return hullIntersects(c2.fPts, c2.kPointCount, isLinear);
203}
204
205bool SkDCubic::hullIntersects(const SkDQuad& quad, bool* isLinear) const {
206    return hullIntersects(quad.fPts, quad.kPointCount, isLinear);
207}
208
209bool SkDCubic::hullIntersects(const SkDConic& conic, bool* isLinear) const {
210
211    return hullIntersects(conic.fPts, isLinear);
212}
213
214bool SkDCubic::isLinear(int startIndex, int endIndex) const {
215    SkLineParameters lineParameters;
216    lineParameters.cubicEndPoints(*this, startIndex, endIndex);
217    // FIXME: maybe it's possible to avoid this and compare non-normalized
218    lineParameters.normalize();
219    double tiniest = SkTMin(SkTMin(SkTMin(SkTMin(SkTMin(SkTMin(SkTMin(fPts[0].fX, fPts[0].fY),
220            fPts[1].fX), fPts[1].fY), fPts[2].fX), fPts[2].fY), fPts[3].fX), fPts[3].fY);
221    double largest = SkTMax(SkTMax(SkTMax(SkTMax(SkTMax(SkTMax(SkTMax(fPts[0].fX, fPts[0].fY),
222            fPts[1].fX), fPts[1].fY), fPts[2].fX), fPts[2].fY), fPts[3].fX), fPts[3].fY);
223    largest = SkTMax(largest, -tiniest);
224    double distance = lineParameters.controlPtDistance(*this, 1);
225    if (!approximately_zero_when_compared_to(distance, largest)) {
226        return false;
227    }
228    distance = lineParameters.controlPtDistance(*this, 2);
229    return approximately_zero_when_compared_to(distance, largest);
230}
231
232bool SkDCubic::ComplexBreak(const SkPoint pointsPtr[4], SkScalar* t, CubicType* resultType) {
233    SkScalar d[3];
234    SkCubicType cubicType = SkClassifyCubic(pointsPtr, d);
235    if (cubicType == kLoop_SkCubicType) {
236        // crib code from gpu path utils that finds t values where loop self-intersects
237        // use it to find mid of t values which should be a friendly place to chop
238        SkScalar tempSqrt = SkScalarSqrt(4.f * d[0] * d[2] - 3.f * d[1] * d[1]);
239        SkScalar ls = d[1] - tempSqrt;
240        SkScalar lt = 2.f * d[0];
241        SkScalar ms = d[1] + tempSqrt;
242        SkScalar mt = 2.f * d[0];
243        if (between(0, ls, lt) || between(0, ms, mt)) {
244            ls = ls / lt;
245            ms = ms / mt;
246            SkScalar smaller = SkTMax(0.f, SkTMin(ls, ms));
247            SkScalar larger = SkTMin(1.f, SkTMax(ls, ms));
248            *t = (smaller + larger) / 2;
249            *resultType = kSplitAtLoop_SkDCubicType;
250            return *t > 0 && *t < 1;
251        }
252    } else if (kSerpentine_SkCubicType == cubicType || kCusp_SkCubicType == cubicType) {
253        SkDCubic cubic;
254        cubic.set(pointsPtr);
255        double inflectionTs[2];
256        int infTCount = cubic.findInflections(inflectionTs);
257        if (infTCount == 2) {
258            double maxCurvature[3];
259            int roots = cubic.findMaxCurvature(maxCurvature);
260#if DEBUG_CUBIC_SPLIT
261            SkDebugf("%s\n", __FUNCTION__);
262            cubic.dump();
263            for (int index = 0; index < infTCount; ++index) {
264                SkDebugf("inflectionsTs[%d]=%1.9g ", index, inflectionTs[index]);
265                SkDPoint pt = cubic.ptAtT(inflectionTs[index]);
266                SkDVector dPt = cubic.dxdyAtT(inflectionTs[index]);
267                SkDLine perp = {{pt - dPt, pt + dPt}};
268                perp.dump();
269            }
270            for (int index = 0; index < roots; ++index) {
271                SkDebugf("maxCurvature[%d]=%1.9g ", index, maxCurvature[index]);
272                SkDPoint pt = cubic.ptAtT(maxCurvature[index]);
273                SkDVector dPt = cubic.dxdyAtT(maxCurvature[index]);
274                SkDLine perp = {{pt - dPt, pt + dPt}};
275                perp.dump();
276            }
277#endif
278            for (int index = 0; index < roots; ++index) {
279                if (between(inflectionTs[0], maxCurvature[index], inflectionTs[1])) {
280                    *t = maxCurvature[index];
281                    *resultType = kSplitAtMaxCurvature_SkDCubicType;
282                    return true;
283                }
284            }
285        } else if (infTCount == 1) {
286            *t = inflectionTs[0];
287            *resultType = kSplitAtInflection_SkDCubicType;
288            return *t > 0 && *t < 1;
289        }
290    }
291    return false;
292}
293
294bool SkDCubic::monotonicInX() const {
295    return precisely_between(fPts[0].fX, fPts[1].fX, fPts[3].fX)
296            && precisely_between(fPts[0].fX, fPts[2].fX, fPts[3].fX);
297}
298
299bool SkDCubic::monotonicInY() const {
300    return precisely_between(fPts[0].fY, fPts[1].fY, fPts[3].fY)
301            && precisely_between(fPts[0].fY, fPts[2].fY, fPts[3].fY);
302}
303
304void SkDCubic::otherPts(int index, const SkDPoint* o1Pts[kPointCount - 1]) const {
305    int offset = (int) !SkToBool(index);
306    o1Pts[0] = &fPts[offset];
307    o1Pts[1] = &fPts[++offset];
308    o1Pts[2] = &fPts[++offset];
309}
310
311int SkDCubic::searchRoots(double extremeTs[6], int extrema, double axisIntercept,
312        SearchAxis xAxis, double* validRoots) const {
313    extrema += findInflections(&extremeTs[extrema]);
314    extremeTs[extrema++] = 0;
315    extremeTs[extrema] = 1;
316    SkTQSort(extremeTs, extremeTs + extrema);
317    int validCount = 0;
318    for (int index = 0; index < extrema; ) {
319        double min = extremeTs[index];
320        double max = extremeTs[++index];
321        if (min == max) {
322            continue;
323        }
324        double newT = binarySearch(min, max, axisIntercept, xAxis);
325        if (newT >= 0) {
326            validRoots[validCount++] = newT;
327        }
328    }
329    return validCount;
330}
331
332// cubic roots
333
334static const double PI = 3.141592653589793;
335
336// from SkGeometry.cpp (and Numeric Solutions, 5.6)
337int SkDCubic::RootsValidT(double A, double B, double C, double D, double t[3]) {
338    double s[3];
339    int realRoots = RootsReal(A, B, C, D, s);
340    int foundRoots = SkDQuad::AddValidTs(s, realRoots, t);
341    for (int index = 0; index < realRoots; ++index) {
342        double tValue = s[index];
343        if (!approximately_one_or_less(tValue) && between(1, tValue, 1.00005)) {
344            for (int idx2 = 0; idx2 < foundRoots; ++idx2) {
345                if (approximately_equal(t[idx2], 1)) {
346                    goto nextRoot;
347                }
348            }
349            SkASSERT(foundRoots < 3);
350            t[foundRoots++] = 1;
351        } else if (!approximately_zero_or_more(tValue) && between(-0.00005, tValue, 0)) {
352            for (int idx2 = 0; idx2 < foundRoots; ++idx2) {
353                if (approximately_equal(t[idx2], 0)) {
354                    goto nextRoot;
355                }
356            }
357            SkASSERT(foundRoots < 3);
358            t[foundRoots++] = 0;
359        }
360nextRoot:
361        ;
362    }
363    return foundRoots;
364}
365
366int SkDCubic::RootsReal(double A, double B, double C, double D, double s[3]) {
367#ifdef SK_DEBUG
368    // create a string mathematica understands
369    // GDB set print repe 15 # if repeated digits is a bother
370    //     set print elements 400 # if line doesn't fit
371    char str[1024];
372    sk_bzero(str, sizeof(str));
373    SK_SNPRINTF(str, sizeof(str), "Solve[%1.19g x^3 + %1.19g x^2 + %1.19g x + %1.19g == 0, x]",
374            A, B, C, D);
375    SkPathOpsDebug::MathematicaIze(str, sizeof(str));
376#if ONE_OFF_DEBUG && ONE_OFF_DEBUG_MATHEMATICA
377    SkDebugf("%s\n", str);
378#endif
379#endif
380    if (approximately_zero(A)
381            && approximately_zero_when_compared_to(A, B)
382            && approximately_zero_when_compared_to(A, C)
383            && approximately_zero_when_compared_to(A, D)) {  // we're just a quadratic
384        return SkDQuad::RootsReal(B, C, D, s);
385    }
386    if (approximately_zero_when_compared_to(D, A)
387            && approximately_zero_when_compared_to(D, B)
388            && approximately_zero_when_compared_to(D, C)) {  // 0 is one root
389        int num = SkDQuad::RootsReal(A, B, C, s);
390        for (int i = 0; i < num; ++i) {
391            if (approximately_zero(s[i])) {
392                return num;
393            }
394        }
395        s[num++] = 0;
396        return num;
397    }
398    if (approximately_zero(A + B + C + D)) {  // 1 is one root
399        int num = SkDQuad::RootsReal(A, A + B, -D, s);
400        for (int i = 0; i < num; ++i) {
401            if (AlmostDequalUlps(s[i], 1)) {
402                return num;
403            }
404        }
405        s[num++] = 1;
406        return num;
407    }
408    double a, b, c;
409    {
410        double invA = 1 / A;
411        a = B * invA;
412        b = C * invA;
413        c = D * invA;
414    }
415    double a2 = a * a;
416    double Q = (a2 - b * 3) / 9;
417    double R = (2 * a2 * a - 9 * a * b + 27 * c) / 54;
418    double R2 = R * R;
419    double Q3 = Q * Q * Q;
420    double R2MinusQ3 = R2 - Q3;
421    double adiv3 = a / 3;
422    double r;
423    double* roots = s;
424    if (R2MinusQ3 < 0) {   // we have 3 real roots
425        double theta = acos(R / sqrt(Q3));
426        double neg2RootQ = -2 * sqrt(Q);
427
428        r = neg2RootQ * cos(theta / 3) - adiv3;
429        *roots++ = r;
430
431        r = neg2RootQ * cos((theta + 2 * PI) / 3) - adiv3;
432        if (!AlmostDequalUlps(s[0], r)) {
433            *roots++ = r;
434        }
435        r = neg2RootQ * cos((theta - 2 * PI) / 3) - adiv3;
436        if (!AlmostDequalUlps(s[0], r) && (roots - s == 1 || !AlmostDequalUlps(s[1], r))) {
437            *roots++ = r;
438        }
439    } else {  // we have 1 real root
440        double sqrtR2MinusQ3 = sqrt(R2MinusQ3);
441        double A = fabs(R) + sqrtR2MinusQ3;
442        A = SkDCubeRoot(A);
443        if (R > 0) {
444            A = -A;
445        }
446        if (A != 0) {
447            A += Q / A;
448        }
449        r = A - adiv3;
450        *roots++ = r;
451        if (AlmostDequalUlps((double) R2, (double) Q3)) {
452            r = -A / 2 - adiv3;
453            if (!AlmostDequalUlps(s[0], r)) {
454                *roots++ = r;
455            }
456        }
457    }
458    return static_cast<int>(roots - s);
459}
460
461// from http://www.cs.sunysb.edu/~qin/courses/geometry/4.pdf
462// c(t)  = a(1-t)^3 + 3bt(1-t)^2 + 3c(1-t)t^2 + dt^3
463// c'(t) = -3a(1-t)^2 + 3b((1-t)^2 - 2t(1-t)) + 3c(2t(1-t) - t^2) + 3dt^2
464//       = 3(b-a)(1-t)^2 + 6(c-b)t(1-t) + 3(d-c)t^2
465static double derivative_at_t(const double* src, double t) {
466    double one_t = 1 - t;
467    double a = src[0];
468    double b = src[2];
469    double c = src[4];
470    double d = src[6];
471    return 3 * ((b - a) * one_t * one_t + 2 * (c - b) * t * one_t + (d - c) * t * t);
472}
473
474// OPTIMIZE? compute t^2, t(1-t), and (1-t)^2 and pass them to another version of derivative at t?
475SkDVector SkDCubic::dxdyAtT(double t) const {
476    SkDVector result = { derivative_at_t(&fPts[0].fX, t), derivative_at_t(&fPts[0].fY, t) };
477    return result;
478}
479
480// OPTIMIZE? share code with formulate_F1DotF2
481int SkDCubic::findInflections(double tValues[]) const {
482    double Ax = fPts[1].fX - fPts[0].fX;
483    double Ay = fPts[1].fY - fPts[0].fY;
484    double Bx = fPts[2].fX - 2 * fPts[1].fX + fPts[0].fX;
485    double By = fPts[2].fY - 2 * fPts[1].fY + fPts[0].fY;
486    double Cx = fPts[3].fX + 3 * (fPts[1].fX - fPts[2].fX) - fPts[0].fX;
487    double Cy = fPts[3].fY + 3 * (fPts[1].fY - fPts[2].fY) - fPts[0].fY;
488    return SkDQuad::RootsValidT(Bx * Cy - By * Cx, Ax * Cy - Ay * Cx, Ax * By - Ay * Bx, tValues);
489}
490
491static void formulate_F1DotF2(const double src[], double coeff[4]) {
492    double a = src[2] - src[0];
493    double b = src[4] - 2 * src[2] + src[0];
494    double c = src[6] + 3 * (src[2] - src[4]) - src[0];
495    coeff[0] = c * c;
496    coeff[1] = 3 * b * c;
497    coeff[2] = 2 * b * b + c * a;
498    coeff[3] = a * b;
499}
500
501/** SkDCubic'(t) = At^2 + Bt + C, where
502    A = 3(-a + 3(b - c) + d)
503    B = 6(a - 2b + c)
504    C = 3(b - a)
505    Solve for t, keeping only those that fit between 0 < t < 1
506*/
507int SkDCubic::FindExtrema(const double src[], double tValues[2]) {
508    // we divide A,B,C by 3 to simplify
509    double a = src[0];
510    double b = src[2];
511    double c = src[4];
512    double d = src[6];
513    double A = d - a + 3 * (b - c);
514    double B = 2 * (a - b - b + c);
515    double C = b - a;
516
517    return SkDQuad::RootsValidT(A, B, C, tValues);
518}
519
520/*  from SkGeometry.cpp
521    Looking for F' dot F'' == 0
522
523    A = b - a
524    B = c - 2b + a
525    C = d - 3c + 3b - a
526
527    F' = 3Ct^2 + 6Bt + 3A
528    F'' = 6Ct + 6B
529
530    F' dot F'' -> CCt^3 + 3BCt^2 + (2BB + CA)t + AB
531*/
532int SkDCubic::findMaxCurvature(double tValues[]) const {
533    double coeffX[4], coeffY[4];
534    int i;
535    formulate_F1DotF2(&fPts[0].fX, coeffX);
536    formulate_F1DotF2(&fPts[0].fY, coeffY);
537    for (i = 0; i < 4; i++) {
538        coeffX[i] = coeffX[i] + coeffY[i];
539    }
540    return RootsValidT(coeffX[0], coeffX[1], coeffX[2], coeffX[3], tValues);
541}
542
543SkDPoint SkDCubic::ptAtT(double t) const {
544    if (0 == t) {
545        return fPts[0];
546    }
547    if (1 == t) {
548        return fPts[3];
549    }
550    double one_t = 1 - t;
551    double one_t2 = one_t * one_t;
552    double a = one_t2 * one_t;
553    double b = 3 * one_t2 * t;
554    double t2 = t * t;
555    double c = 3 * one_t * t2;
556    double d = t2 * t;
557    SkDPoint result = {a * fPts[0].fX + b * fPts[1].fX + c * fPts[2].fX + d * fPts[3].fX,
558            a * fPts[0].fY + b * fPts[1].fY + c * fPts[2].fY + d * fPts[3].fY};
559    return result;
560}
561
562/*
563 Given a cubic c, t1, and t2, find a small cubic segment.
564
565 The new cubic is defined as points A, B, C, and D, where
566 s1 = 1 - t1
567 s2 = 1 - t2
568 A = c[0]*s1*s1*s1 + 3*c[1]*s1*s1*t1 + 3*c[2]*s1*t1*t1 + c[3]*t1*t1*t1
569 D = c[0]*s2*s2*s2 + 3*c[1]*s2*s2*t2 + 3*c[2]*s2*t2*t2 + c[3]*t2*t2*t2
570
571 We don't have B or C. So We define two equations to isolate them.
572 First, compute two reference T values 1/3 and 2/3 from t1 to t2:
573
574 c(at (2*t1 + t2)/3) == E
575 c(at (t1 + 2*t2)/3) == F
576
577 Next, compute where those values must be if we know the values of B and C:
578
579 _12   =  A*2/3 + B*1/3
580 12_   =  A*1/3 + B*2/3
581 _23   =  B*2/3 + C*1/3
582 23_   =  B*1/3 + C*2/3
583 _34   =  C*2/3 + D*1/3
584 34_   =  C*1/3 + D*2/3
585 _123  = (A*2/3 + B*1/3)*2/3 + (B*2/3 + C*1/3)*1/3 = A*4/9 + B*4/9 + C*1/9
586 123_  = (A*1/3 + B*2/3)*1/3 + (B*1/3 + C*2/3)*2/3 = A*1/9 + B*4/9 + C*4/9
587 _234  = (B*2/3 + C*1/3)*2/3 + (C*2/3 + D*1/3)*1/3 = B*4/9 + C*4/9 + D*1/9
588 234_  = (B*1/3 + C*2/3)*1/3 + (C*1/3 + D*2/3)*2/3 = B*1/9 + C*4/9 + D*4/9
589 _1234 = (A*4/9 + B*4/9 + C*1/9)*2/3 + (B*4/9 + C*4/9 + D*1/9)*1/3
590       =  A*8/27 + B*12/27 + C*6/27 + D*1/27
591       =  E
592 1234_ = (A*1/9 + B*4/9 + C*4/9)*1/3 + (B*1/9 + C*4/9 + D*4/9)*2/3
593       =  A*1/27 + B*6/27 + C*12/27 + D*8/27
594       =  F
595 E*27  =  A*8    + B*12   + C*6     + D
596 F*27  =  A      + B*6    + C*12    + D*8
597
598Group the known values on one side:
599
600 M       = E*27 - A*8 - D     = B*12 + C* 6
601 N       = F*27 - A   - D*8   = B* 6 + C*12
602 M*2 - N = B*18
603 N*2 - M = C*18
604 B       = (M*2 - N)/18
605 C       = (N*2 - M)/18
606 */
607
608static double interp_cubic_coords(const double* src, double t) {
609    double ab = SkDInterp(src[0], src[2], t);
610    double bc = SkDInterp(src[2], src[4], t);
611    double cd = SkDInterp(src[4], src[6], t);
612    double abc = SkDInterp(ab, bc, t);
613    double bcd = SkDInterp(bc, cd, t);
614    double abcd = SkDInterp(abc, bcd, t);
615    return abcd;
616}
617
618SkDCubic SkDCubic::subDivide(double t1, double t2) const {
619    if (t1 == 0 || t2 == 1) {
620        if (t1 == 0 && t2 == 1) {
621            return *this;
622        }
623        SkDCubicPair pair = chopAt(t1 == 0 ? t2 : t1);
624        SkDCubic dst = t1 == 0 ? pair.first() : pair.second();
625        return dst;
626    }
627    SkDCubic dst;
628    double ax = dst[0].fX = interp_cubic_coords(&fPts[0].fX, t1);
629    double ay = dst[0].fY = interp_cubic_coords(&fPts[0].fY, t1);
630    double ex = interp_cubic_coords(&fPts[0].fX, (t1*2+t2)/3);
631    double ey = interp_cubic_coords(&fPts[0].fY, (t1*2+t2)/3);
632    double fx = interp_cubic_coords(&fPts[0].fX, (t1+t2*2)/3);
633    double fy = interp_cubic_coords(&fPts[0].fY, (t1+t2*2)/3);
634    double dx = dst[3].fX = interp_cubic_coords(&fPts[0].fX, t2);
635    double dy = dst[3].fY = interp_cubic_coords(&fPts[0].fY, t2);
636    double mx = ex * 27 - ax * 8 - dx;
637    double my = ey * 27 - ay * 8 - dy;
638    double nx = fx * 27 - ax - dx * 8;
639    double ny = fy * 27 - ay - dy * 8;
640    /* bx = */ dst[1].fX = (mx * 2 - nx) / 18;
641    /* by = */ dst[1].fY = (my * 2 - ny) / 18;
642    /* cx = */ dst[2].fX = (nx * 2 - mx) / 18;
643    /* cy = */ dst[2].fY = (ny * 2 - my) / 18;
644    // FIXME: call align() ?
645    return dst;
646}
647
648void SkDCubic::subDivide(const SkDPoint& a, const SkDPoint& d,
649                         double t1, double t2, SkDPoint dst[2]) const {
650    SkASSERT(t1 != t2);
651    // this approach assumes that the control points computed directly are accurate enough
652    SkDCubic sub = subDivide(t1, t2);
653    dst[0] = sub[1] + (a - sub[0]);
654    dst[1] = sub[2] + (d - sub[3]);
655    if (t1 == 0 || t2 == 0) {
656        align(0, 1, t1 == 0 ? &dst[0] : &dst[1]);
657    }
658    if (t1 == 1 || t2 == 1) {
659        align(3, 2, t1 == 1 ? &dst[0] : &dst[1]);
660    }
661    if (AlmostBequalUlps(dst[0].fX, a.fX)) {
662        dst[0].fX = a.fX;
663    }
664    if (AlmostBequalUlps(dst[0].fY, a.fY)) {
665        dst[0].fY = a.fY;
666    }
667    if (AlmostBequalUlps(dst[1].fX, d.fX)) {
668        dst[1].fX = d.fX;
669    }
670    if (AlmostBequalUlps(dst[1].fY, d.fY)) {
671        dst[1].fY = d.fY;
672    }
673}
674
675double SkDCubic::top(const SkDCubic& dCurve, double startT, double endT, SkDPoint*topPt) const {
676    double extremeTs[2];
677    double topT = -1;
678    int roots = SkDCubic::FindExtrema(&fPts[0].fY, extremeTs);
679    for (int index = 0; index < roots; ++index) {
680        double t = startT + (endT - startT) * extremeTs[index];
681        SkDPoint mid = dCurve.ptAtT(t);
682        if (topPt->fY > mid.fY || (topPt->fY == mid.fY && topPt->fX > mid.fX)) {
683            topT = t;
684            *topPt = mid;
685        }
686    }
687    return topT;
688}
689