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
LineCubicIntersections(const SkDCubic & c,const SkDLine & l,SkIntersections * i)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
allowNear(bool allow)92 void allowNear(bool allow) {
93 fAllowNear = allow;
94 }
95
96 // see parallel routine in line quadratic intersections
intersectRay(double roots[3])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 c;
101 for (int n = 0; n < 4; ++n) {
102 c[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(&c[0].fX, &A, &B, &C, &D);
106 int count = SkDCubic::RootsValidT(A, B, C, D, roots);
107 for (int index = 0; index < count; ++index) {
108 SkDPoint calcPt = c.ptAtT(roots[index]);
109 if (!approximately_zero(calcPt.fX)) {
110 for (int n = 0; n < 4; ++n) {
111 c[n].fY = (fCubic[n].fY - fLine[0].fY) * opp
112 + (fCubic[n].fX - fLine[0].fX) * adj;
113 }
114 double extremeTs[6];
115 int extrema = SkDCubic::FindExtrema(c[0].fX, c[1].fX, c[2].fX, c[3].fX, extremeTs);
116 count = c.searchRoots(extremeTs, extrema, 0, SkDCubic::kXAxis, roots);
117 break;
118 }
119 }
120 return count;
121 }
122
intersect()123 int intersect() {
124 addExactEndPoints();
125 if (fAllowNear) {
126 addNearEndPoints();
127 }
128 double rootVals[3];
129 int roots = intersectRay(rootVals);
130 for (int index = 0; index < roots; ++index) {
131 double cubicT = rootVals[index];
132 double lineT = findLineT(cubicT);
133 SkDPoint pt;
134 if (pinTs(&cubicT, &lineT, &pt, kPointUninitialized)) {
135 #if ONE_OFF_DEBUG
136 SkDPoint cPt = fCubic.ptAtT(cubicT);
137 SkDebugf("%s pt=(%1.9g,%1.9g) cPt=(%1.9g,%1.9g)\n", __FUNCTION__, pt.fX, pt.fY,
138 cPt.fX, cPt.fY);
139 #endif
140 for (int inner = 0; inner < fIntersections->used(); ++inner) {
141 if (fIntersections->pt(inner) != pt) {
142 continue;
143 }
144 double existingCubicT = (*fIntersections)[0][inner];
145 if (cubicT == existingCubicT) {
146 goto skipInsert;
147 }
148 // check if midway on cubic is also same point. If so, discard this
149 double cubicMidT = (existingCubicT + cubicT) / 2;
150 SkDPoint cubicMidPt = fCubic.ptAtT(cubicMidT);
151 if (cubicMidPt.approximatelyEqual(pt)) {
152 goto skipInsert;
153 }
154 }
155 fIntersections->insert(cubicT, lineT, pt);
156 skipInsert:
157 ;
158 }
159 }
160 return fIntersections->used();
161 }
162
HorizontalIntersect(const SkDCubic & c,double axisIntercept,double roots[3])163 static int HorizontalIntersect(const SkDCubic& c, double axisIntercept, double roots[3]) {
164 double A, B, C, D;
165 SkDCubic::Coefficients(&c[0].fY, &A, &B, &C, &D);
166 D -= axisIntercept;
167 int count = SkDCubic::RootsValidT(A, B, C, D, roots);
168 for (int index = 0; index < count; ++index) {
169 SkDPoint calcPt = c.ptAtT(roots[index]);
170 if (!approximately_equal(calcPt.fY, axisIntercept)) {
171 double extremeTs[6];
172 int extrema = SkDCubic::FindExtrema(c[0].fY, c[1].fY, c[2].fY, c[3].fY, extremeTs);
173 count = c.searchRoots(extremeTs, extrema, axisIntercept, SkDCubic::kYAxis, roots);
174 break;
175 }
176 }
177 return count;
178 }
179
horizontalIntersect(double axisIntercept,double left,double right,bool flipped)180 int horizontalIntersect(double axisIntercept, double left, double right, bool flipped) {
181 addExactHorizontalEndPoints(left, right, axisIntercept);
182 if (fAllowNear) {
183 addNearHorizontalEndPoints(left, right, axisIntercept);
184 }
185 double roots[3];
186 int count = HorizontalIntersect(fCubic, axisIntercept, roots);
187 for (int index = 0; index < count; ++index) {
188 double cubicT = roots[index];
189 SkDPoint pt;
190 pt.fX = fCubic.ptAtT(cubicT).fX;
191 pt.fY = axisIntercept;
192 double lineT = (pt.fX - left) / (right - left);
193 if (pinTs(&cubicT, &lineT, &pt, kPointInitialized)) {
194 fIntersections->insert(cubicT, lineT, pt);
195 }
196 }
197 if (flipped) {
198 fIntersections->flip();
199 }
200 return fIntersections->used();
201 }
202
VerticalIntersect(const SkDCubic & c,double axisIntercept,double roots[3])203 static int VerticalIntersect(const SkDCubic& c, double axisIntercept, double roots[3]) {
204 double A, B, C, D;
205 SkDCubic::Coefficients(&c[0].fX, &A, &B, &C, &D);
206 D -= axisIntercept;
207 int count = SkDCubic::RootsValidT(A, B, C, D, roots);
208 for (int index = 0; index < count; ++index) {
209 SkDPoint calcPt = c.ptAtT(roots[index]);
210 if (!approximately_equal(calcPt.fX, axisIntercept)) {
211 double extremeTs[6];
212 int extrema = SkDCubic::FindExtrema(c[0].fX, c[1].fX, c[2].fX, c[3].fX, extremeTs);
213 count = c.searchRoots(extremeTs, extrema, axisIntercept, SkDCubic::kXAxis, roots);
214 break;
215 }
216 }
217 return count;
218 }
219
verticalIntersect(double axisIntercept,double top,double bottom,bool flipped)220 int verticalIntersect(double axisIntercept, double top, double bottom, bool flipped) {
221 addExactVerticalEndPoints(top, bottom, axisIntercept);
222 if (fAllowNear) {
223 addNearVerticalEndPoints(top, bottom, axisIntercept);
224 }
225 double roots[3];
226 int count = VerticalIntersect(fCubic, axisIntercept, roots);
227 for (int index = 0; index < count; ++index) {
228 double cubicT = roots[index];
229 SkDPoint pt;
230 pt.fX = axisIntercept;
231 pt.fY = fCubic.ptAtT(cubicT).fY;
232 double lineT = (pt.fY - top) / (bottom - top);
233 if (pinTs(&cubicT, &lineT, &pt, kPointInitialized)) {
234 fIntersections->insert(cubicT, lineT, pt);
235 }
236 }
237 if (flipped) {
238 fIntersections->flip();
239 }
240 return fIntersections->used();
241 }
242
243 protected:
244
addExactEndPoints()245 void addExactEndPoints() {
246 for (int cIndex = 0; cIndex < 4; cIndex += 3) {
247 double lineT = fLine.exactPoint(fCubic[cIndex]);
248 if (lineT < 0) {
249 continue;
250 }
251 double cubicT = (double) (cIndex >> 1);
252 fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
253 }
254 }
255
256 /* Note that this does not look for endpoints of the line that are near the cubic.
257 These points are found later when check ends looks for missing points */
addNearEndPoints()258 void addNearEndPoints() {
259 for (int cIndex = 0; cIndex < 4; cIndex += 3) {
260 double cubicT = (double) (cIndex >> 1);
261 if (fIntersections->hasT(cubicT)) {
262 continue;
263 }
264 double lineT = fLine.nearPoint(fCubic[cIndex], NULL);
265 if (lineT < 0) {
266 continue;
267 }
268 fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
269 }
270 }
271
addExactHorizontalEndPoints(double left,double right,double y)272 void addExactHorizontalEndPoints(double left, double right, double y) {
273 for (int cIndex = 0; cIndex < 4; cIndex += 3) {
274 double lineT = SkDLine::ExactPointH(fCubic[cIndex], left, right, y);
275 if (lineT < 0) {
276 continue;
277 }
278 double cubicT = (double) (cIndex >> 1);
279 fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
280 }
281 }
282
addNearHorizontalEndPoints(double left,double right,double y)283 void addNearHorizontalEndPoints(double left, double right, double y) {
284 for (int cIndex = 0; cIndex < 4; cIndex += 3) {
285 double cubicT = (double) (cIndex >> 1);
286 if (fIntersections->hasT(cubicT)) {
287 continue;
288 }
289 double lineT = SkDLine::NearPointH(fCubic[cIndex], left, right, y);
290 if (lineT < 0) {
291 continue;
292 }
293 fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
294 }
295 // FIXME: see if line end is nearly on cubic
296 }
297
addExactVerticalEndPoints(double top,double bottom,double x)298 void addExactVerticalEndPoints(double top, double bottom, double x) {
299 for (int cIndex = 0; cIndex < 4; cIndex += 3) {
300 double lineT = SkDLine::ExactPointV(fCubic[cIndex], top, bottom, x);
301 if (lineT < 0) {
302 continue;
303 }
304 double cubicT = (double) (cIndex >> 1);
305 fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
306 }
307 }
308
addNearVerticalEndPoints(double top,double bottom,double x)309 void addNearVerticalEndPoints(double top, double bottom, double x) {
310 for (int cIndex = 0; cIndex < 4; cIndex += 3) {
311 double cubicT = (double) (cIndex >> 1);
312 if (fIntersections->hasT(cubicT)) {
313 continue;
314 }
315 double lineT = SkDLine::NearPointV(fCubic[cIndex], top, bottom, x);
316 if (lineT < 0) {
317 continue;
318 }
319 fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
320 }
321 // FIXME: see if line end is nearly on cubic
322 }
323
findLineT(double t)324 double findLineT(double t) {
325 SkDPoint xy = fCubic.ptAtT(t);
326 double dx = fLine[1].fX - fLine[0].fX;
327 double dy = fLine[1].fY - fLine[0].fY;
328 if (fabs(dx) > fabs(dy)) {
329 return (xy.fX - fLine[0].fX) / dx;
330 }
331 return (xy.fY - fLine[0].fY) / dy;
332 }
333
pinTs(double * cubicT,double * lineT,SkDPoint * pt,PinTPoint ptSet)334 bool pinTs(double* cubicT, double* lineT, SkDPoint* pt, PinTPoint ptSet) {
335 if (!approximately_one_or_less(*lineT)) {
336 return false;
337 }
338 if (!approximately_zero_or_more(*lineT)) {
339 return false;
340 }
341 double cT = *cubicT = SkPinT(*cubicT);
342 double lT = *lineT = SkPinT(*lineT);
343 SkDPoint lPt = fLine.ptAtT(lT);
344 SkDPoint cPt = fCubic.ptAtT(cT);
345 if (!lPt.moreRoughlyEqual(cPt)) {
346 return false;
347 }
348 // FIXME: if points are roughly equal but not approximately equal, need to do
349 // a binary search like quad/quad intersection to find more precise t values
350 if (lT == 0 || lT == 1 || (ptSet == kPointUninitialized && cT != 0 && cT != 1)) {
351 *pt = lPt;
352 } else if (ptSet == kPointUninitialized) {
353 *pt = cPt;
354 }
355 SkPoint gridPt = pt->asSkPoint();
356 if (gridPt == fLine[0].asSkPoint()) {
357 *lineT = 0;
358 } else if (gridPt == fLine[1].asSkPoint()) {
359 *lineT = 1;
360 }
361 if (gridPt == fCubic[0].asSkPoint() && approximately_equal(*cubicT, 0)) {
362 *cubicT = 0;
363 } else if (gridPt == fCubic[3].asSkPoint() && approximately_equal(*cubicT, 1)) {
364 *cubicT = 1;
365 }
366 return true;
367 }
368
369 private:
370 const SkDCubic& fCubic;
371 const SkDLine& fLine;
372 SkIntersections* fIntersections;
373 bool fAllowNear;
374 };
375
horizontal(const SkDCubic & cubic,double left,double right,double y,bool flipped)376 int SkIntersections::horizontal(const SkDCubic& cubic, double left, double right, double y,
377 bool flipped) {
378 SkDLine line = {{{ left, y }, { right, y }}};
379 LineCubicIntersections c(cubic, line, this);
380 return c.horizontalIntersect(y, left, right, flipped);
381 }
382
vertical(const SkDCubic & cubic,double top,double bottom,double x,bool flipped)383 int SkIntersections::vertical(const SkDCubic& cubic, double top, double bottom, double x,
384 bool flipped) {
385 SkDLine line = {{{ x, top }, { x, bottom }}};
386 LineCubicIntersections c(cubic, line, this);
387 return c.verticalIntersect(x, top, bottom, flipped);
388 }
389
intersect(const SkDCubic & cubic,const SkDLine & line)390 int SkIntersections::intersect(const SkDCubic& cubic, const SkDLine& line) {
391 LineCubicIntersections c(cubic, line, this);
392 c.allowNear(fAllowNear);
393 return c.intersect();
394 }
395
intersectRay(const SkDCubic & cubic,const SkDLine & line)396 int SkIntersections::intersectRay(const SkDCubic& cubic, const SkDLine& line) {
397 LineCubicIntersections c(cubic, line, this);
398 fUsed = c.intersectRay(fT[0]);
399 for (int index = 0; index < fUsed; ++index) {
400 fPt[index] = cubic.ptAtT(fT[0][index]);
401 }
402 return fUsed;
403 }
404