shape ops work in progress
add quartic solution for intersecting quadratics
git-svn-id: http://skia.googlecode.com/svn/trunk@5541 2bbb7eff-a529-9590-31e7-b0007b416f81
diff --git a/experimental/Intersection/CubeRoot.cpp b/experimental/Intersection/CubeRoot.cpp
index 7610e82..82d2732 100644
--- a/experimental/Intersection/CubeRoot.cpp
+++ b/experimental/Intersection/CubeRoot.cpp
@@ -374,7 +374,14 @@
#endif
double cube_root(double x) {
- return halley_cbrt3d(x);
+ if (approximately_zero(x)) {
+ return 0;
+ }
+ double result = halley_cbrt3d(fabs(x));
+ if (x < 0) {
+ result = -result;
+ }
+ return result;
}
#if TEST_ALTERNATIVES
diff --git a/experimental/Intersection/CubicUtilities.h b/experimental/Intersection/CubicUtilities.h
index 0b014ac..ffaeeb2 100644
--- a/experimental/Intersection/CubicUtilities.h
+++ b/experimental/Intersection/CubicUtilities.h
@@ -4,6 +4,9 @@
* Use of this source code is governed by a BSD-style license that can be
* found in the LICENSE file.
*/
+#if !defined CUBIC_UTILITIES_H
+#define CUBIC_UTILITIES_H
+
#include "DataTypes.h"
double cube_root(double x);
@@ -16,3 +19,5 @@
bool rotate(const Cubic& cubic, int zero, int index, Cubic& rotPath);
double secondDerivativeAtT(const double* cubic, double t);
void xy_at_t(const Cubic& , double t, double& x, double& y);
+
+#endif
diff --git a/experimental/Intersection/CurveIntersection.h b/experimental/Intersection/CurveIntersection.h
index cd2153a..664145b 100644
--- a/experimental/Intersection/CurveIntersection.h
+++ b/experimental/Intersection/CurveIntersection.h
@@ -53,6 +53,11 @@
int intersect(const Cubic& cubic, const _Line& line, double cRange[3], double lRange[3]);
bool intersect(const Quadratic& q1, const Quadratic& q2, Intersections& );
int intersect(const Quadratic& quad, const _Line& line, Intersections& );
+// the following flavor uses the implicit form instead of convex hulls
+bool intersect2(const Quadratic& q1, const Quadratic& q2, Intersections& i);
+int intersectRay(const Quadratic& quad, const _Line& line, Intersections& i);
+
+
bool isLinear(const Quadratic& quad, int startIndex, int endIndex);
bool isLinear(const Cubic& cubic, int startIndex, int endIndex);
double leftMostT(const Cubic& , double startT, double endT);
diff --git a/experimental/Intersection/DataTypes.h b/experimental/Intersection/DataTypes.h
index f4ef1ab..dc61699 100644
--- a/experimental/Intersection/DataTypes.h
+++ b/experimental/Intersection/DataTypes.h
@@ -129,6 +129,17 @@
return x > -FLT_EPSILON;
}
+inline bool approximately_between(double a, double b, double c) {
+ assert(a <= c);
+ return a <= c ? approximately_negative(a - b) && approximately_negative(b - c)
+ : approximately_negative(b - a) && approximately_negative(c - b);
+}
+
+// returns true if (a <= b <= c) || (a >= b >= c)
+inline bool between(double a, double b, double c) {
+ assert(((a <= b && b <= c) || (a >= b && b >= c)) == ((a - b) * (c - b) <= 0));
+ return (a - b) * (c - b) <= 0;
+}
#endif
@@ -179,6 +190,12 @@
bottom = pt.y;
}
}
+
+ // FIXME: used by debugging only ?
+ bool contains(const _Point& pt) {
+ return approximately_between(left, pt.x, right)
+ && approximately_between(top, pt.y, bottom);
+ }
void set(const _Point& pt) {
left = right = pt.x;
diff --git a/experimental/Intersection/EdgeMain.cpp b/experimental/Intersection/EdgeMain.cpp
index c210beb..442c629 100644
--- a/experimental/Intersection/EdgeMain.cpp
+++ b/experimental/Intersection/EdgeMain.cpp
@@ -7,7 +7,7 @@
#include "Intersection_Tests.h"
-int main(int /*argc*/, char* /*argv*/) {
+int main(int /*argc*/, char** /*argv*/) {
Intersection_Tests();
return 0;
}
diff --git a/experimental/Intersection/Intersection_Tests.cpp b/experimental/Intersection/Intersection_Tests.cpp
index 58d7cc4..36582fe 100644
--- a/experimental/Intersection/Intersection_Tests.cpp
+++ b/experimental/Intersection/Intersection_Tests.cpp
@@ -14,6 +14,7 @@
void Intersection_Tests() {
int testsRun = 0;
+ QuarticRoot_Test();
// QuadraticIntersection_Test();
SimplifyNew_Test();
Simplify4x4QuadraticsThreaded_Test(testsRun);
diff --git a/experimental/Intersection/Intersection_Tests.h b/experimental/Intersection/Intersection_Tests.h
index 80c8418..6e6aa89 100644
--- a/experimental/Intersection/Intersection_Tests.h
+++ b/experimental/Intersection/Intersection_Tests.h
@@ -42,3 +42,4 @@
void QuadraticIntersection_Test();
void QuadraticParameterization_Test();
void QuadraticReduceOrder_Test();
+void QuarticRoot_Test();
diff --git a/experimental/Intersection/Intersections.cpp b/experimental/Intersection/Intersections.cpp
new file mode 100644
index 0000000..089a432
--- /dev/null
+++ b/experimental/Intersection/Intersections.cpp
@@ -0,0 +1,16 @@
+/*
+ * Copyright 2012 Google Inc.
+ *
+ * Use of this source code is governed by a BSD-style license that can be
+ * found in the LICENSE file.
+ */
+
+#include "DataTypes.h"
+#include "Intersections.h"
+
+void Intersections::cleanUp() {
+ assert(fCoincidentUsed);
+ assert(fUsed);
+ // find any entries in fT that could be part of the coincident range
+
+}
diff --git a/experimental/Intersection/Intersections.h b/experimental/Intersection/Intersections.h
index 57cebf7..464e1a8 100644
--- a/experimental/Intersection/Intersections.h
+++ b/experimental/Intersection/Intersections.h
@@ -7,10 +7,13 @@
#ifndef Intersections_DEFINE
#define Intersections_DEFINE
+#include <algorithm> // for std::min
+
class Intersections {
public:
Intersections()
: fUsed(0)
+ , fUsed2(0)
, fCoincidentUsed(0)
, fSwap(0)
{
@@ -26,27 +29,83 @@
return;
}
}
+ assert(fUsed < 9);
fT[fSwap][fUsed] = one;
fT[fSwap ^ 1][fUsed] = two;
++fUsed;
}
// start if index == 0 : end if index == 1
- void addCoincident(double one, double two, bool cancel) {
+ void addCoincident(double one, double two) {
for (int index = 0; index < fCoincidentUsed; ++index) {
if (approximately_equal(fCoincidentT[fSwap][index], one)
&& approximately_equal(fCoincidentT[fSwap ^ 1][index], two)) {
- if (cancel) {
- --fCoincidentUsed;
- }
return;
}
}
+ assert(fCoincidentUsed < 9);
fCoincidentT[fSwap][fCoincidentUsed] = one;
fCoincidentT[fSwap ^ 1][fCoincidentUsed] = two;
++fCoincidentUsed;
}
+ void addCoincident(double s1, double e1, double s2, double e2) {
+ assert((fCoincidentUsed & 1) != 1);
+ for (int index = 0; index < fCoincidentUsed; index += 2) {
+ double cs1 = fCoincidentT[fSwap][index];
+ double ce1 = fCoincidentT[fSwap][index + 1];
+ bool s1in = approximately_between(cs1, s1, ce1);
+ bool e1in = approximately_between(cs1, e1, ce1);
+ double cs2 = fCoincidentT[fSwap ^ 1][index];
+ double ce2 = fCoincidentT[fSwap ^ 1][index + 1];
+ bool s2in = approximately_between(cs2, s2, ce2);
+ bool e2in = approximately_between(cs2, e2, ce2);
+ if ((s1in | e1in) & (s2in | e2in)) {
+ double lesser1 = std::min(cs1, ce1);
+ index += cs1 > ce1;
+ if (s1in < lesser1) {
+ fCoincidentT[fSwap][index] = s1in;
+ } else if (e1in < lesser1) {
+ fCoincidentT[fSwap][index] = e1in;
+ }
+ index ^= 1;
+ double greater1 = fCoincidentT[fSwap][index];
+ if (s1in > greater1) {
+ fCoincidentT[fSwap][index] = s1in;
+ } else if (e1in > greater1) {
+ fCoincidentT[fSwap][index] = e1in;
+ }
+ index &= ~1;
+ double lesser2 = std::min(cs2, ce2);
+ index += cs2 > ce2;
+ if (s2in < lesser2) {
+ fCoincidentT[fSwap ^ 1][index] = s2in;
+ } else if (e2in < lesser2) {
+ fCoincidentT[fSwap ^ 1][index] = e2in;
+ }
+ index ^= 1;
+ double greater2 = fCoincidentT[fSwap ^ 1][index];
+ if (s2in > greater2) {
+ fCoincidentT[fSwap ^ 1][index] = s2in;
+ } else if (e2in > greater2) {
+ fCoincidentT[fSwap ^ 1][index] = e2in;
+ }
+ return;
+ }
+ }
+ assert(fCoincidentUsed < 9);
+ fCoincidentT[fSwap][fCoincidentUsed] = s1;
+ fCoincidentT[fSwap ^ 1][fCoincidentUsed] = s2;
+ ++fCoincidentUsed;
+ fCoincidentT[fSwap][fCoincidentUsed] = e1;
+ fCoincidentT[fSwap ^ 1][fCoincidentUsed] = e2;
+ ++fCoincidentUsed;
+ }
+
+ // FIXME: this is necessary because curve/curve intersections are noisy
+ // remove once curve/curve intersections are improved
+ void cleanUp();
+
int coincidentUsed() {
return fCoincidentUsed;
}
@@ -59,10 +118,58 @@
fT[fSwap][index] = val;
}
}
+
+ void insert(double one, double two) {
+ assert(fUsed <= 1 || fT[0][0] < fT[0][1]);
+ int index;
+ for (index = 0; index < fUsed; ++index) {
+ if (approximately_equal(fT[0][index], one)
+ && approximately_equal(fT[1][index], two)) {
+ return;
+ }
+ if (fT[0][index] > one) {
+ break;
+ }
+ }
+ assert(fUsed < 9);
+ int remaining = fUsed - index;
+ if (remaining > 0) {
+ memmove(&fT[0][index + 1], &fT[0][index], sizeof(fT[0][0]) * remaining);
+ memmove(&fT[1][index + 1], &fT[1][index], sizeof(fT[1][0]) * remaining);
+ }
+ fT[0][index] = one;
+ fT[1][index] = two;
+ ++fUsed;
+ }
+
+ void insertOne(double t, int side) {
+ int used = side ? fUsed2 : fUsed;
+ assert(used <= 1 || fT[side][0] < fT[side][1]);
+ int index;
+ for (index = 0; index < used; ++index) {
+ if (approximately_equal(fT[side][index], t)) {
+ return;
+ }
+ if (fT[side][index] > t) {
+ break;
+ }
+ }
+ assert(used < 9);
+ int remaining = used - index;
+ if (remaining > 0) {
+ memmove(&fT[side][index + 1], &fT[side][index], sizeof(fT[side][0]) * remaining);
+ }
+ fT[side][index] = t;
+ side ? ++fUsed2 : ++fUsed;
+ }
- bool intersected() {
+ bool intersected() const {
return fUsed > 0;
}
+
+ bool insertBalanced() const {
+ return fUsed == fUsed2;
+ }
void swap() {
fSwap ^= 1;
@@ -79,6 +186,7 @@
double fT[2][9];
double fCoincidentT[2][9];
int fUsed;
+ int fUsed2;
int fCoincidentUsed;
private:
int fSwap;
diff --git a/experimental/Intersection/LineQuadraticIntersection.cpp b/experimental/Intersection/LineQuadraticIntersection.cpp
index 1ad3735..f730250 100644
--- a/experimental/Intersection/LineQuadraticIntersection.cpp
+++ b/experimental/Intersection/LineQuadraticIntersection.cpp
@@ -97,7 +97,7 @@
, intersections(i) {
}
-int intersect() {
+int intersectRay() {
/*
solve by rotating line+quad so line is horizontal, then finding the roots
set up matrix to rotate quad to x-axis
@@ -124,7 +124,11 @@
double C = r[0];
A += C - 2 * B; // A = a - 2*b + c
B -= C; // B = -(b - c)
- int roots = quadraticRoots(A, B, C, intersections.fT[0]);
+ return quadraticRoots(A, B, C, intersections.fT[0]);
+}
+
+int intersect() {
+ int roots = intersectRay();
for (int index = 0; index < roots; ) {
double lineT = findLineT(intersections.fT[0][index]);
if (lineIntersects(lineT, index, roots)) {
@@ -183,18 +187,20 @@
protected:
double findLineT(double t) {
- const double* qPtr;
- const double* lPtr;
- if (moreHorizontal) {
- qPtr = &quad[0].x;
- lPtr = &line[0].x;
- } else {
- qPtr = &quad[0].y;
- lPtr = &line[0].y;
+ double x, y;
+ xy_at_t(quad, t, x, y);
+ if (approximately_equal(x, line[0].x) && approximately_equal(y, line[0].y)) {
+ return 0;
}
- double s = 1 - t;
- double quadVal = qPtr[0] * s * s + 2 * qPtr[2] * s * t + qPtr[4] * t * t;
- return (quadVal - lPtr[0]) / (lPtr[2] - lPtr[0]);
+ if (approximately_equal(x, line[1].x) && approximately_equal(y, line[1].y)) {
+ return 1;
+ }
+ double dx = line[1].x - line[0].x;
+ double dy = line[1].y - line[0].y;
+ if (fabs(dx) > fabs(dy)) {
+ return (x - line[0].x) / dx;
+ }
+ return (y - line[0].y) / dy;
}
bool lineIntersects(double lineT, const int x, int& roots) {
@@ -218,8 +224,6 @@
const Quadratic& quad;
const _Line& line;
Intersections& intersections;
-bool moreHorizontal;
-
};
// utility for pairs of coincident quads
@@ -307,3 +311,8 @@
LineQuadraticIntersections q(quad, line, i);
return q.intersect();
}
+
+int intersectRay(const Quadratic& quad, const _Line& line, Intersections& i) {
+ LineQuadraticIntersections q(quad, line, i);
+ return q.intersectRay();
+}
diff --git a/experimental/Intersection/QuadraticBezierClip.cpp b/experimental/Intersection/QuadraticBezierClip.cpp
index 6df2d1b..6100914 100644
--- a/experimental/Intersection/QuadraticBezierClip.cpp
+++ b/experimental/Intersection/QuadraticBezierClip.cpp
@@ -9,6 +9,8 @@
#include "LineParameters.h"
#include <algorithm> // used for std::swap
+#define DEBUG_BEZIER_CLIP 1
+
// return false if unable to clip (e.g., unable to create implicit line)
// caller should subdivide, or create degenerate if the values are too small
bool bezier_clip(const Quadratic& q1, const Quadratic& q2, double& minT, double& maxT) {
@@ -65,5 +67,17 @@
x_at(distance2y[idx], distance2y[next], top, bottom, flags, minT, maxT);
idx = next;
} while (idx);
+#if DEBUG_BEZIER_CLIP
+ _Rect r1, r2;
+ r1.setBounds(q1);
+ r2.setBounds(q2);
+ _Point testPt = {0.487, 0.337};
+ if (r1.contains(testPt) && r2.contains(testPt)) {
+ printf("%s q1=(%1.9g,%1.9g %1.9g,%1.9g %1.9g,%1.9g)"
+ " q2=(%1.9g,%1.9g %1.9g,%1.9g %1.9g,%1.9g) minT=%1.9g maxT=%1.9g\n",
+ __FUNCTION__, q1[0].x, q1[0].y, q1[1].x, q1[1].y, q1[2].x, q1[2].y,
+ q2[0].x, q2[0].y, q2[1].x, q2[1].y, q2[2].x, q2[2].y, minT, maxT);
+ }
+#endif
return minT < maxT; // returns false if distance shows no intersection
}
diff --git a/experimental/Intersection/QuadraticImplicit.cpp b/experimental/Intersection/QuadraticImplicit.cpp
new file mode 100644
index 0000000..b5efd03
--- /dev/null
+++ b/experimental/Intersection/QuadraticImplicit.cpp
@@ -0,0 +1,132 @@
+// 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 "CurveIntersection.h"
+#include "Intersections.h"
+#include "QuadraticParameterization.h"
+#include "QuarticRoot.h"
+#include "QuadraticUtilities.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 QuadImplicitForm& i, const Quadratic& q2, double roots[4]) {
+ double a, b, c;
+ set_abc(&q2[0].x, a, b, c);
+ double d, e, f;
+ set_abc(&q2[0].y, 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();
+ return quarticRoots(t4, t3, t2, t1, t0, roots);
+}
+
+static void addValidRoots(const double roots[4], const int count, const int side, Intersections& i) {
+ 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;
+ }
+ i.insertOne(t, side);
+ }
+}
+
+bool intersect2(const Quadratic& q1, const Quadratic& q2, Intersections& i) {
+ QuadImplicitForm i1(q1);
+ QuadImplicitForm i2(q2);
+ if (i1.implicit_match(i2)) {
+ // FIXME: compute T values
+ // compute the intersections of the ends to find the coincident span
+ bool useVertical = fabs(q1[0].x - q1[2].x) < fabs(q1[0].y - q1[2].y);
+ double t;
+ if ((t = axialIntersect(q1, q2[0], useVertical)) >= 0) {
+ i.addCoincident(t, 0);
+ }
+ if ((t = axialIntersect(q1, q2[2], useVertical)) >= 0) {
+ i.addCoincident(t, 1);
+ }
+ useVertical = fabs(q2[0].x - q2[2].x) < fabs(q2[0].y - q2[2].y);
+ if ((t = axialIntersect(q2, q1[0], useVertical)) >= 0) {
+ i.addCoincident(0, t);
+ }
+ if ((t = axialIntersect(q2, q1[2], useVertical)) >= 0) {
+ i.addCoincident(1, t);
+ }
+ assert(i.fCoincidentUsed <= 2);
+ return i.fCoincidentUsed > 0;
+ }
+ double roots1[4], roots2[4];
+ int rootCount = findRoots(i2, q1, roots1);
+ // OPTIMIZATION: could short circuit here if all roots are < 0 or > 1
+ int rootCount2 = findRoots(i1, q2, roots2);
+ assert(rootCount == rootCount2);
+ addValidRoots(roots1, rootCount, 0, i);
+ addValidRoots(roots2, rootCount, 1, i);
+ _Point pts[4];
+ bool matches[4];
+ int index;
+ for (index = 0; index < i.fUsed2; ++index) {
+ xy_at_t(q2, i.fT[1][index], pts[index].x, pts[index].y);
+ matches[index] = false;
+ }
+ for (index = 0; index < i.fUsed; ) {
+ _Point xy;
+ xy_at_t(q1, i.fT[0][index], xy.x, xy.y);
+ for (int inner = 0; inner < i.fUsed2; ++inner) {
+ if (approximately_equal(pts[inner].x, xy.x) && approximately_equal(pts[inner].y, xy.y)) {
+ matches[index] = true;
+ goto next;
+ }
+ }
+ if (--i.fUsed > index) {
+ memmove(&i.fT[0][index], &i.fT[0][index + 1], (i.fUsed - index) * sizeof(i.fT[0][0]));
+ continue;
+ }
+ next:
+ ++index;
+ }
+ for (index = 0; index < i.fUsed2; ) {
+ if (!matches[index]) {
+ if (--i.fUsed2 > index) {
+ memmove(&i.fT[1][index], &i.fT[1][index + 1], (i.fUsed2 - index) * sizeof(i.fT[1][0]));
+ continue;
+ }
+ }
+ ++index;
+ }
+ assert(i.insertBalanced());
+ return i.intersected();
+}
diff --git a/experimental/Intersection/QuadraticIntersection.cpp b/experimental/Intersection/QuadraticIntersection.cpp
index 6d77efc..28299dd 100644
--- a/experimental/Intersection/QuadraticIntersection.cpp
+++ b/experimental/Intersection/QuadraticIntersection.cpp
@@ -13,7 +13,9 @@
#include "QuadraticUtilities.h"
#include <algorithm> // for swap
-class QuadraticIntersections : public Intersections {
+static const double tClipLimit = 0.8; // http://cagd.cs.byu.edu/~tom/papers/bezclip.pdf see Multiple intersections
+
+class QuadraticIntersections {
public:
QuadraticIntersections(const Quadratic& q1, const Quadratic& q2, Intersections& i)
@@ -21,7 +23,8 @@
, quad2(q2)
, intersections(i)
, depth(0)
- , splits(0) {
+ , splits(0)
+ , coinMinT1(-1) {
}
bool intersect() {
@@ -128,6 +131,20 @@
std::swap(minT1, minT2);
std::swap(maxT1, maxT2);
}
+ if (coinMinT1 >= 0) {
+ bool earlyExit;
+ if ((earlyExit = coinMaxT1 == minT1)) {
+ coinMaxT1 = maxT1;
+ }
+ if (coinMaxT2 == minT2) {
+ coinMaxT2 = maxT2;
+ return true;
+ }
+ if (earlyExit) {
+ return true;
+ }
+ coinMinT1 = -1;
+ }
// do line/quadratic or even line/line intersection instead
if (treat1AsLine) {
xy_at_t(quad1, minT1, line1[0].x, line1[0].y);
@@ -138,48 +155,66 @@
xy_at_t(quad2, maxT2, line2[1].x, line2[1].y);
}
int pts;
- double smallT, largeT;
+ double smallT1, largeT1, smallT2, largeT2;
if (treat1AsLine & treat2AsLine) {
double t1[2], t2[2];
pts = ::intersect(line1, line2, t1, t2);
- for (int index = 0; index < pts; ++index) {
- smallT = interp(minT1, maxT1, t1[index]);
- largeT = interp(minT2, maxT2, t2[index]);
- if (pts == 2) {
- intersections.addCoincident(smallT, largeT, true);
- } else {
- intersections.add(smallT, largeT);
- }
+ if (pts == 2) {
+ smallT1 = interp(minT1, maxT1, t1[0]);
+ largeT1 = interp(minT2, maxT2, t2[0]);
+ smallT2 = interp(minT1, maxT1, t1[1]);
+ largeT2 = interp(minT2, maxT2, t2[1]);
+ intersections.addCoincident(smallT1, smallT2, largeT1, largeT2);
+ } else {
+ smallT1 = interp(minT1, maxT1, t1[0]);
+ largeT1 = interp(minT2, maxT2, t2[0]);
+ intersections.add(smallT1, largeT1);
}
} else {
Intersections lq;
pts = ::intersect(treat1AsLine ? quad2 : quad1,
treat1AsLine ? line1 : line2, lq);
- bool coincident = false;
if (pts == 2) { // if the line and edge are coincident treat differently
_Point midQuad, midLine;
double midQuadT = (lq.fT[0][0] + lq.fT[0][1]) / 2;
xy_at_t(treat1AsLine ? quad2 : quad1, midQuadT, midQuad.x, midQuad.y);
double lineT = t_at(treat1AsLine ? line1 : line2, midQuad);
xy_at_t(treat1AsLine ? line1 : line2, lineT, midLine.x, midLine.y);
- coincident = approximately_equal(midQuad.x, midLine.x)
- && approximately_equal(midQuad.y, midLine.y);
+ if (approximately_equal(midQuad.x, midLine.x)
+ && approximately_equal(midQuad.y, midLine.y)) {
+ smallT1 = lq.fT[0][0];
+ largeT1 = lq.fT[1][0];
+ smallT2 = lq.fT[0][1];
+ largeT2 = lq.fT[1][1];
+ if (treat2AsLine) {
+ smallT1 = interp(minT1, maxT1, smallT1);
+ smallT2 = interp(minT1, maxT1, smallT2);
+ } else {
+ largeT1 = interp(minT2, maxT2, largeT1);
+ largeT2 = interp(minT2, maxT2, largeT2);
+ }
+ intersections.addCoincident(smallT1, smallT2, largeT1, largeT2);
+ goto setCoinMinMax;
+ }
}
for (int index = 0; index < pts; ++index) {
- smallT = lq.fT[0][index];
- largeT = lq.fT[1][index];
- if (treat1AsLine) {
- smallT = interp(minT1, maxT1, smallT);
+ smallT1 = lq.fT[0][index];
+ largeT1 = lq.fT[1][index];
+ if (treat2AsLine) {
+ smallT1 = interp(minT1, maxT1, smallT1);
} else {
- largeT = interp(minT2, maxT2, largeT);
+ largeT1 = interp(minT2, maxT2, largeT1);
}
- if (coincident) {
- intersections.addCoincident(smallT, largeT, true);
- } else {
- intersections.add(smallT, largeT);
- }
+ intersections.add(smallT1, largeT1);
}
}
+ if (pts > 0) {
+setCoinMinMax:
+ coinMinT1 = minT1;
+ coinMaxT1 = maxT1;
+ coinMinT2 = minT2;
+ coinMaxT2 = maxT2;
+ }
return pts > 0;
}
@@ -210,7 +245,6 @@
private:
-static const double tClipLimit = 0.8; // http://cagd.cs.byu.edu/~tom/papers/bezclip.pdf see Multiple intersections
const Quadratic& quad1;
const Quadratic& quad2;
Intersections& intersections;
@@ -218,8 +252,123 @@
int splits;
double quad1Divisions; // line segments to approximate original within error
double quad2Divisions;
+double coinMinT1; // range of Ts where approximate line intersected curve
+double coinMaxT1;
+double coinMinT2;
+double coinMaxT2;
};
+#include "LineParameters.h"
+
+static void hackToFixPartialCoincidence(const Quadratic& q1, const Quadratic& q2, Intersections& i) {
+ // look to see if non-coincident data basically has unsortable tangents
+
+ // look to see if a point between non-coincident data is on the curve
+ int cIndex;
+ for (int uIndex = 0; uIndex < i.fUsed; ) {
+ double bestDist1 = 1;
+ double bestDist2 = 1;
+ int closest1 = -1;
+ int closest2 = -1;
+ for (cIndex = 0; cIndex < i.fCoincidentUsed; ++cIndex) {
+ double dist = fabs(i.fT[0][uIndex] - i.fCoincidentT[0][cIndex]);
+ if (bestDist1 > dist) {
+ bestDist1 = dist;
+ closest1 = cIndex;
+ }
+ dist = fabs(i.fT[1][uIndex] - i.fCoincidentT[1][cIndex]);
+ if (bestDist2 > dist) {
+ bestDist2 = dist;
+ closest2 = cIndex;
+ }
+ }
+ _Line ends;
+ _Point mid;
+ double t1 = i.fT[0][uIndex];
+ xy_at_t(q1, t1, ends[0].x, ends[0].y);
+ xy_at_t(q1, i.fCoincidentT[0][closest1], ends[1].x, ends[1].y);
+ double midT = (t1 + i.fCoincidentT[0][closest1]) / 2;
+ xy_at_t(q1, midT, mid.x, mid.y);
+ LineParameters params;
+ params.lineEndPoints(ends);
+ double midDist = params.pointDistance(mid);
+ // Note that we prefer to always measure t error, which does not scale,
+ // instead of point error, which is scale dependent. FIXME
+ if (!approximately_zero(midDist)) {
+ ++uIndex;
+ continue;
+ }
+ double t2 = i.fT[1][uIndex];
+ xy_at_t(q2, t2, ends[0].x, ends[0].y);
+ xy_at_t(q2, i.fCoincidentT[1][closest2], ends[1].x, ends[1].y);
+ midT = (t2 + i.fCoincidentT[1][closest2]) / 2;
+ xy_at_t(q2, midT, mid.x, mid.y);
+ params.lineEndPoints(ends);
+ midDist = params.pointDistance(mid);
+ if (!approximately_zero(midDist)) {
+ ++uIndex;
+ continue;
+ }
+ // if both midpoints are close to the line, lengthen coincident span
+ int cEnd = closest1 ^ 1; // assume coincidence always travels in pairs
+ if (!between(i.fCoincidentT[0][cEnd], t1, i.fCoincidentT[0][closest1])) {
+ i.fCoincidentT[0][closest1] = t1;
+ }
+ cEnd = closest2 ^ 1;
+ if (!between(i.fCoincidentT[0][cEnd], t2, i.fCoincidentT[0][closest2])) {
+ i.fCoincidentT[0][closest2] = t2;
+ }
+ int remaining = --i.fUsed - uIndex;
+ if (remaining > 0) {
+ memmove(&i.fT[0][uIndex], &i.fT[0][uIndex + 1], sizeof(i.fT[0][0]) * remaining);
+ memmove(&i.fT[1][uIndex], &i.fT[1][uIndex + 1], sizeof(i.fT[1][0]) * remaining);
+ }
+ }
+ // if coincident data is subjectively a tiny span, replace it with a single point
+ for (cIndex = 0; cIndex < i.fCoincidentUsed; ) {
+ double start1 = i.fCoincidentT[0][cIndex];
+ double end1 = i.fCoincidentT[0][cIndex + 1];
+ _Line ends1;
+ xy_at_t(q1, start1, ends1[0].x, ends1[0].y);
+ xy_at_t(q1, end1, ends1[1].x, ends1[1].y);
+ if (!approximately_equal(ends1[0].x, ends1[1].x) || approximately_equal(ends1[0].y, ends1[1].y)) {
+ cIndex += 2;
+ continue;
+ }
+ double start2 = i.fCoincidentT[1][cIndex];
+ double end2 = i.fCoincidentT[1][cIndex + 1];
+ _Line ends2;
+ xy_at_t(q2, start2, ends2[0].x, ends2[0].y);
+ xy_at_t(q2, end2, ends2[1].x, ends2[1].y);
+ // again, approximately should be used with T values, not points FIXME
+ if (!approximately_equal(ends2[0].x, ends2[1].x) || approximately_equal(ends2[0].y, ends2[1].y)) {
+ cIndex += 2;
+ continue;
+ }
+ if (approximately_less_than_zero(start1) || approximately_less_than_zero(end1)) {
+ start1 = 0;
+ } else if (approximately_greater_than_one(start1) || approximately_greater_than_one(end1)) {
+ start1 = 1;
+ } else {
+ start1 = (start1 + end1) / 2;
+ }
+ if (approximately_less_than_zero(start2) || approximately_less_than_zero(end2)) {
+ start2 = 0;
+ } else if (approximately_greater_than_one(start2) || approximately_greater_than_one(end2)) {
+ start2 = 1;
+ } else {
+ start2 = (start2 + end2) / 2;
+ }
+ i.insert(start1, start2);
+ i.fCoincidentUsed -= 2;
+ int remaining = i.fCoincidentUsed - cIndex;
+ if (remaining > 0) {
+ memmove(&i.fCoincidentT[0][cIndex], &i.fCoincidentT[0][cIndex + 2], sizeof(i.fCoincidentT[0][0]) * remaining);
+ memmove(&i.fCoincidentT[1][cIndex], &i.fCoincidentT[1][cIndex + 2], sizeof(i.fCoincidentT[1][0]) * remaining);
+ }
+ }
+}
+
bool intersect(const Quadratic& q1, const Quadratic& q2, Intersections& i) {
if (implicit_matches(q1, q2)) {
// FIXME: compute T values
@@ -227,64 +376,28 @@
bool useVertical = fabs(q1[0].x - q1[2].x) < fabs(q1[0].y - q1[2].y);
double t;
if ((t = axialIntersect(q1, q2[0], useVertical)) >= 0) {
- i.addCoincident(t, 0, false);
+ i.addCoincident(t, 0);
}
if ((t = axialIntersect(q1, q2[2], useVertical)) >= 0) {
- i.addCoincident(t, 1, false);
+ i.addCoincident(t, 1);
}
useVertical = fabs(q2[0].x - q2[2].x) < fabs(q2[0].y - q2[2].y);
if ((t = axialIntersect(q2, q1[0], useVertical)) >= 0) {
- i.addCoincident(0, t, false);
+ i.addCoincident(0, t);
}
if ((t = axialIntersect(q2, q1[2], useVertical)) >= 0) {
- i.addCoincident(1, t, false);
+ i.addCoincident(1, t);
}
assert(i.fCoincidentUsed <= 2);
return i.fCoincidentUsed > 0;
}
QuadraticIntersections q(q1, q2, i);
- return q.intersect();
+ bool result = q.intersect();
+ // FIXME: partial coincidence detection is currently poor. For now, try
+ // to fix up the data after the fact. In the future, revisit the error
+ // term to try to avoid this kind of result in the first place.
+ if (i.fUsed && i.fCoincidentUsed) {
+ hackToFixPartialCoincidence(q1, q2, i);
+ }
+ return result;
}
-
-
-// Another approach is to start with the implicit form of one curve and solve
-// 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
-/*
-given x^4 + ax^3 + bx^2 + cx + d
-the resolvent cubic is x^3 - 2bx^2 + (b^2 + ac - 4d)x + (c^2 + a^2d - abc)
-use the cubic formula (CubicRoots.cpp) to find the radical expressions t1, t2, and t3.
-
-(x - r1 r2) (x - r3 r4) = x^2 - (t2 + t3 - t1) / 2 x + d
-s = r1*r2 = ((t2 + t3 - t1) + sqrt((t2 + t3 - t1)^2 - 16*d)) / 4
-t = r3*r4 = ((t2 + t3 - t1) - sqrt((t2 + t3 - t1)^2 - 16*d)) / 4
-
-u = r1+r2 = (-a + sqrt(a^2 - 4*t1)) / 2
-v = r3+r4 = (-a - sqrt(a^2 - 4*t1)) / 2
-
-r1 = (u + sqrt(u^2 - 4*s)) / 2
-r2 = (u - sqrt(u^2 - 4*s)) / 2
-r3 = (v + sqrt(v^2 - 4*t)) / 2
-r4 = (v - sqrt(v^2 - 4*t)) / 2
-*/
-
-
-/* square root of complex number
-http://en.wikipedia.org/wiki/Square_root#Square_roots_of_negative_and_complex_numbers
-Algebraic formula
-When the number is expressed using Cartesian coordinates the following formula
- can be used for the principal square root:[5][6]
-
-sqrt(x + iy) = sqrt((r + x) / 2) +/- i*sqrt((r - x) / 2)
-
-where the sign of the imaginary part of the root is taken to be same as the sign
- of the imaginary part of the original number, and
-
-r = abs(x + iy) = sqrt(x^2 + y^2)
-
-is the absolute value or modulus of the original number. The real part of the
-principal value is always non-negative.
-The other square root is simply –1 times the principal square root; in other
-words, the two square roots of a number sum to 0.
- */
diff --git a/experimental/Intersection/QuadraticIntersection_Test.cpp b/experimental/Intersection/QuadraticIntersection_Test.cpp
index 74037be..08736cc 100644
--- a/experimental/Intersection/QuadraticIntersection_Test.cpp
+++ b/experimental/Intersection/QuadraticIntersection_Test.cpp
@@ -28,8 +28,10 @@
printf("[%d] quad2 order=%d\n", (int) index, order2);
}
if (order1 == 3 && order2 == 3) {
- Intersections intersections;
+ Intersections intersections, intersections2;
intersect(reduce1, reduce2, intersections);
+ intersect2(reduce1, reduce2, intersections2);
+ SkASSERT(intersections.used() == intersections2.used());
if (intersections.intersected()) {
for (int pt = 0; pt < intersections.used(); ++pt) {
double tt1 = intersections.fT[0][pt];
@@ -46,6 +48,10 @@
printf("%s [%d,%d] y!= t1=%g (%g,%g) t2=%g (%g,%g)\n",
__FUNCTION__, (int)index, pt, tt1, tx1, ty1, tt2, tx2, ty2);
}
+ tt1 = intersections2.fT[0][pt];
+ SkASSERT(approximately_equal(intersections.fT[0][pt], tt1));
+ tt2 = intersections2.fT[1][pt];
+ SkASSERT(approximately_equal(intersections.fT[1][pt], tt2));
}
}
}
@@ -70,11 +76,12 @@
if (!intersections.intersected()) {
SkDebugf("%s no intersection!\n", __FUNCTION__);
}
+ double tt1, tt2;
for (int pt = 0; pt < intersections.used(); ++pt) {
- double tt1 = intersections.fT[0][pt];
+ tt1 = intersections.fT[0][pt];
double tx1, ty1;
xy_at_t(quad1, tt1, tx1, ty1);
- double tt2 = intersections.fT[1][pt];
+ tt2 = intersections.fT[1][pt];
double tx2, ty2;
xy_at_t(quad2, tt2, tx2, ty2);
if (!approximately_equal(tx1, tx2)) {
@@ -90,6 +97,15 @@
SkDebugf("%s [%d][%d] t1=%1.9g (%1.9g, %1.9g) t2=%1.9g\n", __FUNCTION__,
outer, inner, tt1, tx1, tx2, tt2);
}
+ Intersections intersections2;
+ intersect2(quad1, quad2, intersections2);
+ SkASSERT(intersections.used() == intersections2.used());
+ for (int pt = 0; pt < intersections2.used(); ++pt) {
+ tt1 = intersections2.fT[0][pt];
+ SkASSERT(approximately_equal(intersections.fT[0][pt], tt1));
+ tt2 = intersections2.fT[1][pt];
+ SkASSERT(approximately_equal(intersections.fT[1][pt], tt2));
+ }
}
}
}
@@ -105,19 +121,29 @@
for (size_t testIndex = 0; testIndex < coincidentTestSetCount - 1; testIndex += 2) {
const Quadratic& quad1 = coincidentTestSet[testIndex];
const Quadratic& quad2 = coincidentTestSet[testIndex + 1];
- Intersections intersections;
+ Intersections intersections, intersections2;
intersect(quad1, quad2, intersections);
SkASSERT(intersections.coincidentUsed() == 2);
- for (int pt = 0; pt < intersections.coincidentUsed(); ++pt) {
- double tt1 = intersections.fT[0][pt];
+ int pt;
+ double tt1, tt2;
+ for (pt = 0; pt < intersections.coincidentUsed(); ++pt) {
+ tt1 = intersections.fT[0][pt];
double tx1, ty1;
xy_at_t(quad1, tt1, tx1, ty1);
- double tt2 = intersections.fT[1][pt];
+ tt2 = intersections.fT[1][pt];
double tx2, ty2;
xy_at_t(quad2, tt2, tx2, ty2);
SkDebugf("%s [%d,%d] t1=%g (%g,%g) t2=%g (%g,%g)\n",
__FUNCTION__, (int)testIndex, pt, tt1, tx1, ty1, tt2, tx2, ty2);
}
+ intersect2(quad1, quad2, intersections2);
+ SkASSERT(intersections2.coincidentUsed() == 2);
+ for (pt = 0; pt < intersections2.coincidentUsed(); ++pt) {
+ tt1 = intersections2.fT[0][pt];
+ SkASSERT(approximately_equal(intersections.fT[0][pt], tt1));
+ tt2 = intersections2.fT[1][pt];
+ SkASSERT(approximately_equal(intersections.fT[1][pt], tt2));
+ }
}
}
diff --git a/experimental/Intersection/QuadraticParameterization.cpp b/experimental/Intersection/QuadraticParameterization.cpp
index 7d86f7f..8e7f1a2 100644
--- a/experimental/Intersection/QuadraticParameterization.cpp
+++ b/experimental/Intersection/QuadraticParameterization.cpp
@@ -5,6 +5,7 @@
* found in the LICENSE file.
*/
#include "CurveIntersection.h"
+#include "QuadraticParameterization.h"
#include "QuadraticUtilities.h"
/* from http://tom.cs.byu.edu/~tom/papers/cvgip84.pdf 4.1
@@ -49,19 +50,10 @@
*
*/
-enum {
- xx_coeff,
- xy_coeff,
- yy_coeff,
- x_coeff,
- y_coeff,
- c_coeff,
- coeff_count
-};
static bool straight_forward = true;
-static void implicit_coefficients(const Quadratic& q, double p[coeff_count]) {
+QuadImplicitForm::QuadImplicitForm(const Quadratic& q) {
double a, b, c;
set_abc(&q[0].x, a, b, c);
double d, e, f;
@@ -103,28 +95,30 @@
* lazily compute the coefficients, comparing the easiest to compute first.
* xx and yy first; then xy; and so on.
*/
-bool implicit_matches(const Quadratic& one, const Quadratic& two) {
- double p1[coeff_count]; // a'xx , b'xy , c'yy , d'x , e'y , f
- double p2[coeff_count];
- implicit_coefficients(one, p1);
- implicit_coefficients(two, p2);
+bool QuadImplicitForm::implicit_match(const QuadImplicitForm& p2) const {
int first = 0;
for (int index = 0; index < coeff_count; ++index) {
- if (approximately_zero(p1[index]) || approximately_zero(p2[index])) {
+ if (approximately_zero(p[index]) && approximately_zero(p2.p[index])) {
first += first == index;
continue;
}
if (first == index) {
continue;
}
- if (!approximately_equal(p1[index] * p2[first],
- p1[first] * p2[index])) {
+ if (!approximately_equal(p[index] * p2.p[first],
+ p[first] * p2.p[index])) {
return false;
}
}
return true;
}
+bool implicit_matches(const Quadratic& quad1, const Quadratic& quad2) {
+ QuadImplicitForm i1(quad1); // a'xx , b'xy , c'yy , d'x , e'y , f
+ QuadImplicitForm i2(quad2);
+ return i1.implicit_match(i2);
+}
+
static double tangent(const double* quadratic, double t) {
double a, b, c;
set_abc(quadratic, a, b, c);
@@ -136,5 +130,7 @@
result.y = tangent(&quadratic[0].y, t);
}
+
+
// unit test to return and validate parametric coefficients
#include "QuadraticParameterization_TestUtility.cpp"
diff --git a/experimental/Intersection/QuadraticParameterization.h b/experimental/Intersection/QuadraticParameterization.h
new file mode 100644
index 0000000..5c651b6
--- /dev/null
+++ b/experimental/Intersection/QuadraticParameterization.h
@@ -0,0 +1,27 @@
+#include "DataTypes.h"
+
+class QuadImplicitForm {
+public:
+ QuadImplicitForm(const Quadratic& q);
+ bool implicit_match(const QuadImplicitForm& two) const;
+
+ double x2() const { return p[xx_coeff]; }
+ double xy() const { return p[xy_coeff]; }
+ double y2() const { return p[yy_coeff]; }
+ double x() const { return p[x_coeff]; }
+ double y() const { return p[y_coeff]; }
+ double c() const { return p[c_coeff]; }
+
+private:
+ enum Coeffs {
+ xx_coeff,
+ xy_coeff,
+ yy_coeff,
+ x_coeff,
+ y_coeff,
+ c_coeff,
+ coeff_count
+ };
+
+ double p[coeff_count];
+};
diff --git a/experimental/Intersection/QuadraticParameterization_TestUtility.cpp b/experimental/Intersection/QuadraticParameterization_TestUtility.cpp
index 08a562b..7c91eb5 100644
--- a/experimental/Intersection/QuadraticParameterization_TestUtility.cpp
+++ b/experimental/Intersection/QuadraticParameterization_TestUtility.cpp
@@ -4,14 +4,13 @@
#include "Parameterization_Test.h"
bool point_on_parameterized_curve(const Quadratic& quad, const _Point& point) {
- double coeffs[coeff_count];
- implicit_coefficients(quad, coeffs);
- double xx = coeffs[ xx_coeff] * point.x * point.x;
- double xy = coeffs[ xy_coeff] * point.x * point.y;
- double yy = coeffs[ yy_coeff] * point.y * point.y;
- double x = coeffs[ x_coeff] * point.x;
- double y = coeffs[ y_coeff] * point.y;
- double c = coeffs[ c_coeff];
+ QuadImplicitForm q(quad);
+ double xx = q.x2() * point.x * point.x;
+ double xy = q.xy() * point.x * point.y;
+ double yy = q.y2() * point.y * point.y;
+ double x = q.x() * point.x;
+ double y = q.y() * point.y;
+ double c = q.c();
double sum = xx + xy + yy + x + y + c;
return approximately_zero(sum);
}
diff --git a/experimental/Intersection/QuadraticUtilities.cpp b/experimental/Intersection/QuadraticUtilities.cpp
index 8971fca..95be90a 100644
--- a/experimental/Intersection/QuadraticUtilities.cpp
+++ b/experimental/Intersection/QuadraticUtilities.cpp
@@ -28,8 +28,11 @@
int quadraticRoots(double A, double B, double C, double t[2]) {
B *= 2;
double square = B * B - 4 * A * C;
- if (square < 0) {
- return 0;
+ if (approximately_negative(square)) {
+ if (!approximately_positive(square)) {
+ return 0;
+ }
+ square = 0;
}
double squareRt = sqrt(square);
double Q = (B + (B < 0 ? -squareRt : squareRt)) / -2;
diff --git a/experimental/Intersection/QuarticRoot.cpp b/experimental/Intersection/QuarticRoot.cpp
new file mode 100644
index 0000000..509a4f0
--- /dev/null
+++ b/experimental/Intersection/QuarticRoot.cpp
@@ -0,0 +1,188 @@
+// from http://tog.acm.org/resources/GraphicsGems/gems/Roots3And4.c
+/*
+ * Roots3And4.c
+ *
+ * Utility functions to find cubic and quartic roots,
+ * coefficients are passed like this:
+ *
+ * c[0] + c[1]*x + c[2]*x^2 + c[3]*x^3 + c[4]*x^4 = 0
+ *
+ * The functions return the number of non-complex roots and
+ * put the values into the s array.
+ *
+ * Author: Jochen Schwarze (schwarze@isa.de)
+ *
+ * Jan 26, 1990 Version for Graphics Gems
+ * Oct 11, 1990 Fixed sign problem for negative q's in SolveQuartic
+ * (reported by Mark Podlipec),
+ * Old-style function definitions,
+ * IsZero() as a macro
+ * Nov 23, 1990 Some systems do not declare acos() and cbrt() in
+ * <math.h>, though the functions exist in the library.
+ * If large coefficients are used, EQN_EPS should be
+ * reduced considerably (e.g. to 1E-30), results will be
+ * correct but multiple roots might be reported more
+ * than once.
+ */
+
+#include <math.h>
+#include "CubicUtilities.h"
+#include "QuarticRoot.h"
+
+const double PI = 4 * atan(1);
+
+// unlike quadraticRoots in QuadraticUtilities.cpp, this does not discard
+// real roots <= 0 or >= 1
+static int quadraticRootsX(const double A, const double B, const double C,
+ double s[2]) {
+ if (approximately_zero(A)) {
+ if (approximately_zero(B)) {
+ s[0] = 0;
+ return C == 0;
+ }
+ s[0] = -C / B;
+ return 1;
+ }
+ /* normal form: x^2 + px + q = 0 */
+ const double p = B / (2 * A);
+ const double q = C / A;
+ const double D = p * p - q;
+ if (approximately_zero(D)) {
+ s[0] = -p;
+ return 1;
+ } else if (D < 0) {
+ return 0;
+ } else {
+ assert(D > 0);
+ double sqrt_D = sqrt(D);
+ s[0] = sqrt_D - p;
+ s[1] = -sqrt_D - p;
+ return 2;
+ }
+}
+
+// unlike cubicRoots in CubicUtilities.cpp, this does not discard
+// real roots <= 0 or >= 1
+static int cubicRootsX(const double A, const double B, const double C,
+ const double D, double s[3]) {
+ int num;
+ /* normal form: x^3 + Ax^2 + Bx + C = 0 */
+ const double invA = 1 / A;
+ const double a = B * invA;
+ const double b = C * invA;
+ const double c = D * invA;
+ /* substitute x = y - a/3 to eliminate quadric term:
+ x^3 +px + q = 0 */
+ const double a2 = a * a;
+ const double Q = (-a2 + b * 3) / 9;
+ const double R = (2 * a2 * a - 9 * a * b + 27 * c) / 54;
+ /* use Cardano's formula */
+ const double Q3 = Q * Q * Q;
+ const double R2plusQ3 = R * R + Q3;
+ if (approximately_zero(R2plusQ3)) {
+ if (approximately_zero(R)) {/* one triple solution */
+ s[0] = 0;
+ num = 1;
+ } else { /* one single and one double solution */
+
+ double u = cube_root(-R);
+ s[0] = 2 * u;
+ s[1] = -u;
+ num = 2;
+ }
+ }
+ else if (R2plusQ3 < 0) { /* Casus irreducibilis: three real solutions */
+ const double theta = 1.0/3 * acos(-R / sqrt(-Q3));
+ const double _2RootQ = 2 * sqrt(-Q);
+ s[0] = _2RootQ * cos(theta);
+ s[1] = -_2RootQ * cos(theta + PI / 3);
+ s[2] = -_2RootQ * cos(theta - PI / 3);
+ num = 3;
+ } else { /* one real solution */
+ const double sqrt_D = sqrt(R2plusQ3);
+ const double u = cube_root(sqrt_D - R);
+ const double v = -cube_root(sqrt_D + R);
+ s[0] = u + v;
+ num = 1;
+ }
+ /* resubstitute */
+ const double sub = 1.0/3 * a;
+ for (int i = 0; i < num; ++i) {
+ s[i] -= sub;
+ }
+ return num;
+}
+
+int quarticRoots(const double A, const double B, const double C, const double D,
+ const double E, double s[4]) {
+ if (approximately_zero(A)) {
+ if (approximately_zero(B)) {
+ return quadraticRootsX(C, D, E, s);
+ }
+ return cubicRootsX(B, C, D, E, s);
+ }
+ double u, v;
+ int num;
+ /* normal form: x^4 + Ax^3 + Bx^2 + Cx + D = 0 */
+ const double invA = 1 / A;
+ const double a = B * invA;
+ const double b = C * invA;
+ const double c = D * invA;
+ const double d = E * invA;
+ /* substitute x = y - a/4 to eliminate cubic term:
+ x^4 + px^2 + qx + r = 0 */
+ const double a2 = a * a;
+ const double p = -3 * a2 / 8 + b;
+ const double q = a2 * a / 8 - a * b / 2 + c;
+ const double r = -3 * a2 * a2 / 256 + a2 * b / 16 - a * c / 4 + d;
+ if (approximately_zero(r)) {
+ /* no absolute term: y(y^3 + py + q) = 0 */
+ num = cubicRootsX(1, 0, p, q, s);
+ s[num++] = 0;
+ } else {
+ /* solve the resolvent cubic ... */
+ (void) cubicRootsX(1, -p / 2, -r, r * p / 2 - q * q / 8, s);
+ /* ... and take the one real solution ... */
+ const double z = s[0];
+ /* ... to build two quadric equations */
+ u = z * z - r;
+ v = 2 * z - p;
+ if (approximately_zero(u)) {
+ u = 0;
+ } else if (u > 0) {
+ u = sqrt(u);
+ } else {
+ return 0;
+ }
+ if (approximately_zero(v)) {
+ v = 0;
+ } else if (v > 0) {
+ v = sqrt(v);
+ } else {
+ return 0;
+ }
+ num = quadraticRootsX(1, q < 0 ? -v : v, z - u, s);
+ num += quadraticRootsX(1, q < 0 ? v : -v, z + u, s + num);
+ }
+ // eliminate duplicates
+ int i;
+ for (i = 0; i < num - 1; ++i) {
+ for (int j = i + 1; j < num; ) {
+ if (approximately_equal(s[i], s[j])) {
+ if (j < --num) {
+ s[j] = s[num];
+ }
+ } else {
+ ++j;
+ }
+ }
+ }
+ /* resubstitute */
+ const double sub = a / 4;
+ for (i = 0; i < num; ++i) {
+ s[i] -= sub;
+ }
+ return num;
+}
+
+
diff --git a/experimental/Intersection/QuarticRoot.h b/experimental/Intersection/QuarticRoot.h
new file mode 100644
index 0000000..8baa842
--- /dev/null
+++ b/experimental/Intersection/QuarticRoot.h
@@ -0,0 +1,2 @@
+int quarticRoots(const double A, const double B, const double C, const double D,
+ const double E, double s[4]);
diff --git a/experimental/Intersection/QuarticRoot_Test.cpp b/experimental/Intersection/QuarticRoot_Test.cpp
new file mode 100644
index 0000000..2fe3e72
--- /dev/null
+++ b/experimental/Intersection/QuarticRoot_Test.cpp
@@ -0,0 +1,148 @@
+#include <assert.h>
+#include <math.h>
+#include "CubicUtilities.h"
+#include "Intersection_Tests.h"
+
+namespace QuarticRootTest {
+
+#include "QuarticRoot.cpp"
+
+}
+
+double mulA[] = {-3, -1, 1, 3};
+size_t mulACount = sizeof(mulA) / sizeof(mulA[0]);
+double rootB[] = {-9, -6, -3, -1, 0, 1, 3, 6, 9};
+size_t rootBCount = sizeof(rootB) / sizeof(rootB[0]);
+double rootC[] = {-8, -6, -2, -1, 0, 1, 2, 6, 8};
+size_t rootCCount = sizeof(rootC) / sizeof(rootC[0]);
+double rootD[] = {-7, -4, -1, 0, 1, 2, 5};
+size_t rootDCount = sizeof(rootD) / sizeof(rootD[0]);
+double rootE[] = {-5, -1, 0, 1, 7};
+size_t rootECount = sizeof(rootE) / sizeof(rootE[0]);
+
+static void quadraticTest() {
+ // (x - a)(x - b) == x^2 - (a + b)x + ab
+ for (size_t aIndex = 0; aIndex < mulACount; ++aIndex) {
+ for (size_t bIndex = 0; bIndex < rootBCount; ++bIndex) {
+ for (size_t cIndex = 0; cIndex < rootCCount; ++cIndex) {
+ const double A = mulA[aIndex];
+ const double B = rootB[bIndex];
+ const double C = rootC[cIndex];
+ const double b = A * (B + C);
+ const double c = A * B * C;
+ double roots[2];
+ const int rootCount = QuarticRootTest::quadraticRootsX(A, b, c, roots);
+ const int expected = 1 + (B != C);
+ assert(rootCount == expected);
+ assert(approximately_equal(roots[0], -B)
+ || approximately_equal(roots[0], -C));
+ if (B != C) {
+ assert(!approximately_equal(roots[0], roots[1]));
+ assert(approximately_equal(roots[1], -B)
+ || approximately_equal(roots[1], -C));
+ }
+ }
+ }
+ }
+}
+
+static void cubicTest() {
+ // (x - a)(x - b)(x - c) == x^3 - (a + b + c)x^2 + (ab + bc + ac)x - abc
+ for (size_t aIndex = 0; aIndex < mulACount; ++aIndex) {
+ for (size_t bIndex = 0; bIndex < rootBCount; ++bIndex) {
+ for (size_t cIndex = 0; cIndex < rootCCount; ++cIndex) {
+ for (size_t dIndex = 0; dIndex < rootDCount; ++dIndex) {
+ const double A = mulA[aIndex];
+ const double B = rootB[bIndex];
+ const double C = rootC[cIndex];
+ const double D = rootD[dIndex];
+ const double b = A * (B + C + D);
+ const double c = A * (B * C + C * D + B * D);
+ const double d = A * B * C * D;
+ double roots[3];
+ const int rootCount = QuarticRootTest::cubicRootsX(A, b, c, d, roots);
+ const int expected = 1 + (B != C) + (B != D && C != D);
+ assert(rootCount == expected);
+ assert(approximately_equal(roots[0], -B)
+ || approximately_equal(roots[0], -C)
+ || approximately_equal(roots[0], -D));
+ if (expected > 1) {
+ assert(!approximately_equal(roots[0], roots[1]));
+ assert(approximately_equal(roots[1], -B)
+ || approximately_equal(roots[1], -C)
+ || approximately_equal(roots[1], -D));
+ if (expected > 2) {
+ assert(!approximately_equal(roots[0], roots[2])
+ && !approximately_equal(roots[1], roots[2]));
+ assert(approximately_equal(roots[2], -B)
+ || approximately_equal(roots[2], -C)
+ || approximately_equal(roots[2], -D));
+ }
+ }
+ }
+ }
+ }
+ }
+}
+
+static void quarticTest() {
+ // (x - a)(x - b)(x - c)(x - d) == x^4 - (a + b + c + d)x^3
+ // + (ab + bc + cd + ac + bd + cd)x^2 - (abc + bcd + abd + acd) * x + abcd
+ for (size_t aIndex = 0; aIndex < mulACount; ++aIndex) {
+ for (size_t bIndex = 0; bIndex < rootBCount; ++bIndex) {
+ for (size_t cIndex = 0; cIndex < rootCCount; ++cIndex) {
+ for (size_t dIndex = 0; dIndex < rootDCount; ++dIndex) {
+ for (size_t eIndex = 0; eIndex < rootECount; ++eIndex) {
+ const double A = mulA[aIndex];
+ const double B = rootB[bIndex];
+ const double C = rootC[cIndex];
+ const double D = rootD[dIndex];
+ const double E = rootE[eIndex];
+ const double b = A * (B + C + D + E);
+ const double c = A * (B * C + C * D + B * D + B * E + C * E + D * E);
+ const double d = A * (B * C * D + B * C * E + B * D * E + C * D * E);
+ const double e = A * B * C * D * E;
+ double roots[4];
+ const int rootCount = QuarticRootTest::quarticRoots(A, b, c, d, e, roots);
+ const int expected = 1 + (B != C) + (B != D && C != D) + (B != E && C != E && D != E);
+ assert(rootCount == expected);
+ assert(approximately_equal(roots[0], -B)
+ || approximately_equal(roots[0], -C)
+ || approximately_equal(roots[0], -D)
+ || approximately_equal(roots[0], -E));
+ if (expected > 1) {
+ assert(!approximately_equal(roots[0], roots[1]));
+ assert(approximately_equal(roots[1], -B)
+ || approximately_equal(roots[1], -C)
+ || approximately_equal(roots[1], -D)
+ || approximately_equal(roots[1], -E));
+ if (expected > 2) {
+ assert(!approximately_equal(roots[0], roots[2])
+ && !approximately_equal(roots[1], roots[2]));
+ assert(approximately_equal(roots[2], -B)
+ || approximately_equal(roots[2], -C)
+ || approximately_equal(roots[2], -D)
+ || approximately_equal(roots[2], -E));
+ if (expected > 3) {
+ assert(!approximately_equal(roots[0], roots[3])
+ && !approximately_equal(roots[1], roots[3])
+ && !approximately_equal(roots[2], roots[3]));
+ assert(approximately_equal(roots[3], -B)
+ || approximately_equal(roots[3], -C)
+ || approximately_equal(roots[3], -D)
+ || approximately_equal(roots[3], -E));
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+}
+
+void QuarticRoot_Test() {
+ quadraticTest();
+ cubicTest();
+ quarticTest();
+}
diff --git a/experimental/Intersection/ShapeOps.cpp b/experimental/Intersection/ShapeOps.cpp
new file mode 100644
index 0000000..74e1626
--- /dev/null
+++ b/experimental/Intersection/ShapeOps.cpp
@@ -0,0 +1,166 @@
+/*
+ * Copyright 2012 Google Inc.
+ *
+ * Use of this source code is governed by a BSD-style license that can be
+ * found in the LICENSE file.
+ */
+#include "Simplify.h"
+
+namespace Op {
+
+#include "Simplify.cpp"
+
+static bool windingIsActive(int winding, int spanWinding, int oppWinding,
+ const ShapeOp op) {
+ return winding * spanWinding <= 0 && abs(winding) <= abs(spanWinding)
+ && (!winding || !spanWinding || winding == -spanWinding);
+}
+
+static void bridgeOp(SkTDArray<Contour*>& contourList, const ShapeOp op,
+ const int aXorMask, const int bXorMask, SkPath& simple) {
+ bool firstContour = true;
+ do {
+ Segment* topStart = findTopContour(contourList);
+ if (!topStart) {
+ break;
+ }
+ // Start at the top. Above the top is outside, below is inside.
+ // follow edges to intersection by changing the index by direction.
+ int index, endIndex;
+ Segment* current = topStart->findTop(index, endIndex);
+ int contourWinding;
+ if (firstContour) {
+ contourWinding = 0;
+ firstContour = false;
+ } else {
+ int sumWinding = current->windSum(SkMin32(index, endIndex));
+ // FIXME: don't I have to adjust windSum to get contourWinding?
+ if (sumWinding == SK_MinS32) {
+ sumWinding = current->computeSum(index, endIndex);
+ }
+ if (sumWinding == SK_MinS32) {
+ contourWinding = innerContourCheck(contourList, current,
+ index, endIndex);
+ } else {
+ contourWinding = sumWinding;
+ int spanWinding = current->spanSign(index, endIndex);
+ bool inner = useInnerWinding(sumWinding - spanWinding, sumWinding);
+ if (inner) {
+ contourWinding -= spanWinding;
+ }
+#if DEBUG_WINDING
+ SkDebugf("%s sumWinding=%d spanWinding=%d sign=%d inner=%d result=%d\n", __FUNCTION__,
+ sumWinding, spanWinding, SkSign32(index - endIndex),
+ inner, contourWinding);
+#endif
+ }
+#if DEBUG_WINDING
+ // SkASSERT(current->debugVerifyWinding(index, endIndex, contourWinding));
+ SkDebugf("%s contourWinding=%d\n", __FUNCTION__, contourWinding);
+#endif
+ }
+ SkPoint lastPt;
+ int winding = contourWinding;
+ int spanWinding = current->spanSign(index, endIndex);
+ int oppWinding = current->oppSign(index, endIndex);
+ bool active = windingIsActive(winding, spanWinding, oppWinding, op);
+ SkTDArray<Span*> chaseArray;
+ do {
+ #if DEBUG_WINDING
+ SkDebugf("%s active=%s winding=%d spanWinding=%d\n",
+ __FUNCTION__, active ? "true" : "false",
+ winding, spanWinding);
+ #endif
+ const SkPoint* firstPt = NULL;
+ do {
+ SkASSERT(!current->done());
+ int nextStart = index;
+ int nextEnd = endIndex;
+ Segment* next = current->findNextOp(chaseArray, active,
+ nextStart, nextEnd, winding, spanWinding, op,
+ aXorMask, bXorMask);
+ if (!next) {
+ if (active && firstPt && current->verb() != SkPath::kLine_Verb && *firstPt != lastPt) {
+ lastPt = current->addCurveTo(index, endIndex, simple, true);
+ SkASSERT(*firstPt == lastPt);
+ }
+ break;
+ }
+ if (!firstPt) {
+ firstPt = ¤t->addMoveTo(index, simple, active);
+ }
+ lastPt = current->addCurveTo(index, endIndex, simple, active);
+ current = next;
+ index = nextStart;
+ endIndex = nextEnd;
+ } while (*firstPt != lastPt && (active || !current->done()));
+ if (firstPt && active) {
+ #if DEBUG_PATH_CONSTRUCTION
+ SkDebugf("%s close\n", __FUNCTION__);
+ #endif
+ simple.close();
+ }
+ current = findChase(chaseArray, index, endIndex, contourWinding);
+ #if DEBUG_ACTIVE_SPANS
+ debugShowActiveSpans(contourList);
+ #endif
+ if (!current) {
+ break;
+ }
+ int lesser = SkMin32(index, endIndex);
+ spanWinding = current->spanSign(index, endIndex);
+ winding = current->windSum(lesser);
+ bool inner = useInnerWinding(winding - spanWinding, winding);
+ #if DEBUG_WINDING
+ SkDebugf("%s id=%d t=%1.9g spanWinding=%d winding=%d sign=%d"
+ " inner=%d result=%d\n",
+ __FUNCTION__, current->debugID(), current->t(lesser),
+ spanWinding, winding, SkSign32(index - endIndex),
+ useInnerWinding(winding - spanWinding, winding),
+ inner ? winding - spanWinding : winding);
+ #endif
+ if (inner) {
+ winding -= spanWinding;
+ }
+ int oppWinding = current->oppSign(index, endIndex);
+ active = windingIsActive(winding, spanWinding, oppWinding, op);
+ } while (true);
+ } while (true);
+}
+
+} // end of Op namespace
+
+
+void operate(const SkPath& one, const SkPath& two, ShapeOp op, SkPath& result) {
+ result.reset();
+ result.setFillType(SkPath::kEvenOdd_FillType);
+ // turn path into list of segments
+ SkTArray<Op::Contour> contours;
+ // FIXME: add self-intersecting cubics' T values to segment
+ Op::EdgeBuilder builder(one, contours);
+ const int aXorMask = builder.xorMask();
+ builder.addOperand(two);
+ const int bXorMask = builder.xorMask();
+ builder.finish();
+ SkTDArray<Op::Contour*> contourList;
+ makeContourList(contours, contourList);
+ Op::Contour** currentPtr = contourList.begin();
+ if (!currentPtr) {
+ return;
+ }
+ Op::Contour** listEnd = contourList.end();
+ // find all intersections between segments
+ do {
+ Op::Contour** nextPtr = currentPtr;
+ Op::Contour* current = *currentPtr++;
+ Op::Contour* next;
+ do {
+ next = *nextPtr++;
+ } while (addIntersectTs(current, next) && nextPtr != listEnd);
+ } while (currentPtr != listEnd);
+ // eat through coincident edges
+ coincidenceCheck(contourList);
+ fixOtherTIndex(contourList);
+ // construct closed contours
+ bridgeOp(contourList, op, aXorMask, bXorMask, result);
+}
diff --git a/experimental/Intersection/ShapeOps.h b/experimental/Intersection/ShapeOps.h
index d39fd0e..8ea3691 100644
--- a/experimental/Intersection/ShapeOps.h
+++ b/experimental/Intersection/ShapeOps.h
@@ -6,8 +6,29 @@
*/
#include "SkPath.h"
+// region-inspired approach
void contourBounds(const SkPath& path, SkTDArray<SkRect>& boundsArray);
void simplify(const SkPath& path, bool asFill, SkPath& simple);
+
+// contour outer edge walking approach
+#ifndef DEFINE_SHAPE_OP
+// FIXME: namespace testing doesn't allow global enums like this
+#define DEFINE_SHAPE_OP
+enum ShapeOp {
+ kDifference_Op,
+ kIntersect_Op,
+ kUnion_Op,
+ kXor_Op
+};
+
+enum ShapeOpMask {
+ kWinding_Mask = -1,
+ kNo_Mask = 0,
+ kEvenOdd_Mask = 1
+};
+#endif
+
+void operate(const SkPath& one, const SkPath& two, ShapeOp op, SkPath& result);
void simplifyx(const SkPath& path, SkPath& simple);
// FIXME: remove this section once debugging is complete
diff --git a/experimental/Intersection/Simplify.cpp b/experimental/Intersection/Simplify.cpp
index 0c81b2a..78c36f6 100644
--- a/experimental/Intersection/Simplify.cpp
+++ b/experimental/Intersection/Simplify.cpp
@@ -121,7 +121,12 @@
Intersections& intersections) {
MAKE_CONST_QUAD(aQuad, a);
MAKE_CONST_QUAD(bQuad, b);
+#define TRY_QUARTIC_SOLUTION 1
+#if TRY_QUARTIC_SOLUTION
+ intersect2(aQuad, bQuad, intersections);
+#else
intersect(aQuad, bQuad, intersections);
+#endif
return intersections.fUsed ? intersections.fUsed : intersections.fCoincidentUsed;
}
@@ -455,6 +460,13 @@
CubicLeftMost
};
+static int QuadRayIntersect(const SkPoint a[3], const SkPoint b[2],
+ Intersections& intersections) {
+ MAKE_CONST_QUAD(aQuad, a);
+ MAKE_CONST_LINE(bLine, b);
+ return intersectRay(aQuad, bLine, intersections);
+}
+
class Segment;
// sorting angles
@@ -481,13 +493,17 @@
maybe I set up LineParameters lazily
*/
bool operator<(const Angle& rh) const {
- if ((fDy < 0) ^ (rh.fDy < 0)) { // OPTIMIZATION: better to use fDy * rh.fDy < 0 ?
- return fDy < 0;
+ double y = dy();
+ double ry = rh.dy();
+ if ((y < 0) ^ (ry < 0)) { // OPTIMIZATION: better to use y * ry < 0 ?
+ return y < 0;
}
- if (fDy == 0 && rh.fDy == 0 && fDx * rh.fDx < 0) {
- return fDx < rh.fDx;
+ double x = dx();
+ double rx = rh.dx();
+ if (y == 0 && ry == 0 && x * rx < 0) {
+ return x < rx;
}
- AngleValue cmp = fDx * rh.fDy - rh.fDx * fDy;
+ AngleValue cmp = x * ry - rx * y;
if (!approximately_zero(cmp)) {
return cmp < 0;
}
@@ -523,25 +539,79 @@
SkASSERT(fSide != rh.fSide);
return fSide < rh.fSide;
}
- if (fDy != rh.fDy) {
- return fabs(fDy) < fabs(rh.fDy);
+ #if 0 // the following code is a bust. In particular, tangent length doesn't provide a sort
+ if (y != ry) {
+ return (fabs(y) < fabs(ry)) ^ (fSide > 0);
}
- if (fDx != rh.fDx) {
- return fabs(fDx) < fabs(rh.fDx);
+ if (x != rx) {
+ return (fabs(x) < fabs(rx)) ^ (fSide > 0);
}
- if (fSide != rh.fSide) {
- return fSide < rh.fSide;
+ // sort by second tangent angle
+ cmp = fD2x * rh.fD2y - rh.fD2x * fD2y;
+ if (!approximately_zero(cmp)) {
+ return cmp < 0;
}
SkASSERT(0); // add code for cubic case
+ #else
+ SkASSERT(fVerb == SkPath::kQuad_Verb); // worry about cubics later
+ SkASSERT(rh.fVerb == SkPath::kQuad_Verb);
+ // FIXME: until I can think of something better, project a ray perpendicular from the
+ // end of the shorter tangent through both curves and use the resulting angle to sort
+ // FIXME: some of this setup can be moved to set() if it works, or cached if it's expensive
+ double len = fTangent1.normalSquared();
+ double rlen = rh.fTangent1.normalSquared();
+ SkPoint ray[2];
+ const Quadratic& q = len < rlen ? fQ : rh.fQ;
+ ray[0].fX = SkDoubleToScalar(q[1].x);
+ ray[0].fY = SkDoubleToScalar(q[1].y);
+ ray[1].fX = SkDoubleToScalar((q[0].x + q[2].x) / 2);
+ ray[1].fY = SkDoubleToScalar((q[0].y + q[2].y) / 2);
+ Intersections i, ri;
+ int roots = QuadRayIntersect(fPts, ray, i);
+ SkASSERT(roots > 0);
+ int rroots = QuadRayIntersect(rh.fPts, ray, ri);
+ SkASSERT(rroots > 0);
+ SkPoint loc;
+ SkScalar best = SK_ScalarInfinity;
+ SkScalar dx, dy, dist;
+ int index;
+ for (index = 0; index < roots; ++index) {
+ QuadXYAtT(fPts, i.fT[0][index], &loc);
+ dx = loc.fX - ray[0].fX;
+ dy = loc.fY - ray[0].fY;
+ dist = dx * dx + dy * dy;
+ if (best > dist) {
+ best = dist;
+ }
+ }
+ for (index = 0; index < rroots; ++index) {
+ QuadXYAtT(rh.fPts, ri.fT[0][index], &loc);
+ dx = loc.fX - ray[0].fX;
+ dy = loc.fY - ray[0].fY;
+ dist = dx * dx + dy * dy;
+ if (best > dist) {
+ return fSide < 0;
+ }
+ }
+ return fSide > 0;
+ #endif
#endif
}
double dx() const {
+#if HIGH_DEF_ANGLES
+ return fTangent1.dx();
+#else
return fDx;
+#endif
}
double dy() const {
+#if HIGH_DEF_ANGLES
+ return fTangent1.dy();
+#else
return fDy;
+#endif
}
int end() const {
@@ -549,7 +619,11 @@
}
bool isHorizontal() const {
+#if HIGH_DEF_ANGLES
+ return dy() == 0 && fVerb == SkPath::kLine_Verb;
+#else
return fDy == 0 && fDDy == 0 && fDDDy == 0;
+#endif
}
// high precision version
@@ -559,38 +633,31 @@
fSegment = segment;
fStart = start;
fEnd = end;
+ fPts = orig;
+ fVerb = verb;
switch (verb) {
case SkPath::kLine_Verb:
_Line l;
LineSubDivideHD(orig, startT, endT, l);
- fDDx = fDDy = fDDDx = fDDDy = 0;
- fSide = 0;
// OPTIMIZATION: for pure line compares, we never need fTangent1.c
fTangent1.lineEndPoints(l);
+ fSide = 0;
break;
case SkPath::kQuad_Verb:
- Quadratic q;
- QuadSubDivideHD(orig, startT, endT, q);
- fDDx = approximately_pin(q[2].x - 2 * q[1].x + q[0].x);
- fDDy = approximately_pin(q[2].y - 2 * q[1].y + q[0].y);
- fDDDx = fDDDy = 0;
- fTangent1.quadEndPoints(q, 0, 1);
- fSide = -fTangent1.pointDistance(q[2]);
+ QuadSubDivideHD(orig, startT, endT, fQ);
+ fTangent1.quadEndPoints(fQ, 0, 1);
+ fSide = -fTangent1.pointDistance(fQ[2]); // not normalized -- compare sign only
break;
case SkPath::kCubic_Verb:
Cubic c;
CubicSubDivideHD(orig, startT, endT, c);
- fDDx = approximately_pin(c[2].x - 2 * c[1].x + c[0].x);
- fDDy = approximately_pin(c[2].y - 2 * c[1].y + c[0].y);
- fDDDx = approximately_pin(c[3].x + 3 * (c[1].x - c[2].x) - c[0].x);
- fDDDy = approximately_pin(c[3].y + 3 * (c[1].y - c[2].y) - c[0].y);
fTangent1.cubicEndPoints(c, 0, 1);
- fSide = -fTangent1.pointDistance(c[2]);
+ fSide = -fTangent1.pointDistance(c[2]); // not normalized -- compare sign only
break;
+ default:
+ SkASSERT(0);
}
// OPTIMIZATION: don't need these, access fTangent1 directly
- fDx = fTangent1.dx();
- fDy = fTangent1.dy();
}
#else
@@ -682,6 +749,10 @@
#if DEBUG_ANGLE
void debugShow(const SkPoint& a) const {
+#if HIGH_DEF_ANGLES
+ SkDebugf(" d=(%1.9g,%1.9g) side=%1.9g\n",
+ dx(), dy(), fSide);
+#else
SkDebugf(" d=(%1.9g,%1.9g) dd=(%1.9g,%1.9g) ddd=(%1.9g,%1.9g)\n",
fDx, fDy, fDDx, fDDy, fDDDx, fDDDy);
AngleValue ax = (AngleValue) a.fX;
@@ -708,18 +779,25 @@
" {SkPath::kCubic_Verb, {{%1.9g, %1.9g}, {%1.9g, %1.9g}, {%1.9g, %1.9g}, {%1.9g, %1.9g}}},\n",
ax, ay, bx, by, cx, cy, dx, dy);
}
+#endif
}
#endif
private:
+#if HIGH_DEF_ANGLES
+ const SkPoint* fPts;
+ Quadratic fQ;
+ SkPath::Verb fVerb;
+ double fSide;
+ LineParameters fTangent1;
+#else
AngleValue fDx;
AngleValue fDy;
AngleValue fDDx;
AngleValue fDDy;
AngleValue fDDDx;
AngleValue fDDDy;
- double fSide;
- LineParameters fTangent1; // FIXME: replaces fDx, fDy
+#endif
const Segment* fSegment;
int fStart;
int fEnd;
@@ -763,7 +841,7 @@
bool isEmpty() {
return fLeft > fRight || fTop > fBottom
- || fLeft == fRight && fTop == fBottom
+ || (fLeft == fRight && fTop == fBottom)
|| isnan(fLeft) || isnan(fRight)
|| isnan(fTop) || isnan(fBottom);
}
@@ -807,9 +885,23 @@
int fOtherIndex; // can't be used during intersection
int fWindSum; // accumulated from contours surrounding this one
int fWindValue; // 0 == canceled; 1 == normal; >1 == coincident
+ int fWindValueOpp; // opposite value, if any (for binary ops with coincidence)
bool fDone; // if set, this span to next higher T has been processed
};
+static const bool opLookup[][2][2] = {
+ // ==0 !=0
+ // b a b a
+ {{true , false}, {false, true }}, // a - b
+ {{false, false}, {true , true }}, // a & b
+ {{true , true }, {false, false}}, // a | b
+ {{true , true }, {true , true }}, // a ^ b
+};
+
+static bool activeOp(bool angleIsOp, int otherNonZero, ShapeOp op) {
+ return opLookup[op][otherNonZero][angleIsOp];
+}
+
class Segment {
public:
Segment() {
@@ -1013,8 +1105,8 @@
} while (tStart < 1 && oStart < 1 && !approximately_negative(oEnd - oStart));
}
- void addCubic(const SkPoint pts[4]) {
- init(pts, SkPath::kCubic_Verb);
+ void addCubic(const SkPoint pts[4], bool operand) {
+ init(pts, SkPath::kCubic_Verb, operand);
fBounds.setCubicBounds(pts);
}
@@ -1046,13 +1138,15 @@
path.cubicTo(edge[1].fX, edge[1].fY, edge[2].fX, edge[2].fY,
edge[3].fX, edge[3].fY);
break;
+ default:
+ SkASSERT(0);
}
}
return edge[fVerb];
}
- void addLine(const SkPoint pts[2]) {
- init(pts, SkPath::kLine_Verb);
+ void addLine(const SkPoint pts[2], bool operand) {
+ init(pts, SkPath::kLine_Verb, operand);
fBounds.set(pts, 2);
}
@@ -1074,8 +1168,8 @@
span.fOtherIndex = otherIndex;
}
- void addQuad(const SkPoint pts[3]) {
- init(pts, SkPath::kQuad_Verb);
+ void addQuad(const SkPoint pts[3], bool operand) {
+ init(pts, SkPath::kQuad_Verb, operand);
fBounds.setQuadBounds(pts);
}
@@ -1122,6 +1216,7 @@
span->fPt.fX = SK_ScalarNaN;
span->fWindSum = SK_MinS32;
span->fWindValue = 1;
+ span->fWindValueOpp = 0;
if ((span->fDone = newT == 1)) {
++fDoneSpans;
}
@@ -1141,6 +1236,7 @@
double oStartT, double oEndT) {
SkASSERT(!approximately_negative(endT - startT));
SkASSERT(!approximately_negative(oEndT - oStartT));
+ bool binary = fOperand != other.fOperand;
int index = 0;
while (!approximately_negative(startT - fTs[index].fT)) {
++index;
@@ -1161,7 +1257,11 @@
Span* span = test;
do {
if (decrement) {
- decrementSpan(span);
+ if (binary) {
+ --(span->fWindValueOpp);
+ } else {
+ decrementSpan(span);
+ }
} else if (track && span->fT < 1 && oTestT < 1) {
TrackOutside(outsideTs, span->fT, oTestT);
}
@@ -1211,10 +1311,11 @@
// set spans from start to end to increment the greater by one and decrement
// the lesser
- void addTCoincident(const int xorMask, double startT, double endT, Segment& other,
- double oStartT, double oEndT) {
+ void addTCoincident(const bool isXor, double startT, double endT,
+ Segment& other, double oStartT, double oEndT) {
SkASSERT(!approximately_negative(endT - startT));
SkASSERT(!approximately_negative(oEndT - oStartT));
+ bool binary = fOperand != other.fOperand;
int index = 0;
while (!approximately_negative(startT - fTs[index].fT)) {
++index;
@@ -1232,9 +1333,8 @@
SkTDArray<double> oxOutsideTs;
do {
bool transfer = test->fWindValue && oTest->fWindValue;
- bool winding = xorMask < 0;
- bool decrementThis = (test->fWindValue < oTest->fWindValue) & winding;
- bool decrementOther = (test->fWindValue >= oTest->fWindValue) & winding;
+ bool decrementThis = (test->fWindValue < oTest->fWindValue) & !isXor;
+ bool decrementOther = (test->fWindValue >= oTest->fWindValue) & !isXor;
Span* end = test;
double startT = end->fT;
int startIndex = index;
@@ -1246,7 +1346,11 @@
#ifdef SK_DEBUG
SkASSERT(abs(end->fWindValue) < gDebugMaxWindValue);
#endif
- ++(end->fWindValue);
+ if (binary) {
+ ++(end->fWindValueOpp);
+ } else {
+ ++(end->fWindValue);
+ }
} else if (decrementSpan(end)) {
TrackOutside(outsideTs, end->fT, oStartT);
}
@@ -1270,7 +1374,11 @@
#ifdef SK_DEBUG
SkASSERT(abs(oEnd->fWindValue) < gDebugMaxWindValue);
#endif
- ++(oEnd->fWindValue);
+ if (binary) {
+ ++(oEnd->fWindValueOpp);
+ } else {
+ ++(oEnd->fWindValue);
+ }
} else if (other.decrementSpan(oEnd)) {
TrackOutside(oOutsideTs, oEnd->fT, startT);
}
@@ -1475,15 +1583,16 @@
// if the intersection is edge on, wait for another one
continue;
}
- SkASSERT(pts == 1); // FIXME: more code required to disambiguate
- SkPoint pt;
- double foundT = intersections.fT[0][0];
- double testT = startT + (endT - startT) * foundT;
- (*SegmentXYAtT[fVerb])(fPts, testT, &pt);
- if (bestY < pt.fY && pt.fY < basePt.fY) {
- bestY = pt.fY;
- bestT = foundT < 1 ? start : end;
- hitT = testT;
+ for (int index = 0; index < pts; ++index) {
+ SkPoint pt;
+ double foundT = intersections.fT[0][index];
+ double testT = startT + (endT - startT) * foundT;
+ (*SegmentXYAtT[fVerb])(fPts, testT, &pt);
+ if (bestY < pt.fY && pt.fY < basePt.fY) {
+ bestY = pt.fY;
+ bestT = foundT < 1 ? start : end;
+ hitT = testT;
+ }
}
} while (fTs[end].fT != 1);
return bestT;
@@ -1544,6 +1653,205 @@
const Span& mSpan = fTs[SkMin32(start, end)];
return mSpan.fDone;
}
+
+ Segment* findNextOp(SkTDArray<Span*>& chase, bool active,
+ int& nextStart, int& nextEnd, int& winding, int& spanWinding, ShapeOp op,
+ const int aXorMask, const int bXorMask) {
+ const int startIndex = nextStart;
+ const int endIndex = nextEnd;
+ int outerWinding = winding;
+ int innerWinding = winding + spanWinding;
+ #if DEBUG_WINDING
+ SkDebugf("%s winding=%d spanWinding=%d outerWinding=%d innerWinding=%d\n",
+ __FUNCTION__, winding, spanWinding, outerWinding, innerWinding);
+ #endif
+ if (useInnerWinding(outerWinding, innerWinding)) {
+ outerWinding = innerWinding;
+ }
+ SkASSERT(startIndex != endIndex);
+ int count = fTs.count();
+ SkASSERT(startIndex < endIndex ? startIndex < count - 1
+ : startIndex > 0);
+ int step = SkSign32(endIndex - startIndex);
+ int end = nextSpan(startIndex, step);
+ SkASSERT(end >= 0);
+ Span* endSpan = &fTs[end];
+ Segment* other;
+ if (isSimple(end)) {
+ // mark the smaller of startIndex, endIndex done, and all adjacent
+ // spans with the same T value (but not 'other' spans)
+ #if DEBUG_WINDING
+ SkDebugf("%s simple\n", __FUNCTION__);
+ #endif
+ markDone(SkMin32(startIndex, endIndex), outerWinding);
+ other = endSpan->fOther;
+ nextStart = endSpan->fOtherIndex;
+ double startT = other->fTs[nextStart].fT;
+ nextEnd = nextStart;
+ do {
+ nextEnd += step;
+ } while (approximately_zero(startT - other->fTs[nextEnd].fT));
+ SkASSERT(step < 0 ? nextEnd >= 0 : nextEnd < other->fTs.count());
+ return other;
+ }
+ // more than one viable candidate -- measure angles to find best
+ SkTDArray<Angle> angles;
+ SkASSERT(startIndex - endIndex != 0);
+ SkASSERT((startIndex - endIndex < 0) ^ (step < 0));
+ addTwoAngles(startIndex, end, angles);
+ buildAngles(end, angles);
+ SkTDArray<Angle*> sorted;
+ sortAngles(angles, sorted);
+ int angleCount = angles.count();
+ int firstIndex = findStartingEdge(sorted, startIndex, end);
+ SkASSERT(firstIndex >= 0);
+ #if DEBUG_SORT
+ debugShowSort(__FUNCTION__, sorted, firstIndex, winding);
+ #endif
+ SkASSERT(sorted[firstIndex]->segment() == this);
+ #if DEBUG_WINDING
+ SkDebugf("%s [%d] sign=%d\n", __FUNCTION__, firstIndex, sorted[firstIndex]->sign());
+ #endif
+ int aSumWinding = winding;
+ int bSumWinding = winding;
+ bool angleIsOp = sorted[firstIndex]->segment()->operand();
+ int angleSpan = spanSign(sorted[firstIndex]);
+ if (angleIsOp) {
+ bSumWinding -= angleSpan;
+ } else {
+ aSumWinding -= angleSpan;
+ }
+ int nextIndex = firstIndex + 1;
+ int lastIndex = firstIndex != 0 ? firstIndex : angleCount;
+ const Angle* foundAngle = NULL;
+ // FIXME: found done logic probably fails if there are more than 4
+ // sorted angles. It should bias towards the first and last undone
+ // edges -- but not sure that it won't choose a middle (incorrect)
+ // edge if one is undone
+ bool foundDone = false;
+ bool foundDone2 = false;
+ // iterate through the angle, and compute everyone's winding
+ bool altFlipped = false;
+ bool foundFlipped = false;
+ int foundMax = SK_MinS32;
+ int foundSum = SK_MinS32;
+ Segment* nextSegment;
+ int lastNonZeroSum = winding;
+ do {
+ if (nextIndex == angleCount) {
+ nextIndex = 0;
+ }
+ const Angle* nextAngle = sorted[nextIndex];
+ nextSegment = nextAngle->segment();
+ angleIsOp = nextSegment->operand();
+ int sumWinding = angleIsOp ? bSumWinding : aSumWinding;
+ int maxWinding = sumWinding;
+ if (sumWinding) {
+ lastNonZeroSum = sumWinding;
+ }
+ sumWinding -= nextSegment->spanSign(nextAngle);
+ int xorMask = nextSegment->operand() ? bXorMask : aXorMask;
+ bool otherNonZero;
+ if (angleIsOp) {
+ bSumWinding = sumWinding;
+ otherNonZero = aSumWinding & aXorMask;
+ } else {
+ aSumWinding = sumWinding;
+ otherNonZero = bSumWinding & bXorMask;
+ }
+ altFlipped ^= lastNonZeroSum * sumWinding < 0; // flip if different signs
+ #if DEBUG_WINDING
+ SkASSERT(abs(sumWinding) <= gDebugMaxWindSum);
+ SkDebugf("%s [%d] maxWinding=%d sumWinding=%d sign=%d altFlipped=%d\n", __FUNCTION__,
+ nextIndex, maxWinding, sumWinding, nextAngle->sign(), altFlipped);
+ #endif
+
+ if (!(sumWinding & xorMask) && activeOp(angleIsOp, otherNonZero, op)) {
+ if (!active) {
+ markDone(SkMin32(startIndex, endIndex), outerWinding);
+ // FIXME: seems like a bug that this isn't calling userInnerWinding
+ nextSegment->markWinding(SkMin32(nextAngle->start(),
+ nextAngle->end()), maxWinding);
+ #if DEBUG_WINDING
+ SkDebugf("%s [%d] inactive\n", __FUNCTION__, nextIndex);
+ #endif
+ return NULL;
+ }
+ if (!foundAngle || foundDone) {
+ foundAngle = nextAngle;
+ foundDone = nextSegment->done(*nextAngle);
+ foundFlipped = altFlipped;
+ foundMax = maxWinding;
+ }
+ continue;
+ }
+ if (!(maxWinding & xorMask) && (!foundAngle || foundDone2)
+ && activeOp(angleIsOp, otherNonZero, op)) {
+ #if DEBUG_WINDING
+ if (foundAngle && foundDone2) {
+ SkDebugf("%s [%d] !foundAngle && foundDone2\n", __FUNCTION__, nextIndex);
+ }
+ #endif
+ foundAngle = nextAngle;
+ foundDone2 = nextSegment->done(*nextAngle);
+ foundFlipped = altFlipped;
+ foundSum = sumWinding;
+ }
+ if (nextSegment->done()) {
+ continue;
+ }
+ // if the winding is non-zero, nextAngle does not connect to
+ // current chain. If we haven't done so already, mark the angle
+ // as done, record the winding value, and mark connected unambiguous
+ // segments as well.
+ if (nextSegment->windSum(nextAngle) == SK_MinS32) {
+ if (useInnerWinding(maxWinding, sumWinding)) {
+ maxWinding = sumWinding;
+ }
+ Span* last;
+ if (foundAngle) {
+ last = nextSegment->markAndChaseWinding(nextAngle, maxWinding);
+ } else {
+ last = nextSegment->markAndChaseDone(nextAngle, maxWinding);
+ }
+ if (last) {
+ *chase.append() = last;
+ }
+ }
+ } while (++nextIndex != lastIndex);
+ markDone(SkMin32(startIndex, endIndex), outerWinding);
+ if (!foundAngle) {
+ return NULL;
+ }
+ nextStart = foundAngle->start();
+ nextEnd = foundAngle->end();
+ nextSegment = foundAngle->segment();
+ int flipped = foundFlipped ? -1 : 1;
+ spanWinding = SkSign32(spanWinding) * flipped * nextSegment->windValue(
+ SkMin32(nextStart, nextEnd));
+ if (winding) {
+ #if DEBUG_WINDING
+ SkDebugf("%s ---6 winding=%d foundSum=", __FUNCTION__, winding);
+ if (foundSum == SK_MinS32) {
+ SkDebugf("?");
+ } else {
+ SkDebugf("%d", foundSum);
+ }
+ SkDebugf(" foundMax=");
+ if (foundMax == SK_MinS32) {
+ SkDebugf("?");
+ } else {
+ SkDebugf("%d", foundMax);
+ }
+ SkDebugf("\n");
+ #endif
+ winding = foundSum;
+ }
+ #if DEBUG_WINDING
+ SkDebugf("%s spanWinding=%d flipped=%d\n", __FUNCTION__, spanWinding, flipped);
+ #endif
+ return nextSegment;
+ }
// so the span needs to contain the pairing info found here
// this should include the winding computed for the edge, and
@@ -1840,7 +2148,7 @@
}
// FIXME: this is tricky code; needs its own unit test
- void findTooCloseToCall(int xorMask) {
+ void findTooCloseToCall(bool isXor) {
int count = fTs.count();
if (count < 3) { // require t=0, x, 1 at minimum
return;
@@ -1963,7 +2271,10 @@
if (flipped) {
mOther->addTCancel(moStartT, moEndT, *tOther, toEndT, toStartT);
} else {
- mOther->addTCoincident(xorMask, moStartT, moEndT, *tOther, toStartT, toEndT);
+ // FIXME: this is bogus for multiple ops
+ // the xorMask needs to be accumulated from the union of the two
+ // edges -- which means that the segment must have its own copy of the mask
+ mOther->addTCoincident(isXor, moStartT, moEndT, *tOther, toStartT, toEndT);
}
}
}
@@ -2086,10 +2397,11 @@
return last;
}
- void init(const SkPoint pts[], SkPath::Verb verb) {
+ void init(const SkPoint pts[], SkPath::Verb verb, bool operand) {
+ fDoneSpans = 0;
+ fOperand = operand;
fPts = pts;
fVerb = verb;
- fDoneSpans = 0;
}
bool intersected() const {
@@ -2101,6 +2413,10 @@
|| fTs[endIndex].fWindSum != SK_MinS32;
}
+ bool isHorizontal() const {
+ return fBounds.fTop == fBounds.fBottom;
+ }
+
bool isLinear(int start, int end) const {
if (fVerb == SkPath::kLine_Verb) {
return true;
@@ -2144,10 +2460,6 @@
return false;
}
- bool isHorizontal() const {
- return fBounds.fTop == fBounds.fBottom;
- }
-
bool isVertical() const {
return fBounds.fLeft == fBounds.fRight;
}
@@ -2291,13 +2603,26 @@
}
return -1;
}
+
+ bool operand() const {
+ return fOperand;
+ }
+
+ int oppSign(int startIndex, int endIndex) const {
+ int result = startIndex < endIndex ? -fTs[startIndex].fWindValueOpp :
+ fTs[endIndex].fWindValueOpp;
+#if DEBUG_WIND_BUMP
+ SkDebugf("%s spanSign=%d\n", __FUNCTION__, result);
+#endif
+ return result;
+ }
const SkPoint* pts() const {
return fPts;
}
void reset() {
- init(NULL, (SkPath::Verb) -1);
+ init(NULL, (SkPath::Verb) -1, false);
fBounds.set(SK_ScalarMax, SK_ScalarMax, SK_ScalarMax, SK_ScalarMax);
fTs.reset();
}
@@ -2307,6 +2632,11 @@
return fTs[tIndex];
}
+ int spanSign(const Angle* angle) const {
+ SkASSERT(angle->segment() == this);
+ return spanSign(angle->start(), angle->end());
+ }
+
int spanSign(int startIndex, int endIndex) const {
int result = startIndex < endIndex ? -fTs[startIndex].fWindValue :
fTs[endIndex].fWindValue;
@@ -2316,11 +2646,6 @@
return result;
}
- int spanSign(const Angle* angle) const {
- SkASSERT(angle->segment() == this);
- return spanSign(angle->start(), angle->end());
- }
-
// OPTIMIZATION: mark as debugging only if used solely by tests
double t(int tIndex) const {
return fTs[tIndex].fT;
@@ -2385,7 +2710,7 @@
SkScalar xAtT(const Span* span) const {
return xyAtT(span).fX;
}
-
+
const SkPoint& xyAtT(int index) const {
return xyAtT(&fTs[index]);
}
@@ -2593,6 +2918,7 @@
Bounds fBounds;
SkTDArray<Span> fTs; // two or more (always includes t=0 t=1)
int fDoneSpans; // quick check that segment is finished
+ bool fOperand;
#if DEBUG_DUMP
int fID;
#endif
@@ -2604,6 +2930,7 @@
Contour* fContours[2];
int fSegments[2];
double fTs[2][2];
+ bool fXor;
};
class Contour {
@@ -2642,6 +2969,7 @@
coincidence.fTs[!swap][0] = ts.fCoincidentT[1][0];
coincidence.fTs[!swap][1] = ts.fCoincidentT[1][1];
}
+ coincidence.fXor = fOperand == other->fOperand ? fXor : true;
}
void addCross(const Contour* crosser) {
@@ -2654,12 +2982,12 @@
}
void addCubic(const SkPoint pts[4]) {
- fSegments.push_back().addCubic(pts);
+ fSegments.push_back().addCubic(pts, fOperand);
fContainsCurves = true;
}
int addLine(const SkPoint pts[2]) {
- fSegments.push_back().addLine(pts);
+ fSegments.push_back().addLine(pts, fOperand);
return fSegments.count();
}
@@ -2668,7 +2996,7 @@
}
int addQuad(const SkPoint pts[3]) {
- fSegments.push_back().addQuad(pts);
+ fSegments.push_back().addQuad(pts, fOperand);
fContainsCurves = true;
return fSegments.count();
}
@@ -2741,10 +3069,10 @@
return false;
}
- void findTooCloseToCall(int xorMask) {
+ void findTooCloseToCall() {
int segmentCount = fSegments.count();
for (int sIndex = 0; sIndex < segmentCount; ++sIndex) {
- fSegments[sIndex].findTooCloseToCall(xorMask);
+ fSegments[sIndex].findTooCloseToCall(fXor);
}
}
@@ -2761,7 +3089,9 @@
fContainsCurves = fContainsIntercepts = false;
}
- void resolveCoincidence(int xorMask) {
+ // FIXME: for binary ops, need to keep both ops winding contributions separately
+ // in edge array
+ void resolveCoincidence() {
int count = fCoincidences.count();
for (int index = 0; index < count; ++index) {
Coincidence& coincidence = fCoincidences[index];
@@ -2811,7 +3141,7 @@
|| thisOne.isMissing(endT) || other.isMissing(oEndT)) {
other.addTPair(oEndT, thisOne, endT, true);
}
- thisOne.addTCoincident(xorMask, startT, endT, other, oStartT, oEndT);
+ thisOne.addTCoincident(coincidence.fXor, startT, endT, other, oStartT, oEndT);
}
#if DEBUG_CONCIDENT
thisOne.debugShowTs();
@@ -2824,6 +3154,14 @@
return fSegments;
}
+ void setOperand(bool isOp) {
+ fOperand = isOp;
+ }
+
+ void setXor(bool isXor) {
+ fXor = isXor;
+ }
+
// OPTIMIZATION: feel pretty uneasy about this. It seems like once again
// we need to sort and walk edges in y, but that on the surface opens the
// same can of worms as before. But then, this is a rough sort based on
@@ -2941,6 +3279,8 @@
Bounds fBounds;
bool fContainsIntercepts;
bool fContainsCurves;
+ bool fOperand; // true for the second argument to a binary operator
+ bool fXor;
#if DEBUG_DUMP
int fID;
#endif
@@ -2950,15 +3290,52 @@
public:
EdgeBuilder(const SkPath& path, SkTArray<Contour>& contours)
- : fPath(path)
- , fCurrentContour(NULL)
+ : fPath(&path)
, fContours(contours)
+ , fCurrentContour(NULL)
+ , fOperand(false)
{
+ fXorMask = (path.getFillType() & 1) ? kEvenOdd_Mask : kWinding_Mask;
#if DEBUG_DUMP
gContourID = 0;
gSegmentID = 0;
#endif
+ fSecondHalf = preFetch();
+}
+
+void addOperand(const SkPath& path) {
+ fPath = &path;
+ fXorMask = (path.getFillType() & 1) ? kEvenOdd_Mask : kWinding_Mask;
+ preFetch();
+}
+
+void finish() {
walk();
+ complete();
+ if (fCurrentContour && !fCurrentContour->segments().count()) {
+ fContours.pop_back();
+ }
+ // correct pointers in contours since fReducePts may have moved as it grew
+ int cIndex = 0;
+ int extraCount = fExtra.count();
+ SkASSERT(extraCount == 0 || fExtra[0] == -1);
+ int eIndex = 0;
+ int rIndex = 0;
+ while (++eIndex < extraCount) {
+ int offset = fExtra[eIndex];
+ if (offset < 0) {
+ ++cIndex;
+ continue;
+ }
+ fCurrentContour = &fContours[cIndex];
+ rIndex += fCurrentContour->updateSegment(offset - 1,
+ &fReducePts[rIndex]);
+ }
+ fExtra.reset(); // we're done with this
+}
+
+ShapeOpMask xorMask() const {
+ return fXorMask;
}
protected:
@@ -2970,9 +3347,9 @@
}
}
-void walk() {
- // FIXME:remove once we can access path pts directly
- SkPath::RawIter iter(fPath); // FIXME: access path directly when allowed
+// FIXME:remove once we can access path pts directly
+int preFetch() {
+ SkPath::RawIter iter(*fPath); // FIXME: access path directly when allowed
SkPoint pts[4];
SkPath::Verb verb;
do {
@@ -2984,19 +3361,25 @@
fPathPts.append(verb, &pts[1]);
}
} while (verb != SkPath::kDone_Verb);
- // FIXME: end of section to remove once path pts are accessed directly
+ return fPathVerbs.count();
+}
+void walk() {
SkPath::Verb reducedVerb;
uint8_t* verbPtr = fPathVerbs.begin();
+ uint8_t* endOfFirstHalf = &verbPtr[fSecondHalf];
const SkPoint* pointsPtr = fPathPts.begin();
const SkPoint* finalCurveStart = NULL;
const SkPoint* finalCurveEnd = NULL;
+ SkPath::Verb verb;
while ((verb = (SkPath::Verb) *verbPtr++) != SkPath::kDone_Verb) {
switch (verb) {
case SkPath::kMove_Verb:
complete();
if (!fCurrentContour) {
fCurrentContour = fContours.push_back_n(1);
+ fCurrentContour->setOperand(fOperand);
+ fCurrentContour->setXor(fXorMask == kEvenOdd_Mask);
*fExtra.append() = -1; // start new contour
}
finalCurveEnd = pointsPtr++;
@@ -3056,38 +3439,23 @@
finalCurveStart = &pointsPtr[verb - 1];
pointsPtr += verb;
SkASSERT(fCurrentContour);
+ if (verbPtr == endOfFirstHalf) {
+ fOperand = true;
+ }
}
- complete();
- if (fCurrentContour && !fCurrentContour->segments().count()) {
- fContours.pop_back();
- }
- // correct pointers in contours since fReducePts may have moved as it grew
- int cIndex = 0;
- int extraCount = fExtra.count();
- SkASSERT(extraCount == 0 || fExtra[0] == -1);
- int eIndex = 0;
- int rIndex = 0;
- while (++eIndex < extraCount) {
- int offset = fExtra[eIndex];
- if (offset < 0) {
- ++cIndex;
- continue;
- }
- fCurrentContour = &fContours[cIndex];
- rIndex += fCurrentContour->updateSegment(offset - 1,
- &fReducePts[rIndex]);
- }
- fExtra.reset(); // we're done with this
}
private:
- const SkPath& fPath;
+ const SkPath* fPath;
SkTDArray<SkPoint> fPathPts; // FIXME: point directly to path pts instead
SkTDArray<uint8_t> fPathVerbs; // FIXME: remove
Contour* fCurrentContour;
SkTArray<Contour>& fContours;
SkTDArray<SkPoint> fReducePts; // segments created on the fly
SkTDArray<int> fExtra; // -1 marks new contour, > 0 offsets into contour
+ ShapeOpMask fXorMask;
+ int fSecondHalf;
+ bool fOperand;
};
class Work {
@@ -3466,15 +3834,15 @@
// resolve any coincident pairs found while intersecting, and
// see if coincidence is formed by clipping non-concident segments
-static void coincidenceCheck(SkTDArray<Contour*>& contourList, int xorMask) {
+static void coincidenceCheck(SkTDArray<Contour*>& contourList) {
int contourCount = contourList.count();
for (int cIndex = 0; cIndex < contourCount; ++cIndex) {
Contour* contour = contourList[cIndex];
- contour->resolveCoincidence(xorMask);
+ contour->resolveCoincidence();
}
for (int cIndex = 0; cIndex < contourCount; ++cIndex) {
Contour* contour = contourList[cIndex];
- contour->findTooCloseToCall(xorMask);
+ contour->findTooCloseToCall();
}
}
@@ -3962,7 +4330,6 @@
void simplifyx(const SkPath& path, SkPath& simple) {
// returns 1 for evenodd, -1 for winding, regardless of inverse-ness
- int xorMask = (path.getFillType() & 1) ? 1 : -1;
simple.reset();
simple.setFillType(SkPath::kEvenOdd_FillType);
@@ -3970,6 +4337,7 @@
SkTArray<Contour> contours;
// FIXME: add self-intersecting cubics' T values to segment
EdgeBuilder builder(path, contours);
+ builder.finish();
SkTDArray<Contour*> contourList;
makeContourList(contours, contourList);
Contour** currentPtr = contourList.begin();
@@ -3987,10 +4355,10 @@
} while (addIntersectTs(current, next) && nextPtr != listEnd);
} while (currentPtr != listEnd);
// eat through coincident edges
- coincidenceCheck(contourList, xorMask);
+ coincidenceCheck(contourList);
fixOtherTIndex(contourList);
// construct closed contours
- if (xorMask < 0) {
+ if (builder.xorMask() == kWinding_Mask) {
bridgeWinding(contourList, simple);
} else {
bridgeXor(contourList, simple);
diff --git a/experimental/Intersection/SimplifyNew_Test.cpp b/experimental/Intersection/SimplifyNew_Test.cpp
index da38135..980baa8 100644
--- a/experimental/Intersection/SimplifyNew_Test.cpp
+++ b/experimental/Intersection/SimplifyNew_Test.cpp
@@ -2538,12 +2538,180 @@
testSimplifyx(path);
}
-static void (*firstTest)() = 0; // testQuadratic20;
+static void testQuadratic21() {
+ SkPath path;
+ path.moveTo(0, 0);
+ path.quadTo(1, 0, 0, 1);
+ path.lineTo(0, 1);
+ path.close();
+ path.moveTo(0, 0);
+ path.lineTo(0, 0);
+ path.quadTo(1, 0, 0, 2);
+ path.close();
+ testSimplifyx(path);
+}
+
+static void testIntersect1() {
+ SkPath one, two, result;
+ one.addRect(0, 0, 6, 6, (SkPath::Direction) 0);
+ two.addRect(3, 3, 9, 9, (SkPath::Direction) 0);
+ operate(one, two, kIntersect_Op, result);
+ SkASSERT(result.countPoints() == 4);
+ SkASSERT(result.countVerbs() == 6); // move, 4 lines, close
+ SkRect bounds = result.getBounds();
+ SkASSERT(bounds.fLeft == 3);
+ SkASSERT(bounds.fTop == 3);
+ SkASSERT(bounds.fRight == 6);
+ SkASSERT(bounds.fBottom == 6);
+}
+
+static void testUnion1() {
+ SkPath one, two, result;
+ one.addRect(0, 0, 6, 6, (SkPath::Direction) 0);
+ two.addRect(3, 3, 9, 9, (SkPath::Direction) 0);
+ operate(one, two, kIntersect_Op, result);
+ SkASSERT(result.countPoints() == 8);
+ SkASSERT(result.countVerbs() == 10); // move, 8 lines, close
+ SkRect bounds = result.getBounds();
+ SkASSERT(bounds.fLeft == 0);
+ SkASSERT(bounds.fTop == 0);
+ SkASSERT(bounds.fRight == 9);
+ SkASSERT(bounds.fBottom == 9);
+}
+
+static void testDiff1() {
+ SkPath one, two, result;
+ one.addRect(0, 0, 6, 6, (SkPath::Direction) 0);
+ two.addRect(3, 3, 9, 9, (SkPath::Direction) 0);
+ operate(one, two, kIntersect_Op, result);
+ SkASSERT(result.countPoints() == 6);
+ SkASSERT(result.countVerbs() == 8); // move, 8 lines, close
+ SkRect bounds = result.getBounds();
+ SkASSERT(bounds.fLeft == 0);
+ SkASSERT(bounds.fTop == 0);
+ SkASSERT(bounds.fRight == 6);
+ SkASSERT(bounds.fBottom == 6);
+}
+
+static void testXor1() {
+ SkPath one, two, result;
+ one.addRect(0, 0, 6, 6, (SkPath::Direction) 0);
+ two.addRect(3, 3, 9, 9, (SkPath::Direction) 0);
+ operate(one, two, kIntersect_Op, result);
+ SkASSERT(result.countPoints() == 10);
+ SkASSERT(result.countVerbs() == 12); // move, 8 lines, close
+ SkRect bounds = result.getBounds();
+ SkASSERT(bounds.fLeft == 0);
+ SkASSERT(bounds.fTop == 0);
+ SkASSERT(bounds.fRight == 12);
+ SkASSERT(bounds.fBottom == 12);
+}
+
+static void testQuadratic22() {
+ SkPath path;
+ path.moveTo(0, 0);
+ path.quadTo(1, 0, 0, 1);
+ path.lineTo(0, 1);
+ path.close();
+ path.moveTo(0, 0);
+ path.lineTo(0, 0);
+ path.quadTo(0, 1, 2, 1);
+ path.close();
+ testSimplifyx(path);
+}
+
+static void testQuadratic23() {
+ SkPath path;
+ path.moveTo(0, 0);
+ path.quadTo(1, 0, 0, 1);
+ path.lineTo(0, 1);
+ path.close();
+ path.moveTo(0, 0);
+ path.lineTo(0, 0);
+ path.quadTo(0, 2, 1, 2);
+ path.close();
+ testSimplifyx(path);
+}
+
+static void testQuadratic24() {
+ SkPath path;
+ path.moveTo(0, 0);
+ path.quadTo(1, 0, 0, 1);
+ path.lineTo(0, 1);
+ path.close();
+ path.moveTo(0, 0);
+ path.lineTo(1, 0);
+ path.quadTo(2, 0, 0, 1);
+ path.close();
+ testSimplifyx(path);
+}
+
+static void testQuadratic25() {
+ SkPath path;
+ path.moveTo(0, 0);
+ path.quadTo(1, 0, 1, 1);
+ path.lineTo(1, 1);
+ path.close();
+ path.moveTo(0, 0);
+ path.lineTo(0, 0);
+ path.quadTo(2, 1, 0, 2);
+ path.close();
+ testSimplifyx(path);
+}
+
+static void testQuadratic26() {
+ SkPath path;
+ path.moveTo(0, 0);
+ path.quadTo(1, 0, 1, 1);
+ path.lineTo(0, 2);
+ path.close();
+ path.moveTo(0, 0);
+ path.lineTo(0, 0);
+ path.quadTo(1, 0, 0, 1);
+ path.close();
+ testSimplifyx(path);
+}
+
+static void testQuadratic27() {
+ SkPath path;
+ path.moveTo(0, 0);
+ path.quadTo(1, 0, 1, 1);
+ path.lineTo(2, 1);
+ path.close();
+ path.moveTo(0, 0);
+ path.lineTo(0, 0);
+ path.quadTo(2, 1, 0, 2);
+ path.close();
+ testSimplifyx(path);
+}
+
+static void testQuadratic28() {
+ SkPath path;
+ path.moveTo(0, 0);
+ path.quadTo(1, 0, 0, 1);
+ path.lineTo(0, 1);
+ path.close();
+ path.moveTo(0, 0);
+ path.lineTo(0, 2);
+ path.quadTo(1, 2, 0, 3);
+ path.close();
+ testSimplifyx(path);
+}
+
+static void (*firstTest)() = testQuadratic28;
static struct {
void (*fun)();
const char* str;
} tests[] = {
+ TEST(testQuadratic28),
+ TEST(testQuadratic27),
+ TEST(testQuadratic26),
+ TEST(testQuadratic25),
+ TEST(testQuadratic24),
+ TEST(testQuadratic23),
+ TEST(testQuadratic22),
+ TEST(testQuadratic21),
TEST(testQuadratic20),
TEST(testQuadratic19),
TEST(testQuadratic18),
@@ -2795,19 +2963,17 @@
void (*fun)();
const char* str;
} subTests[] = {
- TEST(testQuadralateral6),
- TEST(testQuadralateral6a),
- TEST(testFauxQuadralateral6d),
- TEST(testFauxQuadralateral6c),
- TEST(testFauxQuadralateral6b),
- TEST(testFauxQuadralateral6a),
- TEST(testFauxQuadralateral6),
+ TEST(testXor1),
+ TEST(testDiff1),
+ TEST(testUnion1),
+ TEST(testIntersect1),
};
static const size_t subTestCount = sizeof(subTests) / sizeof(subTests[0]);
static bool skipAll = false;
static bool runSubTests = false;
+static bool runReverse = false;
void SimplifyNew_Test() {
if (skipAll) {
@@ -2833,13 +2999,18 @@
SkDebugf(" %s [%s]\n", __FUNCTION__, tests[index].str);
(*tests[index].fun)();
}
- index = testCount - 1;
+ index = runReverse ? testCount - 1 : 0;
+ size_t last = runReverse ? 0 : testCount - 1;
bool firstTestComplete = false;
do {
SkDebugf(" %s [%s]\n", __FUNCTION__, tests[index].str);
(*tests[index].fun)();
firstTestComplete = true;
- } while (index--);
+ if (index == last) {
+ break;
+ }
+ index += runReverse ? -1 : 1;
+ } while (true);
#ifdef SK_DEBUG
gDebugMaxWindSum = SK_MaxS32;
gDebugMaxWindValue = SK_MaxS32;
diff --git a/experimental/Intersection/bc.htm b/experimental/Intersection/bc.htm
index 684f1c0..8a659f7 100644
--- a/experimental/Intersection/bc.htm
+++ b/experimental/Intersection/bc.htm
@@ -39,11 +39,62 @@
}}
</div>
+<div id="quad21a">
+bezier_clip q1=(0,0 1,0 0,2) q2=(0.5,0.25 0.5,0.5 0,1) minT=0 maxT=1
+</div>
+<div id="quad21b">
+bezier_clip q1=(0.5,0.25 0.5,0.375 0.375,0.5625) q2=(0,0 1,0 0,2) minT=0.3 maxT=0.78125
+</div>
+<div id="quad21c">
+bezier_clip q1=(0.42,0.18 0.6125,0.46875 0.341796875,1.22070312) q2=(0.5,0.25 0.5,0.375 0.375,0.5625) minT=0 maxT=0.926710098
+</div>
+<div id="quad21d">
+bezier_clip q1=(0.5,0.25 0.5,0.307919381 0.473162762,0.379257381) q2=(0.42,0.18 0.6125,0.46875 0.341796875,1.22070312) minT=0.187231244 maxT=0.729263299
+</div>
+<div id="quad21e">
+bezier_clip q1=(0.475846194,0.304363878 0.53317904,0.507883959 0.454423387,0.847492538) q2=(0.5,0.25 0.5,0.307919381 0.473162762,0.379257381) minT=0 maxT=1
+</div>
+<div id="quad21f">
+bezier_clip q1=(0.493290691,0.311274036 0.486581381,0.343588381 0.473162762,0.379257381) q2=(0.475846194,0.304363878 0.53317904,0.507883959 0.454423387,0.847492538) minT=0.0828748517 maxT=0.150086861
+</div>
+
+<div id="quad21g">
+(gdb) p smaller
+$1 = {{
+ x = 0.48441440743366754,
+ y = 0.33903196011243797
+ }, {
+ x = 0.48750982503868118,
+ y = 0.35346899178071778
+ }, {
+ x = 0.48999046908865357,
+ y = 0.368520797004039
+ }}
+(gdb) p larger
+$2 = {{
+ x = 0.49329069058425024,
+ y = 0.31127403581536672
+ }, {
+ x = 0.48658138116850047,
+ y = 0.34358838107698753
+ }, {
+ x = 0.47316276233700094,
+ y = 0.37925738104648321
+ }}
+</div>
+
</div>
<script type="text/javascript">
var testDivs = [
+ quad21g,
+ quad21a,
+ quad21b,
+ quad21c,
+ quad21d,
+ quad21e,
+ quad21f,
clip1,
];
@@ -63,7 +114,9 @@
function parse(test, title) {
var curveStrs = test.split("{{");
- var pattern = /\$?-?\d+\.*\d*/g;
+ if (curveStrs.length == 1)
+ curveStrs = test.split("=(");
+ var pattern = /[a-z$=]?-?\d+\.*\d*/g;
var curves = [];
for (var c in curveStrs) {
var curveStr = curveStrs[c];
@@ -77,7 +130,7 @@
if (pts.length > 0)
curves.push(pts);
}
- if (curves.length == 2) {
+ if (curves.length >= 2) {
tests.push(curves);
testTitles.push(title);
}
@@ -104,8 +157,13 @@
}
}
var subscale = 1;
- while ((xmax - xmin) * subscale < 0.1 && (ymax - ymin) * subscale < 0.1) {
- subscale *= 10;
+ if (xmax != xmin && ymax != ymin) {
+ while ((xmax - xmin) * subscale < 0.1 && (ymax - ymin) * subscale < 0.1) {
+ subscale *= 10;
+ if (subscale > 100000) {
+ break;
+ }
+ }
}
columns = Math.ceil(xmax) - Math.floor(xmin) + 1;
rows = Math.ceil(ymax) - Math.floor(ymin) + 1;
@@ -196,7 +254,7 @@
ctx.lineTo(xoffset + curve[6] * unit, yoffset + curve[7] * unit);
ctx.stroke();
}
- // optionally draw fat lines for cruve
+ // optionally draw fat lines for curve
if (fat1)
drawFat(test[0], xoffset, yoffset, unit);
if (fat2)
diff --git a/experimental/Intersection/op.htm b/experimental/Intersection/op.htm
index ae23166..4643ade 100644
--- a/experimental/Intersection/op.htm
+++ b/experimental/Intersection/op.htm
@@ -1322,11 +1322,107 @@
path.close();
</div>
+<div id="testQuadratic21">
+ path.moveTo(0, 0);
+ path.quadTo(1, 0, 0, 1);
+ path.lineTo(0, 1);
+ path.close();
+ path.moveTo(0, 0);
+ path.lineTo(0, 0);
+ path.quadTo(1, 0, 0, 2);
+ path.close();
+</div>
+
+<div id="testQuadratic22">
+ path.moveTo(0, 0);
+ path.quadTo(1, 0, 0, 1);
+ path.lineTo(0, 1);
+ path.close();
+ path.moveTo(0, 0);
+ path.lineTo(0, 0);
+ path.quadTo(0, 1, 2, 1);
+ path.close();
+</div>
+
+<div id="testQuadratic23">
+ path.moveTo(0, 0);
+ path.quadTo(1, 0, 0, 1);
+ path.lineTo(0, 1);
+ path.close();
+ path.moveTo(0, 0);
+ path.lineTo(0, 0);
+ path.quadTo(0, 2, 1, 2);
+ path.close();
+</div>
+
+<div id="testQuadratic24">
+ path.moveTo(0, 0);
+ path.quadTo(1, 0, 0, 1);
+ path.lineTo(0, 1);
+ path.close();
+ path.moveTo(0, 0);
+ path.lineTo(1, 0);
+ path.quadTo(2, 0, 0, 1);
+ path.close();
+</div>
+
+<div id="testQuadratic25">
+ path.moveTo(0, 0);
+ path.quadTo(1, 0, 1, 1);
+ path.lineTo(1, 1);
+ path.close();
+ path.moveTo(0, 0);
+ path.lineTo(0, 0);
+ path.quadTo(2, 1, 0, 2);
+ path.close();
+</div>
+
+<div id="testQuadratic26">
+ path.moveTo(0, 0);
+ path.quadTo(1, 0, 1, 1);
+ path.lineTo(0, 2);
+ path.close();
+ path.moveTo(0, 0);
+ path.lineTo(0, 0);
+ path.quadTo(1, 0, 0, 1);
+ path.close();
+</div>
+
+<div id="testQuadratic27">
+ path.moveTo(0, 0);
+ path.quadTo(1, 0, 1, 1);
+ path.lineTo(2, 1);
+ path.close();
+ path.moveTo(0, 0);
+ path.lineTo(0, 0);
+ path.quadTo(2, 1, 0, 2);
+ path.close();
+</div>
+
+<div id="testQuadratic28">
+ path.moveTo(0, 0);
+ path.quadTo(1, 0, 0, 1);
+ path.lineTo(0, 1);
+ path.close();
+ path.moveTo(0, 0);
+ path.lineTo(0, 2);
+ path.quadTo(1, 2, 0, 3);
+ path.close();
+</div>
+
</div>
<script type="text/javascript">
var testDivs = [
+ testQuadratic28,
+ testQuadratic27,
+ testQuadratic26,
+ testQuadratic25,
+ testQuadratic24,
+ testQuadratic23,
+ testQuadratic22,
+ testQuadratic21,
testQuadratic20,
testQuadratic19,
testQuadratic18,