shape ops work in progress
first 100,000 random cubic/cubic intersections working
git-svn-id: http://skia.googlecode.com/svn/trunk@7380 2bbb7eff-a529-9590-31e7-b0007b416f81
diff --git a/experimental/Intersection/CubeRoot.cpp b/experimental/Intersection/CubeRoot.cpp
index 82d2732..5f785a0 100644
--- a/experimental/Intersection/CubeRoot.cpp
+++ b/experimental/Intersection/CubeRoot.cpp
@@ -374,7 +374,7 @@
#endif
double cube_root(double x) {
- if (approximately_zero(x)) {
+ if (approximately_zero_cubed(x)) {
return 0;
}
double result = halley_cbrt3d(fabs(x));
diff --git a/experimental/Intersection/CubicIntersection.cpp b/experimental/Intersection/CubicIntersection.cpp
index 4ae0d84..a5b4261 100644
--- a/experimental/Intersection/CubicIntersection.cpp
+++ b/experimental/Intersection/CubicIntersection.cpp
@@ -10,6 +10,7 @@
#include "Intersections.h"
#include "IntersectionUtilities.h"
#include "LineIntersection.h"
+#include "LineUtilities.h"
static const double tClipLimit = 0.8; // http://cagd.cs.byu.edu/~tom/papers/bezclip.pdf see Multiple intersections
@@ -163,27 +164,49 @@
#include "CubicUtilities.h"
-// FIXME: ? if needed, compute the error term from the tangent length
-start here;
-// need better delta computation -- assert fails
-static double computeDelta(const Cubic& cubic, double t, double scale) {
- double attempt = scale / precisionUnit * 2;
-#if SK_DEBUG
- double precision = calcPrecision(cubic, t, scale);
- _Point dxy;
+static void cubicTangent(const Cubic& cubic, double t, _Line& tangent, _Point& pt, _Point& dxy) {
+ xy_at_t(cubic, t, tangent[0].x, tangent[0].y);
+ pt = tangent[1] = tangent[0];
dxdy_at_t(cubic, t, dxy);
- _Point p1, p2;
- xy_at_t(cubic, std::max(t - attempt, 0.), p1.x, p1.y);
- xy_at_t(cubic, std::min(t + attempt, 1.), p2.x, p2.y);
- double dx = p1.x - p2.x;
- double dy = p1.y - p2.y;
- double distSq = dx * dx + dy * dy;
- double dist = sqrt(distSq);
- assert(dist > precision);
-#endif
- return attempt;
+ tangent[0] -= dxy;
+ tangent[1] += dxy;
}
+static double cubicDelta(const _Point& dxy, _Line& tangent, double scale) {
+ double tangentLen = dxy.length();
+ tangent[0] -= tangent[1];
+ double intersectLen = tangent[0].length();
+ double result = intersectLen / tangentLen + scale;
+ return result;
+}
+
+// FIXME: after testing, make this static
+void computeDelta(const Cubic& c1, double t1, double scale1, const Cubic& c2, double t2,
+ double scale2, double& delta1, double& delta2) {
+ _Line tangent1, tangent2, line1, line2;
+ _Point dxy1, dxy2;
+ cubicTangent(c1, t1, line1, tangent1[0], dxy1);
+ cubicTangent(c2, t2, line2, tangent2[0], dxy2);
+ double range1[2], range2[2];
+ int found = intersect(line1, line2, range1, range2);
+ if (found == 0) {
+ range1[0] = 0.5;
+ } else {
+ SkASSERT(found == 1);
+ }
+ xy_at_t(line1, range1[0], tangent1[1].x, tangent1[1].y);
+#if SK_DEBUG
+ if (found == 1) {
+ xy_at_t(line2, range2[0], tangent2[1].x, tangent2[1].y);
+ SkASSERT(tangent2[1].approximatelyEqual(tangent1[1]));
+ }
+#endif
+ tangent2[1] = tangent1[1];
+ delta1 = cubicDelta(dxy1, tangent1, scale1 / precisionUnit);
+ delta2 = cubicDelta(dxy2, tangent2, scale2 / precisionUnit);
+}
+
+int debugDepth;
// this flavor approximates the cubics with quads to find the intersecting ts
// OPTIMIZE: if this strategy proves successful, the quad approximations, or the ts used
// to create the approximations, could be stored in the cubic segment
@@ -252,12 +275,16 @@
if (p1.approximatelyEqual(p2)) {
i.insert(i.swapped() ? to2 : to1, i.swapped() ? to1 : to2);
} else {
- double dt1 = computeDelta(cubic1, to1, t1e - t1s);
- double dt2 = computeDelta(cubic2, to2, t2e - t2s);
+ double dt1, dt2;
+ computeDelta(cubic1, to1, (t1e - t1s), cubic2, to2, (t2e - t2s), dt1, dt2);
+ ++debugDepth;
+ assert(debugDepth < 10);
i.swap();
- intersect2(cubic2, std::max(to2 - dt2, 0.), std::min(to2 + dt2, 1.),
- cubic1, std::max(to1 - dt1, 0.), std::min(to1 + dt1, 1.), i);
+ intersect2(cubic2, SkTMax(to2 - dt2, 0.), SkTMin(to2 + dt2, 1.),
+ cubic1, SkTMax(to1 - dt1, 0.), SkTMin(to1 + dt1, 1.), i);
i.swap();
+ --debugDepth;
+
}
}
t2Start = t2;
@@ -309,6 +336,7 @@
tMin = std::min(tMin, local2.fT[0][index]);
tMax = std::max(tMax, local2.fT[0][index]);
}
+ debugDepth = 0;
return intersect2(cubic1, start ? 0 : 1, start ? 1.0 / precisionUnit : 1 - 1.0 / precisionUnit,
cubic2, tMin, tMax, i);
}
@@ -318,6 +346,7 @@
// line segments intersect the cubic, then use the intersections to construct a subdivision for
// quadratic curve fitting.
bool intersect2(const Cubic& c1, const Cubic& c2, Intersections& i) {
+ debugDepth = 0;
bool result = intersect2(c1, 0, 1, c2, 0, 1, i);
// FIXME: pass in cached bounds from caller
_Rect c1Bounds, c2Bounds;
diff --git a/experimental/Intersection/CubicIntersection_Test.cpp b/experimental/Intersection/CubicIntersection_Test.cpp
index 8c2263e..2e738fc 100644
--- a/experimental/Intersection/CubicIntersection_Test.cpp
+++ b/experimental/Intersection/CubicIntersection_Test.cpp
@@ -57,23 +57,29 @@
}
}
+#define ONE_OFF_DEBUG 1
+
static void oneOff(const Cubic& cubic1, const Cubic& cubic2) {
SkTDArray<Quadratic> quads1;
cubic_to_quadratics(cubic1, calcPrecision(cubic1), quads1);
+#if ONE_OFF_DEBUG
for (int index = 0; index < quads1.count(); ++index) {
const Quadratic& q = quads1[index];
- SkDebugf("{{%1.9g,%1.9g}, {%1.9g,%1.9g}, {%1.9g,%1.9g}},\n", q[0].x, q[0].y,
+ SkDebugf(" {{%1.9g,%1.9g}, {%1.9g,%1.9g}, {%1.9g,%1.9g}},\n", q[0].x, q[0].y,
q[1].x, q[1].y, q[2].x, q[2].y);
}
SkDebugf("\n");
+#endif
SkTDArray<Quadratic> quads2;
cubic_to_quadratics(cubic2, calcPrecision(cubic2), quads2);
+#if ONE_OFF_DEBUG
for (int index = 0; index < quads2.count(); ++index) {
const Quadratic& q = quads2[index];
- SkDebugf("{{%1.9g,%1.9g}, {%1.9g,%1.9g}, {%1.9g,%1.9g}},\n", q[0].x, q[0].y,
+ SkDebugf(" {{%1.9g,%1.9g}, {%1.9g,%1.9g}, {%1.9g,%1.9g}},\n", q[0].x, q[0].y,
q[1].x, q[1].y, q[2].x, q[2].y);
}
SkDebugf("\n");
+#endif
Intersections intersections2;
intersect2(cubic1, cubic2, intersections2);
for (int pt = 0; pt < intersections2.used(); ++pt) {
@@ -83,13 +89,30 @@
int pt2 = intersections2.fFlip ? intersections2.used() - pt - 1 : pt;
double tt2 = intersections2.fT[1][pt2];
xy_at_t(cubic2, tt2, xy2.x, xy2.y);
+#if ONE_OFF_DEBUG
SkDebugf("%s t1=%1.9g (%1.9g, %1.9g) (%1.9g, %1.9g) t2=%1.9g\n", __FUNCTION__,
tt1, xy1.x, xy1.y, xy2.x, xy2.y, tt2);
+#endif
assert(xy1.approximatelyEqual(xy2));
}
}
static const Cubic testSet[] = {
+{{65.454505973241524, 93.881892270353575}, {45.867360264932437, 92.723972719499827}, {2.1464054482739447, 74.636369140183717}, {33.774068594804994, 40.770872887582925}},
+{{72.963387832494163, 95.659300729473728}, {11.809496633619768, 82.209921247423594}, {13.456139067865974, 57.329313623406605}, {36.060621606214262, 70.867335643091849}},
+
+{{32.484981432782945, 75.082940782924624}, {42.467313093350882, 48.131159948246157}, {3.5963115764764657, 43.208665839959245}, {79.442476890721579, 89.709102357602262}},
+{{18.98573861410177, 93.308887208490106}, {40.405250173250792, 91.039661826118675}, {8.0467721950480584, 42.100282172719147}, {40.883324221187891, 26.030185504830527}},
+
+{{7.5374809128872498, 82.441702896003477}, {22.444346930107265, 22.138854312775123}, {66.76091829629658, 50.753805856571446}, {78.193478508942519, 97.7932997968948}},
+{{97.700573130371311, 53.53260215070685}, {87.72443481149358, 84.575876772671876}, {19.215031396232092, 47.032676472809484}, {11.989686410869325, 10.659507480757082}},
+
+{{26.192053931854691, 9.8504326817814416}, {10.174241480498686, 98.476562741434464}, {21.177712558385782, 33.814968789841501}, {75.329030899018534, 55.02231980442177}},
+{{56.222082700683771, 24.54395039218662}, {95.589995289030483, 81.050822735322086}, {28.180450866082897, 28.837706255185282}, {60.128952916771617, 87.311672180570511}},
+
+{{42.449716172390481, 52.379709366885805}, {27.896043159019225, 48.797373636065686}, {92.770268299044233, 89.899302036454571}, {12.102066544863426, 99.43241951960718}},
+{{45.77532924980639, 45.958701495993274}, {37.458701356062065, 68.393691335056758}, {37.569326692060258, 27.673713456687381}, {60.674866037757539, 62.47349659096146}},
+
{{67.426548091427676, 37.993772624988935}, {23.483695892376684, 90.476863174921306}, {35.597065061143162, 79.872482633158796}, {75.38634169631932, 18.244890038969412}},
{{61.336508189019057, 82.693132843213675}, {44.639380902349664, 54.074825790745592}, {16.815615499771951, 20.049704667203923}, {41.866884958868326, 56.735503699973002}},
@@ -104,10 +127,14 @@
void CubicIntersection_OneOffTest() {
for (size_t outer = 0; outer < testSetCount - 1; ++outer) {
+#if ONE_OFF_DEBUG
SkDebugf("%s quads1[%d]\n", __FUNCTION__, outer);
+#endif
const Cubic& cubic1 = testSet[outer];
for (size_t inner = outer + 1; inner < testSetCount; ++inner) {
+#if ONE_OFF_DEBUG
SkDebugf("%s quads2[%d]\n", __FUNCTION__, inner);
+#endif
const Cubic& cubic2 = testSet[inner];
oneOff(cubic1, cubic2);
}
@@ -309,9 +336,56 @@
int pt2 = intersections2.fFlip ? intersections2.used() - pt - 1 : pt;
double tt2 = intersections2.fT[1][pt2];
xy_at_t(cubic2, tt2, xy2.x, xy2.y);
+ #if 0
SkDebugf("%s t1=%1.9g (%1.9g, %1.9g) (%1.9g, %1.9g) t2=%1.9g\n", __FUNCTION__,
tt1, xy1.x, xy1.y, xy2.x, xy2.y, tt2);
+ #endif
assert(xy1.approximatelyEqual(xy2));
}
}
}
+
+static Cubic deltaTestSet[] = {
+ {{1, 4}, {1, 4.*2/3}, {1, 4.*1/3}, {1, 0}},
+ {{0, 3}, {1, 2}, {2, 1}, {3, 0}},
+ {{1, 4}, {1, 4.*2/3}, {1, 4.*1/3}, {1, 0}},
+ {{3.5, 1}, {2.5, 2}, {1.5, 3}, {0.5, 4}}
+};
+
+size_t deltaTestSetLen = sizeof(deltaTestSet) / sizeof(deltaTestSet[0]);
+
+static double deltaTestSetT[] = {
+ 3./8,
+ 5./12,
+ 6./8,
+ 9./12
+};
+
+size_t deltaTestSetTLen = sizeof(deltaTestSetT) / sizeof(deltaTestSetT[0]);
+
+static double expectedT[] = {
+ 0.5,
+ 1./3,
+ 1./8,
+ 5./6
+};
+
+size_t expectedTLen = sizeof(expectedT) / sizeof(expectedT[0]);
+
+// FIXME: this test no longer valid -- does not take minimum scale contribution into account
+void CubicIntersection_ComputeDeltaTest() {
+ SkASSERT(deltaTestSetLen == deltaTestSetTLen);
+ SkASSERT(expectedTLen == deltaTestSetTLen);
+ for (size_t index = 0; index < deltaTestSetLen; index += 2) {
+ const Cubic& c1 = deltaTestSet[index];
+ const Cubic& c2 = deltaTestSet[index + 1];
+ double t1 = deltaTestSetT[index];
+ double t2 = deltaTestSetT[index + 1];
+ double d1, d2;
+ computeDelta(c1, t1, 1, c2, t2, 1, d1, d2);
+ SkASSERT(approximately_equal(t1 + d1, expectedT[index])
+ || approximately_equal(t1 - d1, expectedT[index]));
+ SkASSERT(approximately_equal(t2 + d2, expectedT[index + 1])
+ || approximately_equal(t2 - d2, expectedT[index + 1]));
+ }
+}
diff --git a/experimental/Intersection/CubicUtilities.cpp b/experimental/Intersection/CubicUtilities.cpp
index 19f16c6..3e2f474 100644
--- a/experimental/Intersection/CubicUtilities.cpp
+++ b/experimental/Intersection/CubicUtilities.cpp
@@ -21,7 +21,7 @@
#if SK_DEBUG
double calcPrecision(const Cubic& cubic, double t, double scale) {
Cubic part;
- sub_divide(cubic, SkMax32(0, t - scale), SkMin32(1, t + scale), part);
+ sub_divide(cubic, SkTMax(0., t - scale), SkTMin(1., t + scale), part);
return calcPrecision(part);
}
#endif
@@ -41,14 +41,11 @@
const double PI = 4 * atan(1);
-static bool is_unit_interval(double x) {
- return x > 0 && x < 1;
-}
-
// from SkGeometry.cpp (and Numeric Solutions, 5.6)
-int cubicRoots(double A, double B, double C, double D, double t[3]) {
+int cubicRootsValidT(double A, double B, double C, double D, double t[3]) {
+#if 0
if (approximately_zero(A)) { // we're just a quadratic
- return quadraticRoots(B, C, D, t);
+ return quadraticRootsValidT(B, C, D, t);
}
double a, b, c;
{
@@ -98,6 +95,113 @@
*roots++ = r;
}
return (int)(roots - t);
+#else
+ double s[3];
+ int realRoots = cubicRootsReal(A, B, C, D, s);
+ int foundRoots = add_valid_ts(s, realRoots, t);
+ return foundRoots;
+#endif
+}
+
+int cubicRootsReal(double A, double B, double C, double D, double s[3]) {
+#if SK_DEBUG
+ // create a string mathematica understands
+ // GDB set print repe 15 # if repeated digits is a bother
+ // set print elements 400 # if line doesn't fit
+ char str[1024];
+ bzero(str, sizeof(str));
+ sprintf(str, "Solve[%1.19g x^3 + %1.19g x^2 + %1.19g x + %1.19g == 0, x]", A, B, C, D);
+#endif
+ if (approximately_zero(A)) { // we're just a quadratic
+ return quadraticRootsReal(B, C, D, s);
+ }
+ if (approximately_zero(D)) { // 0 is one root
+ int num = quadraticRootsReal(A, B, C, s);
+ for (int i = 0; i < num; ++i) {
+ if (approximately_zero(s[i])) {
+ return num;
+ }
+ }
+ s[num++] = 0;
+ return num;
+ }
+ if (approximately_zero(A + B + C + D)) { // 1 is one root
+ int num = quadraticRootsReal(A, A + B, -D, s);
+ for (int i = 0; i < num; ++i) {
+ if (AlmostEqualUlps(s[i], 1)) {
+ return num;
+ }
+ }
+ s[num++] = 1;
+ return num;
+ }
+ double a, b, c;
+ {
+ double invA = 1 / A;
+ a = B * invA;
+ b = C * invA;
+ c = D * invA;
+ }
+ double a2 = a * a;
+ double Q = (a2 - b * 3) / 9;
+ double R = (2 * a2 * a - 9 * a * b + 27 * c) / 54;
+ double R2 = R * R;
+ double Q3 = Q * Q * Q;
+ double R2MinusQ3 = R2 - Q3;
+ double adiv3 = a / 3;
+ double r;
+ double* roots = s;
+#if 0
+ if (approximately_zero_squared(R2MinusQ3) && AlmostEqualUlps(R2, Q3)) {
+ if (approximately_zero_squared(R)) {/* one triple solution */
+ *roots++ = -adiv3;
+ } else { /* one single and one double solution */
+
+ double u = cube_root(-R);
+ *roots++ = 2 * u - adiv3;
+ *roots++ = -u - adiv3;
+ }
+ }
+ else
+#endif
+ if (R2MinusQ3 < 0) // we have 3 real roots
+ {
+ double theta = acos(R / sqrt(Q3));
+ double neg2RootQ = -2 * sqrt(Q);
+
+ r = neg2RootQ * cos(theta / 3) - adiv3;
+ *roots++ = r;
+
+ r = neg2RootQ * cos((theta + 2 * PI) / 3) - adiv3;
+ if (!AlmostEqualUlps(s[0], r)) {
+ *roots++ = r;
+ }
+ r = neg2RootQ * cos((theta - 2 * PI) / 3) - adiv3;
+ if (!AlmostEqualUlps(s[0], r) && (roots - s == 1 || !AlmostEqualUlps(s[1], r))) {
+ *roots++ = r;
+ }
+ }
+ else // we have 1 real root
+ {
+ double sqrtR2MinusQ3 = sqrt(R2MinusQ3);
+ double A = fabs(R) + sqrtR2MinusQ3;
+ A = cube_root(A);
+ if (R > 0) {
+ A = -A;
+ }
+ if (A != 0) {
+ A += Q / A;
+ }
+ r = A - adiv3;
+ *roots++ = r;
+ if (AlmostEqualUlps(R2, Q3)) {
+ r = -A / 2 - adiv3;
+ if (!AlmostEqualUlps(s[0], r)) {
+ *roots++ = r;
+ }
+ }
+ }
+ return (int)(roots - s);
}
// from http://www.cs.sunysb.edu/~qin/courses/geometry/4.pdf
@@ -136,14 +240,14 @@
double By = src[2].y - 2 * src[1].y + src[0].y;
double Cx = src[3].x + 3 * (src[1].x - src[2].x) - src[0].x;
double Cy = src[3].y + 3 * (src[1].y - src[2].y) - src[0].y;
- return quadraticRoots(Bx * Cy - By * Cx, (Ax * Cy - Ay * Cx) / 2, Ax * By - Ay * Bx, tValues);
+ return quadraticRootsValidT(Bx * Cy - By * Cx, (Ax * Cy - Ay * Cx), Ax * By - Ay * Bx, tValues);
}
bool rotate(const Cubic& cubic, int zero, int index, Cubic& rotPath) {
double dy = cubic[index].y - cubic[zero].y;
double dx = cubic[index].x - cubic[zero].x;
- if (approximately_equal(dy, 0)) {
- if (approximately_equal(dx, 0)) {
+ if (approximately_zero(dy)) {
+ if (approximately_zero(dx)) {
return false;
}
memcpy(rotPath, cubic, sizeof(Cubic));
diff --git a/experimental/Intersection/CubicUtilities.h b/experimental/Intersection/CubicUtilities.h
index 5205574..186ed43 100644
--- a/experimental/Intersection/CubicUtilities.h
+++ b/experimental/Intersection/CubicUtilities.h
@@ -15,13 +15,16 @@
double calcPrecision(const Cubic& cubic, double t, double scale);
#endif
void chop_at(const Cubic& src, CubicPair& dst, double t);
+// FIXME: should be private static but public here for testing
+void computeDelta(const Cubic& c1, double t1, double scale1, const Cubic& c2, double t2,
+ double scale2, double& delta1, double& delta2);
double cube_root(double x);
int cubic_to_quadratics(const Cubic& cubic, double precision,
SkTDArray<Quadratic>& quadratics);
void cubic_to_quadratics(const Cubic& cubic, double precision, SkTDArray<double>& ts);
void coefficients(const double* cubic, double& A, double& B, double& C, double& D);
-int cubicRoots(double A, double B, double C, double D, double t[3]);
-int cubicRootsX(double A, double B, double C, double D, double s[3]);
+int cubicRootsValidT(double A, double B, double C, double D, double t[3]);
+int cubicRootsReal(double A, double B, double C, double D, double s[3]);
void demote_cubic_to_quad(const Cubic& cubic, Quadratic& quad);
double dx_at_t(const Cubic& , double t);
double dy_at_t(const Cubic& , double t);
diff --git a/experimental/Intersection/DataTypes.h b/experimental/Intersection/DataTypes.h
index afdb8e9..de763ed 100644
--- a/experimental/Intersection/DataTypes.h
+++ b/experimental/Intersection/DataTypes.h
@@ -24,6 +24,8 @@
// FIXME: delete
int UlpsDiff(float A, float B);
+// FLT_EPSILON == 1.19209290E-07 == 1 / (2 ^ 23)
+const double FLT_EPSILON_CUBED = FLT_EPSILON * FLT_EPSILON * FLT_EPSILON;
const double FLT_EPSILON_SQUARED = FLT_EPSILON * FLT_EPSILON;
const double FLT_EPSILON_SQRT = sqrt(FLT_EPSILON);
@@ -42,6 +44,10 @@
return fabs(x) < FLT_EPSILON;
}
+inline bool approximately_zero_cubed(double x) {
+ return fabs(x) < FLT_EPSILON_CUBED;
+}
+
inline bool approximately_zero_squared(double x) {
return fabs(x) < FLT_EPSILON_SQUARED;
}
@@ -50,8 +56,28 @@
return fabs(x) < FLT_EPSILON_SQRT;
}
+// Use this for comparing Ts in the range of 0 to 1. For general numbers (larger and smaller) use
+// AlmostEqualUlps instead.
inline bool approximately_equal(double x, double y) {
+#if 1
return approximately_zero(x - y);
+#else
+// see http://visualstudiomagazine.com/blogs/tool-tracker/2011/11/compare-floating-point-numbers.aspx
+// this allows very small (e.g. degenerate) values to compare unequally, but in this case,
+// AlmostEqualUlps should be used instead.
+ if (x == y) {
+ return true;
+ }
+ double absY = fabs(y);
+ if (x == 0) {
+ return absY < FLT_EPSILON;
+ }
+ double absX = fabs(x);
+ if (y == 0) {
+ return absX < FLT_EPSILON;
+ }
+ return fabs(x - y) < (absX > absY ? absX : absY) * FLT_EPSILON;
+#endif
}
inline bool approximately_equal_squared(double x, double y) {
@@ -160,7 +186,7 @@
}
friend bool operator!=(const _Point& a, const _Point& b) {
- return a.x!= b.x || a.y != b.y;
+ return a.x != b.x || a.y != b.y;
}
// note: this can not be implemented with
@@ -171,9 +197,22 @@
&& AlmostEqualUlps((float) y, (float) a.y);
}
- double dot(const _Point& a) {
+ double cross(const _Point& a) const {
+ return x * a.y - y * a.x;
+ }
+
+ double dot(const _Point& a) const {
return x * a.x + y * a.y;
}
+
+ double length() const {
+ return sqrt(lengthSquared());
+ }
+
+ double lengthSquared() const {
+ return x * x + y * y;
+ }
+
};
typedef _Point _Line[2];
@@ -243,4 +282,17 @@
_Point pts[5];
};
+// FIXME: move these into SkTypes.h
+template <typename T> inline T SkTMax(T a, T b) {
+ if (a < b)
+ a = b;
+ return a;
+}
+
+template <typename T> inline T SkTMin(T a, T b) {
+ if (a > b)
+ a = b;
+ return a;
+}
+
#endif // __DataTypes_h__
diff --git a/experimental/Intersection/Intersection_Tests.cpp b/experimental/Intersection/Intersection_Tests.cpp
index cfa340f..bf86c07 100644
--- a/experimental/Intersection/Intersection_Tests.cpp
+++ b/experimental/Intersection/Intersection_Tests.cpp
@@ -15,9 +15,10 @@
void Intersection_Tests() {
int testsRun = 0;
- CubicIntersection_RandTest();
QuadraticIntersection_Test();
CubicIntersection_OneOffTest();
+ QuarticRoot_Test();
+ CubicIntersection_RandTest();
SimplifyNew_Test();
CubicsToQuadratics_RandTest();
CubicToQuadratics_Test();
@@ -32,7 +33,6 @@
LineQuadraticIntersection_Test();
MiniSimplify_Test();
SimplifyAngle_Test();
- QuarticRoot_Test();
QuadraticBezierClip_Test();
SimplifyFindNext_Test();
SimplifyFindTop_Test();
diff --git a/experimental/Intersection/Intersection_Tests.h b/experimental/Intersection/Intersection_Tests.h
index 6fc0c75..4bf3068 100644
--- a/experimental/Intersection/Intersection_Tests.h
+++ b/experimental/Intersection/Intersection_Tests.h
@@ -11,6 +11,7 @@
void ConvexHull_X_Test();
void CubicBezierClip_Test();
void CubicCoincidence_Test();
+void CubicIntersection_ComputeDeltaTest();
void CubicIntersection_OneOffTest();
void CubicIntersection_Test();
void CubicIntersection_RandTest();
diff --git a/experimental/Intersection/LineCubicIntersection.cpp b/experimental/Intersection/LineCubicIntersection.cpp
index 333e30f..cc63099 100644
--- a/experimental/Intersection/LineCubicIntersection.cpp
+++ b/experimental/Intersection/LineCubicIntersection.cpp
@@ -96,7 +96,7 @@
}
double A, B, C, D;
coefficients(&r[0].x, A, B, C, D);
- return cubicRoots(A, B, C, D, roots);
+ return cubicRootsValidT(A, B, C, D, roots);
}
int intersect() {
@@ -117,7 +117,7 @@
double A, B, C, D;
coefficients(&cubic[0].y, A, B, C, D);
D -= axisIntercept;
- return cubicRoots(A, B, C, D, roots);
+ return cubicRootsValidT(A, B, C, D, roots);
}
int horizontalIntersect(double axisIntercept, double left, double right, bool flipped) {
@@ -143,7 +143,7 @@
double A, B, C, D;
coefficients(&cubic[0].x, A, B, C, D);
D -= axisIntercept;
- return cubicRoots(A, B, C, D, roots);
+ return cubicRootsValidT(A, B, C, D, roots);
}
int verticalIntersect(double axisIntercept, double top, double bottom, bool flipped) {
diff --git a/experimental/Intersection/LineQuadraticIntersection.cpp b/experimental/Intersection/LineQuadraticIntersection.cpp
index 00c3554..9ec3fb7 100644
--- a/experimental/Intersection/LineQuadraticIntersection.cpp
+++ b/experimental/Intersection/LineQuadraticIntersection.cpp
@@ -124,7 +124,7 @@
double C = r[0];
A += C - 2 * B; // A = a - 2*b + c
B -= C; // B = -(b - c)
- return quadraticRoots(A, B, C, roots);
+ return quadraticRootsValidT(A, 2 * B, C, roots);
}
int intersect() {
@@ -148,7 +148,7 @@
D += F - 2 * E; // D = d - 2*e + f
E -= F; // E = -(d - e)
F -= axisIntercept;
- return quadraticRoots(D, E, F, roots);
+ return quadraticRootsValidT(D, 2 * E, F, roots);
}
int horizontalIntersect(double axisIntercept, double left, double right, bool flipped) {
@@ -177,7 +177,7 @@
D += F - 2 * E; // D = d - 2*e + f
E -= F; // E = -(d - e)
F -= axisIntercept;
- return quadraticRoots(D, E, F, roots);
+ return quadraticRootsValidT(D, 2 * E, F, roots);
}
int verticalIntersect(double axisIntercept, double top, double bottom, bool flipped) {
diff --git a/experimental/Intersection/QuadraticImplicit.cpp b/experimental/Intersection/QuadraticImplicit.cpp
index 2d9d9b9..9cd24e9 100644
--- a/experimental/Intersection/QuadraticImplicit.cpp
+++ b/experimental/Intersection/QuadraticImplicit.cpp
@@ -13,8 +13,6 @@
#include "QuadraticUtilities.h"
#include "TSearch.h"
-#include <algorithm> // for std::min, max
-
/* 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
@@ -22,14 +20,8 @@
* 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
*/
-#if SK_DEBUG
-#define QUARTIC_DEBUG 1
-#else
-#define QUARTIC_DEBUG 0
-#endif
-
static int findRoots(const QuadImplicitForm& i, const Quadratic& q2, double roots[4],
- bool useCubic, bool& disregardCount) {
+ bool oneHint) {
double a, b, c;
set_abc(&q2[0].x, a, b, c);
double d, e, f;
@@ -56,43 +48,11 @@
+ i.x() * c
+ i.y() * f
+ i.c();
-#if QUARTIC_DEBUG
- // create a string mathematica understands
- char str[1024];
- bzero(str, sizeof(str));
- sprintf(str, "Solve[%1.19g x^4 + %1.19g x^3 + %1.19g x^2 + %1.19g x + %1.19g == 0, x]",
- t4, t3, t2, t1, t0);
-#endif
- if (approximately_zero(t4)) {
- disregardCount = true;
- if (approximately_zero(t3)) {
- return quadraticRootsX(t2, t1, t0, roots);
- }
- return cubicRootsX(t3, t2, t1, t0, roots);
+ int rootCount = reducedQuarticRoots(t4, t3, t2, t1, t0, oneHint, roots);
+ if (rootCount >= 0) {
+ return rootCount;
}
- if (approximately_zero(t0)) { // 0 is one root
- disregardCount = true;
- int num = cubicRootsX(t4, t3, t2, t1, roots);
- for (int i = 0; i < num; ++i) {
- if (approximately_zero(roots[i])) {
- return num;
- }
- }
- roots[num++] = 0;
- return num;
- }
- if (useCubic) {
- assert(approximately_zero(t4 + t3 + t2 + t1 + t0)); // 1 is one root
- int num = cubicRootsX(t4, t4 + t3, -(t1 + t0), -t0, roots); // note that -C==A+B+D+E
- for (int i = 0; i < num; ++i) {
- if (approximately_equal(roots[i], 1)) {
- return num;
- }
- }
- roots[num++] = 1;
- return num;
- }
- return quarticRoots(t4, t3, t2, t1, t0, roots);
+ return quarticRootsReal(t4, t3, t2, t1, t0, roots);
}
static void addValidRoots(const double roots[4], const int count, const int side, Intersections& i) {
@@ -196,7 +156,10 @@
line[1].y += dxdy.y;
Intersections rootTs;
int roots = intersect(q1, line, rootTs);
- assert(roots == 1);
+ if (roots == 2) {
+ return false;
+ }
+ SkASSERT(roots == 1);
_Point pt2;
xy_at_t(q1, rootTs.fT[0][0], pt2.x, pt2.y);
if (!pt2.approximatelyEqual(mid)) {
@@ -288,12 +251,12 @@
}
static double flatMeasure(const Quadratic& q) {
- _Point mid;
- xy_at_t(q, 0.5, mid.x, mid.y);
- double dx = q[2].x - q[0].x;
- double dy = q[2].y - q[0].y;
- double length = sqrt(dx * dx + dy * dy); // OPTIMIZE: get rid of sqrt
- return ((mid.x - q[0].x) * dy - (mid.y - q[0].y) * dx) / length;
+ _Point mid = q[1];
+ mid -= q[0];
+ _Point dxy = q[2];
+ dxy -= q[0];
+ double length = dxy.length(); // OPTIMIZE: get rid of sqrt
+ return fabs(mid.cross(dxy) / length);
}
// FIXME ? should this measure both and then use the quad that is the flattest as the line?
@@ -306,11 +269,18 @@
return isLinearInner(q1, 0, 1, q2, 0, 1, i);
}
+// FIXME: if flat measure is sufficiently large, then probably the quartic solution failed
static bool relaxedIsLinear(const Quadratic& q1, const Quadratic& q2, Intersections& i) {
double m1 = flatMeasure(q1);
double m2 = flatMeasure(q2);
+#if SK_DEBUG
+ double min = SkTMin(m1, m2);
+ if (min > 5) {
+ SkDebugf("%s maybe not flat enough.. %1.9g\n", __FUNCTION__, min);
+ }
+#endif
i.reset();
- if (fabs(m1) < fabs(m2)) {
+ if (m1 < m2) {
isLinearInner(q1, 0, 1, q2, 0, 1, i);
return false;
} else {
@@ -400,18 +370,10 @@
return i.fCoincidentUsed > 0;
}
double roots1[4], roots2[4];
- bool disregardCount1 = false;
- bool disregardCount2 = false;
bool useCubic = q1[0] == q2[0] || q1[0] == q2[2] || q1[2] == q2[0];
- int rootCount = findRoots(i2, q1, roots1, useCubic, disregardCount1);
+ int rootCount = findRoots(i2, q1, roots1, useCubic);
// OPTIMIZATION: could short circuit here if all roots are < 0 or > 1
- int rootCount2 = findRoots(i1, q2, roots2, useCubic, disregardCount2);
- #if 0
- if (rootCount != rootCount2 && !disregardCount1 && !disregardCount2) {
- unsortableExpanse(q1, q2, i);
- return false;
- }
- #endif
+ int rootCount2 = findRoots(i1, q2, roots2, useCubic);
addValidRoots(roots1, rootCount, 0, i);
addValidRoots(roots2, rootCount2, 1, i);
if (i.insertBalanced() && i.fUsed <= 1) {
diff --git a/experimental/Intersection/QuadraticIntersection_Test.cpp b/experimental/Intersection/QuadraticIntersection_Test.cpp
index 432f614..055a18e8 100644
--- a/experimental/Intersection/QuadraticIntersection_Test.cpp
+++ b/experimental/Intersection/QuadraticIntersection_Test.cpp
@@ -59,6 +59,21 @@
}
static const Quadratic testSet[] = {
+ {{41.5072916,87.1234036}, {28.2747836,80.9545395}, {23.5780771,69.3344126}},
+ {{72.9633878,95.6593007}, {42.7738746,88.4730382}, {31.1932785,80.2458029}},
+
+ {{31.1663962,54.7302484}, {31.1662882,54.7301074}, {31.1663969,54.7302485}},
+ {{26.0404936,45.4260361}, {27.7887523,33.1863051}, {40.8833242,26.0301855}},
+
+ {{29.9404074,49.1672596}, {44.3131071,45.3915253}, {58.1067559,59.5061814}},
+ {{72.6510251,64.2972928}, {53.6989659,60.1862397}, {35.2053722,44.8391126}},
+
+{{52.14807018377202, 65.012420045148644}, {44.778669050208237, 66.315562705604378}, {51.619118408823567, 63.787827046262684}},
+{{30.004993234763383, 93.921296668202288}, {53.384822003076991, 60.732180341802753}, {58.652998934338584, 43.111073088306185}},
+
+{{80.897794748143198, 49.236332042718459}, {81.082078218891212, 64.066749904488631}, {69.972305057149981, 72.968595519850993}},
+{{72.503745601281395, 32.952320736577882}, {88.030880716061645, 38.137194847810164}, {73.193774825517906, 67.773492479591397}},
+
{{67.426548091427676, 37.993772624988935}, {51.129513170665035, 57.542281234563646}, {44.594748190899189, 65.644267382683879}},
{{61.336508189019057, 82.693132843213675}, {54.825078921449354, 71.663932799212432}, {47.727444217558926, 61.4049645128392}},
@@ -119,37 +134,47 @@
const size_t testSetCount = sizeof(testSet) / sizeof(testSet[0]);
+#define ONE_OFF_DEBUG 0
+
+static void oneOffTest1(size_t outer, size_t inner) {
+ const Quadratic& quad1 = testSet[outer];
+ const Quadratic& quad2 = testSet[inner];
+ Intersections intersections2;
+ intersect2(quad1, quad2, intersections2);
+ if (intersections2.fUnsortable) {
+ SkASSERT(0);
+ return;
+ }
+ for (int pt = 0; pt < intersections2.used(); ++pt) {
+ double tt1 = intersections2.fT[0][pt];
+ double tx1, ty1;
+ xy_at_t(quad1, tt1, tx1, ty1);
+ int pt2 = intersections2.fFlip ? intersections2.used() - pt - 1 : pt;
+ double tt2 = intersections2.fT[1][pt2];
+ double tx2, ty2;
+ xy_at_t(quad2, tt2, tx2, ty2);
+ if (!AlmostEqualUlps(tx1, tx2)) {
+ SkDebugf("%s [%d,%d] x!= t1=%g (%g,%g) t2=%g (%g,%g)\n",
+ __FUNCTION__, (int)outer, (int)inner, tt1, tx1, ty1, tt2, tx2, ty2);
+ SkASSERT(0);
+ }
+ if (!AlmostEqualUlps(ty1, ty2)) {
+ SkDebugf("%s [%d,%d] y!= t1=%g (%g,%g) t2=%g (%g,%g)\n",
+ __FUNCTION__, (int)outer, (int)inner, tt1, tx1, ty1, tt2, tx2, ty2);
+ SkASSERT(0);
+ }
+#if ONE_OFF_DEBUG
+ SkDebugf("%s [%d][%d] t1=%1.9g (%1.9g, %1.9g) t2=%1.9g\n", __FUNCTION__,
+ outer, inner, tt1, tx1, tx2, tt2);
+#endif
+ }
+}
+
static void oneOffTest() {
+// oneOffTest1(0, 1);
for (size_t outer = 0; outer < testSetCount - 1; ++outer) {
for (size_t inner = outer + 1; inner < testSetCount; ++inner) {
- const Quadratic& quad1 = testSet[outer];
- const Quadratic& quad2 = testSet[inner];
- Intersections intersections2;
- intersect2(quad1, quad2, intersections2);
- if (intersections2.fUnsortable) {
- continue;
- }
- for (int pt = 0; pt < intersections2.used(); ++pt) {
- double tt1 = intersections2.fT[0][pt];
- double tx1, ty1;
- xy_at_t(quad1, tt1, tx1, ty1);
- int pt2 = intersections2.fFlip ? intersections2.used() - pt - 1 : pt;
- double tt2 = intersections2.fT[1][pt2];
- double tx2, ty2;
- xy_at_t(quad2, tt2, tx2, ty2);
- if (!AlmostEqualUlps(tx1, tx2)) {
- SkDebugf("%s [%d,%d] x!= t1=%g (%g,%g) t2=%g (%g,%g)\n",
- __FUNCTION__, (int)outer, (int)inner, tt1, tx1, ty1, tt2, tx2, ty2);
- SkASSERT(0);
- }
- if (!AlmostEqualUlps(ty1, ty2)) {
- SkDebugf("%s [%d,%d] y!= t1=%g (%g,%g) t2=%g (%g,%g)\n",
- __FUNCTION__, (int)outer, (int)inner, tt1, tx1, ty1, tt2, tx2, ty2);
- SkASSERT(0);
- }
- SkDebugf("%s [%d][%d] t1=%1.9g (%1.9g, %1.9g) t2=%1.9g\n", __FUNCTION__,
- outer, inner, tt1, tx1, tx2, tt2);
- }
+ oneOffTest1(outer, inner);
}
}
}
diff --git a/experimental/Intersection/QuadraticReduceOrder.cpp b/experimental/Intersection/QuadraticReduceOrder.cpp
index aa1d058..096f6b6 100644
--- a/experimental/Intersection/QuadraticReduceOrder.cpp
+++ b/experimental/Intersection/QuadraticReduceOrder.cpp
@@ -92,7 +92,7 @@
if (root) {
_Point extrema;
extrema.x = interp_quad_coords(quad[0].x, quad[1].x, quad[2].x, tValue);
- extrema.y = interp_quad_coords(quad[0].x, quad[1].x, quad[2].x, tValue);
+ extrema.y = interp_quad_coords(quad[0].y, quad[1].y, quad[2].y, tValue);
// sameSide > 0 means mid is smaller than either [0] or [2], so replace smaller
int replace;
if (useX) {
diff --git a/experimental/Intersection/QuadraticUtilities.cpp b/experimental/Intersection/QuadraticUtilities.cpp
index 8b02b4c..475ec13 100644
--- a/experimental/Intersection/QuadraticUtilities.cpp
+++ b/experimental/Intersection/QuadraticUtilities.cpp
@@ -5,6 +5,7 @@
* found in the LICENSE file.
*/
#include "QuadraticUtilities.h"
+#include "SkTypes.h"
#include <math.h>
/*
@@ -20,12 +21,37 @@
*/
+
+int add_valid_ts(double s[], int realRoots, double* t) {
+ int foundRoots = 0;
+ for (int index = 0; index < realRoots; ++index) {
+ double tValue = s[index];
+ if (approximately_zero_or_more(tValue) && approximately_one_or_less(tValue)) {
+ if (approximately_less_than_zero(tValue)) {
+ tValue = 0;
+ } else if (approximately_greater_than_one(tValue)) {
+ tValue = 1;
+ }
+ for (int idx2 = 0; idx2 < foundRoots; ++idx2) {
+ if (approximately_equal(t[idx2], tValue)) {
+ goto nextRoot;
+ }
+ }
+ t[foundRoots++] = tValue;
+ }
+nextRoot:
+ ;
+ }
+ return foundRoots;
+}
+
// note: caller expects multiple results to be sorted smaller first
// note: http://en.wikipedia.org/wiki/Loss_of_significance has an interesting
// analysis of the quadratic equation, suggesting why the following looks at
// the sign of B -- and further suggesting that the greatest loss of precision
// is in b squared less two a c
-int quadraticRoots(double A, double B, double C, double t[2]) {
+int quadraticRootsValidT(double A, double B, double C, double t[2]) {
+#if 0
B *= 2;
double square = B * B - 4 * A * C;
if (approximately_negative(square)) {
@@ -61,9 +87,64 @@
t[0] = ratio;
}
}
+#else
+ double s[2];
+ int realRoots = quadraticRootsReal(A, B, C, s);
+ int foundRoots = add_valid_ts(s, realRoots, t);
+#endif
return foundRoots;
}
+// unlike quadratic roots, this does not discard real roots <= 0 or >= 1
+int quadraticRootsReal(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 p2 = p * p;
+#if 0
+ double D = AlmostEqualUlps(p2, q) ? 0 : p2 - q;
+ if (D <= 0) {
+ if (D < 0) {
+ return 0;
+ }
+ s[0] = -p;
+ SkDebugf("[%d] %1.9g\n", 1, s[0]);
+ return 1;
+ }
+ double sqrt_D = sqrt(D);
+ s[0] = sqrt_D - p;
+ s[1] = -sqrt_D - p;
+ SkDebugf("[%d] %1.9g %1.9g\n", 2, s[0], s[1]);
+ return 2;
+#else
+ if (!AlmostEqualUlps(p2, q) && p2 < q) {
+ return 0;
+ }
+ double sqrt_D = 0;
+ if (p2 > q) {
+ sqrt_D = sqrt(p2 - q);
+ }
+ s[0] = sqrt_D - p;
+ s[1] = -sqrt_D - p;
+#if 0
+ if (AlmostEqualUlps(s[0], s[1])) {
+ SkDebugf("[%d] %1.9g\n", 1, s[0]);
+ } else {
+ SkDebugf("[%d] %1.9g %1.9g\n", 2, s[0], s[1]);
+ }
+#endif
+ return 1 + !AlmostEqualUlps(s[0], s[1]);
+#endif
+}
+
static double derivativeAtT(const double* quad, double t) {
double a = t - 1;
double b = 1 - 2 * t;
diff --git a/experimental/Intersection/QuadraticUtilities.h b/experimental/Intersection/QuadraticUtilities.h
index 0fe8be8..f423dc6 100644
--- a/experimental/Intersection/QuadraticUtilities.h
+++ b/experimental/Intersection/QuadraticUtilities.h
@@ -10,6 +10,7 @@
#include "DataTypes.h"
+int add_valid_ts(double s[], int realRoots, double* t);
void chop_at(const Quadratic& src, QuadraticPair& dst, double t);
double dx_at_t(const Quadratic& , double t);
double dy_at_t(const Quadratic& , double t);
@@ -30,8 +31,8 @@
b -= c; // b = 2*B - 2*C
}
-int quadraticRoots(double A, double B, double C, double t[2]);
-int quadraticRootsX(const double A, const double B, const double C, double s[2]);
+int quadraticRootsReal(double A, double B, double C, double t[2]);
+int quadraticRootsValidT(const double A, const double B, const double C, double s[2]);
void sub_divide(const Quadratic& src, double t1, double t2, Quadratic& dst);
void xy_at_t(const Quadratic& , double t, double& x, double& y);
diff --git a/experimental/Intersection/QuarticRoot.cpp b/experimental/Intersection/QuarticRoot.cpp
index 86ea7a6..6941935 100644
--- a/experimental/Intersection/QuarticRoot.cpp
+++ b/experimental/Intersection/QuarticRoot.cpp
@@ -30,190 +30,48 @@
#include "QuadraticUtilities.h"
#include "QuarticRoot.h"
+int reducedQuarticRoots(const double t4, const double t3, const double t2, const double t1,
+ const double t0, const bool oneHint, double roots[4]) {
#if SK_DEBUG
-#define QUARTIC_DEBUG 1
-#else
-#define QUARTIC_DEBUG 0
-#endif
-
-const double PI = 4 * atan(1);
-
-// unlike quadraticRoots in QuadraticUtilities.cpp, this does not discard
-// real roots <= 0 or >= 1
-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;
- double D = p * p - q;
- if (D < 0) {
- if (approximately_positive_squared(D)) {
- D = 0;
- } else {
- return 0;
- }
- }
- double sqrt_D = sqrt(D);
- if (approximately_less_than_zero(sqrt_D)) {
- s[0] = -p;
- return 1;
- }
- s[0] = sqrt_D - p;
- s[1] = -sqrt_D - p;
- return 2;
-}
-
-#define USE_GEMS 0
-#if USE_GEMS
-// unlike cubicRoots in CubicUtilities.cpp, this does not discard
-// real roots <= 0 or >= 1
-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 = acos(-R / sqrt(-Q3)) / 3;
- 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 = a / 3;
- for (int i = 0; i < num; ++i) {
- s[i] -= sub;
- }
- return num;
-}
-#else
-
-int cubicRootsX(double A, double B, double C, double D, double s[3]) {
-#if QUARTIC_DEBUG
// create a string mathematica understands
+ // GDB set print repe 15 # if repeated digits is a bother
+ // set print elements 400 # if line doesn't fit
char str[1024];
bzero(str, sizeof(str));
- sprintf(str, "Solve[%1.19g x^3 + %1.19g x^2 + %1.19g x + %1.19g == 0, x]", A, B, C, D);
+ sprintf(str, "Solve[%1.19g x^4 + %1.19g x^3 + %1.19g x^2 + %1.19g x + %1.19g == 0, x]",
+ t4, t3, t2, t1, t0);
#endif
- if (approximately_zero(A)) { // we're just a quadratic
- return quadraticRootsX(B, C, D, s);
+ if (approximately_zero(t4)) {
+ if (approximately_zero(t3)) {
+ return quadraticRootsReal(t2, t1, t0, roots);
+ }
+ return cubicRootsReal(t3, t2, t1, t0, roots);
}
- if (approximately_zero(D)) { // 0 is one root
- int num = quadraticRootsX(A, B, C, s);
+ if (approximately_zero(t0)) { // 0 is one root
+ int num = cubicRootsReal(t4, t3, t2, t1, roots);
for (int i = 0; i < num; ++i) {
- if (approximately_zero(s[i])) {
+ if (approximately_zero(roots[i])) {
return num;
}
}
- s[num++] = 0;
+ roots[num++] = 0;
return num;
}
- if (approximately_zero(A + B + C + D)) { // 1 is one root
- int num = quadraticRootsX(A, A + B, -D, s);
+ if (oneHint) {
+ assert(approximately_zero(t4 + t3 + t2 + t1 + t0)); // 1 is one root
+ int num = cubicRootsReal(t4, t4 + t3, -(t1 + t0), -t0, roots); // note that -C==A+B+D+E
for (int i = 0; i < num; ++i) {
- if (approximately_equal(s[i], 1)) {
+ if (approximately_equal(roots[i], 1)) {
return num;
}
}
- s[num++] = 1;
+ roots[num++] = 1;
return num;
}
- double a, b, c;
- {
- double invA = 1 / A;
- a = B * invA;
- b = C * invA;
- c = D * invA;
- }
- double a2 = a * a;
- double Q = (a2 - b * 3) / 9;
- double R = (2 * a2 * a - 9 * a * b + 27 * c) / 54;
- double Q3 = Q * Q * Q;
- double R2MinusQ3 = R * R - Q3;
- double adiv3 = a / 3;
- double r;
- double* roots = s;
-
- if (approximately_zero_squared(R2MinusQ3)) {
- if (approximately_zero(R)) {/* one triple solution */
- *roots++ = -adiv3;
- } else { /* one single and one double solution */
-
- double u = cube_root(-R);
- *roots++ = 2 * u - adiv3;
- *roots++ = -u - adiv3;
- }
- }
- else if (R2MinusQ3 < 0) // we have 3 real roots
- {
- double theta = acos(R / sqrt(Q3));
- double neg2RootQ = -2 * sqrt(Q);
-
- r = neg2RootQ * cos(theta / 3) - adiv3;
- *roots++ = r;
-
- r = neg2RootQ * cos((theta + 2 * PI) / 3) - adiv3;
- *roots++ = r;
-
- r = neg2RootQ * cos((theta - 2 * PI) / 3) - adiv3;
- *roots++ = r;
- }
- else // we have 1 real root
- {
- double A = fabs(R) + sqrt(R2MinusQ3);
- A = cube_root(A);
- if (R > 0) {
- A = -A;
- }
- if (A != 0) {
- A += Q / A;
- }
- r = A - adiv3;
- *roots++ = r;
- }
- return (int)(roots - s);
+ return -1;
}
-#endif
-int quarticRoots(const double A, const double B, const double C, const double D,
+int quarticRootsReal(const double A, const double B, const double C, const double D,
const double E, double s[4]) {
double u, v;
/* normal form: x^4 + Ax^3 + Bx^2 + Cx + D = 0 */
@@ -231,37 +89,97 @@
int num;
if (approximately_zero(r)) {
/* no absolute term: y(y^3 + py + q) = 0 */
- num = cubicRootsX(1, 0, p, q, s);
+ num = cubicRootsReal(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);
+ double cubicRoots[3];
+ int roots = cubicRootsReal(1, -p / 2, -r, r * p / 2 - q * q / 8, cubicRoots);
+ int index;
+ #if 0 && SK_DEBUG // enable to verify that any cubic root is as good as any other
+ double tries[3][4];
+ int nums[3];
+ for (index = 0; index < roots; ++index) {
+ /* ... and take one real solution ... */
+ const double z = cubicRoots[index];
+ /* ... to build two quadric equations */
+ u = z * z - r;
+ v = 2 * z - p;
+ if (approximately_zero_squared(u)) {
+ u = 0;
+ } else if (u > 0) {
+ u = sqrt(u);
+ } else {
+ SkDebugf("%s u=%1.9g <0\n", __FUNCTION__, u);
+ continue;
+ }
+ if (approximately_zero_squared(v)) {
+ v = 0;
+ } else if (v > 0) {
+ v = sqrt(v);
+ } else {
+ SkDebugf("%s v=%1.9g <0\n", __FUNCTION__, v);
+ continue;
+ }
+ nums[index] = quadraticRootsReal(1, q < 0 ? -v : v, z - u, tries[index]);
+ nums[index] += quadraticRootsReal(1, q < 0 ? v : -v, z + u, tries[index] + nums[index]);
+ /* resubstitute */
+ const double sub = a / 4;
+ for (int i = 0; i < nums[index]; ++i) {
+ tries[index][i] -= sub;
+ }
+ }
+ for (index = 0; index < roots; ++index) {
+ SkDebugf("%s", __FUNCTION__);
+ for (int idx2 = 0; idx2 < nums[index]; ++idx2) {
+ SkDebugf(" %1.9g", tries[index][idx2]);
+ }
+ SkDebugf("\n");
+ }
+ #endif
/* ... and take one real solution ... */
- const double z = s[0];
- /* ... to build two quadric equations */
- u = z * z - r;
- v = 2 * z - p;
- if (approximately_zero_squared(u)) {
- u = 0;
- } else if (u > 0) {
- u = sqrt(u);
- } else {
- return 0;
+ double z;
+ num = 0;
+ int num2 = 0;
+ for (index = 0; index < roots; ++index) {
+ z = cubicRoots[index];
+ /* ... to build two quadric equations */
+ u = z * z - r;
+ v = 2 * z - p;
+ if (approximately_zero_squared(u)) {
+ u = 0;
+ } else if (u > 0) {
+ u = sqrt(u);
+ } else {
+ continue;
+ }
+ if (approximately_zero_squared(v)) {
+ v = 0;
+ } else if (v > 0) {
+ v = sqrt(v);
+ } else {
+ continue;
+ }
+ num = quadraticRootsReal(1, q < 0 ? -v : v, z - u, s);
+ num2 = quadraticRootsReal(1, q < 0 ? v : -v, z + u, s + num);
+ if (!((num | num2) & 1)) {
+ break; // prefer solutions without single quad roots
+ }
}
- if (approximately_zero_squared(v)) {
- v = 0;
- } else if (v > 0) {
- v = sqrt(v);
- } else {
- return 0;
+ num += num2;
+ if (!num) {
+ return 0; // no valid cubic root
}
- num = quadraticRootsX(1, q < 0 ? -v : v, z - u, s);
- num += quadraticRootsX(1, q < 0 ? v : -v, z + u, s + num);
+ }
+ /* resubstitute */
+ const double sub = a / 4;
+ for (int i = 0; i < num; ++i) {
+ s[i] -= sub;
}
// eliminate duplicates
for (int i = 0; i < num - 1; ++i) {
for (int j = i + 1; j < num; ) {
- if (approximately_equal(s[i], s[j])) {
+ if (AlmostEqualUlps(s[i], s[j])) {
if (j < --num) {
s[j] = s[num];
}
@@ -270,10 +188,5 @@
}
}
}
- /* resubstitute */
- const double sub = a / 4;
- for (int i = 0; i < num; ++i) {
- s[i] -= sub;
- }
return num;
}
diff --git a/experimental/Intersection/QuarticRoot.h b/experimental/Intersection/QuarticRoot.h
index 8baa842..f76509e 100644
--- a/experimental/Intersection/QuarticRoot.h
+++ b/experimental/Intersection/QuarticRoot.h
@@ -1,2 +1,5 @@
-int quarticRoots(const double A, const double B, const double C, const double D,
+int reducedQuarticRoots(const double t4, const double t3, const double t2, const double t1,
+ const double t0, const bool oneHint, double s[4]);
+
+int quarticRootsReal(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
index 0f310bb..027d19b 100644
--- a/experimental/Intersection/QuarticRoot_Test.cpp
+++ b/experimental/Intersection/QuarticRoot_Test.cpp
@@ -17,23 +17,37 @@
size_t rootECount = sizeof(rootE) / sizeof(rootE[0]);
-static void quadraticTest() {
+static void quadraticTest(bool limit) {
// (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];
+ double B = rootB[bIndex];
+ double C = rootC[cIndex];
+ if (limit) {
+ B = (B - 6) / 12;
+ C = (C - 6) / 12;
+ }
const double b = A * (B + C);
const double c = A * B * C;
double roots[2];
- const int rootCount = quadraticRootsX(A, b, c, roots);
- const int expected = 1 + (B != C);
+ const int rootCount = limit ? quadraticRootsValidT(A, b, c, roots)
+ : quadraticRootsReal(A, b, c, roots);
+ int expected;
+ if (limit) {
+ expected = B <= 0 && B >= -1;
+ expected += B != C && C <= 0 && C >= -1;
+ } else {
+ expected = 1 + (B != C);
+ }
assert(rootCount == expected);
+ if (!rootCount) {
+ continue;
+ }
assert(approximately_equal(roots[0], -B)
|| approximately_equal(roots[0], -C));
- if (B != C) {
+ if (expected > 1) {
assert(!approximately_equal(roots[0], roots[1]));
assert(approximately_equal(roots[1], -B)
|| approximately_equal(roots[1], -C));
@@ -43,45 +57,119 @@
}
}
-static void cubicTest() {
+static void testOneCubic(bool limit, size_t aIndex, size_t bIndex, size_t cIndex, size_t dIndex) {
+ const double A = mulA[aIndex];
+ double B = rootB[bIndex];
+ double C = rootC[cIndex];
+ double D = rootD[dIndex];
+ if (limit) {
+ B = (B - 6) / 12;
+ C = (C - 6) / 12;
+ D = (C - 2) / 6;
+ }
+ 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 = limit ? cubicRootsValidT(A, b, c, d, roots)
+ : cubicRootsReal(A, b, c, d, roots);
+ int expected;
+ if (limit) {
+ expected = B <= 0 && B >= -1;
+ expected += B != C && C <= 0 && C >= -1;
+ expected += B != D && C != D && D <= 0 && D >= -1;
+ } else {
+ expected = 1 + (B != C) + (B != D && C != D);
+ }
+ assert(rootCount == expected);
+ if (!rootCount) {
+ return;
+ }
+ assert(approximately_equal(roots[0], -B)
+ || approximately_equal(roots[0], -C)
+ || approximately_equal(roots[0], -D));
+ if (expected <= 1) {
+ return;
+ }
+ 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) {
+ return;
+ }
+ 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 cubicTest(bool limit) {
// (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 = 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));
- }
- }
+ testOneCubic(limit, aIndex, bIndex, cIndex, dIndex);
}
}
}
}
}
+static void testOneQuartic(size_t aIndex, size_t bIndex, size_t cIndex, size_t dIndex,
+ size_t 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];
+ bool oneHint = approximately_zero(A + b + c + d + e);
+ int rootCount = reducedQuarticRoots(A, b, c, d, e, oneHint, roots);
+ if (rootCount < 0) {
+ rootCount = quarticRootsReal(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(AlmostEqualUlps(roots[0], -B)
+ || AlmostEqualUlps(roots[0], -C)
+ || AlmostEqualUlps(roots[0], -D)
+ || AlmostEqualUlps(roots[0], -E));
+ if (expected <= 1) {
+ return;
+ }
+ assert(!AlmostEqualUlps(roots[0], roots[1]));
+ assert(AlmostEqualUlps(roots[1], -B)
+ || AlmostEqualUlps(roots[1], -C)
+ || AlmostEqualUlps(roots[1], -D)
+ || AlmostEqualUlps(roots[1], -E));
+ if (expected <= 2) {
+ return;
+ }
+ assert(!AlmostEqualUlps(roots[0], roots[2])
+ && !AlmostEqualUlps(roots[1], roots[2]));
+ assert(AlmostEqualUlps(roots[2], -B)
+ || AlmostEqualUlps(roots[2], -C)
+ || AlmostEqualUlps(roots[2], -D)
+ || AlmostEqualUlps(roots[2], -E));
+ if (expected <= 3) {
+ return;
+ }
+ assert(!AlmostEqualUlps(roots[0], roots[3])
+ && !AlmostEqualUlps(roots[1], roots[3])
+ && !AlmostEqualUlps(roots[2], roots[3]));
+ assert(AlmostEqualUlps(roots[3], -B)
+ || AlmostEqualUlps(roots[3], -C)
+ || AlmostEqualUlps(roots[3], -D)
+ || AlmostEqualUlps(roots[3], -E));
+}
+
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
@@ -90,47 +178,7 @@
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 = 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));
- }
- }
- }
+ testOneQuartic(aIndex, bIndex, cIndex, dIndex, eIndex);
}
}
}
@@ -139,7 +187,11 @@
}
void QuarticRoot_Test() {
- quadraticTest();
- cubicTest();
+ testOneCubic(false, 0, 5, 5, 4);
+ testOneQuartic(0, 0, 2, 4, 3);
+ quadraticTest(true);
+ quadraticTest(false);
+ cubicTest(true);
+ cubicTest(false);
quarticTest();
}
diff --git a/experimental/Intersection/TestUtilities.cpp b/experimental/Intersection/TestUtilities.cpp
index 0477ed4..3ae10d2 100644
--- a/experimental/Intersection/TestUtilities.cpp
+++ b/experimental/Intersection/TestUtilities.cpp
@@ -63,4 +63,3 @@
&& ((cubic[0].y <= cubic[1].y && cubic[0].y <= cubic[2].y && cubic[1].y <= cubic[3].y && cubic[2].y <= cubic[3].y)
|| (cubic[0].y >= cubic[1].y && cubic[0].y >= cubic[2].y && cubic[1].y >= cubic[3].y && cubic[2].x >= cubic[3].y));
}
-
diff --git a/experimental/Intersection/qc.htm b/experimental/Intersection/qc.htm
index 3fe3998..2cf3592 100644
--- a/experimental/Intersection/qc.htm
+++ b/experimental/Intersection/qc.htm
@@ -1545,11 +1545,119 @@
{{61.3365082,82.6931328}, {54.8250789,71.6639328}, {47.7274442,61.4049645}},
</div>
+<div id="quad15">
+{{x = 80.897794748143198, y = 49.236332042718459}, {x = 81.082078218891212, y = 64.066749904488631}, {x = 69.972305057149981, y = 72.968595519850993}}
+{{x = 72.503745601281395, y = 32.952320736577882}, {x = 88.030880716061645, y = 38.137194847810164}, {x = 73.193774825517906, y = 67.773492479591397}}
+</div>
+
+<div id="cubic19">
+{{x = 34.560092601254624, y = 51.476349286491221}, {x = 27.498466254909744, y = 66.722346267999313}, {x = 42.500359724508769, y = 3.5458898188294325}, {x = 73.37353619438295, y = 89.022818994253328}}
+{{x = 63.002458057833124, y = 82.312578001205154}, {x = 2.4737262644217006, y = 75.917326135522373}, {x = 95.77018506628005, y = 9.5004089686555826}, {x = 6.5188364156143912, y = 62.083637231068508}}
+</div>
+
+<div id="cubic20">
+ {{x = 42.449716172390481, y = 52.379709366885805}, {x = 27.896043159019225, y = 48.797373636065686}, {x = 92.770268299044233, y = 89.899302036454571}, {x = 12.102066544863426, y = 99.43241951960718}}
+{{x = 45.77532924980639, y = 45.958701495993274}, {x = 37.458701356062065, y = 68.393691335056758}, {x = 37.569326692060258, y = 27.673713456687381}, {x = 60.674866037757539, y = 62.47349659096146}}
+</div>
+
+<div id="cubic21">
+{{x = 26.192053931854691, y = 9.8504326817814416}, {x = 10.174241480498686, y = 98.476562741434464}, {x = 21.177712558385782, y = 33.814968789841501}, {x = 75.329030899018534, y = 55.02231980442177}}
+{{x = 56.222082700683771, y = 24.54395039218662}, {x = 95.589995289030483, y = 81.050822735322086}, {x = 28.180450866082897, y = 28.837706255185282}, {x = 60.128952916771617, y = 87.311672180570511}}
+</div>
+
+<div id="quad16">
+{{x = 67.965974918365831, y = 52.573040929556633}, {x = 67.973015821010591, y = 52.57495862082331}, {x = 67.980057838863502, y = 52.576878275262274}}
+{{x = 67.975025709349239, y = 52.572750461020817}, {x = 67.973101328974863, y = 52.57506284863603}, {x = 67.971173663444745, y = 52.577372136133093}}
+</div>
+
+<div id="quad17">
+{{x = 52.14807018377202, y = 65.012420045148644}, {x = 44.778669050208237, y = 66.315562705604378}, {x = 51.619118408823567, y = 63.787827046262684}}
+{{x = 30.004993234763383, y = 93.921296668202288}, {x = 53.384822003076991, y = 60.732180341802753}, {x = 58.652998934338584, y = 43.111073088306185}}
+</div>
+
+<div id="quad18">
+{{x = 369.850525, y = 145.67596399999999}, {x = 382.36291499999999, y = 121.29286999999999}, {x = 406.21127300000001, y = 121.29286999999999}}
+{{x = 369.962311, y = 137.976044}, {x = 383.97189300000002, y = 121.29286999999999}, {x = 406.21612499999998, y = 121.29286999999999}}
+</div>
+
+<div id="quad19">
+{{x = 406.23635899999999, y = 121.254936}, {x = 409.44567899999998, y = 121.254936}, {x = 412.97595200000001, y = 121.789818}}
+{{x = 406.23599200000001, y = 121.254936}, {x = 425.70590199999998, y = 121.254936}, {x = 439.71994000000001, y = 137.087616}}
+</div>
+
+<div id="cubic22">
+{{x = 7.5374809128872498, y = 82.441702896003477}, {x = 22.444346930107265, y = 22.138854312775123}, {x = 66.76091829629658, y = 50.753805856571446}, {x = 78.193478508942519, y = 97.7932997968948}}
+{{x = 97.700573130371311, y = 53.53260215070685}, {x = 87.72443481149358, y = 84.575876772671876}, {x = 19.215031396232092, y = 47.032676472809484}, {x = 11.989686410869325, y = 10.659507480757082}}
+ {{7.53748091,82.4417029}, {15.5677076,52.942994}, {29.9404074,49.1672596}},
+ {{29.9404074,49.1672596}, {44.3131071,45.3915253}, {58.1067559,59.5061814}},
+ {{58.1067559,59.5061814}, {71.9004047,73.6208375}, {78.1934785,97.7932998}},
+
+ {{97.7005731,53.5326022}, {91.6030843,68.4083459}, {72.6510251,64.2972928}},
+ {{72.6510251,64.2972928}, {53.6989659,60.1862397}, {35.2053722,44.8391126}},
+ {{35.2053722,44.8391126}, {16.7117786,29.4919856}, {11.9896864,10.6595075}},
+</div>
+
+<div id="quad20">
+ {{29.9404074,49.1672596}, {44.3131071,45.3915253}, {58.1067559,59.5061814}},
+ {{72.6510251,64.2972928}, {53.6989659,60.1862397}, {35.2053722,44.8391126}},
+</div>
+
+<div id="cubic23">
+{{x = 32.484981432782945, y = 75.082940782924624}, {x = 42.467313093350882, y = 48.131159948246157}, {x = 3.5963115764764657, y = 43.208665839959245}, {x = 79.442476890721579, y = 89.709102357602262}}
+{{x = 18.98573861410177, y = 93.308887208490106}, {x = 40.405250173250792, y = 91.039661826118675}, {x = 8.0467721950480584, y = 42.100282172719147}, {x = 40.883324221187891, y = 26.030185504830527}}
+ {{32.4849814,75.0829408}, {35.4553509,65.5763004}, {33.5767697,60.2097835}},
+ {{33.5767697,60.2097835}, {31.6981886,54.8432666}, {31.1663962,54.7302484}},
+ {{31.1663962,54.7302484}, {31.1662882,54.7301074}, {31.1663969,54.7302485}},
+ {{31.1663969,54.7302485}, {30.4117445,54.6146017}, {40.1631726,62.9428436}},
+ {{40.1631726,62.9428436}, {49.9146008,71.2710854}, {79.4424769,89.7091024}},
+
+ {{18.9857386,93.3088872}, {25.7662938,92.3417699}, {26.5917262,85.8225583}},
+ {{26.5917262,85.8225583}, {27.4171586,79.3033467}, {26.141946,69.8089528}},
+ {{26.141946,69.8089528}, {24.2922348,57.665767}, {26.0404936,45.4260361}},
+ {{26.0404936,45.4260361}, {27.7887523,33.1863051}, {40.8833242,26.0301855}},
+</div>
+
+<div id="quad21">
+ {{31.1663962,54.7302484}, {31.1662882,54.7301074}, {31.1663969,54.7302485}},
+ {{26.0404936,45.4260361}, {27.7887523,33.1863051}, {40.8833242,26.0301855}},
+</div>
+
+<div id="cubic24">
+{{x = 65.454505973241524, y = 93.881892270353575}, {x = 45.867360264932437, y = 92.723972719499827}, {x = 2.1464054482739447, y = 74.636369140183717}, {x = 33.774068594804994, y = 40.770872887582925}}
+{{x = 72.963387832494163, y = 95.659300729473728}, {x = 11.809496633619768, y = 82.209921247423594}, {x = 13.456139067865974, y = 57.329313623406605}, {x = 36.060621606214262, y = 70.867335643091849}}
+ {{65.454506,93.8818923}, {54.7397995,93.2922678}, {41.5072916,87.1234036}},
+ {{41.5072916,87.1234036}, {28.2747836,80.9545395}, {23.5780771,69.3344126}},
+ {{23.5780771,69.3344126}, {18.8813706,57.7142857}, {33.7740686,40.7708729}},
+
+ {{72.9633878,95.6593007}, {42.7738746,88.4730382}, {31.1932785,80.2458029}},
+ {{31.1932785,80.2458029}, {19.6126823,72.0185676}, {21.9918152,68.2892325}},
+ {{21.9918152,68.2892325}, {24.370948,64.5598974}, {36.0606216,70.8673356}},
+</div>
+
+<div id="quad22">
+ {{41.5072916,87.1234036}, {28.2747836,80.9545395}, {23.5780771,69.3344126}},
+ {{72.9633878,95.6593007}, {42.7738746,88.4730382}, {31.1932785,80.2458029}},
+</div>
+
</div>
<script type="text/javascript">
var testDivs = [
+ quad22,
+ cubic24,
+ quad21,
+ cubic23,
+ quad20,
+ cubic22,
+ quad19,
+ quad18,
+ quad17,
+ quad16,
+ cubic21,
+ cubic20,
+ cubic19,
+ quad15,
quad14,
cubic18,
cubic17,