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 "include/core/SkPath.h"
8 #include "include/core/SkPoint.h"
9 #include "include/core/SkTypes.h"
10 #include "include/private/base/SkDebug.h"
11 #include "src/pathops/SkIntersections.h"
12 #include "src/pathops/SkPathOpsCubic.h"
13 #include "src/pathops/SkPathOpsCurve.h"
14 #include "src/pathops/SkPathOpsDebug.h"
15 #include "src/pathops/SkPathOpsLine.h"
16 #include "src/pathops/SkPathOpsPoint.h"
17 #include "src/pathops/SkPathOpsTypes.h"
18
19 #include <cmath>
20
21 /*
22 Find the intersection of a line and cubic by solving for valid t values.
23
24 Analogous to line-quadratic intersection, solve line-cubic intersection by
25 representing the cubic as:
26 x = a(1-t)^3 + 2b(1-t)^2t + c(1-t)t^2 + dt^3
27 y = e(1-t)^3 + 2f(1-t)^2t + g(1-t)t^2 + ht^3
28 and the line as:
29 y = i*x + j (if the line is more horizontal)
30 or:
31 x = i*y + j (if the line is more vertical)
32
33 Then using Mathematica, solve for the values of t where the cubic intersects the
34 line:
35
36 (in) Resultant[
37 a*(1 - t)^3 + 3*b*(1 - t)^2*t + 3*c*(1 - t)*t^2 + d*t^3 - x,
38 e*(1 - t)^3 + 3*f*(1 - t)^2*t + 3*g*(1 - t)*t^2 + h*t^3 - i*x - j, x]
39 (out) -e + j +
40 3 e t - 3 f t -
41 3 e t^2 + 6 f t^2 - 3 g t^2 +
42 e t^3 - 3 f t^3 + 3 g t^3 - h t^3 +
43 i ( a -
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
48 if i goes to infinity, we can rewrite the line in terms of x. Mathematica:
49
50 (in) Resultant[
51 a*(1 - t)^3 + 3*b*(1 - t)^2*t + 3*c*(1 - t)*t^2 + d*t^3 - i*y - j,
52 e*(1 - t)^3 + 3*f*(1 - t)^2*t + 3*g*(1 - t)*t^2 + h*t^3 - y, y]
53 (out) a - j -
54 3 a t + 3 b t +
55 3 a t^2 - 6 b t^2 + 3 c t^2 -
56 a t^3 + 3 b t^3 - 3 c t^3 + d t^3 -
57 i ( e -
58 3 e t + 3 f t +
59 3 e t^2 - 6 f t^2 + 3 g t^2 -
60 e t^3 + 3 f t^3 - 3 g t^3 + h t^3 )
61
62 Solving this with Mathematica produces an expression with hundreds of terms;
63 instead, use Numeric Solutions recipe to solve the cubic.
64
65 The near-horizontal case, in terms of: Ax^3 + Bx^2 + Cx + D == 0
66 A = (-(-e + 3*f - 3*g + h) + i*(-a + 3*b - 3*c + d) )
67 B = 3*(-( e - 2*f + g ) + i*( a - 2*b + c ) )
68 C = 3*(-(-e + f ) + i*(-a + b ) )
69 D = (-( e ) + i*( a ) + j )
70
71 The near-vertical case, in terms of: Ax^3 + Bx^2 + Cx + D == 0
72 A = ( (-a + 3*b - 3*c + d) - i*(-e + 3*f - 3*g + h) )
73 B = 3*( ( a - 2*b + c ) - i*( e - 2*f + g ) )
74 C = 3*( (-a + b ) - i*(-e + f ) )
75 D = ( ( a ) - i*( e ) - j )
76
77 For horizontal lines:
78 (in) Resultant[
79 a*(1 - t)^3 + 3*b*(1 - t)^2*t + 3*c*(1 - t)*t^2 + d*t^3 - j,
80 e*(1 - t)^3 + 3*f*(1 - t)^2*t + 3*g*(1 - t)*t^2 + h*t^3 - y, y]
81 (out) e - j -
82 3 e t + 3 f t +
83 3 e t^2 - 6 f t^2 + 3 g t^2 -
84 e t^3 + 3 f t^3 - 3 g t^3 + h t^3
85 */
86
87 class LineCubicIntersections {
88 public:
89 enum PinTPoint {
90 kPointUninitialized,
91 kPointInitialized
92 };
93
LineCubicIntersections(const SkDCubic & c,const SkDLine & l,SkIntersections * i)94 LineCubicIntersections(const SkDCubic& c, const SkDLine& l, SkIntersections* i)
95 : fCubic(c)
96 , fLine(l)
97 , fIntersections(i)
98 , fAllowNear(true) {
99 i->setMax(4);
100 }
101
allowNear(bool allow)102 void allowNear(bool allow) {
103 fAllowNear = allow;
104 }
105
checkCoincident()106 void checkCoincident() {
107 int last = fIntersections->used() - 1;
108 for (int index = 0; index < last; ) {
109 double cubicMidT = ((*fIntersections)[0][index] + (*fIntersections)[0][index + 1]) / 2;
110 SkDPoint cubicMidPt = fCubic.ptAtT(cubicMidT);
111 double t = fLine.nearPoint(cubicMidPt, nullptr);
112 if (t < 0) {
113 ++index;
114 continue;
115 }
116 if (fIntersections->isCoincident(index)) {
117 fIntersections->removeOne(index);
118 --last;
119 } else if (fIntersections->isCoincident(index + 1)) {
120 fIntersections->removeOne(index + 1);
121 --last;
122 } else {
123 fIntersections->setCoincident(index++);
124 }
125 fIntersections->setCoincident(index);
126 }
127 }
128
129 // see parallel routine in line quadratic intersections
intersectRay(double roots[3])130 int intersectRay(double roots[3]) {
131 double adj = fLine[1].fX - fLine[0].fX;
132 double opp = fLine[1].fY - fLine[0].fY;
133 SkDCubic c;
134 SkDEBUGCODE(c.fDebugGlobalState = fIntersections->globalState());
135 for (int n = 0; n < 4; ++n) {
136 c[n].fX = (fCubic[n].fY - fLine[0].fY) * adj - (fCubic[n].fX - fLine[0].fX) * opp;
137 }
138 double A, B, C, D;
139 SkDCubic::Coefficients(&c[0].fX, &A, &B, &C, &D);
140 int count = SkDCubic::RootsValidT(A, B, C, D, roots);
141 for (int index = 0; index < count; ++index) {
142 SkDPoint calcPt = c.ptAtT(roots[index]);
143 if (!approximately_zero(calcPt.fX)) {
144 for (int n = 0; n < 4; ++n) {
145 c[n].fY = (fCubic[n].fY - fLine[0].fY) * opp
146 + (fCubic[n].fX - fLine[0].fX) * adj;
147 }
148 double extremeTs[6];
149 int extrema = SkDCubic::FindExtrema(&c[0].fX, extremeTs);
150 count = c.searchRoots(extremeTs, extrema, 0, SkDCubic::kXAxis, roots);
151 break;
152 }
153 }
154 return count;
155 }
156
intersect()157 int intersect() {
158 addExactEndPoints();
159 if (fAllowNear) {
160 addNearEndPoints();
161 }
162 double rootVals[3];
163 int roots = intersectRay(rootVals);
164 for (int index = 0; index < roots; ++index) {
165 double cubicT = rootVals[index];
166 double lineT = findLineT(cubicT);
167 SkDPoint pt;
168 if (pinTs(&cubicT, &lineT, &pt, kPointUninitialized) && uniqueAnswer(cubicT, pt)) {
169 fIntersections->insert(cubicT, lineT, pt);
170 }
171 }
172 checkCoincident();
173 return fIntersections->used();
174 }
175
HorizontalIntersect(const SkDCubic & c,double axisIntercept,double roots[3])176 static int HorizontalIntersect(const SkDCubic& c, double axisIntercept, double roots[3]) {
177 double A, B, C, D;
178 SkDCubic::Coefficients(&c[0].fY, &A, &B, &C, &D);
179 D -= axisIntercept;
180 int count = SkDCubic::RootsValidT(A, B, C, D, roots);
181 for (int index = 0; index < count; ++index) {
182 SkDPoint calcPt = c.ptAtT(roots[index]);
183 if (!approximately_equal(calcPt.fY, axisIntercept)) {
184 double extremeTs[6];
185 int extrema = SkDCubic::FindExtrema(&c[0].fY, extremeTs);
186 count = c.searchRoots(extremeTs, extrema, axisIntercept, SkDCubic::kYAxis, roots);
187 break;
188 }
189 }
190 return count;
191 }
192
horizontalIntersect(double axisIntercept,double left,double right,bool flipped)193 int horizontalIntersect(double axisIntercept, double left, double right, bool flipped) {
194 addExactHorizontalEndPoints(left, right, axisIntercept);
195 if (fAllowNear) {
196 addNearHorizontalEndPoints(left, right, axisIntercept);
197 }
198 double roots[3];
199 int count = HorizontalIntersect(fCubic, axisIntercept, roots);
200 for (int index = 0; index < count; ++index) {
201 double cubicT = roots[index];
202 SkDPoint pt = { fCubic.ptAtT(cubicT).fX, axisIntercept };
203 double lineT = (pt.fX - left) / (right - left);
204 if (pinTs(&cubicT, &lineT, &pt, kPointInitialized) && uniqueAnswer(cubicT, pt)) {
205 fIntersections->insert(cubicT, lineT, pt);
206 }
207 }
208 if (flipped) {
209 fIntersections->flip();
210 }
211 checkCoincident();
212 return fIntersections->used();
213 }
214
uniqueAnswer(double cubicT,const SkDPoint & pt)215 bool uniqueAnswer(double cubicT, const SkDPoint& pt) {
216 for (int inner = 0; inner < fIntersections->used(); ++inner) {
217 if (fIntersections->pt(inner) != pt) {
218 continue;
219 }
220 double existingCubicT = (*fIntersections)[0][inner];
221 if (cubicT == existingCubicT) {
222 return false;
223 }
224 // check if midway on cubic is also same point. If so, discard this
225 double cubicMidT = (existingCubicT + cubicT) / 2;
226 SkDPoint cubicMidPt = fCubic.ptAtT(cubicMidT);
227 if (cubicMidPt.approximatelyEqual(pt)) {
228 return false;
229 }
230 }
231 #if ONE_OFF_DEBUG
232 SkDPoint cPt = fCubic.ptAtT(cubicT);
233 SkDebugf("%s pt=(%1.9g,%1.9g) cPt=(%1.9g,%1.9g)\n", __FUNCTION__, pt.fX, pt.fY,
234 cPt.fX, cPt.fY);
235 #endif
236 return true;
237 }
238
VerticalIntersect(const SkDCubic & c,double axisIntercept,double roots[3])239 static int VerticalIntersect(const SkDCubic& c, double axisIntercept, double roots[3]) {
240 double A, B, C, D;
241 SkDCubic::Coefficients(&c[0].fX, &A, &B, &C, &D);
242 D -= axisIntercept;
243 int count = SkDCubic::RootsValidT(A, B, C, D, roots);
244 for (int index = 0; index < count; ++index) {
245 SkDPoint calcPt = c.ptAtT(roots[index]);
246 if (!approximately_equal(calcPt.fX, axisIntercept)) {
247 double extremeTs[6];
248 int extrema = SkDCubic::FindExtrema(&c[0].fX, extremeTs);
249 count = c.searchRoots(extremeTs, extrema, axisIntercept, SkDCubic::kXAxis, roots);
250 break;
251 }
252 }
253 return count;
254 }
255
verticalIntersect(double axisIntercept,double top,double bottom,bool flipped)256 int verticalIntersect(double axisIntercept, double top, double bottom, bool flipped) {
257 addExactVerticalEndPoints(top, bottom, axisIntercept);
258 if (fAllowNear) {
259 addNearVerticalEndPoints(top, bottom, axisIntercept);
260 }
261 double roots[3];
262 int count = VerticalIntersect(fCubic, axisIntercept, roots);
263 for (int index = 0; index < count; ++index) {
264 double cubicT = roots[index];
265 SkDPoint pt = { axisIntercept, fCubic.ptAtT(cubicT).fY };
266 double lineT = (pt.fY - top) / (bottom - top);
267 if (pinTs(&cubicT, &lineT, &pt, kPointInitialized) && uniqueAnswer(cubicT, pt)) {
268 fIntersections->insert(cubicT, lineT, pt);
269 }
270 }
271 if (flipped) {
272 fIntersections->flip();
273 }
274 checkCoincident();
275 return fIntersections->used();
276 }
277
278 protected:
279
addExactEndPoints()280 void addExactEndPoints() {
281 for (int cIndex = 0; cIndex < 4; cIndex += 3) {
282 double lineT = fLine.exactPoint(fCubic[cIndex]);
283 if (lineT < 0) {
284 continue;
285 }
286 double cubicT = (double) (cIndex >> 1);
287 fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
288 }
289 }
290
291 /* Note that this does not look for endpoints of the line that are near the cubic.
292 These points are found later when check ends looks for missing points */
addNearEndPoints()293 void addNearEndPoints() {
294 for (int cIndex = 0; cIndex < 4; cIndex += 3) {
295 double cubicT = (double) (cIndex >> 1);
296 if (fIntersections->hasT(cubicT)) {
297 continue;
298 }
299 double lineT = fLine.nearPoint(fCubic[cIndex], nullptr);
300 if (lineT < 0) {
301 continue;
302 }
303 fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
304 }
305 this->addLineNearEndPoints();
306 }
307
addLineNearEndPoints()308 void addLineNearEndPoints() {
309 for (int lIndex = 0; lIndex < 2; ++lIndex) {
310 double lineT = (double) lIndex;
311 if (fIntersections->hasOppT(lineT)) {
312 continue;
313 }
314 double cubicT = ((SkDCurve*) &fCubic)->nearPoint(SkPath::kCubic_Verb,
315 fLine[lIndex], fLine[!lIndex]);
316 if (cubicT < 0) {
317 continue;
318 }
319 fIntersections->insert(cubicT, lineT, fLine[lIndex]);
320 }
321 }
322
addExactHorizontalEndPoints(double left,double right,double y)323 void addExactHorizontalEndPoints(double left, double right, double y) {
324 for (int cIndex = 0; cIndex < 4; cIndex += 3) {
325 double lineT = SkDLine::ExactPointH(fCubic[cIndex], left, right, y);
326 if (lineT < 0) {
327 continue;
328 }
329 double cubicT = (double) (cIndex >> 1);
330 fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
331 }
332 }
333
addNearHorizontalEndPoints(double left,double right,double y)334 void addNearHorizontalEndPoints(double left, double right, double y) {
335 for (int cIndex = 0; cIndex < 4; cIndex += 3) {
336 double cubicT = (double) (cIndex >> 1);
337 if (fIntersections->hasT(cubicT)) {
338 continue;
339 }
340 double lineT = SkDLine::NearPointH(fCubic[cIndex], left, right, y);
341 if (lineT < 0) {
342 continue;
343 }
344 fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
345 }
346 this->addLineNearEndPoints();
347 }
348
addExactVerticalEndPoints(double top,double bottom,double x)349 void addExactVerticalEndPoints(double top, double bottom, double x) {
350 for (int cIndex = 0; cIndex < 4; cIndex += 3) {
351 double lineT = SkDLine::ExactPointV(fCubic[cIndex], top, bottom, x);
352 if (lineT < 0) {
353 continue;
354 }
355 double cubicT = (double) (cIndex >> 1);
356 fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
357 }
358 }
359
addNearVerticalEndPoints(double top,double bottom,double x)360 void addNearVerticalEndPoints(double top, double bottom, double x) {
361 for (int cIndex = 0; cIndex < 4; cIndex += 3) {
362 double cubicT = (double) (cIndex >> 1);
363 if (fIntersections->hasT(cubicT)) {
364 continue;
365 }
366 double lineT = SkDLine::NearPointV(fCubic[cIndex], top, bottom, x);
367 if (lineT < 0) {
368 continue;
369 }
370 fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
371 }
372 this->addLineNearEndPoints();
373 }
374
findLineT(double t)375 double findLineT(double t) {
376 SkDPoint xy = fCubic.ptAtT(t);
377 double dx = fLine[1].fX - fLine[0].fX;
378 double dy = fLine[1].fY - fLine[0].fY;
379 if (fabs(dx) > fabs(dy)) {
380 return (xy.fX - fLine[0].fX) / dx;
381 }
382 return (xy.fY - fLine[0].fY) / dy;
383 }
384
pinTs(double * cubicT,double * lineT,SkDPoint * pt,PinTPoint ptSet)385 bool pinTs(double* cubicT, double* lineT, SkDPoint* pt, PinTPoint ptSet) {
386 if (!approximately_one_or_less(*lineT)) {
387 return false;
388 }
389 if (!approximately_zero_or_more(*lineT)) {
390 return false;
391 }
392 double cT = *cubicT = SkPinT(*cubicT);
393 double lT = *lineT = SkPinT(*lineT);
394 SkDPoint lPt = fLine.ptAtT(lT);
395 SkDPoint cPt = fCubic.ptAtT(cT);
396 if (!lPt.roughlyEqual(cPt)) {
397 return false;
398 }
399 // FIXME: if points are roughly equal but not approximately equal, need to do
400 // a binary search like quad/quad intersection to find more precise t values
401 if (lT == 0 || lT == 1 || (ptSet == kPointUninitialized && cT != 0 && cT != 1)) {
402 *pt = lPt;
403 } else if (ptSet == kPointUninitialized) {
404 *pt = cPt;
405 }
406 SkPoint gridPt = pt->asSkPoint();
407 if (gridPt == fLine[0].asSkPoint()) {
408 *lineT = 0;
409 } else if (gridPt == fLine[1].asSkPoint()) {
410 *lineT = 1;
411 }
412 if (gridPt == fCubic[0].asSkPoint() && approximately_equal(*cubicT, 0)) {
413 *cubicT = 0;
414 } else if (gridPt == fCubic[3].asSkPoint() && approximately_equal(*cubicT, 1)) {
415 *cubicT = 1;
416 }
417 return true;
418 }
419
420 private:
421 const SkDCubic& fCubic;
422 const SkDLine& fLine;
423 SkIntersections* fIntersections;
424 bool fAllowNear;
425 };
426
horizontal(const SkDCubic & cubic,double left,double right,double y,bool flipped)427 int SkIntersections::horizontal(const SkDCubic& cubic, double left, double right, double y,
428 bool flipped) {
429 SkDLine line = {{{ left, y }, { right, y }}};
430 LineCubicIntersections c(cubic, line, this);
431 return c.horizontalIntersect(y, left, right, flipped);
432 }
433
vertical(const SkDCubic & cubic,double top,double bottom,double x,bool flipped)434 int SkIntersections::vertical(const SkDCubic& cubic, double top, double bottom, double x,
435 bool flipped) {
436 SkDLine line = {{{ x, top }, { x, bottom }}};
437 LineCubicIntersections c(cubic, line, this);
438 return c.verticalIntersect(x, top, bottom, flipped);
439 }
440
intersect(const SkDCubic & cubic,const SkDLine & line)441 int SkIntersections::intersect(const SkDCubic& cubic, const SkDLine& line) {
442 LineCubicIntersections c(cubic, line, this);
443 c.allowNear(fAllowNear);
444 return c.intersect();
445 }
446
intersectRay(const SkDCubic & cubic,const SkDLine & line)447 int SkIntersections::intersectRay(const SkDCubic& cubic, const SkDLine& line) {
448 LineCubicIntersections c(cubic, line, this);
449 fUsed = c.intersectRay(fT[0]);
450 for (int index = 0; index < fUsed; ++index) {
451 fPt[index] = cubic.ptAtT(fT[0][index]);
452 }
453 return fUsed;
454 }
455
456 // SkDCubic accessors to Intersection utilities
457
horizontalIntersect(double yIntercept,double roots[3]) const458 int SkDCubic::horizontalIntersect(double yIntercept, double roots[3]) const {
459 return LineCubicIntersections::HorizontalIntersect(*this, yIntercept, roots);
460 }
461
verticalIntersect(double xIntercept,double roots[3]) const462 int SkDCubic::verticalIntersect(double xIntercept, double roots[3]) const {
463 return LineCubicIntersections::VerticalIntersect(*this, xIntercept, roots);
464 }
465