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/SkScalar.h"
10 #include "src/pathops/SkIntersections.h"
11 #include "src/pathops/SkPathOpsCurve.h"
12 #include "src/pathops/SkPathOpsDebug.h"
13 #include "src/pathops/SkPathOpsLine.h"
14 #include "src/pathops/SkPathOpsPoint.h"
15 #include "src/pathops/SkPathOpsQuad.h"
16 #include "src/pathops/SkPathOpsTypes.h"
17
18 #include <cmath>
19
20 /*
21 Find the intersection of a line and quadratic by solving for valid t values.
22
23 From http://stackoverflow.com/questions/1853637/how-to-find-the-mathematical-function-defining-a-bezier-curve
24
25 "A Bezier curve is a parametric function. A quadratic Bezier curve (i.e. three
26 control points) can be expressed as: F(t) = A(1 - t)^2 + B(1 - t)t + Ct^2 where
27 A, B and C are points and t goes from zero to one.
28
29 This will give you two equations:
30
31 x = a(1 - t)^2 + b(1 - t)t + ct^2
32 y = d(1 - t)^2 + e(1 - t)t + ft^2
33
34 If you add for instance the line equation (y = kx + m) to that, you'll end up
35 with three equations and three unknowns (x, y and t)."
36
37 Similar to above, the quadratic is represented as
38 x = a(1-t)^2 + 2b(1-t)t + ct^2
39 y = d(1-t)^2 + 2e(1-t)t + ft^2
40 and the line as
41 y = g*x + h
42
43 Using Mathematica, solve for the values of t where the quadratic intersects the
44 line:
45
46 (in) t1 = Resultant[a*(1 - t)^2 + 2*b*(1 - t)*t + c*t^2 - x,
47 d*(1 - t)^2 + 2*e*(1 - t)*t + f*t^2 - g*x - h, x]
48 (out) -d + h + 2 d t - 2 e t - d t^2 + 2 e t^2 - f t^2 +
49 g (a - 2 a t + 2 b t + a t^2 - 2 b t^2 + c t^2)
50 (in) Solve[t1 == 0, t]
51 (out) {
52 {t -> (-2 d + 2 e + 2 a g - 2 b g -
53 Sqrt[(2 d - 2 e - 2 a g + 2 b g)^2 -
54 4 (-d + 2 e - f + a g - 2 b g + c g) (-d + a g + h)]) /
55 (2 (-d + 2 e - f + a g - 2 b g + c g))
56 },
57 {t -> (-2 d + 2 e + 2 a g - 2 b g +
58 Sqrt[(2 d - 2 e - 2 a g + 2 b g)^2 -
59 4 (-d + 2 e - f + a g - 2 b g + c g) (-d + a g + h)]) /
60 (2 (-d + 2 e - f + a g - 2 b g + c g))
61 }
62 }
63
64 Using the results above (when the line tends towards horizontal)
65 A = (-(d - 2*e + f) + g*(a - 2*b + c) )
66 B = 2*( (d - e ) - g*(a - b ) )
67 C = (-(d ) + g*(a ) + h )
68
69 If g goes to infinity, we can rewrite the line in terms of x.
70 x = g'*y + h'
71
72 And solve accordingly in Mathematica:
73
74 (in) t2 = Resultant[a*(1 - t)^2 + 2*b*(1 - t)*t + c*t^2 - g'*y - h',
75 d*(1 - t)^2 + 2*e*(1 - t)*t + f*t^2 - y, y]
76 (out) a - h' - 2 a t + 2 b t + a t^2 - 2 b t^2 + c t^2 -
77 g' (d - 2 d t + 2 e t + d t^2 - 2 e t^2 + f t^2)
78 (in) Solve[t2 == 0, t]
79 (out) {
80 {t -> (2 a - 2 b - 2 d g' + 2 e g' -
81 Sqrt[(-2 a + 2 b + 2 d g' - 2 e g')^2 -
82 4 (a - 2 b + c - d g' + 2 e g' - f g') (a - d g' - h')]) /
83 (2 (a - 2 b + c - d g' + 2 e g' - f g'))
84 },
85 {t -> (2 a - 2 b - 2 d g' + 2 e g' +
86 Sqrt[(-2 a + 2 b + 2 d g' - 2 e g')^2 -
87 4 (a - 2 b + c - d g' + 2 e g' - f g') (a - d g' - h')])/
88 (2 (a - 2 b + c - d g' + 2 e g' - f g'))
89 }
90 }
91
92 Thus, if the slope of the line tends towards vertical, we use:
93 A = ( (a - 2*b + c) - g'*(d - 2*e + f) )
94 B = 2*(-(a - b ) + g'*(d - e ) )
95 C = ( (a ) - g'*(d ) - h' )
96 */
97
98 class LineQuadraticIntersections {
99 public:
100 enum PinTPoint {
101 kPointUninitialized,
102 kPointInitialized
103 };
104
LineQuadraticIntersections(const SkDQuad & q,const SkDLine & l,SkIntersections * i)105 LineQuadraticIntersections(const SkDQuad& q, const SkDLine& l, SkIntersections* i)
106 : fQuad(q)
107 , fLine(&l)
108 , fIntersections(i)
109 , fAllowNear(true) {
110 i->setMax(5); // allow short partial coincidence plus discrete intersections
111 }
112
LineQuadraticIntersections(const SkDQuad & q)113 LineQuadraticIntersections(const SkDQuad& q)
114 : fQuad(q)
115 SkDEBUGPARAMS(fLine(nullptr))
116 SkDEBUGPARAMS(fIntersections(nullptr))
117 SkDEBUGPARAMS(fAllowNear(false)) {
118 }
119
allowNear(bool allow)120 void allowNear(bool allow) {
121 fAllowNear = allow;
122 }
123
checkCoincident()124 void checkCoincident() {
125 int last = fIntersections->used() - 1;
126 for (int index = 0; index < last; ) {
127 double quadMidT = ((*fIntersections)[0][index] + (*fIntersections)[0][index + 1]) / 2;
128 SkDPoint quadMidPt = fQuad.ptAtT(quadMidT);
129 double t = fLine->nearPoint(quadMidPt, nullptr);
130 if (t < 0) {
131 ++index;
132 continue;
133 }
134 if (fIntersections->isCoincident(index)) {
135 fIntersections->removeOne(index);
136 --last;
137 } else if (fIntersections->isCoincident(index + 1)) {
138 fIntersections->removeOne(index + 1);
139 --last;
140 } else {
141 fIntersections->setCoincident(index++);
142 }
143 fIntersections->setCoincident(index);
144 }
145 }
146
intersectRay(double roots[2])147 int intersectRay(double roots[2]) {
148 /*
149 solve by rotating line+quad so line is horizontal, then finding the roots
150 set up matrix to rotate quad to x-axis
151 |cos(a) -sin(a)|
152 |sin(a) cos(a)|
153 note that cos(a) = A(djacent) / Hypoteneuse
154 sin(a) = O(pposite) / Hypoteneuse
155 since we are computing Ts, we can ignore hypoteneuse, the scale factor:
156 | A -O |
157 | O A |
158 A = line[1].fX - line[0].fX (adjacent side of the right triangle)
159 O = line[1].fY - line[0].fY (opposite side of the right triangle)
160 for each of the three points (e.g. n = 0 to 2)
161 quad[n].fY' = (quad[n].fY - line[0].fY) * A - (quad[n].fX - line[0].fX) * O
162 */
163 double adj = (*fLine)[1].fX - (*fLine)[0].fX;
164 double opp = (*fLine)[1].fY - (*fLine)[0].fY;
165 double r[3];
166 for (int n = 0; n < 3; ++n) {
167 r[n] = (fQuad[n].fY - (*fLine)[0].fY) * adj - (fQuad[n].fX - (*fLine)[0].fX) * opp;
168 }
169 double A = r[2];
170 double B = r[1];
171 double C = r[0];
172 A += C - 2 * B; // A = a - 2*b + c
173 B -= C; // B = -(b - c)
174 return SkDQuad::RootsValidT(A, 2 * B, C, roots);
175 }
176
intersect()177 int intersect() {
178 addExactEndPoints();
179 if (fAllowNear) {
180 addNearEndPoints();
181 }
182 double rootVals[2];
183 int roots = intersectRay(rootVals);
184 for (int index = 0; index < roots; ++index) {
185 double quadT = rootVals[index];
186 double lineT = findLineT(quadT);
187 SkDPoint pt;
188 if (pinTs(&quadT, &lineT, &pt, kPointUninitialized) && uniqueAnswer(quadT, pt)) {
189 fIntersections->insert(quadT, lineT, pt);
190 }
191 }
192 checkCoincident();
193 return fIntersections->used();
194 }
195
horizontalIntersect(double axisIntercept,double roots[2])196 int horizontalIntersect(double axisIntercept, double roots[2]) {
197 double D = fQuad[2].fY; // f
198 double E = fQuad[1].fY; // e
199 double F = fQuad[0].fY; // d
200 D += F - 2 * E; // D = d - 2*e + f
201 E -= F; // E = -(d - e)
202 F -= axisIntercept;
203 return SkDQuad::RootsValidT(D, 2 * E, F, roots);
204 }
205
horizontalIntersect(double axisIntercept,double left,double right,bool flipped)206 int horizontalIntersect(double axisIntercept, double left, double right, bool flipped) {
207 addExactHorizontalEndPoints(left, right, axisIntercept);
208 if (fAllowNear) {
209 addNearHorizontalEndPoints(left, right, axisIntercept);
210 }
211 double rootVals[2];
212 int roots = horizontalIntersect(axisIntercept, rootVals);
213 for (int index = 0; index < roots; ++index) {
214 double quadT = rootVals[index];
215 SkDPoint pt = fQuad.ptAtT(quadT);
216 double lineT = (pt.fX - left) / (right - left);
217 if (pinTs(&quadT, &lineT, &pt, kPointInitialized) && uniqueAnswer(quadT, pt)) {
218 fIntersections->insert(quadT, lineT, pt);
219 }
220 }
221 if (flipped) {
222 fIntersections->flip();
223 }
224 checkCoincident();
225 return fIntersections->used();
226 }
227
uniqueAnswer(double quadT,const SkDPoint & pt)228 bool uniqueAnswer(double quadT, const SkDPoint& pt) {
229 for (int inner = 0; inner < fIntersections->used(); ++inner) {
230 if (fIntersections->pt(inner) != pt) {
231 continue;
232 }
233 double existingQuadT = (*fIntersections)[0][inner];
234 if (quadT == existingQuadT) {
235 return false;
236 }
237 // check if midway on quad is also same point. If so, discard this
238 double quadMidT = (existingQuadT + quadT) / 2;
239 SkDPoint quadMidPt = fQuad.ptAtT(quadMidT);
240 if (quadMidPt.approximatelyEqual(pt)) {
241 return false;
242 }
243 }
244 #if ONE_OFF_DEBUG
245 SkDPoint qPt = fQuad.ptAtT(quadT);
246 SkDebugf("%s pt=(%1.9g,%1.9g) cPt=(%1.9g,%1.9g)\n", __FUNCTION__, pt.fX, pt.fY,
247 qPt.fX, qPt.fY);
248 #endif
249 return true;
250 }
251
verticalIntersect(double axisIntercept,double roots[2])252 int verticalIntersect(double axisIntercept, double roots[2]) {
253 double D = fQuad[2].fX; // f
254 double E = fQuad[1].fX; // e
255 double F = fQuad[0].fX; // d
256 D += F - 2 * E; // D = d - 2*e + f
257 E -= F; // E = -(d - e)
258 F -= axisIntercept;
259 return SkDQuad::RootsValidT(D, 2 * E, F, roots);
260 }
261
verticalIntersect(double axisIntercept,double top,double bottom,bool flipped)262 int verticalIntersect(double axisIntercept, double top, double bottom, bool flipped) {
263 addExactVerticalEndPoints(top, bottom, axisIntercept);
264 if (fAllowNear) {
265 addNearVerticalEndPoints(top, bottom, axisIntercept);
266 }
267 double rootVals[2];
268 int roots = verticalIntersect(axisIntercept, rootVals);
269 for (int index = 0; index < roots; ++index) {
270 double quadT = rootVals[index];
271 SkDPoint pt = fQuad.ptAtT(quadT);
272 double lineT = (pt.fY - top) / (bottom - top);
273 if (pinTs(&quadT, &lineT, &pt, kPointInitialized) && uniqueAnswer(quadT, pt)) {
274 fIntersections->insert(quadT, lineT, pt);
275 }
276 }
277 if (flipped) {
278 fIntersections->flip();
279 }
280 checkCoincident();
281 return fIntersections->used();
282 }
283
284 protected:
285 // add endpoints first to get zero and one t values exactly
addExactEndPoints()286 void addExactEndPoints() {
287 for (int qIndex = 0; qIndex < 3; qIndex += 2) {
288 double lineT = fLine->exactPoint(fQuad[qIndex]);
289 if (lineT < 0) {
290 continue;
291 }
292 double quadT = (double) (qIndex >> 1);
293 fIntersections->insert(quadT, lineT, fQuad[qIndex]);
294 }
295 }
296
addNearEndPoints()297 void addNearEndPoints() {
298 for (int qIndex = 0; qIndex < 3; qIndex += 2) {
299 double quadT = (double) (qIndex >> 1);
300 if (fIntersections->hasT(quadT)) {
301 continue;
302 }
303 double lineT = fLine->nearPoint(fQuad[qIndex], nullptr);
304 if (lineT < 0) {
305 continue;
306 }
307 fIntersections->insert(quadT, lineT, fQuad[qIndex]);
308 }
309 this->addLineNearEndPoints();
310 }
311
addLineNearEndPoints()312 void addLineNearEndPoints() {
313 for (int lIndex = 0; lIndex < 2; ++lIndex) {
314 double lineT = (double) lIndex;
315 if (fIntersections->hasOppT(lineT)) {
316 continue;
317 }
318 double quadT = ((const SkDCurve*) &fQuad)->nearPoint(SkPath::kQuad_Verb,
319 (*fLine)[lIndex], (*fLine)[!lIndex]);
320 if (quadT < 0) {
321 continue;
322 }
323 fIntersections->insert(quadT, lineT, (*fLine)[lIndex]);
324 }
325 }
326
addExactHorizontalEndPoints(double left,double right,double y)327 void addExactHorizontalEndPoints(double left, double right, double y) {
328 for (int qIndex = 0; qIndex < 3; qIndex += 2) {
329 double lineT = SkDLine::ExactPointH(fQuad[qIndex], left, right, y);
330 if (lineT < 0) {
331 continue;
332 }
333 double quadT = (double) (qIndex >> 1);
334 fIntersections->insert(quadT, lineT, fQuad[qIndex]);
335 }
336 }
337
addNearHorizontalEndPoints(double left,double right,double y)338 void addNearHorizontalEndPoints(double left, double right, double y) {
339 for (int qIndex = 0; qIndex < 3; qIndex += 2) {
340 double quadT = (double) (qIndex >> 1);
341 if (fIntersections->hasT(quadT)) {
342 continue;
343 }
344 double lineT = SkDLine::NearPointH(fQuad[qIndex], left, right, y);
345 if (lineT < 0) {
346 continue;
347 }
348 fIntersections->insert(quadT, lineT, fQuad[qIndex]);
349 }
350 this->addLineNearEndPoints();
351 }
352
addExactVerticalEndPoints(double top,double bottom,double x)353 void addExactVerticalEndPoints(double top, double bottom, double x) {
354 for (int qIndex = 0; qIndex < 3; qIndex += 2) {
355 double lineT = SkDLine::ExactPointV(fQuad[qIndex], top, bottom, x);
356 if (lineT < 0) {
357 continue;
358 }
359 double quadT = (double) (qIndex >> 1);
360 fIntersections->insert(quadT, lineT, fQuad[qIndex]);
361 }
362 }
363
addNearVerticalEndPoints(double top,double bottom,double x)364 void addNearVerticalEndPoints(double top, double bottom, double x) {
365 for (int qIndex = 0; qIndex < 3; qIndex += 2) {
366 double quadT = (double) (qIndex >> 1);
367 if (fIntersections->hasT(quadT)) {
368 continue;
369 }
370 double lineT = SkDLine::NearPointV(fQuad[qIndex], top, bottom, x);
371 if (lineT < 0) {
372 continue;
373 }
374 fIntersections->insert(quadT, lineT, fQuad[qIndex]);
375 }
376 this->addLineNearEndPoints();
377 }
378
findLineT(double t)379 double findLineT(double t) {
380 SkDPoint xy = fQuad.ptAtT(t);
381 double dx = (*fLine)[1].fX - (*fLine)[0].fX;
382 double dy = (*fLine)[1].fY - (*fLine)[0].fY;
383 if (fabs(dx) > fabs(dy)) {
384 return (xy.fX - (*fLine)[0].fX) / dx;
385 }
386 return (xy.fY - (*fLine)[0].fY) / dy;
387 }
388
pinTs(double * quadT,double * lineT,SkDPoint * pt,PinTPoint ptSet)389 bool pinTs(double* quadT, double* lineT, SkDPoint* pt, PinTPoint ptSet) {
390 if (!approximately_one_or_less_double(*lineT)) {
391 return false;
392 }
393 if (!approximately_zero_or_more_double(*lineT)) {
394 return false;
395 }
396 double qT = *quadT = SkPinT(*quadT);
397 double lT = *lineT = SkPinT(*lineT);
398 if (lT == 0 || lT == 1 || (ptSet == kPointUninitialized && qT != 0 && qT != 1)) {
399 *pt = (*fLine).ptAtT(lT);
400 } else if (ptSet == kPointUninitialized) {
401 *pt = fQuad.ptAtT(qT);
402 }
403 SkPoint gridPt = pt->asSkPoint();
404 if (SkDPoint::ApproximatelyEqual(gridPt, (*fLine)[0].asSkPoint())) {
405 *pt = (*fLine)[0];
406 *lineT = 0;
407 } else if (SkDPoint::ApproximatelyEqual(gridPt, (*fLine)[1].asSkPoint())) {
408 *pt = (*fLine)[1];
409 *lineT = 1;
410 }
411 if (fIntersections->used() > 0 && approximately_equal((*fIntersections)[1][0], *lineT)) {
412 return false;
413 }
414 if (gridPt == fQuad[0].asSkPoint()) {
415 *pt = fQuad[0];
416 *quadT = 0;
417 } else if (gridPt == fQuad[2].asSkPoint()) {
418 *pt = fQuad[2];
419 *quadT = 1;
420 }
421 return true;
422 }
423
424 private:
425 const SkDQuad& fQuad;
426 const SkDLine* fLine;
427 SkIntersections* fIntersections;
428 bool fAllowNear;
429 };
430
horizontal(const SkDQuad & quad,double left,double right,double y,bool flipped)431 int SkIntersections::horizontal(const SkDQuad& quad, double left, double right, double y,
432 bool flipped) {
433 SkDLine line = {{{ left, y }, { right, y }}};
434 LineQuadraticIntersections q(quad, line, this);
435 return q.horizontalIntersect(y, left, right, flipped);
436 }
437
vertical(const SkDQuad & quad,double top,double bottom,double x,bool flipped)438 int SkIntersections::vertical(const SkDQuad& quad, double top, double bottom, double x,
439 bool flipped) {
440 SkDLine line = {{{ x, top }, { x, bottom }}};
441 LineQuadraticIntersections q(quad, line, this);
442 return q.verticalIntersect(x, top, bottom, flipped);
443 }
444
intersect(const SkDQuad & quad,const SkDLine & line)445 int SkIntersections::intersect(const SkDQuad& quad, const SkDLine& line) {
446 LineQuadraticIntersections q(quad, line, this);
447 q.allowNear(fAllowNear);
448 return q.intersect();
449 }
450
intersectRay(const SkDQuad & quad,const SkDLine & line)451 int SkIntersections::intersectRay(const SkDQuad& quad, const SkDLine& line) {
452 LineQuadraticIntersections q(quad, line, this);
453 fUsed = q.intersectRay(fT[0]);
454 for (int index = 0; index < fUsed; ++index) {
455 fPt[index] = quad.ptAtT(fT[0][index]);
456 }
457 return fUsed;
458 }
459
HorizontalIntercept(const SkDQuad & quad,SkScalar y,double * roots)460 int SkIntersections::HorizontalIntercept(const SkDQuad& quad, SkScalar y, double* roots) {
461 LineQuadraticIntersections q(quad);
462 return q.horizontalIntersect(y, roots);
463 }
464
VerticalIntercept(const SkDQuad & quad,SkScalar x,double * roots)465 int SkIntersections::VerticalIntercept(const SkDQuad& quad, SkScalar x, double* roots) {
466 LineQuadraticIntersections q(quad);
467 return q.verticalIntersect(x, roots);
468 }
469
470 // SkDQuad accessors to Intersection utilities
471
horizontalIntersect(double yIntercept,double roots[2]) const472 int SkDQuad::horizontalIntersect(double yIntercept, double roots[2]) const {
473 return SkIntersections::HorizontalIntercept(*this, yIntercept, roots);
474 }
475
verticalIntersect(double xIntercept,double roots[2]) const476 int SkDQuad::verticalIntersect(double xIntercept, double roots[2]) const {
477 return SkIntersections::VerticalIntercept(*this, xIntercept, roots);
478 }
479