LineCubicIntersection.cpp revision 9f60291c5375457f8adf228dbe6e8ff1186b13e1 
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 "CubicUtilities.h"
9#include "Intersections.h"
10#include "LineUtilities.h"
11
12/*
13Find the interection of a line and cubic by solving for valid t values.
14
15Analogous to line-quadratic intersection, solve line-cubic intersection by
16representing the cubic as:
17  x = a(1-t)^3 + 2b(1-t)^2t + c(1-t)t^2 + dt^3
18  y = e(1-t)^3 + 2f(1-t)^2t + g(1-t)t^2 + ht^3
19and the line as:
20  y = i*x + j  (if the line is more horizontal)
21or:
22  x = i*y + j  (if the line is more vertical)
23
24Then using Mathematica, solve for the values of t where the cubic intersects the
25line:
26
27  (in) Resultant[
28        a*(1 - t)^3 + 3*b*(1 - t)^2*t + 3*c*(1 - t)*t^2 + d*t^3 - x,
29        e*(1 - t)^3 + 3*f*(1 - t)^2*t + 3*g*(1 - t)*t^2 + h*t^3 - i*x - j, x]
30  (out) -e     +   j     +
31       3 e t   - 3 f t   -
32       3 e t^2 + 6 f t^2 - 3 g t^2 +
33         e t^3 - 3 f t^3 + 3 g t^3 - h t^3 +
34     i ( a     -
35       3 a t + 3 b t +
36       3 a t^2 - 6 b t^2 + 3 c t^2 -
37         a t^3 + 3 b t^3 - 3 c t^3 + d t^3 )
38
39if i goes to infinity, we can rewrite the line in terms of x. Mathematica:
40
41  (in) Resultant[
42        a*(1 - t)^3 + 3*b*(1 - t)^2*t + 3*c*(1 - t)*t^2 + d*t^3 - i*y - j,
43        e*(1 - t)^3 + 3*f*(1 - t)^2*t + 3*g*(1 - t)*t^2 + h*t^3 - y,       y]
44  (out)  a     -   j     -
45       3 a t   + 3 b t   +
46       3 a t^2 - 6 b t^2 + 3 c t^2 -
47         a t^3 + 3 b t^3 - 3 c t^3 + d t^3 -
48     i ( e     -
49       3 e t   + 3 f t   +
50       3 e t^2 - 6 f t^2 + 3 g t^2 -
51         e t^3 + 3 f t^3 - 3 g t^3 + h t^3 )
52
53Solving this with Mathematica produces an expression with hundreds of terms;
54instead, use Numeric Solutions recipe to solve the cubic.
55
56The near-horizontal case, in terms of:  Ax^3 + Bx^2 + Cx + D == 0
57    A =   (-(-e + 3*f - 3*g + h) + i*(-a + 3*b - 3*c + d)     )
58    B = 3*(-( e - 2*f +   g    ) + i*( a - 2*b +   c    )     )
59    C = 3*(-(-e +   f          ) + i*(-a +   b          )     )
60    D =   (-( e                ) + i*( a                ) + j )
61
62The near-vertical case, in terms of:  Ax^3 + Bx^2 + Cx + D == 0
63    A =   ( (-a + 3*b - 3*c + d) - i*(-e + 3*f - 3*g + h)     )
64    B = 3*( ( a - 2*b +   c    ) - i*( e - 2*f +   g    )     )
65    C = 3*( (-a +   b          ) - i*(-e +   f          )     )
66    D =   ( ( a                ) - i*( e                ) - j )
67
68For horizontal lines:
69(in) Resultant[
70      a*(1 - t)^3 + 3*b*(1 - t)^2*t + 3*c*(1 - t)*t^2 + d*t^3 - j,
71      e*(1 - t)^3 + 3*f*(1 - t)^2*t + 3*g*(1 - t)*t^2 + h*t^3 - y, y]
72(out)  e     -   j     -
73     3 e t   + 3 f t   +
74     3 e t^2 - 6 f t^2 + 3 g t^2 -
75       e t^3 + 3 f t^3 - 3 g t^3 + h t^3
76So the cubic coefficients are:
77
78 */
79
80class LineCubicIntersections {
81public:
82
83LineCubicIntersections(const Cubic& c, const _Line& l, Intersections& i)
84    : cubic(c)
85    , line(l)
86    , intersections(i) {
87}
88
89// see parallel routine in line quadratic intersections
90int intersectRay(double roots[3]) {
91    double adj = line[1].x - line[0].x;
92    double opp = line[1].y - line[0].y;
93    Cubic r;
94    for (int n = 0; n < 4; ++n) {
95        r[n].x = (cubic[n].y - line[0].y) * adj - (cubic[n].x - line[0].x) * opp;
96    }
97    double A, B, C, D;
98    coefficients(&r[0].x, A, B, C, D);
99    return cubicRootsValidT(A, B, C, D, roots);
100}
101
102int intersect() {
103    addEndPoints();
104    double rootVals[3];
105    int roots = intersectRay(rootVals);
106    for (int index = 0; index < roots; ++index) {
107        double cubicT = rootVals[index];
108        double lineT = findLineT(cubicT);
109        if (pinTs(cubicT, lineT)) {
110            intersections.insert(cubicT, lineT);
111        }
112    }
113    return intersections.fUsed;
114}
115
116int horizontalIntersect(double axisIntercept, double roots[3]) {
117    double A, B, C, D;
118    coefficients(&cubic[0].y, A, B, C, D);
119    D -= axisIntercept;
120    return cubicRootsValidT(A, B, C, D, roots);
121}
122
123int horizontalIntersect(double axisIntercept, double left, double right, bool flipped) {
124    addHorizontalEndPoints(left, right, axisIntercept);
125    double rootVals[3];
126    int roots = horizontalIntersect(axisIntercept, rootVals);
127    for (int index = 0; index < roots; ++index) {
128        double x;
129        double cubicT = rootVals[index];
130        xy_at_t(cubic, cubicT, x, *(double*) NULL);
131        double lineT = (x - left) / (right - left);
132        if (pinTs(cubicT, lineT)) {
133            intersections.insert(cubicT, lineT);
134        }
135    }
136    if (flipped) {
137        flip();
138    }
139    return intersections.fUsed;
140}
141
142int verticalIntersect(double axisIntercept, double roots[3]) {
143    double A, B, C, D;
144    coefficients(&cubic[0].x, A, B, C, D);
145    D -= axisIntercept;
146    return cubicRootsValidT(A, B, C, D, roots);
147}
148
149int verticalIntersect(double axisIntercept, double top, double bottom, bool flipped) {
150    addVerticalEndPoints(top, bottom, axisIntercept);
151    double rootVals[3];
152    int roots = verticalIntersect(axisIntercept, rootVals);
153    for (int index = 0; index < roots; ++index) {
154        double y;
155        double cubicT = rootVals[index];
156        xy_at_t(cubic, cubicT, *(double*) NULL, y);
157        double lineT = (y - top) / (bottom - top);
158        if (pinTs(cubicT, lineT)) {
159            intersections.insert(cubicT, lineT);
160        }
161    }
162    if (flipped) {
163        flip();
164    }
165    return intersections.fUsed;
166}
167
168protected:
169
170void addEndPoints()
171{
172    for (int cIndex = 0; cIndex < 4; cIndex += 3) {
173        for (int lIndex = 0; lIndex < 2; lIndex++) {
174            if (cubic[cIndex] == line[lIndex]) {
175                intersections.insert(cIndex >> 1, lIndex);
176            }
177        }
178    }
179}
180
181void addHorizontalEndPoints(double left, double right, double y)
182{
183    for (int cIndex = 0; cIndex < 4; cIndex += 3) {
184        if (cubic[cIndex].y != y) {
185            continue;
186        }
187        if (cubic[cIndex].x == left) {
188            intersections.insert(cIndex >> 1, 0);
189        }
190        if (cubic[cIndex].x == right) {
191            intersections.insert(cIndex >> 1, 1);
192        }
193    }
194}
195
196void addVerticalEndPoints(double top, double bottom, double x)
197{
198    for (int cIndex = 0; cIndex < 4; cIndex += 3) {
199        if (cubic[cIndex].x != x) {
200            continue;
201        }
202        if (cubic[cIndex].y == top) {
203            intersections.insert(cIndex >> 1, 0);
204        }
205        if (cubic[cIndex].y == bottom) {
206            intersections.insert(cIndex >> 1, 1);
207        }
208    }
209}
210
211double findLineT(double t) {
212    double x, y;
213    xy_at_t(cubic, t, x, y);
214    double dx = line[1].x - line[0].x;
215    double dy = line[1].y - line[0].y;
216    if (fabs(dx) > fabs(dy)) {
217        return (x - line[0].x) / dx;
218    }
219    return (y - line[0].y) / dy;
220}
221
222void flip() {
223    // OPTIMIZATION: instead of swapping, pass original line, use [1].y - [0].y
224    int roots = intersections.fUsed;
225    for (int index = 0; index < roots; ++index) {
226        intersections.fT[1][index] = 1 - intersections.fT[1][index];
227    }
228}
229
230bool pinTs(double& cubicT, double& lineT) {
231    if (!approximately_one_or_less(lineT)) {
232        return false;
233    }
234    if (!approximately_zero_or_more(lineT)) {
235        return false;
236    }
237    if (cubicT < 0) {
238        cubicT = 0;
239    } else if (cubicT > 1) {
240        cubicT = 1;
241    }
242    if (lineT < 0) {
243        lineT = 0;
244    } else if (lineT > 1) {
245        lineT = 1;
246    }
247    return true;
248}
249
250private:
251
252const Cubic& cubic;
253const _Line& line;
254Intersections& intersections;
255};
256
257int horizontalIntersect(const Cubic& cubic, double left, double right, double y,
258        double tRange[3]) {
259    LineCubicIntersections c(cubic, *((_Line*) 0), *((Intersections*) 0));
260    double rootVals[3];
261    int result = c.horizontalIntersect(y, rootVals);
262    int tCount = 0;
263    for (int index = 0; index < result; ++index) {
264        double x, y;
265        xy_at_t(cubic, rootVals[index], x, y);
266        if (x < left || x > right) {
267            continue;
268        }
269        tRange[tCount++] = rootVals[index];
270    }
271    return result;
272}
273
274int horizontalIntersect(const Cubic& cubic, double left, double right, double y,
275        bool flipped, Intersections& intersections) {
276    LineCubicIntersections c(cubic, *((_Line*) 0), intersections);
277    return c.horizontalIntersect(y, left, right, flipped);
278}
279
280int verticalIntersect(const Cubic& cubic, double top, double bottom, double x,
281        bool flipped, Intersections& intersections) {
282    LineCubicIntersections c(cubic, *((_Line*) 0), intersections);
283    return c.verticalIntersect(x, top, bottom, flipped);
284}
285
286int intersect(const Cubic& cubic, const _Line& line, Intersections& i) {
287    LineCubicIntersections c(cubic, line, i);
288    return c.intersect();
289}
290