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