CubicReduceOrder.cpp revision b45a1b46ee25e9b19800b028bb1ca925212ac7b4 
1#include "CurveIntersection.h"
2#include "Extrema.h"
3#include "IntersectionUtilities.h"
4#include "LineParameters.h"
5
6static double interp_cubic_coords(const double* src, double t)
7{
8    double ab = interp(src[0], src[2], t);
9    double bc = interp(src[2], src[4], t);
10    double cd = interp(src[4], src[6], t);
11    double abc = interp(ab, bc, t);
12    double bcd = interp(bc, cd, t);
13    return interp(abc, bcd, t);
14}
15
16static int coincident_line(const Cubic& cubic, Cubic& reduction) {
17    reduction[0] = reduction[1] = cubic[0];
18    return 1;
19}
20
21static int vertical_line(const Cubic& cubic, Cubic& reduction) {
22    double tValues[2];
23    reduction[0] = cubic[0];
24    reduction[1] = cubic[3];
25    int smaller = reduction[1].y > reduction[0].y;
26    int larger = smaller ^ 1;
27    int roots = findExtrema(cubic[0].y, cubic[1].y, cubic[2].y, cubic[3].y, tValues);
28    for (int index = 0; index < roots; ++index) {
29        double yExtrema = interp_cubic_coords(&cubic[0].y, tValues[index]);
30        if (reduction[smaller].y > yExtrema) {
31            reduction[smaller].y = yExtrema;
32            continue;
33        }
34        if (reduction[larger].y < yExtrema) {
35            reduction[larger].y = yExtrema;
36        }
37    }
38    return 2;
39}
40
41static int horizontal_line(const Cubic& cubic, Cubic& reduction) {
42    double tValues[2];
43    reduction[0] = cubic[0];
44    reduction[1] = cubic[3];
45    int smaller = reduction[1].x > reduction[0].x;
46    int larger = smaller ^ 1;
47    int roots = findExtrema(cubic[0].x, cubic[1].x, cubic[2].x, cubic[3].x, tValues);
48    for (int index = 0; index < roots; ++index) {
49        double xExtrema = interp_cubic_coords(&cubic[0].x, tValues[index]);
50        if (reduction[smaller].x > xExtrema) {
51            reduction[smaller].x = xExtrema;
52            continue;
53        }
54        if (reduction[larger].x < xExtrema) {
55            reduction[larger].x = xExtrema;
56        }
57    }
58    return 2;
59}
60
61// check to see if it is a quadratic or a line
62static int check_quadratic(const Cubic& cubic, Cubic& reduction,
63        int minX, int maxX, int minY, int maxY) {
64    double dx10 = cubic[1].x - cubic[0].x;
65    double dx23 = cubic[2].x - cubic[3].x;
66    double midX = cubic[0].x + dx10 * 3 / 2;
67    if (!approximately_equal(midX - cubic[3].x, dx23 * 3 / 2)) {
68        return 0;
69    }
70    double dy10 = cubic[1].y - cubic[0].y;
71    double dy23 = cubic[2].y - cubic[3].y;
72    double midY = cubic[0].y + dy10 * 3 / 2;
73    if (!approximately_equal(midY - cubic[3].y, dy23 * 3 / 2)) {
74        return 0;
75    }
76    reduction[0] = cubic[0];
77    reduction[1].x = midX;
78    reduction[1].y = midY;
79    reduction[2] = cubic[3];
80    return 3;
81}
82
83static int check_linear(const Cubic& cubic, Cubic& reduction,
84        int minX, int maxX, int minY, int maxY) {
85    int startIndex = 0;
86    int endIndex = 3;
87    while (cubic[startIndex].approximatelyEqual(cubic[endIndex])) {
88        --endIndex;
89        if (endIndex == 0) {
90            printf("%s shouldn't get here if all four points are about equal", __FUNCTION__);
91            assert(0);
92        }
93    }
94    if (!isLinear(cubic, startIndex, endIndex)) {
95        return 0;
96    }
97    // four are colinear: return line formed by outside
98    reduction[0] = cubic[0];
99    reduction[1] = cubic[3];
100    int sameSide1;
101    int sameSide2;
102    bool useX = cubic[maxX].x - cubic[minX].x >= cubic[maxY].y - cubic[minY].y;
103    if (useX) {
104        sameSide1 = sign(cubic[0].x - cubic[1].x) + sign(cubic[3].x - cubic[1].x);
105        sameSide2 = sign(cubic[0].x - cubic[2].x) + sign(cubic[3].x - cubic[2].x);
106    } else {
107        sameSide1 = sign(cubic[0].y - cubic[1].y) + sign(cubic[3].y - cubic[1].y);
108        sameSide2 = sign(cubic[0].y - cubic[2].y) + sign(cubic[3].y - cubic[2].y);
109    }
110    if (sameSide1 == sameSide2 && (sameSide1 & 3) != 2) {
111        return 2;
112    }
113    double tValues[2];
114    int roots;
115    if (useX) {
116        roots = findExtrema(cubic[0].x, cubic[1].x, cubic[2].x, cubic[3].x, tValues);
117    } else {
118        roots = findExtrema(cubic[0].y, cubic[1].y, cubic[2].y, cubic[3].y, tValues);
119    }
120    for (int index = 0; index < roots; ++index) {
121        _Point extrema;
122        extrema.x = interp_cubic_coords(&cubic[0].x, tValues[index]);
123        extrema.y = interp_cubic_coords(&cubic[0].y, tValues[index]);
124        // sameSide > 0 means mid is smaller than either [0] or [3], so replace smaller
125        int replace;
126        if (useX) {
127            if (extrema.x < cubic[0].x ^ extrema.x < cubic[3].x) {
128                continue;
129            }
130            replace = (extrema.x < cubic[0].x | extrema.x < cubic[3].x)
131                    ^ cubic[0].x < cubic[3].x;
132        } else {
133            if (extrema.y < cubic[0].y ^ extrema.y < cubic[3].y) {
134                continue;
135            }
136            replace = (extrema.y < cubic[0].y | extrema.y < cubic[3].y)
137                    ^ cubic[0].y < cubic[3].y;
138        }
139        reduction[replace] = extrema;
140    }
141    return 2;
142}
143
144bool isLinear(const Cubic& cubic, int startIndex, int endIndex) {
145    LineParameters lineParameters;
146    lineParameters.cubicEndPoints(cubic, startIndex, endIndex);
147    double normalSquared = lineParameters.normalSquared();
148    double distance[2]; // distance is not normalized
149    int mask = other_two(startIndex, endIndex);
150    int inner1 = startIndex ^ mask;
151    int inner2 = endIndex ^ mask;
152    lineParameters.controlPtDistance(cubic, inner1, inner2, distance);
153    double limit = normalSquared;
154    int index;
155    for (index = 0; index < 2; ++index) {
156        double distSq = distance[index];
157        distSq *= distSq;
158        if (approximately_greater(distSq, limit)) {
159            return false;
160        }
161    }
162    return true;
163}
164
165/* food for thought:
166http://objectmix.com/graphics/132906-fast-precision-driven-cubic-quadratic-piecewise-degree-reduction-algos-2-a.html
167
168Given points c1, c2, c3 and c4 of a cubic Bezier, the points of the
169corresponding quadratic Bezier are (given in convex combinations of
170points):
171
172q1 = (11/13)c1 + (3/13)c2 -(3/13)c3 + (2/13)c4
173q2 = -c1 + (3/2)c2 + (3/2)c3 - c4
174q3 = (2/13)c1 - (3/13)c2 + (3/13)c3 + (11/13)c4
175
176Of course, this curve does not interpolate the end-points, but it would
177be interesting to see the behaviour of such a curve in an applet.
178
179--
180Kalle Rutanen
181http://kaba.hilvi.org
182
183*/
184
185// reduce to a quadratic or smaller
186// look for identical points
187// look for all four points in a line
188    // note that three points in a line doesn't simplify a cubic
189// look for approximation with single quadratic
190    // save approximation with multiple quadratics for later
191int reduceOrder(const Cubic& cubic, Cubic& reduction, ReduceOrder_Flags allowQuadratics) {
192    int index, minX, maxX, minY, maxY;
193    int minXSet, minYSet;
194    minX = maxX = minY = maxY = 0;
195    minXSet = minYSet = 0;
196    for (index = 1; index < 4; ++index) {
197        if (cubic[minX].x > cubic[index].x) {
198            minX = index;
199        }
200        if (cubic[minY].y > cubic[index].y) {
201            minY = index;
202        }
203        if (cubic[maxX].x < cubic[index].x) {
204            maxX = index;
205        }
206        if (cubic[maxY].y < cubic[index].y) {
207            maxY = index;
208        }
209    }
210    for (index = 0; index < 4; ++index) {
211        if (approximately_equal(cubic[index].x, cubic[minX].x)) {
212            minXSet |= 1 << index;
213        }
214        if (approximately_equal(cubic[index].y, cubic[minY].y)) {
215            minYSet |= 1 << index;
216        }
217    }
218    if (minXSet == 0xF) { // test for vertical line
219        if (minYSet == 0xF) { // return 1 if all four are coincident
220            return coincident_line(cubic, reduction);
221        }
222        return vertical_line(cubic, reduction);
223    }
224    if (minYSet == 0xF) { // test for horizontal line
225        return horizontal_line(cubic, reduction);
226    }
227    int result = check_linear(cubic, reduction, minX, maxX, minY, maxY);
228    if (result) {
229        return result;
230    }
231    if (allowQuadratics && (result = check_quadratic(cubic, reduction, minX, maxX, minY, maxY))) {
232        return result;
233    }
234    memcpy(reduction, cubic, sizeof(Cubic));
235    return 4;
236}
237