• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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