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