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 class LineQuadraticIntersections { 90 public: 91 enum PinTPoint { 92 kPointUninitialized, 93 kPointInitialized 94 }; 95 96 LineQuadraticIntersections(const SkDQuad& q, const SkDLine& l, SkIntersections* i) 97 : fQuad(q) 98 , fLine(&l) 99 , fIntersections(i) 100 , fAllowNear(true) { 101 i->setMax(3); // allow short partial coincidence plus discrete intersection 102 } 103 104 LineQuadraticIntersections(const SkDQuad& q) 105 : fQuad(q) 106 SkDEBUGPARAMS(fLine(NULL)) 107 SkDEBUGPARAMS(fIntersections(NULL)) 108 SkDEBUGPARAMS(fAllowNear(false)) { 109 } 110 111 void allowNear(bool allow) { 112 fAllowNear = allow; 113 } 114 115 void checkCoincident() { 116 int last = fIntersections->used() - 1; 117 for (int index = 0; index < last; ) { 118 double quadMidT = ((*fIntersections)[0][index] + (*fIntersections)[0][index + 1]) / 2; 119 SkDPoint quadMidPt = fQuad.ptAtT(quadMidT); 120 double t = fLine->nearPoint(quadMidPt, NULL); 121 if (t < 0) { 122 ++index; 123 continue; 124 } 125 if (fIntersections->isCoincident(index)) { 126 fIntersections->removeOne(index); 127 --last; 128 } else if (fIntersections->isCoincident(index + 1)) { 129 fIntersections->removeOne(index + 1); 130 --last; 131 } else { 132 fIntersections->setCoincident(index++); 133 } 134 fIntersections->setCoincident(index); 135 } 136 } 137 138 int intersectRay(double roots[2]) { 139 /* 140 solve by rotating line+quad so line is horizontal, then finding the roots 141 set up matrix to rotate quad to x-axis 142 |cos(a) -sin(a)| 143 |sin(a) cos(a)| 144 note that cos(a) = A(djacent) / Hypoteneuse 145 sin(a) = O(pposite) / Hypoteneuse 146 since we are computing Ts, we can ignore hypoteneuse, the scale factor: 147 | A -O | 148 | O A | 149 A = line[1].fX - line[0].fX (adjacent side of the right triangle) 150 O = line[1].fY - line[0].fY (opposite side of the right triangle) 151 for each of the three points (e.g. n = 0 to 2) 152 quad[n].fY' = (quad[n].fY - line[0].fY) * A - (quad[n].fX - line[0].fX) * O 153 */ 154 double adj = (*fLine)[1].fX - (*fLine)[0].fX; 155 double opp = (*fLine)[1].fY - (*fLine)[0].fY; 156 double r[3]; 157 for (int n = 0; n < 3; ++n) { 158 r[n] = (fQuad[n].fY - (*fLine)[0].fY) * adj - (fQuad[n].fX - (*fLine)[0].fX) * opp; 159 } 160 double A = r[2]; 161 double B = r[1]; 162 double C = r[0]; 163 A += C - 2 * B; // A = a - 2*b + c 164 B -= C; // B = -(b - c) 165 return SkDQuad::RootsValidT(A, 2 * B, C, roots); 166 } 167 168 int intersect() { 169 addExactEndPoints(); 170 if (fAllowNear) { 171 addNearEndPoints(); 172 } 173 double rootVals[2]; 174 int roots = intersectRay(rootVals); 175 for (int index = 0; index < roots; ++index) { 176 double quadT = rootVals[index]; 177 double lineT = findLineT(quadT); 178 SkDPoint pt; 179 if (pinTs(&quadT, &lineT, &pt, kPointUninitialized) && uniqueAnswer(quadT, pt)) { 180 fIntersections->insert(quadT, lineT, pt); 181 } 182 } 183 checkCoincident(); 184 return fIntersections->used(); 185 } 186 187 int horizontalIntersect(double axisIntercept, double roots[2]) { 188 double D = fQuad[2].fY; // f 189 double E = fQuad[1].fY; // e 190 double F = fQuad[0].fY; // d 191 D += F - 2 * E; // D = d - 2*e + f 192 E -= F; // E = -(d - e) 193 F -= axisIntercept; 194 return SkDQuad::RootsValidT(D, 2 * E, F, roots); 195 } 196 197 int horizontalIntersect(double axisIntercept, double left, double right, bool flipped) { 198 addExactHorizontalEndPoints(left, right, axisIntercept); 199 if (fAllowNear) { 200 addNearHorizontalEndPoints(left, right, axisIntercept); 201 } 202 double rootVals[2]; 203 int roots = horizontalIntersect(axisIntercept, rootVals); 204 for (int index = 0; index < roots; ++index) { 205 double quadT = rootVals[index]; 206 SkDPoint pt = fQuad.ptAtT(quadT); 207 double lineT = (pt.fX - left) / (right - left); 208 if (pinTs(&quadT, &lineT, &pt, kPointInitialized) && uniqueAnswer(quadT, pt)) { 209 fIntersections->insert(quadT, lineT, pt); 210 } 211 } 212 if (flipped) { 213 fIntersections->flip(); 214 } 215 checkCoincident(); 216 return fIntersections->used(); 217 } 218 219 bool uniqueAnswer(double quadT, const SkDPoint& pt) { 220 for (int inner = 0; inner < fIntersections->used(); ++inner) { 221 if (fIntersections->pt(inner) != pt) { 222 continue; 223 } 224 double existingQuadT = (*fIntersections)[0][inner]; 225 if (quadT == existingQuadT) { 226 return false; 227 } 228 // check if midway on quad is also same point. If so, discard this 229 double quadMidT = (existingQuadT + quadT) / 2; 230 SkDPoint quadMidPt = fQuad.ptAtT(quadMidT); 231 if (quadMidPt.approximatelyEqual(pt)) { 232 return false; 233 } 234 } 235 #if ONE_OFF_DEBUG 236 SkDPoint qPt = fQuad.ptAtT(quadT); 237 SkDebugf("%s pt=(%1.9g,%1.9g) cPt=(%1.9g,%1.9g)\n", __FUNCTION__, pt.fX, pt.fY, 238 qPt.fX, qPt.fY); 239 #endif 240 return true; 241 } 242 243 int verticalIntersect(double axisIntercept, double roots[2]) { 244 double D = fQuad[2].fX; // f 245 double E = fQuad[1].fX; // e 246 double F = fQuad[0].fX; // d 247 D += F - 2 * E; // D = d - 2*e + f 248 E -= F; // E = -(d - e) 249 F -= axisIntercept; 250 return SkDQuad::RootsValidT(D, 2 * E, F, roots); 251 } 252 253 int verticalIntersect(double axisIntercept, double top, double bottom, bool flipped) { 254 addExactVerticalEndPoints(top, bottom, axisIntercept); 255 if (fAllowNear) { 256 addNearVerticalEndPoints(top, bottom, axisIntercept); 257 } 258 double rootVals[2]; 259 int roots = verticalIntersect(axisIntercept, rootVals); 260 for (int index = 0; index < roots; ++index) { 261 double quadT = rootVals[index]; 262 SkDPoint pt = fQuad.ptAtT(quadT); 263 double lineT = (pt.fY - top) / (bottom - top); 264 if (pinTs(&quadT, &lineT, &pt, kPointInitialized) && uniqueAnswer(quadT, pt)) { 265 fIntersections->insert(quadT, lineT, pt); 266 } 267 } 268 if (flipped) { 269 fIntersections->flip(); 270 } 271 checkCoincident(); 272 return fIntersections->used(); 273 } 274 275 protected: 276 // add endpoints first to get zero and one t values exactly 277 void addExactEndPoints() { 278 for (int qIndex = 0; qIndex < 3; qIndex += 2) { 279 double lineT = fLine->exactPoint(fQuad[qIndex]); 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 addNearEndPoints() { 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 = fLine->nearPoint(fQuad[qIndex], NULL); 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 void addExactHorizontalEndPoints(double left, double right, double y) { 304 for (int qIndex = 0; qIndex < 3; qIndex += 2) { 305 double lineT = SkDLine::ExactPointH(fQuad[qIndex], left, right, y); 306 if (lineT < 0) { 307 continue; 308 } 309 double quadT = (double) (qIndex >> 1); 310 fIntersections->insert(quadT, lineT, fQuad[qIndex]); 311 } 312 } 313 314 void addNearHorizontalEndPoints(double left, double right, double y) { 315 for (int qIndex = 0; qIndex < 3; qIndex += 2) { 316 double quadT = (double) (qIndex >> 1); 317 if (fIntersections->hasT(quadT)) { 318 continue; 319 } 320 double lineT = SkDLine::NearPointH(fQuad[qIndex], left, right, y); 321 if (lineT < 0) { 322 continue; 323 } 324 fIntersections->insert(quadT, lineT, fQuad[qIndex]); 325 } 326 // FIXME: see if line end is nearly on quad 327 } 328 329 void addExactVerticalEndPoints(double top, double bottom, double x) { 330 for (int qIndex = 0; qIndex < 3; qIndex += 2) { 331 double lineT = SkDLine::ExactPointV(fQuad[qIndex], top, bottom, x); 332 if (lineT < 0) { 333 continue; 334 } 335 double quadT = (double) (qIndex >> 1); 336 fIntersections->insert(quadT, lineT, fQuad[qIndex]); 337 } 338 } 339 340 void addNearVerticalEndPoints(double top, double bottom, double x) { 341 for (int qIndex = 0; qIndex < 3; qIndex += 2) { 342 double quadT = (double) (qIndex >> 1); 343 if (fIntersections->hasT(quadT)) { 344 continue; 345 } 346 double lineT = SkDLine::NearPointV(fQuad[qIndex], top, bottom, x); 347 if (lineT < 0) { 348 continue; 349 } 350 fIntersections->insert(quadT, lineT, fQuad[qIndex]); 351 } 352 // FIXME: see if line end is nearly on quad 353 } 354 355 double findLineT(double t) { 356 SkDPoint xy = fQuad.ptAtT(t); 357 double dx = (*fLine)[1].fX - (*fLine)[0].fX; 358 double dy = (*fLine)[1].fY - (*fLine)[0].fY; 359 if (fabs(dx) > fabs(dy)) { 360 return (xy.fX - (*fLine)[0].fX) / dx; 361 } 362 return (xy.fY - (*fLine)[0].fY) / dy; 363 } 364 365 bool pinTs(double* quadT, double* lineT, SkDPoint* pt, PinTPoint ptSet) { 366 if (!approximately_one_or_less_double(*lineT)) { 367 return false; 368 } 369 if (!approximately_zero_or_more_double(*lineT)) { 370 return false; 371 } 372 double qT = *quadT = SkPinT(*quadT); 373 double lT = *lineT = SkPinT(*lineT); 374 if (lT == 0 || lT == 1 || (ptSet == kPointUninitialized && qT != 0 && qT != 1)) { 375 *pt = (*fLine).ptAtT(lT); 376 } else if (ptSet == kPointUninitialized) { 377 *pt = fQuad.ptAtT(qT); 378 } 379 SkPoint gridPt = pt->asSkPoint(); 380 if (SkDPoint::ApproximatelyEqual(gridPt, (*fLine)[0].asSkPoint())) { 381 *pt = (*fLine)[0]; 382 *lineT = 0; 383 } else if (SkDPoint::ApproximatelyEqual(gridPt, (*fLine)[1].asSkPoint())) { 384 *pt = (*fLine)[1]; 385 *lineT = 1; 386 } 387 if (fIntersections->used() > 0 && approximately_equal((*fIntersections)[1][0], *lineT)) { 388 return false; 389 } 390 if (gridPt == fQuad[0].asSkPoint()) { 391 *pt = fQuad[0]; 392 *quadT = 0; 393 } else if (gridPt == fQuad[2].asSkPoint()) { 394 *pt = fQuad[2]; 395 *quadT = 1; 396 } 397 return true; 398 } 399 400 private: 401 const SkDQuad& fQuad; 402 const SkDLine* fLine; 403 SkIntersections* fIntersections; 404 bool fAllowNear; 405 }; 406 407 int SkIntersections::horizontal(const SkDQuad& quad, double left, double right, double y, 408 bool flipped) { 409 SkDLine line = {{{ left, y }, { right, y }}}; 410 LineQuadraticIntersections q(quad, line, this); 411 return q.horizontalIntersect(y, left, right, flipped); 412 } 413 414 int SkIntersections::vertical(const SkDQuad& quad, double top, double bottom, double x, 415 bool flipped) { 416 SkDLine line = {{{ x, top }, { x, bottom }}}; 417 LineQuadraticIntersections q(quad, line, this); 418 return q.verticalIntersect(x, top, bottom, flipped); 419 } 420 421 int SkIntersections::intersect(const SkDQuad& quad, const SkDLine& line) { 422 LineQuadraticIntersections q(quad, line, this); 423 q.allowNear(fAllowNear); 424 return q.intersect(); 425 } 426 427 int SkIntersections::intersectRay(const SkDQuad& quad, const SkDLine& line) { 428 LineQuadraticIntersections q(quad, line, this); 429 fUsed = q.intersectRay(fT[0]); 430 for (int index = 0; index < fUsed; ++index) { 431 fPt[index] = quad.ptAtT(fT[0][index]); 432 } 433 return fUsed; 434 } 435 436 int SkIntersections::HorizontalIntercept(const SkDQuad& quad, SkScalar y, double* roots) { 437 LineQuadraticIntersections q(quad); 438 return q.horizontalIntersect(y, roots); 439 } 440 441 int SkIntersections::VerticalIntercept(const SkDQuad& quad, SkScalar x, double* roots) { 442 LineQuadraticIntersections q(quad); 443 return q.verticalIntersect(x, roots); 444 } 445