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
96 // see parallel routine in line quadratic intersections
intersectRay(double roots[3])97 int intersectRay(double roots[3]) {
98 double adj = fLine[1].fX - fLine[0].fX;
99 double opp = fLine[1].fY - fLine[0].fY;
100 SkDCubic r;
101 for (int n = 0; n < 4; ++n) {
102 r[n].fX = (fCubic[n].fY - fLine[0].fY) * adj - (fCubic[n].fX - fLine[0].fX) * opp;
103 }
104 double A, B, C, D;
105 SkDCubic::Coefficients(&r[0].fX, &A, &B, &C, &D);
106 return SkDCubic::RootsValidT(A, B, C, D, roots);
107 }
108
intersect()109 int intersect() {
110 addExactEndPoints();
111 if (fAllowNear) {
112 addNearEndPoints();
113 }
114 double rootVals[3];
115 int roots = intersectRay(rootVals);
116 for (int index = 0; index < roots; ++index) {
117 double cubicT = rootVals[index];
118 double lineT = findLineT(cubicT);
119 SkDPoint pt;
120 if (pinTs(&cubicT, &lineT, &pt, kPointUninitialized)) {
121 #if ONE_OFF_DEBUG
122 SkDPoint cPt = fCubic.ptAtT(cubicT);
123 SkDebugf("%s pt=(%1.9g,%1.9g) cPt=(%1.9g,%1.9g)\n", __FUNCTION__, pt.fX, pt.fY,
124 cPt.fX, cPt.fY);
125 #endif
126 for (int inner = 0; inner < fIntersections->used(); ++inner) {
127 if (fIntersections->pt(inner) != pt) {
128 continue;
129 }
130 double existingCubicT = (*fIntersections)[0][inner];
131 if (cubicT == existingCubicT) {
132 goto skipInsert;
133 }
134 // check if midway on cubic is also same point. If so, discard this
135 double cubicMidT = (existingCubicT + cubicT) / 2;
136 SkDPoint cubicMidPt = fCubic.ptAtT(cubicMidT);
137 if (cubicMidPt.approximatelyEqual(pt)) {
138 goto skipInsert;
139 }
140 }
141 fIntersections->insert(cubicT, lineT, pt);
142 skipInsert:
143 ;
144 }
145 }
146 return fIntersections->used();
147 }
148
horizontalIntersect(double axisIntercept,double roots[3])149 int horizontalIntersect(double axisIntercept, double roots[3]) {
150 double A, B, C, D;
151 SkDCubic::Coefficients(&fCubic[0].fY, &A, &B, &C, &D);
152 D -= axisIntercept;
153 return SkDCubic::RootsValidT(A, B, C, D, roots);
154 }
155
horizontalIntersect(double axisIntercept,double left,double right,bool flipped)156 int horizontalIntersect(double axisIntercept, double left, double right, bool flipped) {
157 addExactHorizontalEndPoints(left, right, axisIntercept);
158 if (fAllowNear) {
159 addNearHorizontalEndPoints(left, right, axisIntercept);
160 }
161 double rootVals[3];
162 int roots = horizontalIntersect(axisIntercept, rootVals);
163 for (int index = 0; index < roots; ++index) {
164 double cubicT = rootVals[index];
165 SkDPoint pt = fCubic.ptAtT(cubicT);
166 double lineT = (pt.fX - left) / (right - left);
167 if (pinTs(&cubicT, &lineT, &pt, kPointInitialized)) {
168 fIntersections->insert(cubicT, lineT, pt);
169 }
170 }
171 if (flipped) {
172 fIntersections->flip();
173 }
174 return fIntersections->used();
175 }
176
verticalIntersect(double axisIntercept,double roots[3])177 int verticalIntersect(double axisIntercept, double roots[3]) {
178 double A, B, C, D;
179 SkDCubic::Coefficients(&fCubic[0].fX, &A, &B, &C, &D);
180 D -= axisIntercept;
181 return SkDCubic::RootsValidT(A, B, C, D, roots);
182 }
183
verticalIntersect(double axisIntercept,double top,double bottom,bool flipped)184 int verticalIntersect(double axisIntercept, double top, double bottom, bool flipped) {
185 addExactVerticalEndPoints(top, bottom, axisIntercept);
186 if (fAllowNear) {
187 addNearVerticalEndPoints(top, bottom, axisIntercept);
188 }
189 double rootVals[3];
190 int roots = verticalIntersect(axisIntercept, rootVals);
191 for (int index = 0; index < roots; ++index) {
192 double cubicT = rootVals[index];
193 SkDPoint pt = fCubic.ptAtT(cubicT);
194 double lineT = (pt.fY - top) / (bottom - top);
195 if (pinTs(&cubicT, &lineT, &pt, kPointInitialized)) {
196 fIntersections->insert(cubicT, lineT, pt);
197 }
198 }
199 if (flipped) {
200 fIntersections->flip();
201 }
202 return fIntersections->used();
203 }
204
205 protected:
206
addExactEndPoints()207 void addExactEndPoints() {
208 for (int cIndex = 0; cIndex < 4; cIndex += 3) {
209 double lineT = fLine.exactPoint(fCubic[cIndex]);
210 if (lineT < 0) {
211 continue;
212 }
213 double cubicT = (double) (cIndex >> 1);
214 fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
215 }
216 }
217
218 /* Note that this does not look for endpoints of the line that are near the cubic.
219 These points are found later when check ends looks for missing points */
addNearEndPoints()220 void addNearEndPoints() {
221 for (int cIndex = 0; cIndex < 4; cIndex += 3) {
222 double cubicT = (double) (cIndex >> 1);
223 if (fIntersections->hasT(cubicT)) {
224 continue;
225 }
226 double lineT = fLine.nearPoint(fCubic[cIndex]);
227 if (lineT < 0) {
228 continue;
229 }
230 fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
231 }
232 }
233
addExactHorizontalEndPoints(double left,double right,double y)234 void addExactHorizontalEndPoints(double left, double right, double y) {
235 for (int cIndex = 0; cIndex < 4; cIndex += 3) {
236 double lineT = SkDLine::ExactPointH(fCubic[cIndex], left, right, y);
237 if (lineT < 0) {
238 continue;
239 }
240 double cubicT = (double) (cIndex >> 1);
241 fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
242 }
243 }
244
addNearHorizontalEndPoints(double left,double right,double y)245 void addNearHorizontalEndPoints(double left, double right, double y) {
246 for (int cIndex = 0; cIndex < 4; cIndex += 3) {
247 double cubicT = (double) (cIndex >> 1);
248 if (fIntersections->hasT(cubicT)) {
249 continue;
250 }
251 double lineT = SkDLine::NearPointH(fCubic[cIndex], left, right, y);
252 if (lineT < 0) {
253 continue;
254 }
255 fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
256 }
257 // FIXME: see if line end is nearly on cubic
258 }
259
addExactVerticalEndPoints(double top,double bottom,double x)260 void addExactVerticalEndPoints(double top, double bottom, double x) {
261 for (int cIndex = 0; cIndex < 4; cIndex += 3) {
262 double lineT = SkDLine::ExactPointV(fCubic[cIndex], top, bottom, x);
263 if (lineT < 0) {
264 continue;
265 }
266 double cubicT = (double) (cIndex >> 1);
267 fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
268 }
269 }
270
addNearVerticalEndPoints(double top,double bottom,double x)271 void addNearVerticalEndPoints(double top, double bottom, double x) {
272 for (int cIndex = 0; cIndex < 4; cIndex += 3) {
273 double cubicT = (double) (cIndex >> 1);
274 if (fIntersections->hasT(cubicT)) {
275 continue;
276 }
277 double lineT = SkDLine::NearPointV(fCubic[cIndex], top, bottom, x);
278 if (lineT < 0) {
279 continue;
280 }
281 fIntersections->insert(cubicT, lineT, fCubic[cIndex]);
282 }
283 // FIXME: see if line end is nearly on cubic
284 }
285
findLineT(double t)286 double findLineT(double t) {
287 SkDPoint xy = fCubic.ptAtT(t);
288 double dx = fLine[1].fX - fLine[0].fX;
289 double dy = fLine[1].fY - fLine[0].fY;
290 if (fabs(dx) > fabs(dy)) {
291 return (xy.fX - fLine[0].fX) / dx;
292 }
293 return (xy.fY - fLine[0].fY) / dy;
294 }
295
pinTs(double * cubicT,double * lineT,SkDPoint * pt,PinTPoint ptSet)296 bool pinTs(double* cubicT, double* lineT, SkDPoint* pt, PinTPoint ptSet) {
297 if (!approximately_one_or_less(*lineT)) {
298 return false;
299 }
300 if (!approximately_zero_or_more(*lineT)) {
301 return false;
302 }
303 double cT = *cubicT = SkPinT(*cubicT);
304 double lT = *lineT = SkPinT(*lineT);
305 if (lT == 0 || lT == 1 || (ptSet == kPointUninitialized && cT != 0 && cT != 1)) {
306 *pt = fLine.ptAtT(lT);
307 } else if (ptSet == kPointUninitialized) {
308 *pt = fCubic.ptAtT(cT);
309 }
310 SkPoint gridPt = pt->asSkPoint();
311 if (gridPt == fLine[0].asSkPoint()) {
312 *lineT = 0;
313 } else if (gridPt == fLine[1].asSkPoint()) {
314 *lineT = 1;
315 }
316 if (gridPt == fCubic[0].asSkPoint() && approximately_equal(*cubicT, 0)) {
317 *cubicT = 0;
318 } else if (gridPt == fCubic[3].asSkPoint() && approximately_equal(*cubicT, 1)) {
319 *cubicT = 1;
320 }
321 return true;
322 }
323
324 private:
325 const SkDCubic& fCubic;
326 const SkDLine& fLine;
327 SkIntersections* fIntersections;
328 bool fAllowNear;
329 };
330
horizontal(const SkDCubic & cubic,double left,double right,double y,bool flipped)331 int SkIntersections::horizontal(const SkDCubic& cubic, double left, double right, double y,
332 bool flipped) {
333 SkDLine line = {{{ left, y }, { right, y }}};
334 LineCubicIntersections c(cubic, line, this);
335 return c.horizontalIntersect(y, left, right, flipped);
336 }
337
vertical(const SkDCubic & cubic,double top,double bottom,double x,bool flipped)338 int SkIntersections::vertical(const SkDCubic& cubic, double top, double bottom, double x,
339 bool flipped) {
340 SkDLine line = {{{ x, top }, { x, bottom }}};
341 LineCubicIntersections c(cubic, line, this);
342 return c.verticalIntersect(x, top, bottom, flipped);
343 }
344
intersect(const SkDCubic & cubic,const SkDLine & line)345 int SkIntersections::intersect(const SkDCubic& cubic, const SkDLine& line) {
346 LineCubicIntersections c(cubic, line, this);
347 c.allowNear(fAllowNear);
348 return c.intersect();
349 }
350
intersectRay(const SkDCubic & cubic,const SkDLine & line)351 int SkIntersections::intersectRay(const SkDCubic& cubic, const SkDLine& line) {
352 LineCubicIntersections c(cubic, line, this);
353 fUsed = c.intersectRay(fT[0]);
354 for (int index = 0; index < fUsed; ++index) {
355 fPt[index] = cubic.ptAtT(fT[0][index]);
356 }
357 return fUsed;
358 }
359