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