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