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