caryclark@google.com | c682590 | 2012-02-03 22:07:47 +0000 | [diff] [blame] | 1 | #include "CubicUtilities.h" |
| 2 | #include "DataTypes.h" |
| 3 | #include "QuadraticUtilities.h" |
| 4 | |
| 5 | void coefficients(const double* cubic, double& A, double& B, double& C, double& D) { |
| 6 | A = cubic[6]; // d |
| 7 | B = cubic[4] * 3; // 3*c |
| 8 | C = cubic[2] * 3; // 3*b |
| 9 | D = cubic[0]; // a |
| 10 | A -= D - C + B; // A = -a + 3*b - 3*c + d |
| 11 | B += 3 * D - 2 * C; // B = 3*a - 6*b + 3*c |
| 12 | C -= 3 * D; // C = -3*a + 3*b |
| 13 | } |
| 14 | |
| 15 | // cubic roots |
| 16 | |
| 17 | const double PI = 4 * atan(1); |
| 18 | |
| 19 | static bool is_unit_interval(double x) { |
| 20 | return x > 0 && x < 1; |
| 21 | } |
| 22 | |
| 23 | // from SkGeometry.cpp (and Numeric Solutions, 5.6) |
| 24 | int cubicRoots(double A, double B, double C, double D, double t[3]) { |
| 25 | if (approximately_zero(A)) { // we're just a quadratic |
| 26 | return quadraticRoots(B, C, D, t); |
| 27 | } |
| 28 | double a, b, c; |
| 29 | { |
| 30 | double invA = 1 / A; |
| 31 | a = B * invA; |
| 32 | b = C * invA; |
| 33 | c = D * invA; |
| 34 | } |
| 35 | double a2 = a * a; |
| 36 | double Q = (a2 - b * 3) / 9; |
| 37 | double R = (2 * a2 * a - 9 * a * b + 27 * c) / 54; |
| 38 | double Q3 = Q * Q * Q; |
| 39 | double R2MinusQ3 = R * R - Q3; |
| 40 | double adiv3 = a / 3; |
| 41 | double* roots = t; |
| 42 | double r; |
| 43 | |
| 44 | if (R2MinusQ3 < 0) // we have 3 real roots |
| 45 | { |
| 46 | double theta = acos(R / sqrt(Q3)); |
| 47 | double neg2RootQ = -2 * sqrt(Q); |
| 48 | |
| 49 | r = neg2RootQ * cos(theta / 3) - adiv3; |
| 50 | if (is_unit_interval(r)) |
| 51 | *roots++ = r; |
| 52 | |
| 53 | r = neg2RootQ * cos((theta + 2 * PI) / 3) - adiv3; |
| 54 | if (is_unit_interval(r)) |
| 55 | *roots++ = r; |
| 56 | |
| 57 | r = neg2RootQ * cos((theta - 2 * PI) / 3) - adiv3; |
| 58 | if (is_unit_interval(r)) |
| 59 | *roots++ = r; |
| 60 | } |
| 61 | else // we have 1 real root |
| 62 | { |
| 63 | double A = fabs(R) + sqrt(R2MinusQ3); |
| 64 | A = cube_root(A); |
| 65 | if (R > 0) { |
| 66 | A = -A; |
| 67 | } |
| 68 | if (A != 0) { |
| 69 | A += Q / A; |
| 70 | } |
| 71 | r = A - adiv3; |
| 72 | if (is_unit_interval(r)) |
| 73 | *roots++ = r; |
| 74 | } |
| 75 | return (int)(roots - t); |
| 76 | } |
caryclark@google.com | 8dcf114 | 2012-07-02 20:27:02 +0000 | [diff] [blame^] | 77 | |
| 78 | // from http://www.cs.sunysb.edu/~qin/courses/geometry/4.pdf |
| 79 | // c(t) = a(1-t)^3 + 3bt(1-t)^2 + 3c(1-t)t^2 + dt^3 |
| 80 | // c'(t) = -3a(1-t)^2 + 3b((1-t)^2 - 2t(1-t)) + 3c(2t(1-t) - t^2) + 3dt^2 |
| 81 | // = 3(b-a)(1-t)^2 + 6(c-b)t(1-t) + 3(d-c)t^2 |
| 82 | double derivativeAtT(const double* cubic, double t) { |
| 83 | double one_t = 1 - t; |
| 84 | double a = cubic[0]; |
| 85 | double b = cubic[2]; |
| 86 | double c = cubic[4]; |
| 87 | double d = cubic[6]; |
| 88 | return (b - a) * one_t * one_t + 2 * (c - b) * t * one_t + (d - c) * t * t; |
| 89 | } |
| 90 | |
| 91 | // same as derivativeAtT |
| 92 | // which is more accurate? which is faster? |
| 93 | double derivativeAtT_2(const double* cubic, double t) { |
| 94 | double a = cubic[2] - cubic[0]; |
| 95 | double b = cubic[4] - 2 * cubic[2] + cubic[0]; |
| 96 | double c = cubic[6] + 3 * (cubic[2] - cubic[4]) - cubic[0]; |
| 97 | return c * c * t * t + 2 * b * t + a; |
| 98 | } |
| 99 | |
| 100 | void dxdy_at_t(const Cubic& cubic, double t, double& dx, double& dy) { |
| 101 | if (&dx) { |
| 102 | dx = derivativeAtT(&cubic[0].x, t); |
| 103 | } |
| 104 | if (&dy) { |
| 105 | dy = derivativeAtT(&cubic[0].y, t); |
| 106 | } |
| 107 | } |
| 108 | |
| 109 | bool rotate(const Cubic& cubic, int zero, int index, Cubic& rotPath) { |
| 110 | double dy = cubic[index].y - cubic[zero].y; |
| 111 | double dx = cubic[index].x - cubic[zero].x; |
| 112 | if (approximately_equal(dy, 0)) { |
| 113 | if (approximately_equal(dx, 0)) { |
| 114 | return false; |
| 115 | } |
| 116 | memcpy(rotPath, cubic, sizeof(Cubic)); |
| 117 | return true; |
| 118 | } |
| 119 | for (int index = 0; index < 4; ++index) { |
| 120 | rotPath[index].x = cubic[index].x * dx + cubic[index].y * dy; |
| 121 | rotPath[index].y = cubic[index].y * dx - cubic[index].x * dy; |
| 122 | } |
| 123 | return true; |
| 124 | } |
| 125 | |
| 126 | double secondDerivativeAtT(const double* cubic, double t) { |
| 127 | double a = cubic[0]; |
| 128 | double b = cubic[2]; |
| 129 | double c = cubic[4]; |
| 130 | double d = cubic[6]; |
| 131 | return (c - 2 * b + a) * (1 - t) + (d - 2 * c + b) * t; |
| 132 | } |
| 133 | |
| 134 | void xy_at_t(const Cubic& cubic, double t, double& x, double& y) { |
| 135 | double one_t = 1 - t; |
| 136 | double one_t2 = one_t * one_t; |
| 137 | double a = one_t2 * one_t; |
| 138 | double b = 3 * one_t2 * t; |
| 139 | double t2 = t * t; |
| 140 | double c = 3 * one_t * t2; |
| 141 | double d = t2 * t; |
| 142 | if (&x) { |
| 143 | x = a * cubic[0].x + b * cubic[1].x + c * cubic[2].x + d * cubic[3].x; |
| 144 | } |
| 145 | if (&y) { |
| 146 | y = a * cubic[0].y + b * cubic[1].y + c * cubic[2].y + d * cubic[3].y; |
| 147 | } |
| 148 | } |