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