From 9f60291c5375457f8adf228dbe6e8ff1186b13e1 Mon Sep 17 00:00:00 2001 From: "caryclark@google.com" Date: Thu, 24 Jan 2013 21:47:16 +0000 Subject: [PATCH] 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 --- experimental/Intersection/CubeRoot.cpp | 2 +- .../Intersection/CubicIntersection.cpp | 73 +++-- .../Intersection/CubicIntersection_Test.cpp | 78 ++++- experimental/Intersection/CubicUtilities.cpp | 124 +++++++- experimental/Intersection/CubicUtilities.h | 7 +- experimental/Intersection/DataTypes.h | 56 +++- .../Intersection/Intersection_Tests.cpp | 4 +- .../Intersection/Intersection_Tests.h | 1 + .../Intersection/LineCubicIntersection.cpp | 6 +- .../LineQuadraticIntersection.cpp | 6 +- .../Intersection/QuadraticImplicit.cpp | 88 ++---- .../QuadraticIntersection_Test.cpp | 81 +++-- .../Intersection/QuadraticReduceOrder.cpp | 2 +- .../Intersection/QuadraticUtilities.cpp | 83 ++++- .../Intersection/QuadraticUtilities.h | 5 +- experimental/Intersection/QuarticRoot.cpp | 293 ++++++------------ experimental/Intersection/QuarticRoot.h | 5 +- .../Intersection/QuarticRoot_Test.cpp | 206 +++++++----- experimental/Intersection/TestUtilities.cpp | 1 - experimental/Intersection/qc.htm | 108 +++++++ 20 files changed, 818 insertions(+), 411 deletions(-) diff --git a/experimental/Intersection/CubeRoot.cpp b/experimental/Intersection/CubeRoot.cpp index 82d2732f6..5f785a035 100644 --- a/experimental/Intersection/CubeRoot.cpp +++ b/experimental/Intersection/CubeRoot.cpp @@ -374,7 +374,7 @@ static int _tmain() #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 4ae0d84e5..a5b4261df 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 @@ bool intersect(const Cubic& c1, const Cubic& c2, Intersections& i) { #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 @@ static bool intersect2(const Cubic& cubic1, double t1s, double t1e, const Cubic& 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 @@ static bool intersectEnd(const Cubic& cubic1, bool start, const Cubic& cubic2, c 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 @@ static bool intersectEnd(const Cubic& cubic1, bool start, const Cubic& cubic2, c // 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 8c2263e3a..2e738fc0e 100644 --- a/experimental/Intersection/CubicIntersection_Test.cpp +++ b/experimental/Intersection/CubicIntersection_Test.cpp @@ -57,23 +57,29 @@ void CubicIntersection_Test() { } } +#define ONE_OFF_DEBUG 1 + static void oneOff(const Cubic& cubic1, const Cubic& cubic2) { SkTDArray 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 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 @@ static void oneOff(const Cubic& cubic1, const Cubic& cubic2) { 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 @@ const size_t testSetCount = sizeof(testSet) / sizeof(testSet[0]); 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 @@ void CubicIntersection_RandTest() { 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 19f16c6aa..3e2f474d7 100644 --- a/experimental/Intersection/CubicUtilities.cpp +++ b/experimental/Intersection/CubicUtilities.cpp @@ -21,7 +21,7 @@ double calcPrecision(const Cubic& cubic) { #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 @@ void coefficients(const double* cubic, double& A, double& B, double& C, double& 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 @@ int cubicRoots(double A, double B, double C, double D, double t[3]) { *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 @@ int find_cubic_inflections(const Cubic& src, double tValues[]) 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 5205574d2..186ed43b6 100644 --- a/experimental/Intersection/CubicUtilities.h +++ b/experimental/Intersection/CubicUtilities.h @@ -15,13 +15,16 @@ double calcPrecision(const Cubic& cubic); 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& quadratics); void cubic_to_quadratics(const Cubic& cubic, double precision, SkTDArray& 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 afdb8e918..de763edcc 100644 --- a/experimental/Intersection/DataTypes.h +++ b/experimental/Intersection/DataTypes.h @@ -24,6 +24,8 @@ inline bool AlmostEqualUlps(double A, double B) { return AlmostEqualUlps((float) // 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 @@ inline bool approximately_zero(float x) { 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 @@ inline bool approximately_zero_sqrt(double x) { 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 @@ struct _Point { } 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 @@ struct _Point { && 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 @@ struct QuadraticPair { _Point pts[5]; }; +// FIXME: move these into SkTypes.h +template inline T SkTMax(T a, T b) { + if (a < b) + a = b; + return a; +} + +template 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 cfa340f59..bf86c07c7 100644 --- a/experimental/Intersection/Intersection_Tests.cpp +++ b/experimental/Intersection/Intersection_Tests.cpp @@ -15,9 +15,10 @@ void cubecode_test(int test); 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 @@ void Intersection_Tests() { 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 6fc0c754a..4bf3068ef 100644 --- a/experimental/Intersection/Intersection_Tests.h +++ b/experimental/Intersection/Intersection_Tests.h @@ -11,6 +11,7 @@ void ConvexHull_Test(); 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 333e30f0e..cc63099fd 100644 --- a/experimental/Intersection/LineCubicIntersection.cpp +++ b/experimental/Intersection/LineCubicIntersection.cpp @@ -96,7 +96,7 @@ int intersectRay(double roots[3]) { } 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 @@ int horizontalIntersect(double axisIntercept, double roots[3]) { 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 @@ int verticalIntersect(double axisIntercept, double roots[3]) { 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 00c355432..9ec3fb75e 100644 --- a/experimental/Intersection/LineQuadraticIntersection.cpp +++ b/experimental/Intersection/LineQuadraticIntersection.cpp @@ -124,7 +124,7 @@ int intersectRay(double roots[2]) { 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 @@ int horizontalIntersect(double axisIntercept, double roots[2]) { 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 @@ int verticalIntersect(double axisIntercept, double roots[2]) { 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 2d9d9b9a8..9cd24e9f5 100644 --- a/experimental/Intersection/QuadraticImplicit.cpp +++ b/experimental/Intersection/QuadraticImplicit.cpp @@ -13,8 +13,6 @@ #include "QuadraticUtilities.h" #include "TSearch.h" -#include // 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 @@ static int findRoots(const QuadImplicitForm& i, const Quadratic& q2, double root + 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 @@ static bool addIntercept(const Quadratic& q1, const Quadratic& q2, double tMin, 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 bool isLinearInner(const Quadratic& q1, double t1s, double t1e, const Qua } 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 @@ static bool isLinear(const Quadratic& q1, const Quadratic& q2, Intersections& i) 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 @@ bool intersect2(const Quadratic& q1, const Quadratic& q2, Intersections& i) { 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 432f61428..055a18e82 100644 --- a/experimental/Intersection/QuadraticIntersection_Test.cpp +++ b/experimental/Intersection/QuadraticIntersection_Test.cpp @@ -59,6 +59,21 @@ static void standardTestCases() { } 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 @@ static const Quadratic testSet[] = { 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 aa1d0585d..096f6b664 100644 --- a/experimental/Intersection/QuadraticReduceOrder.cpp +++ b/experimental/Intersection/QuadraticReduceOrder.cpp @@ -92,7 +92,7 @@ static int check_linear(const Quadratic& quad, Quadratic& reduction, 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 8b02b4ce6..475ec13b7 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 /* @@ -20,12 +21,37 @@ and using the roots */ + +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 @@ int quadraticRoots(double A, double B, double C, double t[2]) { 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 0fe8be837..f423dc6ec 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 @@ inline void set_abc(const double* quad, double& a, double& b, double& c) { 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 86ea7a63f..6941935f4 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 quarticRoots(const double A, const double B, const double C, const double D, 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 @@ int quarticRoots(const double A, const double B, const double C, const double D, } } } - /* 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 8baa842f3..f76509e2f 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 0f310bb08..027d19bd7 100644 --- a/experimental/Intersection/QuarticRoot_Test.cpp +++ b/experimental/Intersection/QuarticRoot_Test.cpp @@ -17,23 +17,37 @@ double rootE[] = {-5, -1, 0, 1, 7}; 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 quadraticTest() { } } -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 @@ static void quarticTest() { 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 @@ static void quarticTest() { } 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 0477ed48e..3ae10d216 100644 --- a/experimental/Intersection/TestUtilities.cpp +++ b/experimental/Intersection/TestUtilities.cpp @@ -63,4 +63,3 @@ bool controls_inside(const Cubic& cubic) { && ((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 3fe3998b0..2cf3592e5 100644 --- a/experimental/Intersection/qc.htm +++ b/experimental/Intersection/qc.htm @@ -1545,11 +1545,119 @@ $7 = {{x = 24.006224853920855, y = 72.621119847810419}, {x = 29.758671200376888, {{61.3365082,82.6931328}, {54.8250789,71.6639328}, {47.7274442,61.4049645}}, +
+{{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}} +
+ +
+{{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}} +
+ +
+ {{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}} +
+ +
+{{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}} +
+ +
+{{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}} +
+ +
+{{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}} +
+ +
+{{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}} +
+ +
+{{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}} +
+ +
+{{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}}, +
+ +
+ {{29.9404074,49.1672596}, {44.3131071,45.3915253}, {58.1067559,59.5061814}}, + {{72.6510251,64.2972928}, {53.6989659,60.1862397}, {35.2053722,44.8391126}}, +
+ +
+{{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}}, +
+ +
+ {{31.1663962,54.7302484}, {31.1662882,54.7301074}, {31.1663969,54.7302485}}, + {{26.0404936,45.4260361}, {27.7887523,33.1863051}, {40.8833242,26.0301855}}, +
+ +
+{{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}}, +
+ +
+ {{41.5072916,87.1234036}, {28.2747836,80.9545395}, {23.5780771,69.3344126}}, + {{72.9633878,95.6593007}, {42.7738746,88.4730382}, {31.1932785,80.2458029}}, +
+