c682590538
git-svn-id: http://skia.googlecode.com/svn/trunk@3141 2bbb7eff-a529-9590-31e7-b0007b416f81
77 lines
1.9 KiB
C++
77 lines
1.9 KiB
C++
#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);
|
|
}
|