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