The Tor Browser / file comparison

comparison: gfx/skia/trunk/src/pathops/SkDQuadIntersection.cpp

gfx/skia/trunk/src/pathops/SkDQuadIntersection.cpp

branch
TOR_BUG_3246
changeset 7
129ffea94266
equal deleted inserted replaced
-1:000000000000 0:d0fe6163e3a6
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.approximatelyEqual(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.approximatelyPEqual(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 // avoid imprecision incurred with chopAt
253 static void relaxed_is_linear(const SkDQuad* q1, double s1, double e1, const SkDQuad* q2,
254 double s2, double e2, SkIntersections* i) {
255 double m1 = flat_measure(*q1);
256 double m2 = flat_measure(*q2);
257 i->reset();
258 const SkDQuad* rounder, *flatter;
259 double sf, midf, ef, sr, er;
260 if (m2 < m1) {
261 rounder = q1;
262 sr = s1;
263 er = e1;
264 flatter = q2;
265 sf = s2;
266 midf = (s2 + e2) / 2;
267 ef = e2;
268 } else {
269 rounder = q2;
270 sr = s2;
271 er = e2;
272 flatter = q1;
273 sf = s1;
274 midf = (s1 + e1) / 2;
275 ef = e1;
276 }
277 bool subDivide = false;
278 is_linear_inner(*flatter, sf, ef, *rounder, sr, er, i, &subDivide);
279 if (subDivide) {
280 relaxed_is_linear(flatter, sf, midf, rounder, sr, er, i);
281 relaxed_is_linear(flatter, midf, ef, rounder, sr, er, i);
282 }
283 if (m2 < m1) {
284 i->swapPts();
285 }
286 }
287
288 // each time through the loop, this computes values it had from the last loop
289 // if i == j == 1, the center values are still good
290 // otherwise, for i != 1 or j != 1, four of the values are still good
291 // and if i == 1 ^ j == 1, an additional value is good
292 static bool binary_search(const SkDQuad& quad1, const SkDQuad& quad2, double* t1Seed,
293 double* t2Seed, SkDPoint* pt) {
294 double tStep = ROUGH_EPSILON;
295 SkDPoint t1[3], t2[3];
296 int calcMask = ~0;
297 do {
298 if (calcMask & (1 << 1)) t1[1] = quad1.ptAtT(*t1Seed);
299 if (calcMask & (1 << 4)) t2[1] = quad2.ptAtT(*t2Seed);
300 if (t1[1].approximatelyEqual(t2[1])) {
301 *pt = t1[1];
302 #if ONE_OFF_DEBUG
303 SkDebugf("%s t1=%1.9g t2=%1.9g (%1.9g,%1.9g) == (%1.9g,%1.9g)\n", __FUNCTION__,
304 t1Seed, t2Seed, t1[1].fX, t1[1].fY, t2[1].fX, t2[1].fY);
305 #endif
306 return true;
307 }
308 if (calcMask & (1 << 0)) t1[0] = quad1.ptAtT(*t1Seed - tStep);
309 if (calcMask & (1 << 2)) t1[2] = quad1.ptAtT(*t1Seed + tStep);
310 if (calcMask & (1 << 3)) t2[0] = quad2.ptAtT(*t2Seed - tStep);
311 if (calcMask & (1 << 5)) t2[2] = quad2.ptAtT(*t2Seed + tStep);
312 double dist[3][3];
313 // OPTIMIZE: using calcMask value permits skipping some distance calcuations
314 // if prior loop's results are moved to correct slot for reuse
315 dist[1][1] = t1[1].distanceSquared(t2[1]);
316 int best_i = 1, best_j = 1;
317 for (int i = 0; i < 3; ++i) {
318 for (int j = 0; j < 3; ++j) {
319 if (i == 1 && j == 1) {
320 continue;
321 }
322 dist[i][j] = t1[i].distanceSquared(t2[j]);
323 if (dist[best_i][best_j] > dist[i][j]) {
324 best_i = i;
325 best_j = j;
326 }
327 }
328 }
329 if (best_i == 1 && best_j == 1) {
330 tStep /= 2;
331 if (tStep < FLT_EPSILON_HALF) {
332 break;
333 }
334 calcMask = (1 << 0) | (1 << 2) | (1 << 3) | (1 << 5);
335 continue;
336 }
337 if (best_i == 0) {
338 *t1Seed -= tStep;
339 t1[2] = t1[1];
340 t1[1] = t1[0];
341 calcMask = 1 << 0;
342 } else if (best_i == 2) {
343 *t1Seed += tStep;
344 t1[0] = t1[1];
345 t1[1] = t1[2];
346 calcMask = 1 << 2;
347 } else {
348 calcMask = 0;
349 }
350 if (best_j == 0) {
351 *t2Seed -= tStep;
352 t2[2] = t2[1];
353 t2[1] = t2[0];
354 calcMask |= 1 << 3;
355 } else if (best_j == 2) {
356 *t2Seed += tStep;
357 t2[0] = t2[1];
358 t2[1] = t2[2];
359 calcMask |= 1 << 5;
360 }
361 } while (true);
362 #if ONE_OFF_DEBUG
363 SkDebugf("%s t1=%1.9g t2=%1.9g (%1.9g,%1.9g) != (%1.9g,%1.9g) %s\n", __FUNCTION__,
364 t1Seed, t2Seed, t1[1].fX, t1[1].fY, t1[2].fX, t1[2].fY);
365 #endif
366 return false;
367 }
368
369 static void lookNearEnd(const SkDQuad& q1, const SkDQuad& q2, int testT,
370 const SkIntersections& orig, bool swap, SkIntersections* i) {
371 if (orig.used() == 1 && orig[!swap][0] == testT) {
372 return;
373 }
374 if (orig.used() == 2 && orig[!swap][1] == testT) {
375 return;
376 }
377 SkDLine tmpLine;
378 int testTIndex = testT << 1;
379 tmpLine[0] = tmpLine[1] = q2[testTIndex];
380 tmpLine[1].fX += q2[1].fY - q2[testTIndex].fY;
381 tmpLine[1].fY -= q2[1].fX - q2[testTIndex].fX;
382 SkIntersections impTs;
383 impTs.intersectRay(q1, tmpLine);
384 for (int index = 0; index < impTs.used(); ++index) {
385 SkDPoint realPt = impTs.pt(index);
386 if (!tmpLine[0].approximatelyEqual(realPt)) {
387 continue;
388 }
389 if (swap) {
390 i->insert(testT, impTs[0][index], tmpLine[0]);
391 } else {
392 i->insert(impTs[0][index], testT, tmpLine[0]);
393 }
394 }
395 }
396
397 int SkIntersections::intersect(const SkDQuad& q1, const SkDQuad& q2) {
398 fMax = 4;
399 // if the quads share an end point, check to see if they overlap
400 for (int i1 = 0; i1 < 3; i1 += 2) {
401 for (int i2 = 0; i2 < 3; i2 += 2) {
402 if (q1[i1].asSkPoint() == q2[i2].asSkPoint()) {
403 insert(i1 >> 1, i2 >> 1, q1[i1]);
404 }
405 }
406 }
407 SkASSERT(fUsed < 3);
408 if (only_end_pts_in_common(q1, q2)) {
409 return fUsed;
410 }
411 if (only_end_pts_in_common(q2, q1)) {
412 return fUsed;
413 }
414 // see if either quad is really a line
415 // FIXME: figure out why reduce step didn't find this earlier
416 if (is_linear(q1, q2, this)) {
417 return fUsed;
418 }
419 SkIntersections swapped;
420 swapped.setMax(fMax);
421 if (is_linear(q2, q1, &swapped)) {
422 swapped.swapPts();
423 set(swapped);
424 return fUsed;
425 }
426 SkIntersections copyI(*this);
427 lookNearEnd(q1, q2, 0, *this, false, &copyI);
428 lookNearEnd(q1, q2, 1, *this, false, &copyI);
429 lookNearEnd(q2, q1, 0, *this, true, &copyI);
430 lookNearEnd(q2, q1, 1, *this, true, &copyI);
431 int innerEqual = 0;
432 if (copyI.fUsed >= 2) {
433 SkASSERT(copyI.fUsed <= 4);
434 double width = copyI[0][1] - copyI[0][0];
435 int midEnd = 1;
436 for (int index = 2; index < copyI.fUsed; ++index) {
437 double testWidth = copyI[0][index] - copyI[0][index - 1];
438 if (testWidth <= width) {
439 continue;
440 }
441 midEnd = index;
442 }
443 for (int index = 0; index < 2; ++index) {
444 double testT = (copyI[0][midEnd] * (index + 1)
445 + copyI[0][midEnd - 1] * (2 - index)) / 3;
446 SkDPoint testPt1 = q1.ptAtT(testT);
447 testT = (copyI[1][midEnd] * (index + 1) + copyI[1][midEnd - 1] * (2 - index)) / 3;
448 SkDPoint testPt2 = q2.ptAtT(testT);
449 innerEqual += testPt1.approximatelyEqual(testPt2);
450 }
451 }
452 bool expectCoincident = copyI.fUsed >= 2 && innerEqual == 2;
453 if (expectCoincident) {
454 reset();
455 insertCoincident(copyI[0][0], copyI[1][0], copyI.fPt[0]);
456 int last = copyI.fUsed - 1;
457 insertCoincident(copyI[0][last], copyI[1][last], copyI.fPt[last]);
458 return fUsed;
459 }
460 SkDQuadImplicit i1(q1);
461 SkDQuadImplicit i2(q2);
462 int index;
463 bool flip1 = q1[2] == q2[0];
464 bool flip2 = q1[0] == q2[2];
465 bool useCubic = q1[0] == q2[0];
466 double roots1[4];
467 int rootCount = findRoots(i2, q1, roots1, useCubic, flip1, 0);
468 // OPTIMIZATION: could short circuit here if all roots are < 0 or > 1
469 double roots1Copy[4];
470 int r1Count = addValidRoots(roots1, rootCount, roots1Copy);
471 SkDPoint pts1[4];
472 for (index = 0; index < r1Count; ++index) {
473 pts1[index] = q1.ptAtT(roots1Copy[index]);
474 }
475 double roots2[4];
476 int rootCount2 = findRoots(i1, q2, roots2, useCubic, flip2, 0);
477 double roots2Copy[4];
478 int r2Count = addValidRoots(roots2, rootCount2, roots2Copy);
479 SkDPoint pts2[4];
480 for (index = 0; index < r2Count; ++index) {
481 pts2[index] = q2.ptAtT(roots2Copy[index]);
482 }
483 if (r1Count == r2Count && r1Count <= 1) {
484 if (r1Count == 1 && used() == 0) {
485 if (pts1[0].approximatelyEqual(pts2[0])) {
486 insert(roots1Copy[0], roots2Copy[0], pts1[0]);
487 } else if (pts1[0].moreRoughlyEqual(pts2[0])) {
488 // experiment: try to find intersection by chasing t
489 if (binary_search(q1, q2, roots1Copy, roots2Copy, pts1)) {
490 insert(roots1Copy[0], roots2Copy[0], pts1[0]);
491 }
492 }
493 }
494 return fUsed;
495 }
496 int closest[4];
497 double dist[4];
498 bool foundSomething = false;
499 for (index = 0; index < r1Count; ++index) {
500 dist[index] = DBL_MAX;
501 closest[index] = -1;
502 for (int ndex2 = 0; ndex2 < r2Count; ++ndex2) {
503 if (!pts2[ndex2].approximatelyEqual(pts1[index])) {
504 continue;
505 }
506 double dx = pts2[ndex2].fX - pts1[index].fX;
507 double dy = pts2[ndex2].fY - pts1[index].fY;
508 double distance = dx * dx + dy * dy;
509 if (dist[index] <= distance) {
510 continue;
511 }
512 for (int outer = 0; outer < index; ++outer) {
513 if (closest[outer] != ndex2) {
514 continue;
515 }
516 if (dist[outer] < distance) {
517 goto next;
518 }
519 closest[outer] = -1;
520 }
521 dist[index] = distance;
522 closest[index] = ndex2;
523 foundSomething = true;
524 next:
525 ;
526 }
527 }
528 if (r1Count && r2Count && !foundSomething) {
529 relaxed_is_linear(&q1, 0, 1, &q2, 0, 1, this);
530 return fUsed;
531 }
532 int used = 0;
533 do {
534 double lowest = DBL_MAX;
535 int lowestIndex = -1;
536 for (index = 0; index < r1Count; ++index) {
537 if (closest[index] < 0) {
538 continue;
539 }
540 if (roots1Copy[index] < lowest) {
541 lowestIndex = index;
542 lowest = roots1Copy[index];
543 }
544 }
545 if (lowestIndex < 0) {
546 break;
547 }
548 insert(roots1Copy[lowestIndex], roots2Copy[closest[lowestIndex]],
549 pts1[lowestIndex]);
550 closest[lowestIndex] = -1;
551 } while (++used < r1Count);
552 return fUsed;
553 }

Mercurial Repository: The Tor Browser

Europalab SCM Repository Server

mercurial