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