Add base types for path ops
Paths contain lines, quads, and cubics, which are
collectively curves.
To work with path intersections, intermediary curves
are constructed. For now, those intermediates use
doubles to guarantee sufficient precision.
The DVector, DPoint, DLine, DQuad, and DCubic
structs encapsulate these intermediate curves.
The DRect and DTriangle structs are created to
describe intersectable areas of interest.
The Bounds struct inherits from SkRect to create
a SkScalar-based rectangle that intersects shared
edges.
This also includes common math equalities and
debugging that the remainder of path ops builds on,
as well as a temporary top-level interface in
include/pathops/SkPathOps.h.
Review URL: https://codereview.chromium.org/12827020
git-svn-id: http://skia.googlecode.com/svn/trunk@8551 2bbb7eff-a529-9590-31e7-b0007b416f81
diff --git a/src/pathops/SkDQuadIntersection.cpp b/src/pathops/SkDQuadIntersection.cpp
new file mode 100644
index 0000000..d8e3f20
--- /dev/null
+++ b/src/pathops/SkDQuadIntersection.cpp
@@ -0,0 +1,496 @@
+// Another approach is to start with the implicit form of one curve and solve
+// (seek implicit coefficients in QuadraticParameter.cpp
+// by substituting in the parametric form of the other.
+// The downside of this approach is that early rejects are difficult to come by.
+// http://planetmath.org/encyclopedia/GaloisTheoreticDerivationOfTheQuarticFormula.html#step
+
+
+#include "SkDQuadImplicit.h"
+#include "SkIntersections.h"
+#include "SkPathOpsLine.h"
+#include "SkQuarticRoot.h"
+#include "SkTDArray.h"
+#include "TSearch.h"
+
+/* given the implicit form 0 = Ax^2 + Bxy + Cy^2 + Dx + Ey + F
+ * and given x = at^2 + bt + c (the parameterized form)
+ * y = dt^2 + et + f
+ * then
+ * 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
+ */
+
+static int findRoots(const SkDQuadImplicit& i, const SkDQuad& q2, double roots[4],
+ bool oneHint, int firstCubicRoot) {
+ double a, b, c;
+ SkDQuad::SetABC(&q2[0].fX, &a, &b, &c);
+ double d, e, f;
+ SkDQuad::SetABC(&q2[0].fY, &d, &e, &f);
+ const double t4 = i.x2() * a * a
+ + i.xy() * a * d
+ + i.y2() * d * d;
+ const double t3 = 2 * i.x2() * a * b
+ + i.xy() * (a * e + b * d)
+ + 2 * i.y2() * d * e;
+ const double t2 = i.x2() * (b * b + 2 * a * c)
+ + i.xy() * (c * d + b * e + a * f)
+ + i.y2() * (e * e + 2 * d * f)
+ + i.x() * a
+ + i.y() * d;
+ const double t1 = 2 * i.x2() * b * c
+ + i.xy() * (c * e + b * f)
+ + 2 * i.y2() * e * f
+ + i.x() * b
+ + i.y() * e;
+ const double t0 = i.x2() * c * c
+ + i.xy() * c * f
+ + i.y2() * f * f
+ + i.x() * c
+ + i.y() * f
+ + i.c();
+ int rootCount = SkReducedQuarticRoots(t4, t3, t2, t1, t0, oneHint, roots);
+ if (rootCount >= 0) {
+ return rootCount;
+ }
+ return SkQuarticRootsReal(firstCubicRoot, t4, t3, t2, t1, t0, roots);
+}
+
+static int addValidRoots(const double roots[4], const int count, double valid[4]) {
+ int result = 0;
+ int index;
+ for (index = 0; index < count; ++index) {
+ if (!approximately_zero_or_more(roots[index]) || !approximately_one_or_less(roots[index])) {
+ continue;
+ }
+ double t = 1 - roots[index];
+ if (approximately_less_than_zero(t)) {
+ t = 0;
+ } else if (approximately_greater_than_one(t)) {
+ t = 1;
+ }
+ valid[result++] = t;
+ }
+ return result;
+}
+
+static bool only_end_pts_in_common(const SkDQuad& q1, const SkDQuad& q2, SkIntersections* i) {
+// the idea here is to see at minimum do a quick reject by rotating all points
+// to either side of the line formed by connecting the endpoints
+// if the opposite curves points are on the line or on the other side, the
+// curves at most intersect at the endpoints
+ for (int oddMan = 0; oddMan < 3; ++oddMan) {
+ const SkDPoint* endPt[2];
+ for (int opp = 1; opp < 3; ++opp) {
+ int end = oddMan ^ opp;
+ if (end == 3) {
+ end = opp;
+ }
+ endPt[opp - 1] = &q1[end];
+ }
+ double origX = endPt[0]->fX;
+ double origY = endPt[0]->fY;
+ double adj = endPt[1]->fX - origX;
+ double opp = endPt[1]->fY - origY;
+ double sign = (q1[oddMan].fY - origY) * adj - (q1[oddMan].fX - origX) * opp;
+ if (approximately_zero(sign)) {
+ goto tryNextHalfPlane;
+ }
+ for (int n = 0; n < 3; ++n) {
+ double test = (q2[n].fY - origY) * adj - (q2[n].fX - origX) * opp;
+ if (test * sign > 0) {
+ goto tryNextHalfPlane;
+ }
+ }
+ for (int i1 = 0; i1 < 3; i1 += 2) {
+ for (int i2 = 0; i2 < 3; i2 += 2) {
+ if (q1[i1] == q2[i2]) {
+ i->insert(i1 >> 1, i2 >> 1, q1[i1]);
+ }
+ }
+ }
+ SkASSERT(i->used() < 3);
+ return true;
+tryNextHalfPlane:
+ ;
+ }
+ return false;
+}
+
+// returns false if there's more than one intercept or the intercept doesn't match the point
+// returns true if the intercept was successfully added or if the
+// original quads need to be subdivided
+static bool add_intercept(const SkDQuad& q1, const SkDQuad& q2, double tMin, double tMax,
+ SkIntersections* i, bool* subDivide) {
+ double tMid = (tMin + tMax) / 2;
+ SkDPoint mid = q2.xyAtT(tMid);
+ SkDLine line;
+ line[0] = line[1] = mid;
+ SkDVector dxdy = q2.dxdyAtT(tMid);
+ line[0] -= dxdy;
+ line[1] += dxdy;
+ SkIntersections rootTs;
+ int roots = rootTs.intersect(q1, line);
+ if (roots == 0) {
+ if (subDivide) {
+ *subDivide = true;
+ }
+ return true;
+ }
+ if (roots == 2) {
+ return false;
+ }
+ SkDPoint pt2 = q1.xyAtT(rootTs[0][0]);
+ if (!pt2.approximatelyEqualHalf(mid)) {
+ return false;
+ }
+ i->insertSwap(rootTs[0][0], tMid, pt2);
+ return true;
+}
+
+static bool is_linear_inner(const SkDQuad& q1, double t1s, double t1e, const SkDQuad& q2,
+ double t2s, double t2e, SkIntersections* i, bool* subDivide) {
+ SkDQuad hull = q1.subDivide(t1s, t1e);
+ SkDLine line = {{hull[2], hull[0]}};
+ const SkDLine* testLines[] = { &line, (const SkDLine*) &hull[0], (const SkDLine*) &hull[1] };
+ size_t testCount = sizeof(testLines) / sizeof(testLines[0]);
+ SkTDArray<double> tsFound;
+ for (size_t index = 0; index < testCount; ++index) {
+ SkIntersections rootTs;
+ int roots = rootTs.intersect(q2, *testLines[index]);
+ for (int idx2 = 0; idx2 < roots; ++idx2) {
+ double t = rootTs[0][idx2];
+#ifdef SK_DEBUG
+ SkDPoint qPt = q2.xyAtT(t);
+ SkDPoint lPt = testLines[index]->xyAtT(rootTs[1][idx2]);
+ SkASSERT(qPt.approximatelyEqual(lPt));
+#endif
+ if (approximately_negative(t - t2s) || approximately_positive(t - t2e)) {
+ continue;
+ }
+ *tsFound.append() = rootTs[0][idx2];
+ }
+ }
+ int tCount = tsFound.count();
+ if (tCount <= 0) {
+ return true;
+ }
+ double tMin, tMax;
+ if (tCount == 1) {
+ tMin = tMax = tsFound[0];
+ } else if (tCount > 1) {
+ QSort<double>(tsFound.begin(), tsFound.end() - 1);
+ tMin = tsFound[0];
+ tMax = tsFound[tsFound.count() - 1];
+ }
+ SkDPoint end = q2.xyAtT(t2s);
+ bool startInTriangle = hull.pointInHull(end);
+ if (startInTriangle) {
+ tMin = t2s;
+ }
+ end = q2.xyAtT(t2e);
+ bool endInTriangle = hull.pointInHull(end);
+ if (endInTriangle) {
+ tMax = t2e;
+ }
+ int split = 0;
+ SkDVector dxy1, dxy2;
+ if (tMin != tMax || tCount > 2) {
+ dxy2 = q2.dxdyAtT(tMin);
+ for (int index = 1; index < tCount; ++index) {
+ dxy1 = dxy2;
+ dxy2 = q2.dxdyAtT(tsFound[index]);
+ double dot = dxy1.dot(dxy2);
+ if (dot < 0) {
+ split = index - 1;
+ break;
+ }
+ }
+ }
+ if (split == 0) { // there's one point
+ if (add_intercept(q1, q2, tMin, tMax, i, subDivide)) {
+ return true;
+ }
+ i->swap();
+ return is_linear_inner(q2, tMin, tMax, q1, t1s, t1e, i, subDivide);
+ }
+ // At this point, we have two ranges of t values -- treat each separately at the split
+ bool result;
+ if (add_intercept(q1, q2, tMin, tsFound[split - 1], i, subDivide)) {
+ result = true;
+ } else {
+ i->swap();
+ result = is_linear_inner(q2, tMin, tsFound[split - 1], q1, t1s, t1e, i, subDivide);
+ }
+ if (add_intercept(q1, q2, tsFound[split], tMax, i, subDivide)) {
+ result = true;
+ } else {
+ i->swap();
+ result |= is_linear_inner(q2, tsFound[split], tMax, q1, t1s, t1e, i, subDivide);
+ }
+ return result;
+}
+
+static double flat_measure(const SkDQuad& q) {
+ SkDVector mid = q[1] - q[0];
+ SkDVector dxy = q[2] - q[0];
+ double length = dxy.length(); // OPTIMIZE: get rid of sqrt
+ return fabs(mid.cross(dxy) / length);
+}
+
+// FIXME ? should this measure both and then use the quad that is the flattest as the line?
+static bool is_linear(const SkDQuad& q1, const SkDQuad& q2, SkIntersections* i) {
+ double measure = flat_measure(q1);
+ // OPTIMIZE: (get rid of sqrt) use approximately_zero
+ if (!approximately_zero_sqrt(measure)) {
+ return false;
+ }
+ return is_linear_inner(q1, 0, 1, q2, 0, 1, i, NULL);
+}
+
+// FIXME: if flat measure is sufficiently large, then probably the quartic solution failed
+static void relaxed_is_linear(const SkDQuad& q1, const SkDQuad& q2, SkIntersections* i) {
+ double m1 = flat_measure(q1);
+ double m2 = flat_measure(q2);
+#if DEBUG_FLAT_QUADS
+ double min = SkTMin<double>(m1, m2);
+ if (min > 5) {
+ SkDebugf("%s maybe not flat enough.. %1.9g\n", __FUNCTION__, min);
+ }
+#endif
+ i->reset();
+ const SkDQuad& rounder = m2 < m1 ? q1 : q2;
+ const SkDQuad& flatter = m2 < m1 ? q2 : q1;
+ bool subDivide = false;
+ is_linear_inner(flatter, 0, 1, rounder, 0, 1, i, &subDivide);
+ if (subDivide) {
+ SkDQuadPair pair = flatter.chopAt(0.5);
+ SkIntersections firstI, secondI;
+ relaxed_is_linear(pair.first(), rounder, &firstI);
+ for (int index = 0; index < firstI.used(); ++index) {
+ i->insert(firstI[0][index] * 0.5, firstI[1][index], firstI.pt(index));
+ }
+ relaxed_is_linear(pair.second(), rounder, &secondI);
+ for (int index = 0; index < secondI.used(); ++index) {
+ i->insert(0.5 + secondI[0][index] * 0.5, secondI[1][index], secondI.pt(index));
+ }
+ }
+ if (m2 < m1) {
+ i->swapPts();
+ }
+}
+
+// each time through the loop, this computes values it had from the last loop
+// if i == j == 1, the center values are still good
+// otherwise, for i != 1 or j != 1, four of the values are still good
+// and if i == 1 ^ j == 1, an additional value is good
+static bool binary_search(const SkDQuad& quad1, const SkDQuad& quad2, double* t1Seed,
+ double* t2Seed, SkDPoint* pt) {
+ double tStep = ROUGH_EPSILON;
+ SkDPoint t1[3], t2[3];
+ int calcMask = ~0;
+ do {
+ if (calcMask & (1 << 1)) t1[1] = quad1.xyAtT(*t1Seed);
+ if (calcMask & (1 << 4)) t2[1] = quad2.xyAtT(*t2Seed);
+ if (t1[1].approximatelyEqual(t2[1])) {
+ *pt = t1[1];
+ #if ONE_OFF_DEBUG
+ SkDebugf("%s t1=%1.9g t2=%1.9g (%1.9g,%1.9g) == (%1.9g,%1.9g)\n", __FUNCTION__,
+ t1Seed, t2Seed, t1[1].fX, t1[1].fY, t1[2].fX, t1[2].fY);
+ #endif
+ return true;
+ }
+ if (calcMask & (1 << 0)) t1[0] = quad1.xyAtT(*t1Seed - tStep);
+ if (calcMask & (1 << 2)) t1[2] = quad1.xyAtT(*t1Seed + tStep);
+ if (calcMask & (1 << 3)) t2[0] = quad2.xyAtT(*t2Seed - tStep);
+ if (calcMask & (1 << 5)) t2[2] = quad2.xyAtT(*t2Seed + tStep);
+ double dist[3][3];
+ // OPTIMIZE: using calcMask value permits skipping some distance calcuations
+ // if prior loop's results are moved to correct slot for reuse
+ dist[1][1] = t1[1].distanceSquared(t2[1]);
+ int best_i = 1, best_j = 1;
+ for (int i = 0; i < 3; ++i) {
+ for (int j = 0; j < 3; ++j) {
+ if (i == 1 && j == 1) {
+ continue;
+ }
+ dist[i][j] = t1[i].distanceSquared(t2[j]);
+ if (dist[best_i][best_j] > dist[i][j]) {
+ best_i = i;
+ best_j = j;
+ }
+ }
+ }
+ if (best_i == 1 && best_j == 1) {
+ tStep /= 2;
+ if (tStep < FLT_EPSILON_HALF) {
+ break;
+ }
+ calcMask = (1 << 0) | (1 << 2) | (1 << 3) | (1 << 5);
+ continue;
+ }
+ if (best_i == 0) {
+ *t1Seed -= tStep;
+ t1[2] = t1[1];
+ t1[1] = t1[0];
+ calcMask = 1 << 0;
+ } else if (best_i == 2) {
+ *t1Seed += tStep;
+ t1[0] = t1[1];
+ t1[1] = t1[2];
+ calcMask = 1 << 2;
+ } else {
+ calcMask = 0;
+ }
+ if (best_j == 0) {
+ *t2Seed -= tStep;
+ t2[2] = t2[1];
+ t2[1] = t2[0];
+ calcMask |= 1 << 3;
+ } else if (best_j == 2) {
+ *t2Seed += tStep;
+ t2[0] = t2[1];
+ t2[1] = t2[2];
+ calcMask |= 1 << 5;
+ }
+ } while (true);
+#if ONE_OFF_DEBUG
+ SkDebugf("%s t1=%1.9g t2=%1.9g (%1.9g,%1.9g) != (%1.9g,%1.9g) %s\n", __FUNCTION__,
+ t1Seed, t2Seed, t1[1].fX, t1[1].fY, t1[2].fX, t1[2].fY);
+#endif
+ return false;
+}
+
+int SkIntersections::intersect(const SkDQuad& q1, const SkDQuad& q2) {
+ // if the quads share an end point, check to see if they overlap
+
+ if (only_end_pts_in_common(q1, q2, this)) {
+ return fUsed;
+ }
+ if (only_end_pts_in_common(q2, q1, this)) {
+ swapPts();
+ return fUsed;
+ }
+ // see if either quad is really a line
+ if (is_linear(q1, q2, this)) {
+ return fUsed;
+ }
+ if (is_linear(q2, q1, this)) {
+ swapPts();
+ return fUsed;
+ }
+ SkDQuadImplicit i1(q1);
+ SkDQuadImplicit i2(q2);
+ if (i1.match(i2)) {
+ // FIXME: compute T values
+ // compute the intersections of the ends to find the coincident span
+ bool useVertical = fabs(q1[0].fX - q1[2].fX) < fabs(q1[0].fY - q1[2].fY);
+ double t;
+ if ((t = SkIntersections::Axial(q1, q2[0], useVertical)) >= 0) {
+ insertCoincident(t, 0, q2[0]);
+ }
+ if ((t = SkIntersections::Axial(q1, q2[2], useVertical)) >= 0) {
+ insertCoincident(t, 1, q2[2]);
+ }
+ useVertical = fabs(q2[0].fX - q2[2].fX) < fabs(q2[0].fY - q2[2].fY);
+ if ((t = SkIntersections::Axial(q2, q1[0], useVertical)) >= 0) {
+ insertCoincident(0, t, q1[0]);
+ }
+ if ((t = SkIntersections::Axial(q2, q1[2], useVertical)) >= 0) {
+ insertCoincident(1, t, q1[2]);
+ }
+ SkASSERT(coincidentUsed() <= 2);
+ return fUsed;
+ }
+ int index;
+ bool useCubic = q1[0] == q2[0] || q1[0] == q2[2] || q1[2] == q2[0];
+ double roots1[4];
+ int rootCount = findRoots(i2, q1, roots1, useCubic, 0);
+ // OPTIMIZATION: could short circuit here if all roots are < 0 or > 1
+ double roots1Copy[4];
+ int r1Count = addValidRoots(roots1, rootCount, roots1Copy);
+ SkDPoint pts1[4];
+ for (index = 0; index < r1Count; ++index) {
+ pts1[index] = q1.xyAtT(roots1Copy[index]);
+ }
+ double roots2[4];
+ int rootCount2 = findRoots(i1, q2, roots2, useCubic, 0);
+ double roots2Copy[4];
+ int r2Count = addValidRoots(roots2, rootCount2, roots2Copy);
+ SkDPoint pts2[4];
+ for (index = 0; index < r2Count; ++index) {
+ pts2[index] = q2.xyAtT(roots2Copy[index]);
+ }
+ if (r1Count == r2Count && r1Count <= 1) {
+ if (r1Count == 1) {
+ if (pts1[0].approximatelyEqualHalf(pts2[0])) {
+ insert(roots1Copy[0], roots2Copy[0], pts1[0]);
+ } else if (pts1[0].moreRoughlyEqual(pts2[0])) {
+ // experiment: try to find intersection by chasing t
+ rootCount = findRoots(i2, q1, roots1, useCubic, 0);
+ (void) addValidRoots(roots1, rootCount, roots1Copy);
+ rootCount2 = findRoots(i1, q2, roots2, useCubic, 0);
+ (void) addValidRoots(roots2, rootCount2, roots2Copy);
+ if (binary_search(q1, q2, roots1Copy, roots2Copy, pts1)) {
+ insert(roots1Copy[0], roots2Copy[0], pts1[0]);
+ }
+ }
+ }
+ return fUsed;
+ }
+ int closest[4];
+ double dist[4];
+ bool foundSomething = false;
+ for (index = 0; index < r1Count; ++index) {
+ dist[index] = DBL_MAX;
+ closest[index] = -1;
+ for (int ndex2 = 0; ndex2 < r2Count; ++ndex2) {
+ if (!pts2[ndex2].approximatelyEqualHalf(pts1[index])) {
+ continue;
+ }
+ double dx = pts2[ndex2].fX - pts1[index].fX;
+ double dy = pts2[ndex2].fY - pts1[index].fY;
+ double distance = dx * dx + dy * dy;
+ if (dist[index] <= distance) {
+ continue;
+ }
+ for (int outer = 0; outer < index; ++outer) {
+ if (closest[outer] != ndex2) {
+ continue;
+ }
+ if (dist[outer] < distance) {
+ goto next;
+ }
+ closest[outer] = -1;
+ }
+ dist[index] = distance;
+ closest[index] = ndex2;
+ foundSomething = true;
+ next:
+ ;
+ }
+ }
+ if (r1Count && r2Count && !foundSomething) {
+ relaxed_is_linear(q1, q2, this);
+ return fUsed;
+ }
+ int used = 0;
+ do {
+ double lowest = DBL_MAX;
+ int lowestIndex = -1;
+ for (index = 0; index < r1Count; ++index) {
+ if (closest[index] < 0) {
+ continue;
+ }
+ if (roots1Copy[index] < lowest) {
+ lowestIndex = index;
+ lowest = roots1Copy[index];
+ }
+ }
+ if (lowestIndex < 0) {
+ break;
+ }
+ insert(roots1Copy[lowestIndex], roots2Copy[closest[lowestIndex]],
+ pts1[lowestIndex]);
+ closest[lowestIndex] = -1;
+ } while (++used < r1Count);
+ return fUsed;
+}