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        SkOPASSERT(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// get the rough scale of the cubic; used to determine if curvature is extreme
79double SkDCubic::calcPrecision() const {
80    return ((fPts[1] - fPts[0]).length()
81            + (fPts[2] - fPts[1]).length()
82            + (fPts[3] - fPts[2]).length()) / gPrecisionUnit;
83}
84
85/* classic one t subdivision */
86static void interp_cubic_coords(const double* src, double* dst, double t) {
87    double ab = SkDInterp(src[0], src[2], t);
88    double bc = SkDInterp(src[2], src[4], t);
89    double cd = SkDInterp(src[4], src[6], t);
90    double abc = SkDInterp(ab, bc, t);
91    double bcd = SkDInterp(bc, cd, t);
92    double abcd = SkDInterp(abc, bcd, t);
93
94    dst[0] = src[0];
95    dst[2] = ab;
96    dst[4] = abc;
97    dst[6] = abcd;
98    dst[8] = bcd;
99    dst[10] = cd;
100    dst[12] = src[6];
101}
102
103SkDCubicPair SkDCubic::chopAt(double t) const {
104    SkDCubicPair dst;
105    if (t == 0.5) {
106        dst.pts[0] = fPts[0];
107        dst.pts[1].fX = (fPts[0].fX + fPts[1].fX) / 2;
108        dst.pts[1].fY = (fPts[0].fY + fPts[1].fY) / 2;
109        dst.pts[2].fX = (fPts[0].fX + 2 * fPts[1].fX + fPts[2].fX) / 4;
110        dst.pts[2].fY = (fPts[0].fY + 2 * fPts[1].fY + fPts[2].fY) / 4;
111        dst.pts[3].fX = (fPts[0].fX + 3 * (fPts[1].fX + fPts[2].fX) + fPts[3].fX) / 8;
112        dst.pts[3].fY = (fPts[0].fY + 3 * (fPts[1].fY + fPts[2].fY) + fPts[3].fY) / 8;
113        dst.pts[4].fX = (fPts[1].fX + 2 * fPts[2].fX + fPts[3].fX) / 4;
114        dst.pts[4].fY = (fPts[1].fY + 2 * fPts[2].fY + fPts[3].fY) / 4;
115        dst.pts[5].fX = (fPts[2].fX + fPts[3].fX) / 2;
116        dst.pts[5].fY = (fPts[2].fY + fPts[3].fY) / 2;
117        dst.pts[6] = fPts[3];
118        return dst;
119    }
120    interp_cubic_coords(&fPts[0].fX, &dst.pts[0].fX, t);
121    interp_cubic_coords(&fPts[0].fY, &dst.pts[0].fY, t);
122    return dst;
123}
124
125void SkDCubic::Coefficients(const double* src, double* A, double* B, double* C, double* D) {
126    *A = src[6];  // d
127    *B = src[4] * 3;  // 3*c
128    *C = src[2] * 3;  // 3*b
129    *D = src[0];  // a
130    *A -= *D - *C + *B;     // A =   -a + 3*b - 3*c + d
131    *B += 3 * *D - 2 * *C;  // B =  3*a - 6*b + 3*c
132    *C -= 3 * *D;           // C = -3*a + 3*b
133}
134
135bool SkDCubic::endsAreExtremaInXOrY() const {
136    return (between(fPts[0].fX, fPts[1].fX, fPts[3].fX)
137            && between(fPts[0].fX, fPts[2].fX, fPts[3].fX))
138            || (between(fPts[0].fY, fPts[1].fY, fPts[3].fY)
139            && between(fPts[0].fY, fPts[2].fY, fPts[3].fY));
140}
141
142// Do a quick reject by rotating all points relative to a line formed by
143// a pair of one cubic's points. If the 2nd cubic's points
144// are on the line or on the opposite side from the 1st cubic's 'odd man', the
145// curves at most intersect at the endpoints.
146/* if returning true, check contains true if cubic's hull collapsed, making the cubic linear
147   if returning false, check contains true if the the cubic pair have only the end point in common
148*/
149bool SkDCubic::hullIntersects(const SkDPoint* pts, int ptCount, bool* isLinear) const {
150    bool linear = true;
151    char hullOrder[4];
152    int hullCount = convexHull(hullOrder);
153    int end1 = hullOrder[0];
154    int hullIndex = 0;
155    const SkDPoint* endPt[2];
156    endPt[0] = &fPts[end1];
157    do {
158        hullIndex = (hullIndex + 1) % hullCount;
159        int end2 = hullOrder[hullIndex];
160        endPt[1] = &fPts[end2];
161        double origX = endPt[0]->fX;
162        double origY = endPt[0]->fY;
163        double adj = endPt[1]->fX - origX;
164        double opp = endPt[1]->fY - origY;
165        int oddManMask = other_two(end1, end2);
166        int oddMan = end1 ^ oddManMask;
167        double sign = (fPts[oddMan].fY - origY) * adj - (fPts[oddMan].fX - origX) * opp;
168        int oddMan2 = end2 ^ oddManMask;
169        double sign2 = (fPts[oddMan2].fY - origY) * adj - (fPts[oddMan2].fX - origX) * opp;
170        if (sign * sign2 < 0) {
171            continue;
172        }
173        if (approximately_zero(sign)) {
174            sign = sign2;
175            if (approximately_zero(sign)) {
176                continue;
177            }
178        }
179        linear = false;
180        bool foundOutlier = false;
181        for (int n = 0; n < ptCount; ++n) {
182            double test = (pts[n].fY - origY) * adj - (pts[n].fX - origX) * opp;
183            if (test * sign > 0 && !precisely_zero(test)) {
184                foundOutlier = true;
185                break;
186            }
187        }
188        if (!foundOutlier) {
189            return false;
190        }
191        endPt[0] = endPt[1];
192        end1 = end2;
193    } while (hullIndex);
194    *isLinear = linear;
195    return true;
196}
197
198bool SkDCubic::hullIntersects(const SkDCubic& c2, bool* isLinear) const {
199    return hullIntersects(c2.fPts, c2.kPointCount, isLinear);
200}
201
202bool SkDCubic::hullIntersects(const SkDQuad& quad, bool* isLinear) const {
203    return hullIntersects(quad.fPts, quad.kPointCount, isLinear);
204}
205
206bool SkDCubic::hullIntersects(const SkDConic& conic, bool* isLinear) const {
207
208    return hullIntersects(conic.fPts, isLinear);
209}
210
211bool SkDCubic::isLinear(int startIndex, int endIndex) const {
212    if (fPts[0].approximatelyDEqual(fPts[3]))  {
213        return ((const SkDQuad *) this)->isLinear(0, 2);
214    }
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
232// from http://www.cs.sunysb.edu/~qin/courses/geometry/4.pdf
233// c(t)  = a(1-t)^3 + 3bt(1-t)^2 + 3c(1-t)t^2 + dt^3
234// c'(t) = -3a(1-t)^2 + 3b((1-t)^2 - 2t(1-t)) + 3c(2t(1-t) - t^2) + 3dt^2
235//       = 3(b-a)(1-t)^2 + 6(c-b)t(1-t) + 3(d-c)t^2
236static double derivative_at_t(const double* src, double t) {
237    double one_t = 1 - t;
238    double a = src[0];
239    double b = src[2];
240    double c = src[4];
241    double d = src[6];
242    return 3 * ((b - a) * one_t * one_t + 2 * (c - b) * t * one_t + (d - c) * t * t);
243}
244
245int SkDCubic::ComplexBreak(const SkPoint pointsPtr[4], SkScalar* t) {
246    SkDCubic cubic;
247    cubic.set(pointsPtr);
248    if (cubic.monotonicInX() && cubic.monotonicInY()) {
249        return 0;
250    }
251    double tt[2], ss[2];
252    SkCubicType cubicType = SkClassifyCubic(pointsPtr, tt, ss);
253    switch (cubicType) {
254        case SkCubicType::kLoop: {
255            const double &td = tt[0], &te = tt[1], &sd = ss[0], &se = ss[1];
256            if (roughly_between(0, td, sd) && roughly_between(0, te, se)) {
257                SkASSERT(roughly_between(0, td/sd, 1) && roughly_between(0, te/se, 1));
258                t[0] = static_cast<SkScalar>((td * se + te * sd) / (2 * sd * se));
259                SkASSERT(roughly_between(0, *t, 1));
260                return (int) (t[0] > 0 && t[0] < 1);
261            }
262        }
263        // fall through if no t value found
264        case SkCubicType::kSerpentine:
265        case SkCubicType::kLocalCusp:
266        case SkCubicType::kCuspAtInfinity: {
267            double inflectionTs[2];
268            int infTCount = cubic.findInflections(inflectionTs);
269            double maxCurvature[3];
270            int roots = cubic.findMaxCurvature(maxCurvature);
271    #if DEBUG_CUBIC_SPLIT
272            SkDebugf("%s\n", __FUNCTION__);
273            cubic.dump();
274            for (int index = 0; index < infTCount; ++index) {
275                SkDebugf("inflectionsTs[%d]=%1.9g ", index, inflectionTs[index]);
276                SkDPoint pt = cubic.ptAtT(inflectionTs[index]);
277                SkDVector dPt = cubic.dxdyAtT(inflectionTs[index]);
278                SkDLine perp = {{pt - dPt, pt + dPt}};
279                perp.dump();
280            }
281            for (int index = 0; index < roots; ++index) {
282                SkDebugf("maxCurvature[%d]=%1.9g ", index, maxCurvature[index]);
283                SkDPoint pt = cubic.ptAtT(maxCurvature[index]);
284                SkDVector dPt = cubic.dxdyAtT(maxCurvature[index]);
285                SkDLine perp = {{pt - dPt, pt + dPt}};
286                perp.dump();
287            }
288    #endif
289            if (infTCount == 2) {
290                for (int index = 0; index < roots; ++index) {
291                    if (between(inflectionTs[0], maxCurvature[index], inflectionTs[1])) {
292                        t[0] = maxCurvature[index];
293                        return (int) (t[0] > 0 && t[0] < 1);
294                    }
295                }
296            } else {
297                int resultCount = 0;
298                // FIXME: constant found through experimentation -- maybe there's a better way....
299                double precision = cubic.calcPrecision() * 2;
300                for (int index = 0; index < roots; ++index) {
301                    double testT = maxCurvature[index];
302                    if (0 >= testT || testT >= 1) {
303                        continue;
304                    }
305                    // don't call dxdyAtT since we want (0,0) results
306                    SkDVector dPt = { derivative_at_t(&cubic.fPts[0].fX, testT),
307                            derivative_at_t(&cubic.fPts[0].fY, testT) };
308                    double dPtLen = dPt.length();
309                    if (dPtLen < precision) {
310                        t[resultCount++] = testT;
311                    }
312                }
313                if (!resultCount && infTCount == 1) {
314                    t[0] = inflectionTs[0];
315                    resultCount = (int) (t[0] > 0 && t[0] < 1);
316                }
317                return resultCount;
318            }
319        }
320        default:
321            ;
322    }
323    return 0;
324}
325
326bool SkDCubic::monotonicInX() const {
327    return precisely_between(fPts[0].fX, fPts[1].fX, fPts[3].fX)
328            && precisely_between(fPts[0].fX, fPts[2].fX, fPts[3].fX);
329}
330
331bool SkDCubic::monotonicInY() const {
332    return precisely_between(fPts[0].fY, fPts[1].fY, fPts[3].fY)
333            && precisely_between(fPts[0].fY, fPts[2].fY, fPts[3].fY);
334}
335
336void SkDCubic::otherPts(int index, const SkDPoint* o1Pts[kPointCount - 1]) const {
337    int offset = (int) !SkToBool(index);
338    o1Pts[0] = &fPts[offset];
339    o1Pts[1] = &fPts[++offset];
340    o1Pts[2] = &fPts[++offset];
341}
342
343int SkDCubic::searchRoots(double extremeTs[6], int extrema, double axisIntercept,
344        SearchAxis xAxis, double* validRoots) const {
345    extrema += findInflections(&extremeTs[extrema]);
346    extremeTs[extrema++] = 0;
347    extremeTs[extrema] = 1;
348    SkASSERT(extrema < 6);
349    SkTQSort(extremeTs, extremeTs + extrema);
350    int validCount = 0;
351    for (int index = 0; index < extrema; ) {
352        double min = extremeTs[index];
353        double max = extremeTs[++index];
354        if (min == max) {
355            continue;
356        }
357        double newT = binarySearch(min, max, axisIntercept, xAxis);
358        if (newT >= 0) {
359            if (validCount >= 3) {
360                return 0;
361            }
362            validRoots[validCount++] = newT;
363        }
364    }
365    return validCount;
366}
367
368// cubic roots
369
370static const double PI = 3.141592653589793;
371
372// from SkGeometry.cpp (and Numeric Solutions, 5.6)
373int SkDCubic::RootsValidT(double A, double B, double C, double D, double t[3]) {
374    double s[3];
375    int realRoots = RootsReal(A, B, C, D, s);
376    int foundRoots = SkDQuad::AddValidTs(s, realRoots, t);
377    for (int index = 0; index < realRoots; ++index) {
378        double tValue = s[index];
379        if (!approximately_one_or_less(tValue) && between(1, tValue, 1.00005)) {
380            for (int idx2 = 0; idx2 < foundRoots; ++idx2) {
381                if (approximately_equal(t[idx2], 1)) {
382                    goto nextRoot;
383                }
384            }
385            SkASSERT(foundRoots < 3);
386            t[foundRoots++] = 1;
387        } else if (!approximately_zero_or_more(tValue) && between(-0.00005, tValue, 0)) {
388            for (int idx2 = 0; idx2 < foundRoots; ++idx2) {
389                if (approximately_equal(t[idx2], 0)) {
390                    goto nextRoot;
391                }
392            }
393            SkASSERT(foundRoots < 3);
394            t[foundRoots++] = 0;
395        }
396nextRoot:
397        ;
398    }
399    return foundRoots;
400}
401
402int SkDCubic::RootsReal(double A, double B, double C, double D, double s[3]) {
403#ifdef SK_DEBUG
404    // create a string mathematica understands
405    // GDB set print repe 15 # if repeated digits is a bother
406    //     set print elements 400 # if line doesn't fit
407    char str[1024];
408    sk_bzero(str, sizeof(str));
409    SK_SNPRINTF(str, sizeof(str), "Solve[%1.19g x^3 + %1.19g x^2 + %1.19g x + %1.19g == 0, x]",
410            A, B, C, D);
411    SkPathOpsDebug::MathematicaIze(str, sizeof(str));
412#if ONE_OFF_DEBUG && ONE_OFF_DEBUG_MATHEMATICA
413    SkDebugf("%s\n", str);
414#endif
415#endif
416    if (approximately_zero(A)
417            && approximately_zero_when_compared_to(A, B)
418            && approximately_zero_when_compared_to(A, C)
419            && approximately_zero_when_compared_to(A, D)) {  // we're just a quadratic
420        return SkDQuad::RootsReal(B, C, D, s);
421    }
422    if (approximately_zero_when_compared_to(D, A)
423            && approximately_zero_when_compared_to(D, B)
424            && approximately_zero_when_compared_to(D, C)) {  // 0 is one root
425        int num = SkDQuad::RootsReal(A, B, C, s);
426        for (int i = 0; i < num; ++i) {
427            if (approximately_zero(s[i])) {
428                return num;
429            }
430        }
431        s[num++] = 0;
432        return num;
433    }
434    if (approximately_zero(A + B + C + D)) {  // 1 is one root
435        int num = SkDQuad::RootsReal(A, A + B, -D, s);
436        for (int i = 0; i < num; ++i) {
437            if (AlmostDequalUlps(s[i], 1)) {
438                return num;
439            }
440        }
441        s[num++] = 1;
442        return num;
443    }
444    double a, b, c;
445    {
446        double invA = 1 / A;
447        a = B * invA;
448        b = C * invA;
449        c = D * invA;
450    }
451    double a2 = a * a;
452    double Q = (a2 - b * 3) / 9;
453    double R = (2 * a2 * a - 9 * a * b + 27 * c) / 54;
454    double R2 = R * R;
455    double Q3 = Q * Q * Q;
456    double R2MinusQ3 = R2 - Q3;
457    double adiv3 = a / 3;
458    double r;
459    double* roots = s;
460    if (R2MinusQ3 < 0) {   // we have 3 real roots
461        // the divide/root can, due to finite precisions, be slightly outside of -1...1
462        double theta = acos(SkTPin(R / sqrt(Q3), -1., 1.));
463        double neg2RootQ = -2 * sqrt(Q);
464
465        r = neg2RootQ * cos(theta / 3) - adiv3;
466        *roots++ = r;
467
468        r = neg2RootQ * cos((theta + 2 * PI) / 3) - adiv3;
469        if (!AlmostDequalUlps(s[0], r)) {
470            *roots++ = r;
471        }
472        r = neg2RootQ * cos((theta - 2 * PI) / 3) - adiv3;
473        if (!AlmostDequalUlps(s[0], r) && (roots - s == 1 || !AlmostDequalUlps(s[1], r))) {
474            *roots++ = r;
475        }
476    } else {  // we have 1 real root
477        double sqrtR2MinusQ3 = sqrt(R2MinusQ3);
478        double A = fabs(R) + sqrtR2MinusQ3;
479        A = SkDCubeRoot(A);
480        if (R > 0) {
481            A = -A;
482        }
483        if (A != 0) {
484            A += Q / A;
485        }
486        r = A - adiv3;
487        *roots++ = r;
488        if (AlmostDequalUlps((double) R2, (double) Q3)) {
489            r = -A / 2 - adiv3;
490            if (!AlmostDequalUlps(s[0], r)) {
491                *roots++ = r;
492            }
493        }
494    }
495    return static_cast<int>(roots - s);
496}
497
498// OPTIMIZE? compute t^2, t(1-t), and (1-t)^2 and pass them to another version of derivative at t?
499SkDVector SkDCubic::dxdyAtT(double t) const {
500    SkDVector result = { derivative_at_t(&fPts[0].fX, t), derivative_at_t(&fPts[0].fY, t) };
501    if (result.fX == 0 && result.fY == 0) {
502        if (t == 0) {
503            result = fPts[2] - fPts[0];
504        } else if (t == 1) {
505            result = fPts[3] - fPts[1];
506        } else {
507            // incomplete
508            SkDebugf("!c");
509        }
510        if (result.fX == 0 && result.fY == 0 && zero_or_one(t)) {
511            result = fPts[3] - fPts[0];
512        }
513    }
514    return result;
515}
516
517// OPTIMIZE? share code with formulate_F1DotF2
518int SkDCubic::findInflections(double tValues[]) const {
519    double Ax = fPts[1].fX - fPts[0].fX;
520    double Ay = fPts[1].fY - fPts[0].fY;
521    double Bx = fPts[2].fX - 2 * fPts[1].fX + fPts[0].fX;
522    double By = fPts[2].fY - 2 * fPts[1].fY + fPts[0].fY;
523    double Cx = fPts[3].fX + 3 * (fPts[1].fX - fPts[2].fX) - fPts[0].fX;
524    double Cy = fPts[3].fY + 3 * (fPts[1].fY - fPts[2].fY) - fPts[0].fY;
525    return SkDQuad::RootsValidT(Bx * Cy - By * Cx, Ax * Cy - Ay * Cx, Ax * By - Ay * Bx, tValues);
526}
527
528static void formulate_F1DotF2(const double src[], double coeff[4]) {
529    double a = src[2] - src[0];
530    double b = src[4] - 2 * src[2] + src[0];
531    double c = src[6] + 3 * (src[2] - src[4]) - src[0];
532    coeff[0] = c * c;
533    coeff[1] = 3 * b * c;
534    coeff[2] = 2 * b * b + c * a;
535    coeff[3] = a * b;
536}
537
538/** SkDCubic'(t) = At^2 + Bt + C, where
539    A = 3(-a + 3(b - c) + d)
540    B = 6(a - 2b + c)
541    C = 3(b - a)
542    Solve for t, keeping only those that fit between 0 < t < 1
543*/
544int SkDCubic::FindExtrema(const double src[], double tValues[2]) {
545    // we divide A,B,C by 3 to simplify
546    double a = src[0];
547    double b = src[2];
548    double c = src[4];
549    double d = src[6];
550    double A = d - a + 3 * (b - c);
551    double B = 2 * (a - b - b + c);
552    double C = b - a;
553
554    return SkDQuad::RootsValidT(A, B, C, tValues);
555}
556
557/*  from SkGeometry.cpp
558    Looking for F' dot F'' == 0
559
560    A = b - a
561    B = c - 2b + a
562    C = d - 3c + 3b - a
563
564    F' = 3Ct^2 + 6Bt + 3A
565    F'' = 6Ct + 6B
566
567    F' dot F'' -> CCt^3 + 3BCt^2 + (2BB + CA)t + AB
568*/
569int SkDCubic::findMaxCurvature(double tValues[]) const {
570    double coeffX[4], coeffY[4];
571    int i;
572    formulate_F1DotF2(&fPts[0].fX, coeffX);
573    formulate_F1DotF2(&fPts[0].fY, coeffY);
574    for (i = 0; i < 4; i++) {
575        coeffX[i] = coeffX[i] + coeffY[i];
576    }
577    return RootsValidT(coeffX[0], coeffX[1], coeffX[2], coeffX[3], tValues);
578}
579
580SkDPoint SkDCubic::ptAtT(double t) const {
581    if (0 == t) {
582        return fPts[0];
583    }
584    if (1 == t) {
585        return fPts[3];
586    }
587    double one_t = 1 - t;
588    double one_t2 = one_t * one_t;
589    double a = one_t2 * one_t;
590    double b = 3 * one_t2 * t;
591    double t2 = t * t;
592    double c = 3 * one_t * t2;
593    double d = t2 * t;
594    SkDPoint result = {a * fPts[0].fX + b * fPts[1].fX + c * fPts[2].fX + d * fPts[3].fX,
595            a * fPts[0].fY + b * fPts[1].fY + c * fPts[2].fY + d * fPts[3].fY};
596    return result;
597}
598
599/*
600 Given a cubic c, t1, and t2, find a small cubic segment.
601
602 The new cubic is defined as points A, B, C, and D, where
603 s1 = 1 - t1
604 s2 = 1 - t2
605 A = c[0]*s1*s1*s1 + 3*c[1]*s1*s1*t1 + 3*c[2]*s1*t1*t1 + c[3]*t1*t1*t1
606 D = c[0]*s2*s2*s2 + 3*c[1]*s2*s2*t2 + 3*c[2]*s2*t2*t2 + c[3]*t2*t2*t2
607
608 We don't have B or C. So We define two equations to isolate them.
609 First, compute two reference T values 1/3 and 2/3 from t1 to t2:
610
611 c(at (2*t1 + t2)/3) == E
612 c(at (t1 + 2*t2)/3) == F
613
614 Next, compute where those values must be if we know the values of B and C:
615
616 _12   =  A*2/3 + B*1/3
617 12_   =  A*1/3 + B*2/3
618 _23   =  B*2/3 + C*1/3
619 23_   =  B*1/3 + C*2/3
620 _34   =  C*2/3 + D*1/3
621 34_   =  C*1/3 + D*2/3
622 _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
623 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
624 _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
625 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
626 _1234 = (A*4/9 + B*4/9 + C*1/9)*2/3 + (B*4/9 + C*4/9 + D*1/9)*1/3
627       =  A*8/27 + B*12/27 + C*6/27 + D*1/27
628       =  E
629 1234_ = (A*1/9 + B*4/9 + C*4/9)*1/3 + (B*1/9 + C*4/9 + D*4/9)*2/3
630       =  A*1/27 + B*6/27 + C*12/27 + D*8/27
631       =  F
632 E*27  =  A*8    + B*12   + C*6     + D
633 F*27  =  A      + B*6    + C*12    + D*8
634
635Group the known values on one side:
636
637 M       = E*27 - A*8 - D     = B*12 + C* 6
638 N       = F*27 - A   - D*8   = B* 6 + C*12
639 M*2 - N = B*18
640 N*2 - M = C*18
641 B       = (M*2 - N)/18
642 C       = (N*2 - M)/18
643 */
644
645static double interp_cubic_coords(const double* src, double t) {
646    double ab = SkDInterp(src[0], src[2], t);
647    double bc = SkDInterp(src[2], src[4], t);
648    double cd = SkDInterp(src[4], src[6], t);
649    double abc = SkDInterp(ab, bc, t);
650    double bcd = SkDInterp(bc, cd, t);
651    double abcd = SkDInterp(abc, bcd, t);
652    return abcd;
653}
654
655SkDCubic SkDCubic::subDivide(double t1, double t2) const {
656    if (t1 == 0 || t2 == 1) {
657        if (t1 == 0 && t2 == 1) {
658            return *this;
659        }
660        SkDCubicPair pair = chopAt(t1 == 0 ? t2 : t1);
661        SkDCubic dst = t1 == 0 ? pair.first() : pair.second();
662        return dst;
663    }
664    SkDCubic dst;
665    double ax = dst[0].fX = interp_cubic_coords(&fPts[0].fX, t1);
666    double ay = dst[0].fY = interp_cubic_coords(&fPts[0].fY, t1);
667    double ex = interp_cubic_coords(&fPts[0].fX, (t1*2+t2)/3);
668    double ey = interp_cubic_coords(&fPts[0].fY, (t1*2+t2)/3);
669    double fx = interp_cubic_coords(&fPts[0].fX, (t1+t2*2)/3);
670    double fy = interp_cubic_coords(&fPts[0].fY, (t1+t2*2)/3);
671    double dx = dst[3].fX = interp_cubic_coords(&fPts[0].fX, t2);
672    double dy = dst[3].fY = interp_cubic_coords(&fPts[0].fY, t2);
673    double mx = ex * 27 - ax * 8 - dx;
674    double my = ey * 27 - ay * 8 - dy;
675    double nx = fx * 27 - ax - dx * 8;
676    double ny = fy * 27 - ay - dy * 8;
677    /* bx = */ dst[1].fX = (mx * 2 - nx) / 18;
678    /* by = */ dst[1].fY = (my * 2 - ny) / 18;
679    /* cx = */ dst[2].fX = (nx * 2 - mx) / 18;
680    /* cy = */ dst[2].fY = (ny * 2 - my) / 18;
681    // FIXME: call align() ?
682    return dst;
683}
684
685void SkDCubic::subDivide(const SkDPoint& a, const SkDPoint& d,
686                         double t1, double t2, SkDPoint dst[2]) const {
687    SkASSERT(t1 != t2);
688    // this approach assumes that the control points computed directly are accurate enough
689    SkDCubic sub = subDivide(t1, t2);
690    dst[0] = sub[1] + (a - sub[0]);
691    dst[1] = sub[2] + (d - sub[3]);
692    if (t1 == 0 || t2 == 0) {
693        align(0, 1, t1 == 0 ? &dst[0] : &dst[1]);
694    }
695    if (t1 == 1 || t2 == 1) {
696        align(3, 2, t1 == 1 ? &dst[0] : &dst[1]);
697    }
698    if (AlmostBequalUlps(dst[0].fX, a.fX)) {
699        dst[0].fX = a.fX;
700    }
701    if (AlmostBequalUlps(dst[0].fY, a.fY)) {
702        dst[0].fY = a.fY;
703    }
704    if (AlmostBequalUlps(dst[1].fX, d.fX)) {
705        dst[1].fX = d.fX;
706    }
707    if (AlmostBequalUlps(dst[1].fY, d.fY)) {
708        dst[1].fY = d.fY;
709    }
710}
711
712bool SkDCubic::toFloatPoints(SkPoint* pts) const {
713    const double* dCubic = &fPts[0].fX;
714    SkScalar* cubic = &pts[0].fX;
715    for (int index = 0; index < kPointCount * 2; ++index) {
716        cubic[index] = SkDoubleToScalar(dCubic[index]);
717        if (SkScalarAbs(cubic[index]) < FLT_EPSILON_ORDERABLE_ERR) {
718            cubic[index] = 0;
719        }
720    }
721    return SkScalarsAreFinite(&pts->fX, kPointCount * 2);
722}
723
724double SkDCubic::top(const SkDCubic& dCurve, double startT, double endT, SkDPoint*topPt) const {
725    double extremeTs[2];
726    double topT = -1;
727    int roots = SkDCubic::FindExtrema(&fPts[0].fY, extremeTs);
728    for (int index = 0; index < roots; ++index) {
729        double t = startT + (endT - startT) * extremeTs[index];
730        SkDPoint mid = dCurve.ptAtT(t);
731        if (topPt->fY > mid.fY || (topPt->fY == mid.fY && topPt->fX > mid.fX)) {
732            topT = t;
733            *topPt = mid;
734        }
735    }
736    return topT;
737}
738