diff options
author | 2012-02-03 22:07:47 +0000 | |
---|---|---|
committer | 2012-02-03 22:07:47 +0000 | |
commit | c682590538a27d73489bc91c098e000fdfb07ccf (patch) | |
tree | 90b03195e3a74cf47f4fa3a385e99575c46da37b /experimental/Intersection/CubicUtilities.cpp | |
parent | 2c23708e4478a83dcded2e9d5672bc57ee016919 (diff) |
save work in progress
git-svn-id: http://skia.googlecode.com/svn/trunk@3141 2bbb7eff-a529-9590-31e7-b0007b416f81
Diffstat (limited to 'experimental/Intersection/CubicUtilities.cpp')
-rw-r--r-- | experimental/Intersection/CubicUtilities.cpp | 76 |
1 files changed, 76 insertions, 0 deletions
diff --git a/experimental/Intersection/CubicUtilities.cpp b/experimental/Intersection/CubicUtilities.cpp new file mode 100644 index 0000000000..3fab29ec80 --- /dev/null +++ b/experimental/Intersection/CubicUtilities.cpp @@ -0,0 +1,76 @@ +#include "CubicUtilities.h" +#include "DataTypes.h" +#include "QuadraticUtilities.h" + +void coefficients(const double* cubic, double& A, double& B, double& C, double& D) { + A = cubic[6]; // d + B = cubic[4] * 3; // 3*c + C = cubic[2] * 3; // 3*b + D = cubic[0]; // a + A -= D - C + B; // A = -a + 3*b - 3*c + d + B += 3 * D - 2 * C; // B = 3*a - 6*b + 3*c + C -= 3 * D; // C = -3*a + 3*b +} + +// cubic roots + +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]) { + if (approximately_zero(A)) { // we're just a quadratic + return quadraticRoots(B, C, D, t); + } + 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* roots = t; + double r; + + 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; + if (is_unit_interval(r)) + *roots++ = r; + + r = neg2RootQ * cos((theta + 2 * PI) / 3) - adiv3; + if (is_unit_interval(r)) + *roots++ = r; + + r = neg2RootQ * cos((theta - 2 * PI) / 3) - adiv3; + if (is_unit_interval(r)) + *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; + if (is_unit_interval(r)) + *roots++ = r; + } + return (int)(roots - t); +} |