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