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