1 // Another approach is to start with the implicit form of one curve and solve 2 // (seek implicit coefficients in QuadraticParameter.cpp 3 // by substituting in the parametric form of the other. 4 // The downside of this approach is that early rejects are difficult to come by. 5 // http://planetmath.org/encyclopedia/GaloisTheoreticDerivationOfTheQuarticFormula.html#step 6 7 #include "SkDQuadImplicit.h" 8 #include "SkIntersections.h" 9 #include "SkPathOpsLine.h" 10 #include "SkQuarticRoot.h" 11 #include "SkTArray.h" 12 #include "SkTSort.h" 13 14 /* given the implicit form 0 = Ax^2 + Bxy + Cy^2 + Dx + Ey + F 15 * and given x = at^2 + bt + c (the parameterized form) 16 * y = dt^2 + et + f 17 * then 18 * 0 = A(at^2+bt+c)(at^2+bt+c)+B(at^2+bt+c)(dt^2+et+f)+C(dt^2+et+f)(dt^2+et+f)+D(at^2+bt+c)+E(dt^2+et+f)+F 19 */ 20 21 static int findRoots(const SkDQuadImplicit& i, const SkDQuad& quad, double roots[4], 22 bool oneHint, bool flip, int firstCubicRoot) { 23 SkDQuad flipped; 24 const SkDQuad& q = flip ? (flipped = quad.flip()) : quad; 25 double a, b, c; 26 SkDQuad::SetABC(&q[0].fX, &a, &b, &c); 27 double d, e, f; 28 SkDQuad::SetABC(&q[0].fY, &d, &e, &f); 29 const double t4 = i.x2() * a * a 30 + i.xy() * a * d 31 + i.y2() * d * d; 32 const double t3 = 2 * i.x2() * a * b 33 + i.xy() * (a * e + b * d) 34 + 2 * i.y2() * d * e; 35 const double t2 = i.x2() * (b * b + 2 * a * c) 36 + i.xy() * (c * d + b * e + a * f) 37 + i.y2() * (e * e + 2 * d * f) 38 + i.x() * a 39 + i.y() * d; 40 const double t1 = 2 * i.x2() * b * c 41 + i.xy() * (c * e + b * f) 42 + 2 * i.y2() * e * f 43 + i.x() * b 44 + i.y() * e; 45 const double t0 = i.x2() * c * c 46 + i.xy() * c * f 47 + i.y2() * f * f 48 + i.x() * c 49 + i.y() * f 50 + i.c(); 51 int rootCount = SkReducedQuarticRoots(t4, t3, t2, t1, t0, oneHint, roots); 52 if (rootCount < 0) { 53 rootCount = SkQuarticRootsReal(firstCubicRoot, t4, t3, t2, t1, t0, roots); 54 } 55 if (flip) { 56 for (int index = 0; index < rootCount; ++index) { 57 roots[index] = 1 - roots[index]; 58 } 59 } 60 return rootCount; 61 } 62 63 static int addValidRoots(const double roots[4], const int count, double valid[4]) { 64 int result = 0; 65 int index; 66 for (index = 0; index < count; ++index) { 67 if (!approximately_zero_or_more(roots[index]) || !approximately_one_or_less(roots[index])) { 68 continue; 69 } 70 double t = 1 - roots[index]; 71 if (approximately_less_than_zero(t)) { 72 t = 0; 73 } else if (approximately_greater_than_one(t)) { 74 t = 1; 75 } 76 valid[result++] = t; 77 } 78 return result; 79 } 80 81 static bool only_end_pts_in_common(const SkDQuad& q1, const SkDQuad& q2) { 82 // the idea here is to see at minimum do a quick reject by rotating all points 83 // to either side of the line formed by connecting the endpoints 84 // if the opposite curves points are on the line or on the other side, the 85 // curves at most intersect at the endpoints 86 for (int oddMan = 0; oddMan < 3; ++oddMan) { 87 const SkDPoint* endPt[2]; 88 for (int opp = 1; opp < 3; ++opp) { 89 int end = oddMan ^ opp; // choose a value not equal to oddMan 90 if (3 == end) { // and correct so that largest value is 1 or 2 91 end = opp; 92 } 93 endPt[opp - 1] = &q1[end]; 94 } 95 double origX = endPt[0]->fX; 96 double origY = endPt[0]->fY; 97 double adj = endPt[1]->fX - origX; 98 double opp = endPt[1]->fY - origY; 99 double sign = (q1[oddMan].fY - origY) * adj - (q1[oddMan].fX - origX) * opp; 100 if (approximately_zero(sign)) { 101 goto tryNextHalfPlane; 102 } 103 for (int n = 0; n < 3; ++n) { 104 double test = (q2[n].fY - origY) * adj - (q2[n].fX - origX) * opp; 105 if (test * sign > 0 && !precisely_zero(test)) { 106 goto tryNextHalfPlane; 107 } 108 } 109 return true; 110 tryNextHalfPlane: 111 ; 112 } 113 return false; 114 } 115 116 // returns false if there's more than one intercept or the intercept doesn't match the point 117 // returns true if the intercept was successfully added or if the 118 // original quads need to be subdivided 119 static bool add_intercept(const SkDQuad& q1, const SkDQuad& q2, double tMin, double tMax, 120 SkIntersections* i, bool* subDivide) { 121 double tMid = (tMin + tMax) / 2; 122 SkDPoint mid = q2.ptAtT(tMid); 123 SkDLine line; 124 line[0] = line[1] = mid; 125 SkDVector dxdy = q2.dxdyAtT(tMid); 126 line[0] -= dxdy; 127 line[1] += dxdy; 128 SkIntersections rootTs; 129 rootTs.allowNear(false); 130 int roots = rootTs.intersect(q1, line); 131 if (roots == 0) { 132 if (subDivide) { 133 *subDivide = true; 134 } 135 return true; 136 } 137 if (roots == 2) { 138 return false; 139 } 140 SkDPoint pt2 = q1.ptAtT(rootTs[0][0]); 141 if (!pt2.approximatelyEqual(mid)) { 142 return false; 143 } 144 i->insertSwap(rootTs[0][0], tMid, pt2); 145 return true; 146 } 147 148 static bool is_linear_inner(const SkDQuad& q1, double t1s, double t1e, const SkDQuad& q2, 149 double t2s, double t2e, SkIntersections* i, bool* subDivide) { 150 SkDQuad hull = q1.subDivide(t1s, t1e); 151 SkDLine line = {{hull[2], hull[0]}}; 152 const SkDLine* testLines[] = { &line, (const SkDLine*) &hull[0], (const SkDLine*) &hull[1] }; 153 const size_t kTestCount = SK_ARRAY_COUNT(testLines); 154 SkSTArray<kTestCount * 2, double, true> tsFound; 155 for (size_t index = 0; index < kTestCount; ++index) { 156 SkIntersections rootTs; 157 rootTs.allowNear(false); 158 int roots = rootTs.intersect(q2, *testLines[index]); 159 for (int idx2 = 0; idx2 < roots; ++idx2) { 160 double t = rootTs[0][idx2]; 161 #if 0 // def SK_DEBUG // FIXME : accurate for error = 16, error of 17.5 seen 162 // {{{136.08723965397621, 1648.2814535211637}, {593.49031197259478, 1190.8784277439891}, {593.49031197259478, 544.0128173828125}}} 163 // {{{-968.181396484375, 544.0128173828125}, {592.2825927734375, 870.552490234375}, {593.435302734375, 557.8828125}}} 164 165 SkDPoint qPt = q2.ptAtT(t); 166 SkDPoint lPt = testLines[index]->ptAtT(rootTs[1][idx2]); 167 SkASSERT(qPt.approximatelyDEqual(lPt)); 168 #endif 169 if (approximately_negative(t - t2s) || approximately_positive(t - t2e)) { 170 continue; 171 } 172 tsFound.push_back(rootTs[0][idx2]); 173 } 174 } 175 int tCount = tsFound.count(); 176 if (tCount <= 0) { 177 return true; 178 } 179 double tMin, tMax; 180 if (tCount == 1) { 181 tMin = tMax = tsFound[0]; 182 } else { 183 SkASSERT(tCount > 1); 184 SkTQSort<double>(tsFound.begin(), tsFound.end() - 1); 185 tMin = tsFound[0]; 186 tMax = tsFound[tsFound.count() - 1]; 187 } 188 SkDPoint end = q2.ptAtT(t2s); 189 bool startInTriangle = hull.pointInHull(end); 190 if (startInTriangle) { 191 tMin = t2s; 192 } 193 end = q2.ptAtT(t2e); 194 bool endInTriangle = hull.pointInHull(end); 195 if (endInTriangle) { 196 tMax = t2e; 197 } 198 int split = 0; 199 SkDVector dxy1, dxy2; 200 if (tMin != tMax || tCount > 2) { 201 dxy2 = q2.dxdyAtT(tMin); 202 for (int index = 1; index < tCount; ++index) { 203 dxy1 = dxy2; 204 dxy2 = q2.dxdyAtT(tsFound[index]); 205 double dot = dxy1.dot(dxy2); 206 if (dot < 0) { 207 split = index - 1; 208 break; 209 } 210 } 211 } 212 if (split == 0) { // there's one point 213 if (add_intercept(q1, q2, tMin, tMax, i, subDivide)) { 214 return true; 215 } 216 i->swap(); 217 return is_linear_inner(q2, tMin, tMax, q1, t1s, t1e, i, subDivide); 218 } 219 // At this point, we have two ranges of t values -- treat each separately at the split 220 bool result; 221 if (add_intercept(q1, q2, tMin, tsFound[split - 1], i, subDivide)) { 222 result = true; 223 } else { 224 i->swap(); 225 result = is_linear_inner(q2, tMin, tsFound[split - 1], q1, t1s, t1e, i, subDivide); 226 } 227 if (add_intercept(q1, q2, tsFound[split], tMax, i, subDivide)) { 228 result = true; 229 } else { 230 i->swap(); 231 result |= is_linear_inner(q2, tsFound[split], tMax, q1, t1s, t1e, i, subDivide); 232 } 233 return result; 234 } 235 236 static double flat_measure(const SkDQuad& q) { 237 SkDVector mid = q[1] - q[0]; 238 SkDVector dxy = q[2] - q[0]; 239 double length = dxy.length(); // OPTIMIZE: get rid of sqrt 240 return fabs(mid.cross(dxy) / length); 241 } 242 243 // FIXME ? should this measure both and then use the quad that is the flattest as the line? 244 static bool is_linear(const SkDQuad& q1, const SkDQuad& q2, SkIntersections* i) { 245 double measure = flat_measure(q1); 246 // OPTIMIZE: (get rid of sqrt) use approximately_zero 247 if (!approximately_zero_sqrt(measure)) { 248 return false; 249 } 250 return is_linear_inner(q1, 0, 1, q2, 0, 1, i, NULL); 251 } 252 253 // FIXME: if flat measure is sufficiently large, then probably the quartic solution failed 254 // avoid imprecision incurred with chopAt 255 static void relaxed_is_linear(const SkDQuad* q1, double s1, double e1, const SkDQuad* q2, 256 double s2, double e2, SkIntersections* i) { 257 double m1 = flat_measure(*q1); 258 double m2 = flat_measure(*q2); 259 i->reset(); 260 const SkDQuad* rounder, *flatter; 261 double sf, midf, ef, sr, er; 262 if (m2 < m1) { 263 rounder = q1; 264 sr = s1; 265 er = e1; 266 flatter = q2; 267 sf = s2; 268 midf = (s2 + e2) / 2; 269 ef = e2; 270 } else { 271 rounder = q2; 272 sr = s2; 273 er = e2; 274 flatter = q1; 275 sf = s1; 276 midf = (s1 + e1) / 2; 277 ef = e1; 278 } 279 bool subDivide = false; 280 is_linear_inner(*flatter, sf, ef, *rounder, sr, er, i, &subDivide); 281 if (subDivide) { 282 relaxed_is_linear(flatter, sf, midf, rounder, sr, er, i); 283 relaxed_is_linear(flatter, midf, ef, rounder, sr, er, i); 284 } 285 if (m2 < m1) { 286 i->swapPts(); 287 } 288 } 289 290 // each time through the loop, this computes values it had from the last loop 291 // if i == j == 1, the center values are still good 292 // otherwise, for i != 1 or j != 1, four of the values are still good 293 // and if i == 1 ^ j == 1, an additional value is good 294 static bool binary_search(const SkDQuad& quad1, const SkDQuad& quad2, double* t1Seed, 295 double* t2Seed, SkDPoint* pt) { 296 double tStep = ROUGH_EPSILON; 297 SkDPoint t1[3], t2[3]; 298 int calcMask = ~0; 299 do { 300 if (calcMask & (1 << 1)) t1[1] = quad1.ptAtT(*t1Seed); 301 if (calcMask & (1 << 4)) t2[1] = quad2.ptAtT(*t2Seed); 302 if (t1[1].approximatelyEqual(t2[1])) { 303 *pt = t1[1]; 304 #if ONE_OFF_DEBUG 305 SkDebugf("%s t1=%1.9g t2=%1.9g (%1.9g,%1.9g) == (%1.9g,%1.9g)\n", __FUNCTION__, 306 t1Seed, t2Seed, t1[1].fX, t1[1].fY, t2[1].fX, t2[1].fY); 307 #endif 308 return true; 309 } 310 if (calcMask & (1 << 0)) t1[0] = quad1.ptAtT(SkTMax(0., *t1Seed - tStep)); 311 if (calcMask & (1 << 2)) t1[2] = quad1.ptAtT(SkTMin(1., *t1Seed + tStep)); 312 if (calcMask & (1 << 3)) t2[0] = quad2.ptAtT(SkTMax(0., *t2Seed - tStep)); 313 if (calcMask & (1 << 5)) t2[2] = quad2.ptAtT(SkTMin(1., *t2Seed + tStep)); 314 double dist[3][3]; 315 // OPTIMIZE: using calcMask value permits skipping some distance calcuations 316 // if prior loop's results are moved to correct slot for reuse 317 dist[1][1] = t1[1].distanceSquared(t2[1]); 318 int best_i = 1, best_j = 1; 319 for (int i = 0; i < 3; ++i) { 320 for (int j = 0; j < 3; ++j) { 321 if (i == 1 && j == 1) { 322 continue; 323 } 324 dist[i][j] = t1[i].distanceSquared(t2[j]); 325 if (dist[best_i][best_j] > dist[i][j]) { 326 best_i = i; 327 best_j = j; 328 } 329 } 330 } 331 if (best_i == 1 && best_j == 1) { 332 tStep /= 2; 333 if (tStep < FLT_EPSILON_HALF) { 334 break; 335 } 336 calcMask = (1 << 0) | (1 << 2) | (1 << 3) | (1 << 5); 337 continue; 338 } 339 if (best_i == 0) { 340 *t1Seed -= tStep; 341 t1[2] = t1[1]; 342 t1[1] = t1[0]; 343 calcMask = 1 << 0; 344 } else if (best_i == 2) { 345 *t1Seed += tStep; 346 t1[0] = t1[1]; 347 t1[1] = t1[2]; 348 calcMask = 1 << 2; 349 } else { 350 calcMask = 0; 351 } 352 if (best_j == 0) { 353 *t2Seed -= tStep; 354 t2[2] = t2[1]; 355 t2[1] = t2[0]; 356 calcMask |= 1 << 3; 357 } else if (best_j == 2) { 358 *t2Seed += tStep; 359 t2[0] = t2[1]; 360 t2[1] = t2[2]; 361 calcMask |= 1 << 5; 362 } 363 } while (true); 364 #if ONE_OFF_DEBUG 365 SkDebugf("%s t1=%1.9g t2=%1.9g (%1.9g,%1.9g) != (%1.9g,%1.9g) %s\n", __FUNCTION__, 366 t1Seed, t2Seed, t1[1].fX, t1[1].fY, t1[2].fX, t1[2].fY); 367 #endif 368 return false; 369 } 370 371 static void lookNearEnd(const SkDQuad& q1, const SkDQuad& q2, int testT, 372 const SkIntersections& orig, bool swap, SkIntersections* i) { 373 if (orig.used() == 1 && orig[!swap][0] == testT) { 374 return; 375 } 376 if (orig.used() == 2 && orig[!swap][1] == testT) { 377 return; 378 } 379 SkDLine tmpLine; 380 int testTIndex = testT << 1; 381 tmpLine[0] = tmpLine[1] = q2[testTIndex]; 382 tmpLine[1].fX += q2[1].fY - q2[testTIndex].fY; 383 tmpLine[1].fY -= q2[1].fX - q2[testTIndex].fX; 384 SkIntersections impTs; 385 impTs.intersectRay(q1, tmpLine); 386 for (int index = 0; index < impTs.used(); ++index) { 387 SkDPoint realPt = impTs.pt(index); 388 if (!tmpLine[0].approximatelyPEqual(realPt)) { 389 continue; 390 } 391 if (swap) { 392 i->insert(testT, impTs[0][index], tmpLine[0]); 393 } else { 394 i->insert(impTs[0][index], testT, tmpLine[0]); 395 } 396 } 397 } 398 399 int SkIntersections::intersect(const SkDQuad& q1, const SkDQuad& q2) { 400 fMax = 4; 401 // if the quads share an end point, check to see if they overlap 402 for (int i1 = 0; i1 < 3; i1 += 2) { 403 for (int i2 = 0; i2 < 3; i2 += 2) { 404 if (q1[i1].asSkPoint() == q2[i2].asSkPoint()) { 405 insert(i1 >> 1, i2 >> 1, q1[i1]); 406 } 407 } 408 } 409 SkASSERT(fUsed < 3); 410 if (only_end_pts_in_common(q1, q2)) { 411 return fUsed; 412 } 413 if (only_end_pts_in_common(q2, q1)) { 414 return fUsed; 415 } 416 // see if either quad is really a line 417 // FIXME: figure out why reduce step didn't find this earlier 418 if (is_linear(q1, q2, this)) { 419 return fUsed; 420 } 421 SkIntersections swapped; 422 swapped.setMax(fMax); 423 if (is_linear(q2, q1, &swapped)) { 424 swapped.swapPts(); 425 *this = swapped; 426 return fUsed; 427 } 428 SkIntersections copyI(*this); 429 lookNearEnd(q1, q2, 0, *this, false, ©I); 430 lookNearEnd(q1, q2, 1, *this, false, ©I); 431 lookNearEnd(q2, q1, 0, *this, true, ©I); 432 lookNearEnd(q2, q1, 1, *this, true, ©I); 433 int innerEqual = 0; 434 if (copyI.fUsed >= 2) { 435 SkASSERT(copyI.fUsed <= 4); 436 double width = copyI[0][1] - copyI[0][0]; 437 int midEnd = 1; 438 for (int index = 2; index < copyI.fUsed; ++index) { 439 double testWidth = copyI[0][index] - copyI[0][index - 1]; 440 if (testWidth <= width) { 441 continue; 442 } 443 midEnd = index; 444 } 445 for (int index = 0; index < 2; ++index) { 446 double testT = (copyI[0][midEnd] * (index + 1) 447 + copyI[0][midEnd - 1] * (2 - index)) / 3; 448 SkDPoint testPt1 = q1.ptAtT(testT); 449 testT = (copyI[1][midEnd] * (index + 1) + copyI[1][midEnd - 1] * (2 - index)) / 3; 450 SkDPoint testPt2 = q2.ptAtT(testT); 451 innerEqual += testPt1.approximatelyEqual(testPt2); 452 } 453 } 454 bool expectCoincident = copyI.fUsed >= 2 && innerEqual == 2; 455 if (expectCoincident) { 456 reset(); 457 insertCoincident(copyI[0][0], copyI[1][0], copyI.fPt[0]); 458 int last = copyI.fUsed - 1; 459 insertCoincident(copyI[0][last], copyI[1][last], copyI.fPt[last]); 460 return fUsed; 461 } 462 SkDQuadImplicit i1(q1); 463 SkDQuadImplicit i2(q2); 464 int index; 465 bool flip1 = q1[2] == q2[0]; 466 bool flip2 = q1[0] == q2[2]; 467 bool useCubic = q1[0] == q2[0]; 468 double roots1[4]; 469 int rootCount = findRoots(i2, q1, roots1, useCubic, flip1, 0); 470 // OPTIMIZATION: could short circuit here if all roots are < 0 or > 1 471 double roots1Copy[4]; 472 int r1Count = addValidRoots(roots1, rootCount, roots1Copy); 473 SkDPoint pts1[4]; 474 for (index = 0; index < r1Count; ++index) { 475 pts1[index] = q1.ptAtT(roots1Copy[index]); 476 } 477 double roots2[4]; 478 int rootCount2 = findRoots(i1, q2, roots2, useCubic, flip2, 0); 479 double roots2Copy[4]; 480 int r2Count = addValidRoots(roots2, rootCount2, roots2Copy); 481 SkDPoint pts2[4]; 482 for (index = 0; index < r2Count; ++index) { 483 pts2[index] = q2.ptAtT(roots2Copy[index]); 484 } 485 if (r1Count == r2Count && r1Count <= 1) { 486 if (r1Count == 1 && used() == 0) { 487 if (pts1[0].approximatelyEqual(pts2[0])) { 488 insert(roots1Copy[0], roots2Copy[0], pts1[0]); 489 } else { 490 // find intersection by chasing t 491 if (binary_search(q1, q2, roots1Copy, roots2Copy, pts1)) { 492 insert(roots1Copy[0], roots2Copy[0], pts1[0]); 493 } 494 } 495 } 496 return fUsed; 497 } 498 int closest[4]; 499 double dist[4]; 500 bool foundSomething = false; 501 for (index = 0; index < r1Count; ++index) { 502 dist[index] = DBL_MAX; 503 closest[index] = -1; 504 for (int ndex2 = 0; ndex2 < r2Count; ++ndex2) { 505 if (!pts2[ndex2].approximatelyEqual(pts1[index])) { 506 continue; 507 } 508 double dx = pts2[ndex2].fX - pts1[index].fX; 509 double dy = pts2[ndex2].fY - pts1[index].fY; 510 double distance = dx * dx + dy * dy; 511 if (dist[index] <= distance) { 512 continue; 513 } 514 for (int outer = 0; outer < index; ++outer) { 515 if (closest[outer] != ndex2) { 516 continue; 517 } 518 if (dist[outer] < distance) { 519 goto next; 520 } 521 closest[outer] = -1; 522 } 523 dist[index] = distance; 524 closest[index] = ndex2; 525 foundSomething = true; 526 next: 527 ; 528 } 529 } 530 if (r1Count && r2Count && !foundSomething) { 531 relaxed_is_linear(&q1, 0, 1, &q2, 0, 1, this); 532 return fUsed; 533 } 534 int used = 0; 535 do { 536 double lowest = DBL_MAX; 537 int lowestIndex = -1; 538 for (index = 0; index < r1Count; ++index) { 539 if (closest[index] < 0) { 540 continue; 541 } 542 if (roots1Copy[index] < lowest) { 543 lowestIndex = index; 544 lowest = roots1Copy[index]; 545 } 546 } 547 if (lowestIndex < 0) { 548 break; 549 } 550 insert(roots1Copy[lowestIndex], roots2Copy[closest[lowestIndex]], 551 pts1[lowestIndex]); 552 closest[lowestIndex] = -1; 553 } while (++used < r1Count); 554 return fUsed; 555 } 556
