caryclark@google.com | 9e49fb6 | 2012-08-27 14:11:33 +0000 | [diff] [blame] | 1 | /* |
| 2 | * Copyright 2012 Google Inc. |
| 3 | * |
| 4 | * Use of this source code is governed by a BSD-style license that can be |
| 5 | * found in the LICENSE file. |
| 6 | */ |
caryclark@google.com | c682590 | 2012-02-03 22:07:47 +0000 | [diff] [blame] | 7 | #include "CubicUtilities.h" |
caryclark@google.com | 45a8fc6 | 2013-02-14 15:29:11 +0000 | [diff] [blame] | 8 | #include "Extrema.h" |
caryclark@google.com | 1304bb2 | 2013-03-13 20:29:41 +0000 | [diff] [blame] | 9 | #include "LineUtilities.h" |
caryclark@google.com | c682590 | 2012-02-03 22:07:47 +0000 | [diff] [blame] | 10 | #include "QuadraticUtilities.h" |
| 11 | |
caryclark@google.com | 7ff5c84 | 2013-02-26 15:56:05 +0000 | [diff] [blame] | 12 | const int gPrecisionUnit = 256; // FIXME: arbitrary -- should try different values in test framework |
caryclark@google.com | 05c4bad | 2013-01-19 13:22:39 +0000 | [diff] [blame] | 13 | |
| 14 | // FIXME: cache keep the bounds and/or precision with the caller? |
caryclark@google.com | 73ca624 | 2013-01-17 21:02:47 +0000 | [diff] [blame] | 15 | double calcPrecision(const Cubic& cubic) { |
| 16 | _Rect dRect; |
caryclark@google.com | 05c4bad | 2013-01-19 13:22:39 +0000 | [diff] [blame] | 17 | dRect.setBounds(cubic); // OPTIMIZATION: just use setRawBounds ? |
caryclark@google.com | 73ca624 | 2013-01-17 21:02:47 +0000 | [diff] [blame] | 18 | double width = dRect.right - dRect.left; |
| 19 | double height = dRect.bottom - dRect.top; |
caryclark@google.com | 7ff5c84 | 2013-02-26 15:56:05 +0000 | [diff] [blame] | 20 | return (width > height ? width : height) / gPrecisionUnit; |
caryclark@google.com | 73ca624 | 2013-01-17 21:02:47 +0000 | [diff] [blame] | 21 | } |
| 22 | |
mtklein | 753b870 | 2014-08-20 07:38:46 -0700 | [diff] [blame] | 23 | #ifdef SK_DEBUG |
caryclark@google.com | 9d5f99b | 2013-01-22 12:55:54 +0000 | [diff] [blame] | 24 | double calcPrecision(const Cubic& cubic, double t, double scale) { |
| 25 | Cubic part; |
caryclark@google.com | 9f60291 | 2013-01-24 21:47:16 +0000 | [diff] [blame] | 26 | sub_divide(cubic, SkTMax(0., t - scale), SkTMin(1., t + scale), part); |
caryclark@google.com | 9d5f99b | 2013-01-22 12:55:54 +0000 | [diff] [blame] | 27 | return calcPrecision(part); |
| 28 | } |
| 29 | #endif |
| 30 | |
caryclark@google.com | 1304bb2 | 2013-03-13 20:29:41 +0000 | [diff] [blame] | 31 | bool clockwise(const Cubic& c) { |
| 32 | double sum = (c[0].x - c[3].x) * (c[0].y + c[3].y); |
| 33 | for (int idx = 0; idx < 3; ++idx){ |
| 34 | sum += (c[idx + 1].x - c[idx].x) * (c[idx + 1].y + c[idx].y); |
| 35 | } |
| 36 | return sum <= 0; |
| 37 | } |
caryclark@google.com | 9d5f99b | 2013-01-22 12:55:54 +0000 | [diff] [blame] | 38 | |
caryclark@google.com | c682590 | 2012-02-03 22:07:47 +0000 | [diff] [blame] | 39 | void coefficients(const double* cubic, double& A, double& B, double& C, double& D) { |
| 40 | A = cubic[6]; // d |
| 41 | B = cubic[4] * 3; // 3*c |
| 42 | C = cubic[2] * 3; // 3*b |
| 43 | D = cubic[0]; // a |
| 44 | A -= D - C + B; // A = -a + 3*b - 3*c + d |
| 45 | B += 3 * D - 2 * C; // B = 3*a - 6*b + 3*c |
| 46 | C -= 3 * D; // C = -3*a + 3*b |
| 47 | } |
| 48 | |
caryclark@google.com | 1304bb2 | 2013-03-13 20:29:41 +0000 | [diff] [blame] | 49 | bool controls_contained_by_ends(const Cubic& c) { |
| 50 | _Vector startTan = c[1] - c[0]; |
| 51 | if (startTan.x == 0 && startTan.y == 0) { |
| 52 | startTan = c[2] - c[0]; |
| 53 | } |
| 54 | _Vector endTan = c[2] - c[3]; |
| 55 | if (endTan.x == 0 && endTan.y == 0) { |
| 56 | endTan = c[1] - c[3]; |
| 57 | } |
| 58 | if (startTan.dot(endTan) >= 0) { |
| 59 | return false; |
| 60 | } |
| 61 | _Line startEdge = {c[0], c[0]}; |
| 62 | startEdge[1].x -= startTan.y; |
| 63 | startEdge[1].y += startTan.x; |
| 64 | _Line endEdge = {c[3], c[3]}; |
| 65 | endEdge[1].x -= endTan.y; |
| 66 | endEdge[1].y += endTan.x; |
| 67 | double leftStart1 = is_left(startEdge, c[1]); |
| 68 | if (leftStart1 * is_left(startEdge, c[2]) < 0) { |
| 69 | return false; |
| 70 | } |
| 71 | double leftEnd1 = is_left(endEdge, c[1]); |
| 72 | if (leftEnd1 * is_left(endEdge, c[2]) < 0) { |
| 73 | return false; |
| 74 | } |
| 75 | return leftStart1 * leftEnd1 >= 0; |
| 76 | } |
| 77 | |
| 78 | bool ends_are_extrema_in_x_or_y(const Cubic& c) { |
| 79 | return (between(c[0].x, c[1].x, c[3].x) && between(c[0].x, c[2].x, c[3].x)) |
| 80 | || (between(c[0].y, c[1].y, c[3].y) && between(c[0].y, c[2].y, c[3].y)); |
| 81 | } |
| 82 | |
| 83 | bool monotonic_in_y(const Cubic& c) { |
| 84 | return between(c[0].y, c[1].y, c[3].y) && between(c[0].y, c[2].y, c[3].y); |
| 85 | } |
| 86 | |
| 87 | bool serpentine(const Cubic& c) { |
| 88 | if (!controls_contained_by_ends(c)) { |
| 89 | return false; |
| 90 | } |
| 91 | double wiggle = (c[0].x - c[2].x) * (c[0].y + c[2].y); |
| 92 | for (int idx = 0; idx < 2; ++idx){ |
| 93 | wiggle += (c[idx + 1].x - c[idx].x) * (c[idx + 1].y + c[idx].y); |
| 94 | } |
| 95 | double waggle = (c[1].x - c[3].x) * (c[1].y + c[3].y); |
| 96 | for (int idx = 1; idx < 3; ++idx){ |
| 97 | waggle += (c[idx + 1].x - c[idx].x) * (c[idx + 1].y + c[idx].y); |
| 98 | } |
| 99 | return wiggle * waggle < 0; |
| 100 | } |
| 101 | |
caryclark@google.com | c682590 | 2012-02-03 22:07:47 +0000 | [diff] [blame] | 102 | // cubic roots |
| 103 | |
| 104 | const double PI = 4 * atan(1); |
| 105 | |
caryclark@google.com | c682590 | 2012-02-03 22:07:47 +0000 | [diff] [blame] | 106 | // from SkGeometry.cpp (and Numeric Solutions, 5.6) |
caryclark@google.com | 9f60291 | 2013-01-24 21:47:16 +0000 | [diff] [blame] | 107 | int cubicRootsValidT(double A, double B, double C, double D, double t[3]) { |
| 108 | #if 0 |
caryclark@google.com | c682590 | 2012-02-03 22:07:47 +0000 | [diff] [blame] | 109 | if (approximately_zero(A)) { // we're just a quadratic |
caryclark@google.com | 9f60291 | 2013-01-24 21:47:16 +0000 | [diff] [blame] | 110 | return quadraticRootsValidT(B, C, D, t); |
caryclark@google.com | c682590 | 2012-02-03 22:07:47 +0000 | [diff] [blame] | 111 | } |
| 112 | double a, b, c; |
| 113 | { |
| 114 | double invA = 1 / A; |
| 115 | a = B * invA; |
| 116 | b = C * invA; |
| 117 | c = D * invA; |
| 118 | } |
| 119 | double a2 = a * a; |
| 120 | double Q = (a2 - b * 3) / 9; |
| 121 | double R = (2 * a2 * a - 9 * a * b + 27 * c) / 54; |
| 122 | double Q3 = Q * Q * Q; |
| 123 | double R2MinusQ3 = R * R - Q3; |
| 124 | double adiv3 = a / 3; |
| 125 | double* roots = t; |
| 126 | double r; |
| 127 | |
| 128 | if (R2MinusQ3 < 0) // we have 3 real roots |
| 129 | { |
| 130 | double theta = acos(R / sqrt(Q3)); |
| 131 | double neg2RootQ = -2 * sqrt(Q); |
| 132 | |
| 133 | r = neg2RootQ * cos(theta / 3) - adiv3; |
| 134 | if (is_unit_interval(r)) |
| 135 | *roots++ = r; |
| 136 | |
| 137 | r = neg2RootQ * cos((theta + 2 * PI) / 3) - adiv3; |
| 138 | if (is_unit_interval(r)) |
| 139 | *roots++ = r; |
| 140 | |
| 141 | r = neg2RootQ * cos((theta - 2 * PI) / 3) - adiv3; |
| 142 | if (is_unit_interval(r)) |
| 143 | *roots++ = r; |
| 144 | } |
| 145 | else // we have 1 real root |
| 146 | { |
| 147 | double A = fabs(R) + sqrt(R2MinusQ3); |
| 148 | A = cube_root(A); |
| 149 | if (R > 0) { |
| 150 | A = -A; |
| 151 | } |
| 152 | if (A != 0) { |
| 153 | A += Q / A; |
| 154 | } |
| 155 | r = A - adiv3; |
| 156 | if (is_unit_interval(r)) |
| 157 | *roots++ = r; |
| 158 | } |
| 159 | return (int)(roots - t); |
caryclark@google.com | 9f60291 | 2013-01-24 21:47:16 +0000 | [diff] [blame] | 160 | #else |
| 161 | double s[3]; |
| 162 | int realRoots = cubicRootsReal(A, B, C, D, s); |
| 163 | int foundRoots = add_valid_ts(s, realRoots, t); |
| 164 | return foundRoots; |
| 165 | #endif |
| 166 | } |
| 167 | |
| 168 | int cubicRootsReal(double A, double B, double C, double D, double s[3]) { |
mtklein | 753b870 | 2014-08-20 07:38:46 -0700 | [diff] [blame] | 169 | #ifdef SK_DEBUG |
caryclark@google.com | 9f60291 | 2013-01-24 21:47:16 +0000 | [diff] [blame] | 170 | // create a string mathematica understands |
| 171 | // GDB set print repe 15 # if repeated digits is a bother |
| 172 | // set print elements 400 # if line doesn't fit |
| 173 | char str[1024]; |
| 174 | bzero(str, sizeof(str)); |
| 175 | sprintf(str, "Solve[%1.19g x^3 + %1.19g x^2 + %1.19g x + %1.19g == 0, x]", A, B, C, D); |
caryclark@google.com | 5e0500f | 2013-02-20 12:51:37 +0000 | [diff] [blame] | 176 | mathematica_ize(str, sizeof(str)); |
caryclark@google.com | 4aaaaea | 2013-02-28 16:12:39 +0000 | [diff] [blame] | 177 | #if ONE_OFF_DEBUG && ONE_OFF_DEBUG_MATHEMATICA |
caryclark@google.com | 5e0500f | 2013-02-20 12:51:37 +0000 | [diff] [blame] | 178 | SkDebugf("%s\n", str); |
| 179 | #endif |
caryclark@google.com | 9f60291 | 2013-01-24 21:47:16 +0000 | [diff] [blame] | 180 | #endif |
caryclark@google.com | 45a8fc6 | 2013-02-14 15:29:11 +0000 | [diff] [blame] | 181 | if (approximately_zero(A) |
| 182 | && approximately_zero_when_compared_to(A, B) |
caryclark@google.com | beda389 | 2013-02-07 13:13:41 +0000 | [diff] [blame] | 183 | && approximately_zero_when_compared_to(A, C) |
| 184 | && approximately_zero_when_compared_to(A, D)) { // we're just a quadratic |
caryclark@google.com | 9f60291 | 2013-01-24 21:47:16 +0000 | [diff] [blame] | 185 | return quadraticRootsReal(B, C, D, s); |
| 186 | } |
caryclark@google.com | f9502d7 | 2013-02-04 14:06:49 +0000 | [diff] [blame] | 187 | if (approximately_zero_when_compared_to(D, A) |
| 188 | && approximately_zero_when_compared_to(D, B) |
| 189 | && approximately_zero_when_compared_to(D, C)) { // 0 is one root |
caryclark@google.com | 9f60291 | 2013-01-24 21:47:16 +0000 | [diff] [blame] | 190 | int num = quadraticRootsReal(A, B, C, s); |
| 191 | for (int i = 0; i < num; ++i) { |
| 192 | if (approximately_zero(s[i])) { |
| 193 | return num; |
| 194 | } |
| 195 | } |
| 196 | s[num++] = 0; |
| 197 | return num; |
| 198 | } |
| 199 | if (approximately_zero(A + B + C + D)) { // 1 is one root |
| 200 | int num = quadraticRootsReal(A, A + B, -D, s); |
| 201 | for (int i = 0; i < num; ++i) { |
| 202 | if (AlmostEqualUlps(s[i], 1)) { |
| 203 | return num; |
| 204 | } |
| 205 | } |
| 206 | s[num++] = 1; |
| 207 | return num; |
| 208 | } |
| 209 | double a, b, c; |
| 210 | { |
| 211 | double invA = 1 / A; |
| 212 | a = B * invA; |
| 213 | b = C * invA; |
| 214 | c = D * invA; |
| 215 | } |
| 216 | double a2 = a * a; |
| 217 | double Q = (a2 - b * 3) / 9; |
| 218 | double R = (2 * a2 * a - 9 * a * b + 27 * c) / 54; |
| 219 | double R2 = R * R; |
| 220 | double Q3 = Q * Q * Q; |
| 221 | double R2MinusQ3 = R2 - Q3; |
| 222 | double adiv3 = a / 3; |
| 223 | double r; |
| 224 | double* roots = s; |
| 225 | #if 0 |
| 226 | if (approximately_zero_squared(R2MinusQ3) && AlmostEqualUlps(R2, Q3)) { |
| 227 | if (approximately_zero_squared(R)) {/* one triple solution */ |
| 228 | *roots++ = -adiv3; |
| 229 | } else { /* one single and one double solution */ |
| 230 | |
| 231 | double u = cube_root(-R); |
| 232 | *roots++ = 2 * u - adiv3; |
| 233 | *roots++ = -u - adiv3; |
| 234 | } |
| 235 | } |
skia.committer@gmail.com | 4024f32 | 2013-01-25 07:06:46 +0000 | [diff] [blame] | 236 | else |
caryclark@google.com | 9f60291 | 2013-01-24 21:47:16 +0000 | [diff] [blame] | 237 | #endif |
| 238 | if (R2MinusQ3 < 0) // we have 3 real roots |
| 239 | { |
| 240 | double theta = acos(R / sqrt(Q3)); |
| 241 | double neg2RootQ = -2 * sqrt(Q); |
| 242 | |
| 243 | r = neg2RootQ * cos(theta / 3) - adiv3; |
| 244 | *roots++ = r; |
| 245 | |
| 246 | r = neg2RootQ * cos((theta + 2 * PI) / 3) - adiv3; |
| 247 | if (!AlmostEqualUlps(s[0], r)) { |
| 248 | *roots++ = r; |
| 249 | } |
| 250 | r = neg2RootQ * cos((theta - 2 * PI) / 3) - adiv3; |
| 251 | if (!AlmostEqualUlps(s[0], r) && (roots - s == 1 || !AlmostEqualUlps(s[1], r))) { |
| 252 | *roots++ = r; |
| 253 | } |
| 254 | } |
| 255 | else // we have 1 real root |
| 256 | { |
| 257 | double sqrtR2MinusQ3 = sqrt(R2MinusQ3); |
| 258 | double A = fabs(R) + sqrtR2MinusQ3; |
| 259 | A = cube_root(A); |
| 260 | if (R > 0) { |
| 261 | A = -A; |
| 262 | } |
| 263 | if (A != 0) { |
| 264 | A += Q / A; |
| 265 | } |
| 266 | r = A - adiv3; |
| 267 | *roots++ = r; |
| 268 | if (AlmostEqualUlps(R2, Q3)) { |
| 269 | r = -A / 2 - adiv3; |
| 270 | if (!AlmostEqualUlps(s[0], r)) { |
| 271 | *roots++ = r; |
| 272 | } |
| 273 | } |
| 274 | } |
| 275 | return (int)(roots - s); |
caryclark@google.com | c682590 | 2012-02-03 22:07:47 +0000 | [diff] [blame] | 276 | } |
caryclark@google.com | 8dcf114 | 2012-07-02 20:27:02 +0000 | [diff] [blame] | 277 | |
| 278 | // from http://www.cs.sunysb.edu/~qin/courses/geometry/4.pdf |
| 279 | // c(t) = a(1-t)^3 + 3bt(1-t)^2 + 3c(1-t)t^2 + dt^3 |
| 280 | // c'(t) = -3a(1-t)^2 + 3b((1-t)^2 - 2t(1-t)) + 3c(2t(1-t) - t^2) + 3dt^2 |
| 281 | // = 3(b-a)(1-t)^2 + 6(c-b)t(1-t) + 3(d-c)t^2 |
caryclark@google.com | 05c4bad | 2013-01-19 13:22:39 +0000 | [diff] [blame] | 282 | static double derivativeAtT(const double* cubic, double t) { |
caryclark@google.com | 8dcf114 | 2012-07-02 20:27:02 +0000 | [diff] [blame] | 283 | double one_t = 1 - t; |
| 284 | double a = cubic[0]; |
| 285 | double b = cubic[2]; |
| 286 | double c = cubic[4]; |
| 287 | double d = cubic[6]; |
caryclark@google.com | 05c4bad | 2013-01-19 13:22:39 +0000 | [diff] [blame] | 288 | return 3 * ((b - a) * one_t * one_t + 2 * (c - b) * t * one_t + (d - c) * t * t); |
caryclark@google.com | 8dcf114 | 2012-07-02 20:27:02 +0000 | [diff] [blame] | 289 | } |
| 290 | |
caryclark@google.com | 05c4bad | 2013-01-19 13:22:39 +0000 | [diff] [blame] | 291 | double dx_at_t(const Cubic& cubic, double t) { |
| 292 | return derivativeAtT(&cubic[0].x, t); |
caryclark@google.com | 8dcf114 | 2012-07-02 20:27:02 +0000 | [diff] [blame] | 293 | } |
| 294 | |
caryclark@google.com | 05c4bad | 2013-01-19 13:22:39 +0000 | [diff] [blame] | 295 | double dy_at_t(const Cubic& cubic, double t) { |
| 296 | return derivativeAtT(&cubic[0].y, t); |
| 297 | } |
| 298 | |
| 299 | // OPTIMIZE? compute t^2, t(1-t), and (1-t)^2 and pass them to another version of derivative at t? |
caryclark@google.com | 7ff5c84 | 2013-02-26 15:56:05 +0000 | [diff] [blame] | 300 | _Vector dxdy_at_t(const Cubic& cubic, double t) { |
| 301 | _Vector result = { derivativeAtT(&cubic[0].x, t), derivativeAtT(&cubic[0].y, t) }; |
caryclark@google.com | 45a8fc6 | 2013-02-14 15:29:11 +0000 | [diff] [blame] | 302 | return result; |
| 303 | } |
| 304 | |
caryclark@google.com | 1304bb2 | 2013-03-13 20:29:41 +0000 | [diff] [blame] | 305 | // OPTIMIZE? share code with formulate_F1DotF2 |
caryclark@google.com | 73ca624 | 2013-01-17 21:02:47 +0000 | [diff] [blame] | 306 | int find_cubic_inflections(const Cubic& src, double tValues[]) |
| 307 | { |
| 308 | double Ax = src[1].x - src[0].x; |
| 309 | double Ay = src[1].y - src[0].y; |
| 310 | double Bx = src[2].x - 2 * src[1].x + src[0].x; |
| 311 | double By = src[2].y - 2 * src[1].y + src[0].y; |
| 312 | double Cx = src[3].x + 3 * (src[1].x - src[2].x) - src[0].x; |
| 313 | double Cy = src[3].y + 3 * (src[1].y - src[2].y) - src[0].y; |
caryclark@google.com | 1304bb2 | 2013-03-13 20:29:41 +0000 | [diff] [blame] | 314 | return quadraticRootsValidT(Bx * Cy - By * Cx, Ax * Cy - Ay * Cx, Ax * By - Ay * Bx, tValues); |
caryclark@google.com | 73ca624 | 2013-01-17 21:02:47 +0000 | [diff] [blame] | 315 | } |
| 316 | |
caryclark@google.com | 1304bb2 | 2013-03-13 20:29:41 +0000 | [diff] [blame] | 317 | static void formulate_F1DotF2(const double src[], double coeff[4]) |
| 318 | { |
| 319 | double a = src[2] - src[0]; |
| 320 | double b = src[4] - 2 * src[2] + src[0]; |
| 321 | double c = src[6] + 3 * (src[2] - src[4]) - src[0]; |
| 322 | coeff[0] = c * c; |
| 323 | coeff[1] = 3 * b * c; |
| 324 | coeff[2] = 2 * b * b + c * a; |
| 325 | coeff[3] = a * b; |
| 326 | } |
| 327 | |
| 328 | /* from SkGeometry.cpp |
| 329 | Looking for F' dot F'' == 0 |
| 330 | |
| 331 | A = b - a |
| 332 | B = c - 2b + a |
| 333 | C = d - 3c + 3b - a |
| 334 | |
| 335 | F' = 3Ct^2 + 6Bt + 3A |
| 336 | F'' = 6Ct + 6B |
| 337 | |
| 338 | F' dot F'' -> CCt^3 + 3BCt^2 + (2BB + CA)t + AB |
| 339 | */ |
| 340 | int find_cubic_max_curvature(const Cubic& src, double tValues[]) |
| 341 | { |
| 342 | double coeffX[4], coeffY[4]; |
| 343 | int i; |
| 344 | formulate_F1DotF2(&src[0].x, coeffX); |
| 345 | formulate_F1DotF2(&src[0].y, coeffY); |
| 346 | for (i = 0; i < 4; i++) { |
| 347 | coeffX[i] = coeffX[i] + coeffY[i]; |
| 348 | } |
| 349 | return cubicRootsValidT(coeffX[0], coeffX[1], coeffX[2], coeffX[3], tValues); |
| 350 | } |
| 351 | |
| 352 | |
caryclark@google.com | 8dcf114 | 2012-07-02 20:27:02 +0000 | [diff] [blame] | 353 | bool rotate(const Cubic& cubic, int zero, int index, Cubic& rotPath) { |
| 354 | double dy = cubic[index].y - cubic[zero].y; |
| 355 | double dx = cubic[index].x - cubic[zero].x; |
caryclark@google.com | 9f60291 | 2013-01-24 21:47:16 +0000 | [diff] [blame] | 356 | if (approximately_zero(dy)) { |
| 357 | if (approximately_zero(dx)) { |
caryclark@google.com | 8dcf114 | 2012-07-02 20:27:02 +0000 | [diff] [blame] | 358 | return false; |
| 359 | } |
| 360 | memcpy(rotPath, cubic, sizeof(Cubic)); |
| 361 | return true; |
| 362 | } |
| 363 | for (int index = 0; index < 4; ++index) { |
| 364 | rotPath[index].x = cubic[index].x * dx + cubic[index].y * dy; |
| 365 | rotPath[index].y = cubic[index].y * dx - cubic[index].x * dy; |
| 366 | } |
| 367 | return true; |
| 368 | } |
| 369 | |
caryclark@google.com | 05c4bad | 2013-01-19 13:22:39 +0000 | [diff] [blame] | 370 | #if 0 // unused for now |
caryclark@google.com | 8dcf114 | 2012-07-02 20:27:02 +0000 | [diff] [blame] | 371 | double secondDerivativeAtT(const double* cubic, double t) { |
| 372 | double a = cubic[0]; |
| 373 | double b = cubic[2]; |
| 374 | double c = cubic[4]; |
| 375 | double d = cubic[6]; |
| 376 | return (c - 2 * b + a) * (1 - t) + (d - 2 * c + b) * t; |
| 377 | } |
caryclark@google.com | 05c4bad | 2013-01-19 13:22:39 +0000 | [diff] [blame] | 378 | #endif |
caryclark@google.com | 8dcf114 | 2012-07-02 20:27:02 +0000 | [diff] [blame] | 379 | |
caryclark@google.com | 45a8fc6 | 2013-02-14 15:29:11 +0000 | [diff] [blame] | 380 | _Point top(const Cubic& cubic, double startT, double endT) { |
| 381 | Cubic sub; |
| 382 | sub_divide(cubic, startT, endT, sub); |
| 383 | _Point topPt = sub[0]; |
| 384 | if (topPt.y > sub[3].y || (topPt.y == sub[3].y && topPt.x > sub[3].x)) { |
| 385 | topPt = sub[3]; |
| 386 | } |
| 387 | double extremeTs[2]; |
caryclark@google.com | 1304bb2 | 2013-03-13 20:29:41 +0000 | [diff] [blame] | 388 | if (!monotonic_in_y(sub)) { |
caryclark@google.com | 45a8fc6 | 2013-02-14 15:29:11 +0000 | [diff] [blame] | 389 | int roots = findExtrema(sub[0].y, sub[1].y, sub[2].y, sub[3].y, extremeTs); |
| 390 | for (int index = 0; index < roots; ++index) { |
| 391 | _Point mid; |
| 392 | double t = startT + (endT - startT) * extremeTs[index]; |
| 393 | xy_at_t(cubic, t, mid.x, mid.y); |
| 394 | if (topPt.y > mid.y || (topPt.y == mid.y && topPt.x > mid.x)) { |
| 395 | topPt = mid; |
| 396 | } |
| 397 | } |
| 398 | } |
| 399 | return topPt; |
| 400 | } |
| 401 | |
caryclark@google.com | 7ff5c84 | 2013-02-26 15:56:05 +0000 | [diff] [blame] | 402 | // OPTIMIZE: avoid computing the unused half |
| 403 | void xy_at_t(const Cubic& cubic, double t, double& x, double& y) { |
| 404 | _Point xy = xy_at_t(cubic, t); |
| 405 | if (&x) { |
| 406 | x = xy.x; |
| 407 | } |
| 408 | if (&y) { |
| 409 | y = xy.y; |
| 410 | } |
| 411 | } |
| 412 | |
caryclark@google.com | 45a8fc6 | 2013-02-14 15:29:11 +0000 | [diff] [blame] | 413 | _Point xy_at_t(const Cubic& cubic, double t) { |
caryclark@google.com | 8dcf114 | 2012-07-02 20:27:02 +0000 | [diff] [blame] | 414 | double one_t = 1 - t; |
| 415 | double one_t2 = one_t * one_t; |
| 416 | double a = one_t2 * one_t; |
| 417 | double b = 3 * one_t2 * t; |
| 418 | double t2 = t * t; |
| 419 | double c = 3 * one_t * t2; |
| 420 | double d = t2 * t; |
caryclark@google.com | 45a8fc6 | 2013-02-14 15:29:11 +0000 | [diff] [blame] | 421 | _Point result = {a * cubic[0].x + b * cubic[1].x + c * cubic[2].x + d * cubic[3].x, |
| 422 | a * cubic[0].y + b * cubic[1].y + c * cubic[2].y + d * cubic[3].y}; |
| 423 | return result; |
| 424 | } |