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 "SkPathOpsCurve.h" 10 #include "SkPathOpsLine.h" 11 12 /* 13 Find the interection of a line and cubic by solving for valid t values. 14 15 Analogous to line-quadratic intersection, solve line-cubic intersection by 16 representing the cubic as: 17 x = a(1-t)^3 + 2b(1-t)^2t + c(1-t)t^2 + dt^3 18 y = e(1-t)^3 + 2f(1-t)^2t + g(1-t)t^2 + ht^3 19 and the line as: 20 y = i*x + j (if the line is more horizontal) 21 or: 22 x = i*y + j (if the line is more vertical) 23 24 Then using Mathematica, solve for the values of t where the cubic intersects the 25 line: 26 27 (in) Resultant[ 28 a*(1 - t)^3 + 3*b*(1 - t)^2*t + 3*c*(1 - t)*t^2 + d*t^3 - x, 29 e*(1 - t)^3 + 3*f*(1 - t)^2*t + 3*g*(1 - t)*t^2 + h*t^3 - i*x - j, x] 30 (out) -e + j + 31 3 e t - 3 f t - 32 3 e t^2 + 6 f t^2 - 3 g t^2 + 33 e t^3 - 3 f t^3 + 3 g t^3 - h t^3 + 34 i ( a - 35 3 a t + 3 b t + 36 3 a t^2 - 6 b t^2 + 3 c t^2 - 37 a t^3 + 3 b t^3 - 3 c t^3 + d t^3 ) 38 39 if i goes to infinity, we can rewrite the line in terms of x. Mathematica: 40 41 (in) Resultant[ 42 a*(1 - t)^3 + 3*b*(1 - t)^2*t + 3*c*(1 - t)*t^2 + d*t^3 - i*y - j, 43 e*(1 - t)^3 + 3*f*(1 - t)^2*t + 3*g*(1 - t)*t^2 + h*t^3 - y, y] 44 (out) a - j - 45 3 a t + 3 b t + 46 3 a t^2 - 6 b t^2 + 3 c t^2 - 47 a t^3 + 3 b t^3 - 3 c t^3 + d t^3 - 48 i ( e - 49 3 e t + 3 f t + 50 3 e t^2 - 6 f t^2 + 3 g t^2 - 51 e t^3 + 3 f t^3 - 3 g t^3 + h t^3 ) 52 53 Solving this with Mathematica produces an expression with hundreds of terms; 54 instead, use Numeric Solutions recipe to solve the cubic. 55 56 The near-horizontal case, in terms of: Ax^3 + Bx^2 + Cx + D == 0 57 A = (-(-e + 3*f - 3*g + h) + i*(-a + 3*b - 3*c + d) ) 58 B = 3*(-( e - 2*f + g ) + i*( a - 2*b + c ) ) 59 C = 3*(-(-e + f ) + i*(-a + b ) ) 60 D = (-( e ) + i*( a ) + j ) 61 62 The near-vertical case, in terms of: Ax^3 + Bx^2 + Cx + D == 0 63 A = ( (-a + 3*b - 3*c + d) - i*(-e + 3*f - 3*g + h) ) 64 B = 3*( ( a - 2*b + c ) - i*( e - 2*f + g ) ) 65 C = 3*( (-a + b ) - i*(-e + f ) ) 66 D = ( ( a ) - i*( e ) - j ) 67 68 For horizontal lines: 69 (in) Resultant[ 70 a*(1 - t)^3 + 3*b*(1 - t)^2*t + 3*c*(1 - t)*t^2 + d*t^3 - j, 71 e*(1 - t)^3 + 3*f*(1 - t)^2*t + 3*g*(1 - t)*t^2 + h*t^3 - y, y] 72 (out) e - j - 73 3 e t + 3 f t + 74 3 e t^2 - 6 f t^2 + 3 g t^2 - 75 e t^3 + 3 f t^3 - 3 g t^3 + h t^3 76 */ 77 78 class LineCubicIntersections { 79 public: 80 enum PinTPoint { 81 kPointUninitialized, 82 kPointInitialized 83 }; 84 85 LineCubicIntersections(const SkDCubic& c, const SkDLine& l, SkIntersections* i) 86 : fCubic(c) 87 , fLine(l) 88 , fIntersections(i) 89 , fAllowNear(true) { 90 i->setMax(4); 91 } 92 93 void allowNear(bool allow) { 94 fAllowNear = allow; 95 } 96 97 void checkCoincident() { 98 int last = fIntersections->used() - 1; 99 for (int index = 0; index < last; ) { 100 double cubicMidT = ((*fIntersections)[0][index] + (*fIntersections)[0][index + 1]) / 2; 101 SkDPoint cubicMidPt = fCubic.ptAtT(cubicMidT); 102 double t = fLine.nearPoint(cubicMidPt, nullptr); 103 if (t < 0) { 104 ++index; 105 continue; 106 } 107 if (fIntersections->isCoincident(index)) { 108 fIntersections->removeOne(index); 109 --last; 110 } else if (fIntersections->isCoincident(index + 1)) { 111 fIntersections->removeOne(index + 1); 112 --last; 113 } else { 114 fIntersections->setCoincident(index++); 115 } 116 fIntersections->setCoincident(index); 117 } 118 } 119 120 // see parallel routine in line quadratic intersections 121 int intersectRay(double roots[3]) { 122 double adj = fLine[1].fX - fLine[0].fX; 123 double opp = fLine[1].fY - fLine[0].fY; 124 SkDCubic c; 125 SkDEBUGCODE(c.fDebugGlobalState = fIntersections->globalState()); 126 for (int n = 0; n < 4; ++n) { 127 c[n].fX = (fCubic[n].fY - fLine[0].fY) * adj - (fCubic[n].fX - fLine[0].fX) * opp; 128 } 129 double A, B, C, D; 130 SkDCubic::Coefficients(&c[0].fX, &A, &B, &C, &D); 131 int count = SkDCubic::RootsValidT(A, B, C, D, roots); 132 for (int index = 0; index < count; ++index) { 133 SkDPoint calcPt = c.ptAtT(roots[index]); 134 if (!approximately_zero(calcPt.fX)) { 135 for (int n = 0; n < 4; ++n) { 136 c[n].fY = (fCubic[n].fY - fLine[0].fY) * opp 137 + (fCubic[n].fX - fLine[0].fX) * adj; 138 } 139 double extremeTs[6]; 140 int extrema = SkDCubic::FindExtrema(&c[0].fX, extremeTs); 141 count = c.searchRoots(extremeTs, extrema, 0, SkDCubic::kXAxis, roots); 142 break; 143 } 144 } 145 return count; 146 } 147 148 int intersect() { 149 addExactEndPoints(); 150 if (fAllowNear) { 151 addNearEndPoints(); 152 } 153 double rootVals[3]; 154 int roots = intersectRay(rootVals); 155 for (int index = 0; index < roots; ++index) { 156 double cubicT = rootVals[index]; 157 double lineT = findLineT(cubicT); 158 SkDPoint pt; 159 if (pinTs(&cubicT, &lineT, &pt, kPointUninitialized) && uniqueAnswer(cubicT, pt)) { 160 fIntersections->insert(cubicT, lineT, pt); 161 } 162 } 163 checkCoincident(); 164 return fIntersections->used(); 165 } 166 167 static int HorizontalIntersect(const SkDCubic& c, double axisIntercept, double roots[3]) { 168 double A, B, C, D; 169 SkDCubic::Coefficients(&c[0].fY, &A, &B, &C, &D); 170 D -= axisIntercept; 171 int count = SkDCubic::RootsValidT(A, B, C, D, roots); 172 for (int index = 0; index < count; ++index) { 173 SkDPoint calcPt = c.ptAtT(roots[index]); 174 if (!approximately_equal(calcPt.fY, axisIntercept)) { 175 double extremeTs[6]; 176 int extrema = SkDCubic::FindExtrema(&c[0].fY, extremeTs); 177 count = c.searchRoots(extremeTs, extrema, axisIntercept, SkDCubic::kYAxis, roots); 178 break; 179 } 180 } 181 return count; 182 } 183 184 int horizontalIntersect(double axisIntercept, double left, double right, bool flipped) { 185 addExactHorizontalEndPoints(left, right, axisIntercept); 186 if (fAllowNear) { 187 addNearHorizontalEndPoints(left, right, axisIntercept); 188 } 189 double roots[3]; 190 int count = HorizontalIntersect(fCubic, axisIntercept, roots); 191 for (int index = 0; index < count; ++index) { 192 double cubicT = roots[index]; 193 SkDPoint pt = { fCubic.ptAtT(cubicT).fX, axisIntercept }; 194 double lineT = (pt.fX - left) / (right - left); 195 if (pinTs(&cubicT, &lineT, &pt, kPointInitialized) && uniqueAnswer(cubicT, pt)) { 196 fIntersections->insert(cubicT, lineT, pt); 197 } 198 } 199 if (flipped) { 200 fIntersections->flip(); 201 } 202 checkCoincident(); 203 return fIntersections->used(); 204 } 205 206 bool uniqueAnswer(double cubicT, const SkDPoint& pt) { 207 for (int inner = 0; inner < fIntersections->used(); ++inner) { 208 if (fIntersections->pt(inner) != pt) { 209 continue; 210 } 211 double existingCubicT = (*fIntersections)[0][inner]; 212 if (cubicT == existingCubicT) { 213 return false; 214 } 215 // check if midway on cubic is also same point. If so, discard this 216 double cubicMidT = (existingCubicT + cubicT) / 2; 217 SkDPoint cubicMidPt = fCubic.ptAtT(cubicMidT); 218 if (cubicMidPt.approximatelyEqual(pt)) { 219 return false; 220 } 221 } 222 #if ONE_OFF_DEBUG 223 SkDPoint cPt = fCubic.ptAtT(cubicT); 224 SkDebugf("%s pt=(%1.9g,%1.9g) cPt=(%1.9g,%1.9g)\n", __FUNCTION__, pt.fX, pt.fY, 225 cPt.fX, cPt.fY); 226 #endif 227 return true; 228 } 229 230 static int VerticalIntersect(const SkDCubic& c, double axisIntercept, double roots[3]) { 231 double A, B, C, D; 232 SkDCubic::Coefficients(&c[0].fX, &A, &B, &C, &D); 233 D -= axisIntercept; 234 int count = SkDCubic::RootsValidT(A, B, C, D, roots); 235 for (int index = 0; index < count; ++index) { 236 SkDPoint calcPt = c.ptAtT(roots[index]); 237 if (!approximately_equal(calcPt.fX, axisIntercept)) { 238 double extremeTs[6]; 239 int extrema = SkDCubic::FindExtrema(&c[0].fX, extremeTs); 240 count = c.searchRoots(extremeTs, extrema, axisIntercept, SkDCubic::kXAxis, roots); 241 break; 242 } 243 } 244 return count; 245 } 246 247 int verticalIntersect(double axisIntercept, double top, double bottom, bool flipped) { 248 addExactVerticalEndPoints(top, bottom, axisIntercept); 249 if (fAllowNear) { 250 addNearVerticalEndPoints(top, bottom, axisIntercept); 251 } 252 double roots[3]; 253 int count = VerticalIntersect(fCubic, axisIntercept, roots); 254 for (int index = 0; index < count; ++index) { 255 double cubicT = roots[index]; 256 SkDPoint pt = { axisIntercept, fCubic.ptAtT(cubicT).fY }; 257 double lineT = (pt.fY - top) / (bottom - top); 258 if (pinTs(&cubicT, &lineT, &pt, kPointInitialized) && uniqueAnswer(cubicT, pt)) { 259 fIntersections->insert(cubicT, lineT, pt); 260 } 261 } 262 if (flipped) { 263 fIntersections->flip(); 264 } 265 checkCoincident(); 266 return fIntersections->used(); 267 } 268 269 protected: 270 271 void addExactEndPoints() { 272 for (int cIndex = 0; cIndex < 4; cIndex += 3) { 273 double lineT = fLine.exactPoint(fCubic[cIndex]); 274 if (lineT < 0) { 275 continue; 276 } 277 double cubicT = (double) (cIndex >> 1); 278 fIntersections->insert(cubicT, lineT, fCubic[cIndex]); 279 } 280 } 281 282 /* Note that this does not look for endpoints of the line that are near the cubic. 283 These points are found later when check ends looks for missing points */ 284 void addNearEndPoints() { 285 for (int cIndex = 0; cIndex < 4; cIndex += 3) { 286 double cubicT = (double) (cIndex >> 1); 287 if (fIntersections->hasT(cubicT)) { 288 continue; 289 } 290 double lineT = fLine.nearPoint(fCubic[cIndex], nullptr); 291 if (lineT < 0) { 292 continue; 293 } 294 fIntersections->insert(cubicT, lineT, fCubic[cIndex]); 295 } 296 this->addLineNearEndPoints(); 297 } 298 299 void addLineNearEndPoints() { 300 for (int lIndex = 0; lIndex < 2; ++lIndex) { 301 double lineT = (double) lIndex; 302 if (fIntersections->hasOppT(lineT)) { 303 continue; 304 } 305 double cubicT = ((SkDCurve*) &fCubic)->nearPoint(SkPath::kCubic_Verb, 306 fLine[lIndex], fLine[!lIndex]); 307 if (cubicT < 0) { 308 continue; 309 } 310 fIntersections->insert(cubicT, lineT, fLine[lIndex]); 311 } 312 } 313 314 void addExactHorizontalEndPoints(double left, double right, double y) { 315 for (int cIndex = 0; cIndex < 4; cIndex += 3) { 316 double lineT = SkDLine::ExactPointH(fCubic[cIndex], left, right, y); 317 if (lineT < 0) { 318 continue; 319 } 320 double cubicT = (double) (cIndex >> 1); 321 fIntersections->insert(cubicT, lineT, fCubic[cIndex]); 322 } 323 } 324 325 void addNearHorizontalEndPoints(double left, double right, double y) { 326 for (int cIndex = 0; cIndex < 4; cIndex += 3) { 327 double cubicT = (double) (cIndex >> 1); 328 if (fIntersections->hasT(cubicT)) { 329 continue; 330 } 331 double lineT = SkDLine::NearPointH(fCubic[cIndex], left, right, y); 332 if (lineT < 0) { 333 continue; 334 } 335 fIntersections->insert(cubicT, lineT, fCubic[cIndex]); 336 } 337 this->addLineNearEndPoints(); 338 } 339 340 void addExactVerticalEndPoints(double top, double bottom, double x) { 341 for (int cIndex = 0; cIndex < 4; cIndex += 3) { 342 double lineT = SkDLine::ExactPointV(fCubic[cIndex], top, bottom, x); 343 if (lineT < 0) { 344 continue; 345 } 346 double cubicT = (double) (cIndex >> 1); 347 fIntersections->insert(cubicT, lineT, fCubic[cIndex]); 348 } 349 } 350 351 void addNearVerticalEndPoints(double top, double bottom, double x) { 352 for (int cIndex = 0; cIndex < 4; cIndex += 3) { 353 double cubicT = (double) (cIndex >> 1); 354 if (fIntersections->hasT(cubicT)) { 355 continue; 356 } 357 double lineT = SkDLine::NearPointV(fCubic[cIndex], top, bottom, x); 358 if (lineT < 0) { 359 continue; 360 } 361 fIntersections->insert(cubicT, lineT, fCubic[cIndex]); 362 } 363 this->addLineNearEndPoints(); 364 } 365 366 double findLineT(double t) { 367 SkDPoint xy = fCubic.ptAtT(t); 368 double dx = fLine[1].fX - fLine[0].fX; 369 double dy = fLine[1].fY - fLine[0].fY; 370 if (fabs(dx) > fabs(dy)) { 371 return (xy.fX - fLine[0].fX) / dx; 372 } 373 return (xy.fY - fLine[0].fY) / dy; 374 } 375 376 bool pinTs(double* cubicT, double* lineT, SkDPoint* pt, PinTPoint ptSet) { 377 if (!approximately_one_or_less(*lineT)) { 378 return false; 379 } 380 if (!approximately_zero_or_more(*lineT)) { 381 return false; 382 } 383 double cT = *cubicT = SkPinT(*cubicT); 384 double lT = *lineT = SkPinT(*lineT); 385 SkDPoint lPt = fLine.ptAtT(lT); 386 SkDPoint cPt = fCubic.ptAtT(cT); 387 if (!lPt.roughlyEqual(cPt)) { 388 return false; 389 } 390 // FIXME: if points are roughly equal but not approximately equal, need to do 391 // a binary search like quad/quad intersection to find more precise t values 392 if (lT == 0 || lT == 1 || (ptSet == kPointUninitialized && cT != 0 && cT != 1)) { 393 *pt = lPt; 394 } else if (ptSet == kPointUninitialized) { 395 *pt = cPt; 396 } 397 SkPoint gridPt = pt->asSkPoint(); 398 if (gridPt == fLine[0].asSkPoint()) { 399 *lineT = 0; 400 } else if (gridPt == fLine[1].asSkPoint()) { 401 *lineT = 1; 402 } 403 if (gridPt == fCubic[0].asSkPoint() && approximately_equal(*cubicT, 0)) { 404 *cubicT = 0; 405 } else if (gridPt == fCubic[3].asSkPoint() && approximately_equal(*cubicT, 1)) { 406 *cubicT = 1; 407 } 408 return true; 409 } 410 411 private: 412 const SkDCubic& fCubic; 413 const SkDLine& fLine; 414 SkIntersections* fIntersections; 415 bool fAllowNear; 416 }; 417 418 int SkIntersections::horizontal(const SkDCubic& cubic, double left, double right, double y, 419 bool flipped) { 420 SkDLine line = {{{ left, y }, { right, y }}}; 421 LineCubicIntersections c(cubic, line, this); 422 return c.horizontalIntersect(y, left, right, flipped); 423 } 424 425 int SkIntersections::vertical(const SkDCubic& cubic, double top, double bottom, double x, 426 bool flipped) { 427 SkDLine line = {{{ x, top }, { x, bottom }}}; 428 LineCubicIntersections c(cubic, line, this); 429 return c.verticalIntersect(x, top, bottom, flipped); 430 } 431 432 int SkIntersections::intersect(const SkDCubic& cubic, const SkDLine& line) { 433 LineCubicIntersections c(cubic, line, this); 434 c.allowNear(fAllowNear); 435 return c.intersect(); 436 } 437 438 int SkIntersections::intersectRay(const SkDCubic& cubic, const SkDLine& line) { 439 LineCubicIntersections c(cubic, line, this); 440 fUsed = c.intersectRay(fT[0]); 441 for (int index = 0; index < fUsed; ++index) { 442 fPt[index] = cubic.ptAtT(fT[0][index]); 443 } 444 return fUsed; 445 } 446 447 // SkDCubic accessors to Intersection utilities 448 449 int SkDCubic::horizontalIntersect(double yIntercept, double roots[3]) const { 450 return LineCubicIntersections::HorizontalIntersect(*this, yIntercept, roots); 451 } 452 453 int SkDCubic::verticalIntersect(double xIntercept, double roots[3]) const { 454 return LineCubicIntersections::VerticalIntersect(*this, xIntercept, roots); 455 } 456