|
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 "SkIntersections.h" |
|
8 #include "SkPathOpsCubic.h" |
|
9 #include "SkPathOpsLine.h" |
|
10 |
|
11 /* |
|
12 Find the interection of a line and cubic by solving for valid t values. |
|
13 |
|
14 Analogous to line-quadratic intersection, solve line-cubic intersection by |
|
15 representing 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 |
|
18 and the line as: |
|
19 y = i*x + j (if the line is more horizontal) |
|
20 or: |
|
21 x = i*y + j (if the line is more vertical) |
|
22 |
|
23 Then using Mathematica, solve for the values of t where the cubic intersects the |
|
24 line: |
|
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 |
|
38 if 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 |
|
52 Solving this with Mathematica produces an expression with hundreds of terms; |
|
53 instead, use Numeric Solutions recipe to solve the cubic. |
|
54 |
|
55 The 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 |
|
61 The 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 |
|
67 For 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 |
|
77 class LineCubicIntersections { |
|
78 public: |
|
79 enum PinTPoint { |
|
80 kPointUninitialized, |
|
81 kPointInitialized |
|
82 }; |
|
83 |
|
84 LineCubicIntersections(const SkDCubic& c, const SkDLine& l, SkIntersections* i) |
|
85 : fCubic(c) |
|
86 , fLine(l) |
|
87 , fIntersections(i) |
|
88 , fAllowNear(true) { |
|
89 i->setMax(3); |
|
90 } |
|
91 |
|
92 void allowNear(bool allow) { |
|
93 fAllowNear = allow; |
|
94 } |
|
95 |
|
96 // see parallel routine in line quadratic intersections |
|
97 int intersectRay(double roots[3]) { |
|
98 double adj = fLine[1].fX - fLine[0].fX; |
|
99 double opp = fLine[1].fY - fLine[0].fY; |
|
100 SkDCubic r; |
|
101 for (int n = 0; n < 4; ++n) { |
|
102 r[n].fX = (fCubic[n].fY - fLine[0].fY) * adj - (fCubic[n].fX - fLine[0].fX) * opp; |
|
103 } |
|
104 double A, B, C, D; |
|
105 SkDCubic::Coefficients(&r[0].fX, &A, &B, &C, &D); |
|
106 return SkDCubic::RootsValidT(A, B, C, D, roots); |
|
107 } |
|
108 |
|
109 int intersect() { |
|
110 addExactEndPoints(); |
|
111 if (fAllowNear) { |
|
112 addNearEndPoints(); |
|
113 } |
|
114 double rootVals[3]; |
|
115 int roots = intersectRay(rootVals); |
|
116 for (int index = 0; index < roots; ++index) { |
|
117 double cubicT = rootVals[index]; |
|
118 double lineT = findLineT(cubicT); |
|
119 SkDPoint pt; |
|
120 if (pinTs(&cubicT, &lineT, &pt, kPointUninitialized)) { |
|
121 #if ONE_OFF_DEBUG |
|
122 SkDPoint cPt = fCubic.ptAtT(cubicT); |
|
123 SkDebugf("%s pt=(%1.9g,%1.9g) cPt=(%1.9g,%1.9g)\n", __FUNCTION__, pt.fX, pt.fY, |
|
124 cPt.fX, cPt.fY); |
|
125 #endif |
|
126 for (int inner = 0; inner < fIntersections->used(); ++inner) { |
|
127 if (fIntersections->pt(inner) != pt) { |
|
128 continue; |
|
129 } |
|
130 double existingCubicT = (*fIntersections)[0][inner]; |
|
131 if (cubicT == existingCubicT) { |
|
132 goto skipInsert; |
|
133 } |
|
134 // check if midway on cubic is also same point. If so, discard this |
|
135 double cubicMidT = (existingCubicT + cubicT) / 2; |
|
136 SkDPoint cubicMidPt = fCubic.ptAtT(cubicMidT); |
|
137 if (cubicMidPt.approximatelyEqual(pt)) { |
|
138 goto skipInsert; |
|
139 } |
|
140 } |
|
141 fIntersections->insert(cubicT, lineT, pt); |
|
142 skipInsert: |
|
143 ; |
|
144 } |
|
145 } |
|
146 return fIntersections->used(); |
|
147 } |
|
148 |
|
149 int horizontalIntersect(double axisIntercept, double roots[3]) { |
|
150 double A, B, C, D; |
|
151 SkDCubic::Coefficients(&fCubic[0].fY, &A, &B, &C, &D); |
|
152 D -= axisIntercept; |
|
153 return SkDCubic::RootsValidT(A, B, C, D, roots); |
|
154 } |
|
155 |
|
156 int horizontalIntersect(double axisIntercept, double left, double right, bool flipped) { |
|
157 addExactHorizontalEndPoints(left, right, axisIntercept); |
|
158 if (fAllowNear) { |
|
159 addNearHorizontalEndPoints(left, right, axisIntercept); |
|
160 } |
|
161 double rootVals[3]; |
|
162 int roots = horizontalIntersect(axisIntercept, rootVals); |
|
163 for (int index = 0; index < roots; ++index) { |
|
164 double cubicT = rootVals[index]; |
|
165 SkDPoint pt = fCubic.ptAtT(cubicT); |
|
166 double lineT = (pt.fX - left) / (right - left); |
|
167 if (pinTs(&cubicT, &lineT, &pt, kPointInitialized)) { |
|
168 fIntersections->insert(cubicT, lineT, pt); |
|
169 } |
|
170 } |
|
171 if (flipped) { |
|
172 fIntersections->flip(); |
|
173 } |
|
174 return fIntersections->used(); |
|
175 } |
|
176 |
|
177 int verticalIntersect(double axisIntercept, double roots[3]) { |
|
178 double A, B, C, D; |
|
179 SkDCubic::Coefficients(&fCubic[0].fX, &A, &B, &C, &D); |
|
180 D -= axisIntercept; |
|
181 return SkDCubic::RootsValidT(A, B, C, D, roots); |
|
182 } |
|
183 |
|
184 int verticalIntersect(double axisIntercept, double top, double bottom, bool flipped) { |
|
185 addExactVerticalEndPoints(top, bottom, axisIntercept); |
|
186 if (fAllowNear) { |
|
187 addNearVerticalEndPoints(top, bottom, axisIntercept); |
|
188 } |
|
189 double rootVals[3]; |
|
190 int roots = verticalIntersect(axisIntercept, rootVals); |
|
191 for (int index = 0; index < roots; ++index) { |
|
192 double cubicT = rootVals[index]; |
|
193 SkDPoint pt = fCubic.ptAtT(cubicT); |
|
194 double lineT = (pt.fY - top) / (bottom - top); |
|
195 if (pinTs(&cubicT, &lineT, &pt, kPointInitialized)) { |
|
196 fIntersections->insert(cubicT, lineT, pt); |
|
197 } |
|
198 } |
|
199 if (flipped) { |
|
200 fIntersections->flip(); |
|
201 } |
|
202 return fIntersections->used(); |
|
203 } |
|
204 |
|
205 protected: |
|
206 |
|
207 void addExactEndPoints() { |
|
208 for (int cIndex = 0; cIndex < 4; cIndex += 3) { |
|
209 double lineT = fLine.exactPoint(fCubic[cIndex]); |
|
210 if (lineT < 0) { |
|
211 continue; |
|
212 } |
|
213 double cubicT = (double) (cIndex >> 1); |
|
214 fIntersections->insert(cubicT, lineT, fCubic[cIndex]); |
|
215 } |
|
216 } |
|
217 |
|
218 /* Note that this does not look for endpoints of the line that are near the cubic. |
|
219 These points are found later when check ends looks for missing points */ |
|
220 void addNearEndPoints() { |
|
221 for (int cIndex = 0; cIndex < 4; cIndex += 3) { |
|
222 double cubicT = (double) (cIndex >> 1); |
|
223 if (fIntersections->hasT(cubicT)) { |
|
224 continue; |
|
225 } |
|
226 double lineT = fLine.nearPoint(fCubic[cIndex]); |
|
227 if (lineT < 0) { |
|
228 continue; |
|
229 } |
|
230 fIntersections->insert(cubicT, lineT, fCubic[cIndex]); |
|
231 } |
|
232 } |
|
233 |
|
234 void addExactHorizontalEndPoints(double left, double right, double y) { |
|
235 for (int cIndex = 0; cIndex < 4; cIndex += 3) { |
|
236 double lineT = SkDLine::ExactPointH(fCubic[cIndex], left, right, y); |
|
237 if (lineT < 0) { |
|
238 continue; |
|
239 } |
|
240 double cubicT = (double) (cIndex >> 1); |
|
241 fIntersections->insert(cubicT, lineT, fCubic[cIndex]); |
|
242 } |
|
243 } |
|
244 |
|
245 void addNearHorizontalEndPoints(double left, double right, double y) { |
|
246 for (int cIndex = 0; cIndex < 4; cIndex += 3) { |
|
247 double cubicT = (double) (cIndex >> 1); |
|
248 if (fIntersections->hasT(cubicT)) { |
|
249 continue; |
|
250 } |
|
251 double lineT = SkDLine::NearPointH(fCubic[cIndex], left, right, y); |
|
252 if (lineT < 0) { |
|
253 continue; |
|
254 } |
|
255 fIntersections->insert(cubicT, lineT, fCubic[cIndex]); |
|
256 } |
|
257 // FIXME: see if line end is nearly on cubic |
|
258 } |
|
259 |
|
260 void addExactVerticalEndPoints(double top, double bottom, double x) { |
|
261 for (int cIndex = 0; cIndex < 4; cIndex += 3) { |
|
262 double lineT = SkDLine::ExactPointV(fCubic[cIndex], top, bottom, x); |
|
263 if (lineT < 0) { |
|
264 continue; |
|
265 } |
|
266 double cubicT = (double) (cIndex >> 1); |
|
267 fIntersections->insert(cubicT, lineT, fCubic[cIndex]); |
|
268 } |
|
269 } |
|
270 |
|
271 void addNearVerticalEndPoints(double top, double bottom, double x) { |
|
272 for (int cIndex = 0; cIndex < 4; cIndex += 3) { |
|
273 double cubicT = (double) (cIndex >> 1); |
|
274 if (fIntersections->hasT(cubicT)) { |
|
275 continue; |
|
276 } |
|
277 double lineT = SkDLine::NearPointV(fCubic[cIndex], top, bottom, x); |
|
278 if (lineT < 0) { |
|
279 continue; |
|
280 } |
|
281 fIntersections->insert(cubicT, lineT, fCubic[cIndex]); |
|
282 } |
|
283 // FIXME: see if line end is nearly on cubic |
|
284 } |
|
285 |
|
286 double findLineT(double t) { |
|
287 SkDPoint xy = fCubic.ptAtT(t); |
|
288 double dx = fLine[1].fX - fLine[0].fX; |
|
289 double dy = fLine[1].fY - fLine[0].fY; |
|
290 if (fabs(dx) > fabs(dy)) { |
|
291 return (xy.fX - fLine[0].fX) / dx; |
|
292 } |
|
293 return (xy.fY - fLine[0].fY) / dy; |
|
294 } |
|
295 |
|
296 bool pinTs(double* cubicT, double* lineT, SkDPoint* pt, PinTPoint ptSet) { |
|
297 if (!approximately_one_or_less(*lineT)) { |
|
298 return false; |
|
299 } |
|
300 if (!approximately_zero_or_more(*lineT)) { |
|
301 return false; |
|
302 } |
|
303 double cT = *cubicT = SkPinT(*cubicT); |
|
304 double lT = *lineT = SkPinT(*lineT); |
|
305 if (lT == 0 || lT == 1 || (ptSet == kPointUninitialized && cT != 0 && cT != 1)) { |
|
306 *pt = fLine.ptAtT(lT); |
|
307 } else if (ptSet == kPointUninitialized) { |
|
308 *pt = fCubic.ptAtT(cT); |
|
309 } |
|
310 SkPoint gridPt = pt->asSkPoint(); |
|
311 if (gridPt == fLine[0].asSkPoint()) { |
|
312 *lineT = 0; |
|
313 } else if (gridPt == fLine[1].asSkPoint()) { |
|
314 *lineT = 1; |
|
315 } |
|
316 if (gridPt == fCubic[0].asSkPoint() && approximately_equal(*cubicT, 0)) { |
|
317 *cubicT = 0; |
|
318 } else if (gridPt == fCubic[3].asSkPoint() && approximately_equal(*cubicT, 1)) { |
|
319 *cubicT = 1; |
|
320 } |
|
321 return true; |
|
322 } |
|
323 |
|
324 private: |
|
325 const SkDCubic& fCubic; |
|
326 const SkDLine& fLine; |
|
327 SkIntersections* fIntersections; |
|
328 bool fAllowNear; |
|
329 }; |
|
330 |
|
331 int SkIntersections::horizontal(const SkDCubic& cubic, double left, double right, double y, |
|
332 bool flipped) { |
|
333 SkDLine line = {{{ left, y }, { right, y }}}; |
|
334 LineCubicIntersections c(cubic, line, this); |
|
335 return c.horizontalIntersect(y, left, right, flipped); |
|
336 } |
|
337 |
|
338 int SkIntersections::vertical(const SkDCubic& cubic, double top, double bottom, double x, |
|
339 bool flipped) { |
|
340 SkDLine line = {{{ x, top }, { x, bottom }}}; |
|
341 LineCubicIntersections c(cubic, line, this); |
|
342 return c.verticalIntersect(x, top, bottom, flipped); |
|
343 } |
|
344 |
|
345 int SkIntersections::intersect(const SkDCubic& cubic, const SkDLine& line) { |
|
346 LineCubicIntersections c(cubic, line, this); |
|
347 c.allowNear(fAllowNear); |
|
348 return c.intersect(); |
|
349 } |
|
350 |
|
351 int SkIntersections::intersectRay(const SkDCubic& cubic, const SkDLine& line) { |
|
352 LineCubicIntersections c(cubic, line, this); |
|
353 fUsed = c.intersectRay(fT[0]); |
|
354 for (int index = 0; index < fUsed; ++index) { |
|
355 fPt[index] = cubic.ptAtT(fT[0][index]); |
|
356 } |
|
357 return fUsed; |
|
358 } |