1/*
2http://stackoverflow.com/questions/2009160/how-do-i-convert-the-2-control-points-of-a-cubic-curve-to-the-single-control-poi
3*/
4
5/*
6Let's call the control points of the cubic Q0..Q3 and the control points of the quadratic P0..P2.
7Then for degree elevation, the equations are:
8
9Q0 = P0
10Q1 = 1/3 P0 + 2/3 P1
11Q2 = 2/3 P1 + 1/3 P2
12Q3 = P2
13In your case you have Q0..Q3 and you're solving for P0..P2. There are two ways to compute P1 from
14 the equations above:
15
16P1 = 3/2 Q1 - 1/2 Q0
17P1 = 3/2 Q2 - 1/2 Q3
18If this is a degree-elevated cubic, then both equations will give the same answer for P1. Since
19 it's likely not, your best bet is to average them. So,
20
21P1 = -1/4 Q0 + 3/4 Q1 + 3/4 Q2 - 1/4 Q3
22
23
24SkDCubic defined by: P1/2 - anchor points, C1/C2 control points
25|x| is the euclidean norm of x
26mid-point approx of cubic: a quad that shares the same anchors with the cubic and has the
27 control point at C = (3·C2 - P2 + 3·C1 - P1)/4
28
29Algorithm
30
31pick an absolute precision (prec)
32Compute the Tdiv as the root of (cubic) equation
33sqrt(3)/18 · |P2 - 3·C2 + 3·C1 - P1|/2 · Tdiv ^ 3 = prec
34if Tdiv < 0.5 divide the cubic at Tdiv. First segment [0..Tdiv] can be approximated with by a
35 quadratic, with a defect less than prec, by the mid-point approximation.
36 Repeat from step 2 with the second resulted segment (corresponding to 1-Tdiv)
370.5<=Tdiv<1 - simply divide the cubic in two. The two halves can be approximated by the mid-point
38 approximation
39Tdiv>=1 - the entire cubic can be approximated by the mid-point approximation
40
41confirmed by (maybe stolen from)
42http://www.caffeineowl.com/graphics/2d/vectorial/cubic2quad01.html
43// maybe in turn derived from  http://www.cccg.ca/proceedings/2004/36.pdf
44// also stored at http://www.cis.usouthal.edu/~hain/general/Publications/Bezier/bezier%20cccg04%20paper.pdf
45
46*/
47
48#include "SkPathOpsCubic.h"
49#include "SkPathOpsLine.h"
50#include "SkPathOpsQuad.h"
51#include "SkReduceOrder.h"
52#include "SkTArray.h"
53#include "SkTSort.h"
54
55#define USE_CUBIC_END_POINTS 1
56
57static double calc_t_div(const SkDCubic& cubic, double precision, double start) {
58    const double adjust = sqrt(3.) / 36;
59    SkDCubic sub;
60    const SkDCubic* cPtr;
61    if (start == 0) {
62        cPtr = &cubic;
63    } else {
64        // OPTIMIZE: special-case half-split ?
65        sub = cubic.subDivide(start, 1);
66        cPtr = &sub;
67    }
68    const SkDCubic& c = *cPtr;
69    double dx = c[3].fX - 3 * (c[2].fX - c[1].fX) - c[0].fX;
70    double dy = c[3].fY - 3 * (c[2].fY - c[1].fY) - c[0].fY;
71    double dist = sqrt(dx * dx + dy * dy);
72    double tDiv3 = precision / (adjust * dist);
73    double t = SkDCubeRoot(tDiv3);
74    if (start > 0) {
75        t = start + (1 - start) * t;
76    }
77    return t;
78}
79
80SkDQuad SkDCubic::toQuad() const {
81    SkDQuad quad;
82    quad[0] = fPts[0];
83    const SkDPoint fromC1 = {(3 * fPts[1].fX - fPts[0].fX) / 2, (3 * fPts[1].fY - fPts[0].fY) / 2};
84    const SkDPoint fromC2 = {(3 * fPts[2].fX - fPts[3].fX) / 2, (3 * fPts[2].fY - fPts[3].fY) / 2};
85    quad[1].fX = (fromC1.fX + fromC2.fX) / 2;
86    quad[1].fY = (fromC1.fY + fromC2.fY) / 2;
87    quad[2] = fPts[3];
88    return quad;
89}
90
91static bool add_simple_ts(const SkDCubic& cubic, double precision, SkTArray<double, true>* ts) {
92    double tDiv = calc_t_div(cubic, precision, 0);
93    if (tDiv >= 1) {
94        return true;
95    }
96    if (tDiv >= 0.5) {
97        ts->push_back(0.5);
98        return true;
99    }
100    return false;
101}
102
103static void addTs(const SkDCubic& cubic, double precision, double start, double end,
104        SkTArray<double, true>* ts) {
105    double tDiv = calc_t_div(cubic, precision, 0);
106    double parts = ceil(1.0 / tDiv);
107    for (double index = 0; index < parts; ++index) {
108        double newT = start + (index / parts) * (end - start);
109        if (newT > 0 && newT < 1) {
110            ts->push_back(newT);
111        }
112    }
113}
114
115// flavor that returns T values only, deferring computing the quads until they are needed
116// FIXME: when called from recursive intersect 2, this could take the original cubic
117// and do a more precise job when calling chop at and sub divide by computing the fractional ts.
118// it would still take the prechopped cubic for reduce order and find cubic inflections
119void SkDCubic::toQuadraticTs(double precision, SkTArray<double, true>* ts) const {
120    SkReduceOrder reducer;
121    int order = reducer.reduce(*this, SkReduceOrder::kAllow_Quadratics, SkReduceOrder::kFill_Style);
122    if (order < 3) {
123        return;
124    }
125    double inflectT[5];
126    int inflections = findInflections(inflectT);
127    SkASSERT(inflections <= 2);
128    if (!endsAreExtremaInXOrY()) {
129        inflections += findMaxCurvature(&inflectT[inflections]);
130        SkASSERT(inflections <= 5);
131    }
132    SkTQSort<double>(inflectT, &inflectT[inflections - 1]);
133    // OPTIMIZATION: is this filtering common enough that it needs to be pulled out into its
134    // own subroutine?
135    while (inflections && approximately_less_than_zero(inflectT[0])) {
136        memmove(inflectT, &inflectT[1], sizeof(inflectT[0]) * --inflections);
137    }
138    int start = 0;
139    do {
140        int next = start + 1;
141        if (next >= inflections) {
142            break;
143        }
144        if (!approximately_equal(inflectT[start], inflectT[next])) {
145            ++start;
146            continue;
147        }
148        memmove(&inflectT[start], &inflectT[next], sizeof(inflectT[0]) * (--inflections - start));
149    } while (true);
150    while (inflections && approximately_greater_than_one(inflectT[inflections - 1])) {
151        --inflections;
152    }
153    SkDCubicPair pair;
154    if (inflections == 1) {
155        pair = chopAt(inflectT[0]);
156        int orderP1 = reducer.reduce(pair.first(), SkReduceOrder::kNo_Quadratics,
157                SkReduceOrder::kFill_Style);
158        if (orderP1 < 2) {
159            --inflections;
160        } else {
161            int orderP2 = reducer.reduce(pair.second(), SkReduceOrder::kNo_Quadratics,
162                    SkReduceOrder::kFill_Style);
163            if (orderP2 < 2) {
164                --inflections;
165            }
166        }
167    }
168    if (inflections == 0 && add_simple_ts(*this, precision, ts)) {
169        return;
170    }
171    if (inflections == 1) {
172        pair = chopAt(inflectT[0]);
173        addTs(pair.first(), precision, 0, inflectT[0], ts);
174        addTs(pair.second(), precision, inflectT[0], 1, ts);
175        return;
176    }
177    if (inflections > 1) {
178        SkDCubic part = subDivide(0, inflectT[0]);
179        addTs(part, precision, 0, inflectT[0], ts);
180        int last = inflections - 1;
181        for (int idx = 0; idx < last; ++idx) {
182            part = subDivide(inflectT[idx], inflectT[idx + 1]);
183            addTs(part, precision, inflectT[idx], inflectT[idx + 1], ts);
184        }
185        part = subDivide(inflectT[last], 1);
186        addTs(part, precision, inflectT[last], 1, ts);
187        return;
188    }
189    addTs(*this, precision, 0, 1, ts);
190}
191