AuroraMiddleware/skia2
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

shape ops work in progress

Browse Source 
first 100,000 random cubic/cubic intersections working

git-svn-id: http://skia.googlecode.com/svn/trunk@7380 2bbb7eff-a529-9590-31e7-b0007b416f81
This commit is contained in:
caryclark@google.com 2013-01-24 21:47:16 +00:00
parent 6c55b513bf
commit 9f60291c53
20 changed files with 818 additions and 411 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
experimental/Intersection/CubeRoot.cpp
View File

@ -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));

73
experimental/Intersection/CubicIntersection.cpp
View File

@ -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;

78
experimental/Intersection/CubicIntersection_Test.cpp
View File

@ -57,23 +57,29 @@ void CubicIntersection_Test() {
}
}
#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 @@ 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]));
}
}

124
experimental/Intersection/CubicUtilities.cpp
View File

@ -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));

7
experimental/Intersection/CubicUtilities.h
View File

@ -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<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);

56
experimental/Intersection/DataTypes.h
View File

@ -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 <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__

4
experimental/Intersection/Intersection_Tests.cpp
View File

@ -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();

1
experimental/Intersection/Intersection_Tests.h
View File

@ -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();

6
experimental/Intersection/LineCubicIntersection.cpp
View File

@ -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) {

6
experimental/Intersection/LineQuadraticIntersection.cpp
View File

@ -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) {

88
experimental/Intersection/QuadraticImplicit.cpp
View File

@ -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 @@ 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) {

81
experimental/Intersection/QuadraticIntersection_Test.cpp
View File

@ -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);
}
}
}

2
experimental/Intersection/QuadraticReduceOrder.cpp
View File

@ -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) {

83
experimental/Intersection/QuadraticUtilities.cpp
View File

@ -5,6 +5,7 @@
* found in the LICENSE file.
*/
#include "QuadraticUtilities.h"
#include "SkTypes.h"
#include <math.h>
/*
@ -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;

5
experimental/Intersection/QuadraticUtilities.h
View File

@ -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);

293
experimental/Intersection/QuarticRoot.cpp
View File

@ -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;
}

5
experimental/Intersection/QuarticRoot.h
View File

@ -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]);

206
experimental/Intersection/QuarticRoot_Test.cpp
View File

@ -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();
}

1
experimental/Intersection/TestUtilities.cpp
View File

@ -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));
}

108
experimental/Intersection/qc.htm
View File

@ -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}},
</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,