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