CubicToQuadratics.cpp revision 73ca6243b31e225e9fd5b75a96cbc82d62557de6 
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
24Cubic 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 "CubicUtilities.h"
49#include "CurveIntersection.h"
50#include "LineIntersection.h"
51
52const bool AVERAGE_END_POINTS = true; // results in better fitting curves
53
54#define USE_CUBIC_END_POINTS 1
55
56static double calcTDiv(const Cubic& cubic, double precision, double start) {
57    const double adjust = sqrt(3) / 36;
58    Cubic sub;
59    const Cubic* cPtr;
60    if (start == 0) {
61        cPtr = &cubic;
62    } else {
63        // OPTIMIZE: special-case half-split ?
64        sub_divide(cubic, start, 1, sub);
65        cPtr = &sub;
66    }
67    const Cubic& c = *cPtr;
68    double dx = c[3].x - 3 * (c[2].x - c[1].x) - c[0].x;
69    double dy = c[3].y - 3 * (c[2].y - c[1].y) - c[0].y;
70    double dist = sqrt(dx * dx + dy * dy);
71    double tDiv3 = precision / (adjust * dist);
72    double t = cube_root(tDiv3);
73    if (start > 0) {
74        t = start + (1 - start) * t;
75    }
76    return t;
77}
78
79void demote_cubic_to_quad(const Cubic& cubic, Quadratic& quad) {
80    quad[0] = cubic[0];
81if (AVERAGE_END_POINTS) {
82    const _Point fromC1 = { (3 * cubic[1].x - cubic[0].x) / 2, (3 * cubic[1].y - cubic[0].y) / 2 };
83    const _Point fromC2 = { (3 * cubic[2].x - cubic[3].x) / 2, (3 * cubic[2].y - cubic[3].y) / 2 };
84    quad[1].x = (fromC1.x + fromC2.x) / 2;
85    quad[1].y = (fromC1.y + fromC2.y) / 2;
86} else {
87    lineIntersect((const _Line&) cubic[0], (const _Line&) cubic[2], quad[1]);
88}
89    quad[2] = cubic[3];
90}
91
92int cubic_to_quadratics(const Cubic& cubic, double precision, SkTDArray<Quadratic>& quadratics) {
93    SkTDArray<double> ts;
94    cubic_to_quadratics(cubic, precision, ts);
95    int tsCount = ts.count();
96    double t1Start = 0;
97    int order = 0;
98    for (int idx = 0; idx <= tsCount; ++idx) {
99        double t1 = idx < tsCount ? ts[idx] : 1;
100        Cubic part;
101        sub_divide(cubic, t1Start, t1, part);
102        Quadratic q1;
103        demote_cubic_to_quad(part, q1);
104        Quadratic s1;
105        int o1 = reduceOrder(q1, s1);
106        if (order < o1) {
107            order = o1;
108        }
109        memcpy(quadratics.append(), o1 < 2 ? s1 : q1, sizeof(Quadratic));
110        t1Start = t1;
111    }
112    return order;
113}
114
115static bool addSimpleTs(const Cubic& cubic, double precision, SkTDArray<double>& ts) {
116    double tDiv = calcTDiv(cubic, precision, 0);
117    if (tDiv >= 1) {
118        return true;
119    }
120    if (tDiv >= 0.5) {
121        *ts.append() = 0.5;
122        return true;
123    }
124    return false;
125}
126
127static void addTs(const Cubic& cubic, double precision, double start, double end,
128        SkTDArray<double>& ts) {
129    double tDiv = calcTDiv(cubic, precision, 0);
130    double parts = ceil(1.0 / tDiv);
131    for (double index = 0; index < parts; ++index) {
132        double newT = start + (index / parts) * (end - start);
133        if (newT > 0 && newT < 1) {
134            *ts.append() = newT;
135        }
136    }
137}
138
139// flavor that returns T values only, deferring computing the quads until they are needed
140void cubic_to_quadratics(const Cubic& cubic, double precision, SkTDArray<double>& ts) {
141    Cubic reduced;
142    int order = reduceOrder(cubic, reduced, kReduceOrder_QuadraticsAllowed);
143    if (order < 3) {
144        return;
145    }
146    double inflectT[2];
147    int inflections = find_cubic_inflections(cubic, inflectT);
148    SkASSERT(inflections <= 2);
149    if (inflections == 0 && addSimpleTs(cubic, precision, ts)) {
150        return;
151    }
152    if (inflections == 1) {
153        CubicPair pair;
154        chop_at(cubic, pair, inflectT[0]);
155        addTs(pair.first(), precision, 0, inflectT[0], ts);
156        addTs(pair.second(), precision, inflectT[0], 1, ts);
157        return;
158    }
159    if (inflections == 2) {
160        if (inflectT[0] > inflectT[1]) {
161            SkTSwap(inflectT[0], inflectT[1]);
162        }
163        Cubic part;
164        sub_divide(cubic, 0, inflectT[0], part);
165        addTs(part, precision, 0, inflectT[0], ts);
166        sub_divide(cubic, inflectT[0], inflectT[1], part);
167        addTs(part, precision, inflectT[0], inflectT[1], ts);
168        sub_divide(cubic, inflectT[1], 1, part);
169        addTs(part, precision, inflectT[1], 1, ts);
170        return;
171    }
172    addTs(cubic, precision, 0, 1, ts);
173}
174