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