epoger@google.com | ec3ed6a | 2011-07-28 14:26:00 +0000 | [diff] [blame^] | 1 | |
| 2 | /* |
| 3 | * Copyright 2006 The Android Open Source Project |
| 4 | * |
| 5 | * Use of this source code is governed by a BSD-style license that can be |
| 6 | * found in the LICENSE file. |
| 7 | */ |
| 8 | |
reed@android.com | 8a1c16f | 2008-12-17 15:59:43 +0000 | [diff] [blame] | 9 | |
| 10 | #include "SkGeometry.h" |
| 11 | #include "Sk64.h" |
| 12 | #include "SkMatrix.h" |
| 13 | |
kbr@chromium.org | 2e08619 | 2010-07-07 22:20:35 +0000 | [diff] [blame] | 14 | bool SkXRayCrossesLine(const SkXRay& pt, const SkPoint pts[2], bool* ambiguous) { |
| 15 | if (ambiguous) { |
| 16 | *ambiguous = false; |
| 17 | } |
reed@android.com | 945a139 | 2010-02-05 20:41:02 +0000 | [diff] [blame] | 18 | // Determine quick discards. |
| 19 | // Consider query line going exactly through point 0 to not |
| 20 | // intersect, for symmetry with SkXRayCrossesMonotonicCubic. |
kbr@chromium.org | 2e08619 | 2010-07-07 22:20:35 +0000 | [diff] [blame] | 21 | if (pt.fY == pts[0].fY) { |
| 22 | if (ambiguous) { |
| 23 | *ambiguous = true; |
| 24 | } |
reed@android.com | 945a139 | 2010-02-05 20:41:02 +0000 | [diff] [blame] | 25 | return false; |
kbr@chromium.org | 2e08619 | 2010-07-07 22:20:35 +0000 | [diff] [blame] | 26 | } |
reed@android.com | 945a139 | 2010-02-05 20:41:02 +0000 | [diff] [blame] | 27 | if (pt.fY < pts[0].fY && pt.fY < pts[1].fY) |
| 28 | return false; |
| 29 | if (pt.fY > pts[0].fY && pt.fY > pts[1].fY) |
| 30 | return false; |
| 31 | if (pt.fX > pts[0].fX && pt.fX > pts[1].fX) |
| 32 | return false; |
| 33 | // Determine degenerate cases |
| 34 | if (SkScalarNearlyZero(pts[0].fY - pts[1].fY)) |
| 35 | return false; |
kbr@chromium.org | 2e08619 | 2010-07-07 22:20:35 +0000 | [diff] [blame] | 36 | if (SkScalarNearlyZero(pts[0].fX - pts[1].fX)) { |
reed@android.com | 945a139 | 2010-02-05 20:41:02 +0000 | [diff] [blame] | 37 | // We've already determined the query point lies within the |
| 38 | // vertical range of the line segment. |
kbr@chromium.org | 2e08619 | 2010-07-07 22:20:35 +0000 | [diff] [blame] | 39 | if (pt.fX <= pts[0].fX) { |
| 40 | if (ambiguous) { |
| 41 | *ambiguous = (pt.fY == pts[1].fY); |
| 42 | } |
| 43 | return true; |
| 44 | } |
| 45 | return false; |
| 46 | } |
| 47 | // Ambiguity check |
| 48 | if (pt.fY == pts[1].fY) { |
| 49 | if (pt.fX <= pts[1].fX) { |
| 50 | if (ambiguous) { |
| 51 | *ambiguous = true; |
| 52 | } |
| 53 | return true; |
| 54 | } |
| 55 | return false; |
| 56 | } |
reed@android.com | 945a139 | 2010-02-05 20:41:02 +0000 | [diff] [blame] | 57 | // Full line segment evaluation |
| 58 | SkScalar delta_y = pts[1].fY - pts[0].fY; |
| 59 | SkScalar delta_x = pts[1].fX - pts[0].fX; |
| 60 | SkScalar slope = SkScalarDiv(delta_y, delta_x); |
| 61 | SkScalar b = pts[0].fY - SkScalarMul(slope, pts[0].fX); |
| 62 | // Solve for x coordinate at y = pt.fY |
| 63 | SkScalar x = SkScalarDiv(pt.fY - b, slope); |
| 64 | return pt.fX <= x; |
| 65 | } |
| 66 | |
reed@android.com | 8a1c16f | 2008-12-17 15:59:43 +0000 | [diff] [blame] | 67 | /** If defined, this makes eval_quad and eval_cubic do more setup (sometimes |
| 68 | involving integer multiplies by 2 or 3, but fewer calls to SkScalarMul. |
| 69 | May also introduce overflow of fixed when we compute our setup. |
| 70 | */ |
| 71 | #ifdef SK_SCALAR_IS_FIXED |
| 72 | #define DIRECT_EVAL_OF_POLYNOMIALS |
| 73 | #endif |
| 74 | |
| 75 | //////////////////////////////////////////////////////////////////////// |
| 76 | |
| 77 | #ifdef SK_SCALAR_IS_FIXED |
| 78 | static int is_not_monotonic(int a, int b, int c, int d) |
| 79 | { |
| 80 | return (((a - b) | (b - c) | (c - d)) & ((b - a) | (c - b) | (d - c))) >> 31; |
| 81 | } |
| 82 | |
| 83 | static int is_not_monotonic(int a, int b, int c) |
| 84 | { |
| 85 | return (((a - b) | (b - c)) & ((b - a) | (c - b))) >> 31; |
| 86 | } |
| 87 | #else |
| 88 | static int is_not_monotonic(float a, float b, float c) |
| 89 | { |
| 90 | float ab = a - b; |
| 91 | float bc = b - c; |
| 92 | if (ab < 0) |
| 93 | bc = -bc; |
| 94 | return ab == 0 || bc < 0; |
| 95 | } |
| 96 | #endif |
| 97 | |
| 98 | //////////////////////////////////////////////////////////////////////// |
| 99 | |
| 100 | static bool is_unit_interval(SkScalar x) |
| 101 | { |
| 102 | return x > 0 && x < SK_Scalar1; |
| 103 | } |
| 104 | |
| 105 | static int valid_unit_divide(SkScalar numer, SkScalar denom, SkScalar* ratio) |
| 106 | { |
| 107 | SkASSERT(ratio); |
| 108 | |
| 109 | if (numer < 0) |
| 110 | { |
| 111 | numer = -numer; |
| 112 | denom = -denom; |
| 113 | } |
| 114 | |
| 115 | if (denom == 0 || numer == 0 || numer >= denom) |
| 116 | return 0; |
| 117 | |
| 118 | SkScalar r = SkScalarDiv(numer, denom); |
reed@android.com | 1516162 | 2010-03-08 17:44:42 +0000 | [diff] [blame] | 119 | if (SkScalarIsNaN(r)) { |
| 120 | return 0; |
| 121 | } |
reed@android.com | 8a1c16f | 2008-12-17 15:59:43 +0000 | [diff] [blame] | 122 | SkASSERT(r >= 0 && r < SK_Scalar1); |
| 123 | if (r == 0) // catch underflow if numer <<<< denom |
| 124 | return 0; |
| 125 | *ratio = r; |
| 126 | return 1; |
| 127 | } |
| 128 | |
| 129 | /** From Numerical Recipes in C. |
| 130 | |
| 131 | Q = -1/2 (B + sign(B) sqrt[B*B - 4*A*C]) |
| 132 | x1 = Q / A |
| 133 | x2 = C / Q |
| 134 | */ |
| 135 | int SkFindUnitQuadRoots(SkScalar A, SkScalar B, SkScalar C, SkScalar roots[2]) |
| 136 | { |
| 137 | SkASSERT(roots); |
| 138 | |
| 139 | if (A == 0) |
| 140 | return valid_unit_divide(-C, B, roots); |
| 141 | |
| 142 | SkScalar* r = roots; |
| 143 | |
| 144 | #ifdef SK_SCALAR_IS_FLOAT |
| 145 | float R = B*B - 4*A*C; |
reed@android.com | 1516162 | 2010-03-08 17:44:42 +0000 | [diff] [blame] | 146 | if (R < 0 || SkScalarIsNaN(R)) { // complex roots |
reed@android.com | 8a1c16f | 2008-12-17 15:59:43 +0000 | [diff] [blame] | 147 | return 0; |
reed@android.com | 1516162 | 2010-03-08 17:44:42 +0000 | [diff] [blame] | 148 | } |
reed@android.com | 8a1c16f | 2008-12-17 15:59:43 +0000 | [diff] [blame] | 149 | R = sk_float_sqrt(R); |
| 150 | #else |
| 151 | Sk64 RR, tmp; |
| 152 | |
| 153 | RR.setMul(B,B); |
| 154 | tmp.setMul(A,C); |
| 155 | tmp.shiftLeft(2); |
| 156 | RR.sub(tmp); |
| 157 | if (RR.isNeg()) |
| 158 | return 0; |
| 159 | SkFixed R = RR.getSqrt(); |
| 160 | #endif |
| 161 | |
| 162 | SkScalar Q = (B < 0) ? -(B-R)/2 : -(B+R)/2; |
| 163 | r += valid_unit_divide(Q, A, r); |
| 164 | r += valid_unit_divide(C, Q, r); |
| 165 | if (r - roots == 2) |
| 166 | { |
| 167 | if (roots[0] > roots[1]) |
| 168 | SkTSwap<SkScalar>(roots[0], roots[1]); |
| 169 | else if (roots[0] == roots[1]) // nearly-equal? |
| 170 | r -= 1; // skip the double root |
| 171 | } |
| 172 | return (int)(r - roots); |
| 173 | } |
| 174 | |
| 175 | #ifdef SK_SCALAR_IS_FIXED |
| 176 | /** Trim A/B/C down so that they are all <= 32bits |
| 177 | and then call SkFindUnitQuadRoots() |
| 178 | */ |
| 179 | static int Sk64FindFixedQuadRoots(const Sk64& A, const Sk64& B, const Sk64& C, SkFixed roots[2]) |
| 180 | { |
| 181 | int na = A.shiftToMake32(); |
| 182 | int nb = B.shiftToMake32(); |
| 183 | int nc = C.shiftToMake32(); |
| 184 | |
| 185 | int shift = SkMax32(na, SkMax32(nb, nc)); |
| 186 | SkASSERT(shift >= 0); |
| 187 | |
| 188 | return SkFindUnitQuadRoots(A.getShiftRight(shift), B.getShiftRight(shift), C.getShiftRight(shift), roots); |
| 189 | } |
| 190 | #endif |
| 191 | |
| 192 | ///////////////////////////////////////////////////////////////////////////////////// |
| 193 | ///////////////////////////////////////////////////////////////////////////////////// |
| 194 | |
| 195 | static SkScalar eval_quad(const SkScalar src[], SkScalar t) |
| 196 | { |
| 197 | SkASSERT(src); |
| 198 | SkASSERT(t >= 0 && t <= SK_Scalar1); |
| 199 | |
| 200 | #ifdef DIRECT_EVAL_OF_POLYNOMIALS |
| 201 | SkScalar C = src[0]; |
| 202 | SkScalar A = src[4] - 2 * src[2] + C; |
| 203 | SkScalar B = 2 * (src[2] - C); |
| 204 | return SkScalarMulAdd(SkScalarMulAdd(A, t, B), t, C); |
| 205 | #else |
| 206 | SkScalar ab = SkScalarInterp(src[0], src[2], t); |
| 207 | SkScalar bc = SkScalarInterp(src[2], src[4], t); |
| 208 | return SkScalarInterp(ab, bc, t); |
| 209 | #endif |
| 210 | } |
| 211 | |
| 212 | static SkScalar eval_quad_derivative(const SkScalar src[], SkScalar t) |
| 213 | { |
| 214 | SkScalar A = src[4] - 2 * src[2] + src[0]; |
| 215 | SkScalar B = src[2] - src[0]; |
| 216 | |
| 217 | return 2 * SkScalarMulAdd(A, t, B); |
| 218 | } |
| 219 | |
| 220 | static SkScalar eval_quad_derivative_at_half(const SkScalar src[]) |
| 221 | { |
| 222 | SkScalar A = src[4] - 2 * src[2] + src[0]; |
| 223 | SkScalar B = src[2] - src[0]; |
| 224 | return A + 2 * B; |
| 225 | } |
| 226 | |
| 227 | void SkEvalQuadAt(const SkPoint src[3], SkScalar t, SkPoint* pt, SkVector* tangent) |
| 228 | { |
| 229 | SkASSERT(src); |
| 230 | SkASSERT(t >= 0 && t <= SK_Scalar1); |
| 231 | |
| 232 | if (pt) |
| 233 | pt->set(eval_quad(&src[0].fX, t), eval_quad(&src[0].fY, t)); |
| 234 | if (tangent) |
| 235 | tangent->set(eval_quad_derivative(&src[0].fX, t), |
| 236 | eval_quad_derivative(&src[0].fY, t)); |
| 237 | } |
| 238 | |
| 239 | void SkEvalQuadAtHalf(const SkPoint src[3], SkPoint* pt, SkVector* tangent) |
| 240 | { |
| 241 | SkASSERT(src); |
| 242 | |
| 243 | if (pt) |
| 244 | { |
| 245 | SkScalar x01 = SkScalarAve(src[0].fX, src[1].fX); |
| 246 | SkScalar y01 = SkScalarAve(src[0].fY, src[1].fY); |
| 247 | SkScalar x12 = SkScalarAve(src[1].fX, src[2].fX); |
| 248 | SkScalar y12 = SkScalarAve(src[1].fY, src[2].fY); |
| 249 | pt->set(SkScalarAve(x01, x12), SkScalarAve(y01, y12)); |
| 250 | } |
| 251 | if (tangent) |
| 252 | tangent->set(eval_quad_derivative_at_half(&src[0].fX), |
| 253 | eval_quad_derivative_at_half(&src[0].fY)); |
| 254 | } |
| 255 | |
| 256 | static void interp_quad_coords(const SkScalar* src, SkScalar* dst, SkScalar t) |
| 257 | { |
| 258 | SkScalar ab = SkScalarInterp(src[0], src[2], t); |
| 259 | SkScalar bc = SkScalarInterp(src[2], src[4], t); |
| 260 | |
| 261 | dst[0] = src[0]; |
| 262 | dst[2] = ab; |
| 263 | dst[4] = SkScalarInterp(ab, bc, t); |
| 264 | dst[6] = bc; |
| 265 | dst[8] = src[4]; |
| 266 | } |
| 267 | |
| 268 | void SkChopQuadAt(const SkPoint src[3], SkPoint dst[5], SkScalar t) |
| 269 | { |
| 270 | SkASSERT(t > 0 && t < SK_Scalar1); |
| 271 | |
| 272 | interp_quad_coords(&src[0].fX, &dst[0].fX, t); |
| 273 | interp_quad_coords(&src[0].fY, &dst[0].fY, t); |
| 274 | } |
| 275 | |
| 276 | void SkChopQuadAtHalf(const SkPoint src[3], SkPoint dst[5]) |
| 277 | { |
| 278 | SkScalar x01 = SkScalarAve(src[0].fX, src[1].fX); |
| 279 | SkScalar y01 = SkScalarAve(src[0].fY, src[1].fY); |
| 280 | SkScalar x12 = SkScalarAve(src[1].fX, src[2].fX); |
| 281 | SkScalar y12 = SkScalarAve(src[1].fY, src[2].fY); |
| 282 | |
| 283 | dst[0] = src[0]; |
| 284 | dst[1].set(x01, y01); |
| 285 | dst[2].set(SkScalarAve(x01, x12), SkScalarAve(y01, y12)); |
| 286 | dst[3].set(x12, y12); |
| 287 | dst[4] = src[2]; |
| 288 | } |
| 289 | |
| 290 | /** Quad'(t) = At + B, where |
| 291 | A = 2(a - 2b + c) |
| 292 | B = 2(b - a) |
| 293 | Solve for t, only if it fits between 0 < t < 1 |
| 294 | */ |
| 295 | int SkFindQuadExtrema(SkScalar a, SkScalar b, SkScalar c, SkScalar tValue[1]) |
| 296 | { |
| 297 | /* At + B == 0 |
| 298 | t = -B / A |
| 299 | */ |
| 300 | #ifdef SK_SCALAR_IS_FIXED |
| 301 | return is_not_monotonic(a, b, c) && valid_unit_divide(a - b, a - b - b + c, tValue); |
| 302 | #else |
| 303 | return valid_unit_divide(a - b, a - b - b + c, tValue); |
| 304 | #endif |
| 305 | } |
| 306 | |
reed@android.com | fc25abd | 2009-01-15 14:38:33 +0000 | [diff] [blame] | 307 | static inline void flatten_double_quad_extrema(SkScalar coords[14]) |
reed@android.com | 8a1c16f | 2008-12-17 15:59:43 +0000 | [diff] [blame] | 308 | { |
| 309 | coords[2] = coords[6] = coords[4]; |
| 310 | } |
| 311 | |
reed@android.com | 8a1c16f | 2008-12-17 15:59:43 +0000 | [diff] [blame] | 312 | /* Returns 0 for 1 quad, and 1 for two quads, either way the answer is |
reed@android.com | 77f0ef7 | 2009-11-17 18:47:52 +0000 | [diff] [blame] | 313 | stored in dst[]. Guarantees that the 1/2 quads will be monotonic. |
| 314 | */ |
reed@android.com | 8a1c16f | 2008-12-17 15:59:43 +0000 | [diff] [blame] | 315 | int SkChopQuadAtYExtrema(const SkPoint src[3], SkPoint dst[5]) |
| 316 | { |
| 317 | SkASSERT(src); |
| 318 | SkASSERT(dst); |
reed@android.com | 77f0ef7 | 2009-11-17 18:47:52 +0000 | [diff] [blame] | 319 | |
reed@android.com | 8a1c16f | 2008-12-17 15:59:43 +0000 | [diff] [blame] | 320 | #if 0 |
| 321 | static bool once = true; |
| 322 | if (once) |
| 323 | { |
| 324 | once = false; |
| 325 | SkPoint s[3] = { 0, 26398, 0, 26331, 0, 20621428 }; |
| 326 | SkPoint d[6]; |
| 327 | |
| 328 | int n = SkChopQuadAtYExtrema(s, d); |
| 329 | SkDebugf("chop=%d, Y=[%x %x %x %x %x %x]\n", n, d[0].fY, d[1].fY, d[2].fY, d[3].fY, d[4].fY, d[5].fY); |
| 330 | } |
| 331 | #endif |
reed@android.com | 77f0ef7 | 2009-11-17 18:47:52 +0000 | [diff] [blame] | 332 | |
reed@android.com | 8a1c16f | 2008-12-17 15:59:43 +0000 | [diff] [blame] | 333 | SkScalar a = src[0].fY; |
| 334 | SkScalar b = src[1].fY; |
| 335 | SkScalar c = src[2].fY; |
reed@android.com | 77f0ef7 | 2009-11-17 18:47:52 +0000 | [diff] [blame] | 336 | |
reed@android.com | 8a1c16f | 2008-12-17 15:59:43 +0000 | [diff] [blame] | 337 | if (is_not_monotonic(a, b, c)) |
| 338 | { |
| 339 | SkScalar tValue; |
| 340 | if (valid_unit_divide(a - b, a - b - b + c, &tValue)) |
| 341 | { |
| 342 | SkChopQuadAt(src, dst, tValue); |
| 343 | flatten_double_quad_extrema(&dst[0].fY); |
| 344 | return 1; |
| 345 | } |
| 346 | // if we get here, we need to force dst to be monotonic, even though |
| 347 | // we couldn't compute a unit_divide value (probably underflow). |
| 348 | b = SkScalarAbs(a - b) < SkScalarAbs(b - c) ? a : c; |
| 349 | } |
| 350 | dst[0].set(src[0].fX, a); |
| 351 | dst[1].set(src[1].fX, b); |
| 352 | dst[2].set(src[2].fX, c); |
| 353 | return 0; |
| 354 | } |
| 355 | |
reed@android.com | 77f0ef7 | 2009-11-17 18:47:52 +0000 | [diff] [blame] | 356 | /* Returns 0 for 1 quad, and 1 for two quads, either way the answer is |
| 357 | stored in dst[]. Guarantees that the 1/2 quads will be monotonic. |
| 358 | */ |
| 359 | int SkChopQuadAtXExtrema(const SkPoint src[3], SkPoint dst[5]) |
| 360 | { |
| 361 | SkASSERT(src); |
| 362 | SkASSERT(dst); |
| 363 | |
| 364 | SkScalar a = src[0].fX; |
| 365 | SkScalar b = src[1].fX; |
| 366 | SkScalar c = src[2].fX; |
| 367 | |
| 368 | if (is_not_monotonic(a, b, c)) { |
| 369 | SkScalar tValue; |
| 370 | if (valid_unit_divide(a - b, a - b - b + c, &tValue)) { |
| 371 | SkChopQuadAt(src, dst, tValue); |
| 372 | flatten_double_quad_extrema(&dst[0].fX); |
| 373 | return 1; |
| 374 | } |
| 375 | // if we get here, we need to force dst to be monotonic, even though |
| 376 | // we couldn't compute a unit_divide value (probably underflow). |
| 377 | b = SkScalarAbs(a - b) < SkScalarAbs(b - c) ? a : c; |
| 378 | } |
| 379 | dst[0].set(a, src[0].fY); |
| 380 | dst[1].set(b, src[1].fY); |
| 381 | dst[2].set(c, src[2].fY); |
| 382 | return 0; |
| 383 | } |
| 384 | |
reed@android.com | 8a1c16f | 2008-12-17 15:59:43 +0000 | [diff] [blame] | 385 | // F(t) = a (1 - t) ^ 2 + 2 b t (1 - t) + c t ^ 2 |
| 386 | // F'(t) = 2 (b - a) + 2 (a - 2b + c) t |
| 387 | // F''(t) = 2 (a - 2b + c) |
| 388 | // |
| 389 | // A = 2 (b - a) |
| 390 | // B = 2 (a - 2b + c) |
| 391 | // |
| 392 | // Maximum curvature for a quadratic means solving |
| 393 | // Fx' Fx'' + Fy' Fy'' = 0 |
| 394 | // |
| 395 | // t = - (Ax Bx + Ay By) / (Bx ^ 2 + By ^ 2) |
| 396 | // |
| 397 | int SkChopQuadAtMaxCurvature(const SkPoint src[3], SkPoint dst[5]) |
| 398 | { |
| 399 | SkScalar Ax = src[1].fX - src[0].fX; |
| 400 | SkScalar Ay = src[1].fY - src[0].fY; |
| 401 | SkScalar Bx = src[0].fX - src[1].fX - src[1].fX + src[2].fX; |
| 402 | SkScalar By = src[0].fY - src[1].fY - src[1].fY + src[2].fY; |
| 403 | SkScalar t = 0; // 0 means don't chop |
| 404 | |
| 405 | #ifdef SK_SCALAR_IS_FLOAT |
| 406 | (void)valid_unit_divide(-(Ax * Bx + Ay * By), Bx * Bx + By * By, &t); |
| 407 | #else |
| 408 | // !!! should I use SkFloat here? seems like it |
| 409 | Sk64 numer, denom, tmp; |
| 410 | |
| 411 | numer.setMul(Ax, -Bx); |
| 412 | tmp.setMul(Ay, -By); |
| 413 | numer.add(tmp); |
| 414 | |
| 415 | if (numer.isPos()) // do nothing if numer <= 0 |
| 416 | { |
| 417 | denom.setMul(Bx, Bx); |
| 418 | tmp.setMul(By, By); |
| 419 | denom.add(tmp); |
| 420 | SkASSERT(!denom.isNeg()); |
| 421 | if (numer < denom) |
| 422 | { |
| 423 | t = numer.getFixedDiv(denom); |
| 424 | SkASSERT(t >= 0 && t <= SK_Fixed1); // assert that we're numerically stable (ha!) |
| 425 | if ((unsigned)t >= SK_Fixed1) // runtime check for numerical stability |
| 426 | t = 0; // ignore the chop |
| 427 | } |
| 428 | } |
| 429 | #endif |
| 430 | |
| 431 | if (t == 0) |
| 432 | { |
| 433 | memcpy(dst, src, 3 * sizeof(SkPoint)); |
| 434 | return 1; |
| 435 | } |
| 436 | else |
| 437 | { |
| 438 | SkChopQuadAt(src, dst, t); |
| 439 | return 2; |
| 440 | } |
| 441 | } |
| 442 | |
reed@google.com | 6fc321a | 2011-07-27 13:54:36 +0000 | [diff] [blame] | 443 | #ifdef SK_SCALAR_IS_FLOAT |
| 444 | #define SK_ScalarTwoThirds (0.666666666f) |
| 445 | #else |
| 446 | #define SK_ScalarTwoThirds ((SkFixed)(43691)) |
| 447 | #endif |
| 448 | |
| 449 | void SkConvertQuadToCubic(const SkPoint src[3], SkPoint dst[4]) { |
| 450 | const SkScalar scale = SK_ScalarTwoThirds; |
| 451 | dst[0] = src[0]; |
| 452 | dst[1].set(src[0].fX + SkScalarMul(src[1].fX - src[0].fX, scale), |
| 453 | src[0].fY + SkScalarMul(src[1].fY - src[0].fY, scale)); |
| 454 | dst[2].set(src[2].fX + SkScalarMul(src[1].fX - src[2].fX, scale), |
| 455 | src[2].fY + SkScalarMul(src[1].fY - src[2].fY, scale)); |
| 456 | dst[3] = src[2]; |
reed@android.com | 945a139 | 2010-02-05 20:41:02 +0000 | [diff] [blame] | 457 | } |
| 458 | |
reed@android.com | 8a1c16f | 2008-12-17 15:59:43 +0000 | [diff] [blame] | 459 | //////////////////////////////////////////////////////////////////////////////////////// |
| 460 | ///// CUBICS // CUBICS // CUBICS // CUBICS // CUBICS // CUBICS // CUBICS // CUBICS ///// |
| 461 | //////////////////////////////////////////////////////////////////////////////////////// |
| 462 | |
| 463 | static void get_cubic_coeff(const SkScalar pt[], SkScalar coeff[4]) |
| 464 | { |
| 465 | coeff[0] = pt[6] + 3*(pt[2] - pt[4]) - pt[0]; |
| 466 | coeff[1] = 3*(pt[4] - pt[2] - pt[2] + pt[0]); |
| 467 | coeff[2] = 3*(pt[2] - pt[0]); |
| 468 | coeff[3] = pt[0]; |
| 469 | } |
| 470 | |
| 471 | void SkGetCubicCoeff(const SkPoint pts[4], SkScalar cx[4], SkScalar cy[4]) |
| 472 | { |
| 473 | SkASSERT(pts); |
| 474 | |
| 475 | if (cx) |
| 476 | get_cubic_coeff(&pts[0].fX, cx); |
| 477 | if (cy) |
| 478 | get_cubic_coeff(&pts[0].fY, cy); |
| 479 | } |
| 480 | |
| 481 | static SkScalar eval_cubic(const SkScalar src[], SkScalar t) |
| 482 | { |
| 483 | SkASSERT(src); |
| 484 | SkASSERT(t >= 0 && t <= SK_Scalar1); |
| 485 | |
| 486 | if (t == 0) |
| 487 | return src[0]; |
| 488 | |
| 489 | #ifdef DIRECT_EVAL_OF_POLYNOMIALS |
| 490 | SkScalar D = src[0]; |
| 491 | SkScalar A = src[6] + 3*(src[2] - src[4]) - D; |
| 492 | SkScalar B = 3*(src[4] - src[2] - src[2] + D); |
| 493 | SkScalar C = 3*(src[2] - D); |
| 494 | |
| 495 | return SkScalarMulAdd(SkScalarMulAdd(SkScalarMulAdd(A, t, B), t, C), t, D); |
| 496 | #else |
| 497 | SkScalar ab = SkScalarInterp(src[0], src[2], t); |
| 498 | SkScalar bc = SkScalarInterp(src[2], src[4], t); |
| 499 | SkScalar cd = SkScalarInterp(src[4], src[6], t); |
| 500 | SkScalar abc = SkScalarInterp(ab, bc, t); |
| 501 | SkScalar bcd = SkScalarInterp(bc, cd, t); |
| 502 | return SkScalarInterp(abc, bcd, t); |
| 503 | #endif |
| 504 | } |
| 505 | |
| 506 | /** return At^2 + Bt + C |
| 507 | */ |
| 508 | static SkScalar eval_quadratic(SkScalar A, SkScalar B, SkScalar C, SkScalar t) |
| 509 | { |
| 510 | SkASSERT(t >= 0 && t <= SK_Scalar1); |
| 511 | |
| 512 | return SkScalarMulAdd(SkScalarMulAdd(A, t, B), t, C); |
| 513 | } |
| 514 | |
| 515 | static SkScalar eval_cubic_derivative(const SkScalar src[], SkScalar t) |
| 516 | { |
| 517 | SkScalar A = src[6] + 3*(src[2] - src[4]) - src[0]; |
| 518 | SkScalar B = 2*(src[4] - 2 * src[2] + src[0]); |
| 519 | SkScalar C = src[2] - src[0]; |
| 520 | |
| 521 | return eval_quadratic(A, B, C, t); |
| 522 | } |
| 523 | |
| 524 | static SkScalar eval_cubic_2ndDerivative(const SkScalar src[], SkScalar t) |
| 525 | { |
| 526 | SkScalar A = src[6] + 3*(src[2] - src[4]) - src[0]; |
| 527 | SkScalar B = src[4] - 2 * src[2] + src[0]; |
| 528 | |
| 529 | return SkScalarMulAdd(A, t, B); |
| 530 | } |
| 531 | |
| 532 | void SkEvalCubicAt(const SkPoint src[4], SkScalar t, SkPoint* loc, SkVector* tangent, SkVector* curvature) |
| 533 | { |
| 534 | SkASSERT(src); |
| 535 | SkASSERT(t >= 0 && t <= SK_Scalar1); |
| 536 | |
| 537 | if (loc) |
| 538 | loc->set(eval_cubic(&src[0].fX, t), eval_cubic(&src[0].fY, t)); |
| 539 | if (tangent) |
| 540 | tangent->set(eval_cubic_derivative(&src[0].fX, t), |
| 541 | eval_cubic_derivative(&src[0].fY, t)); |
| 542 | if (curvature) |
| 543 | curvature->set(eval_cubic_2ndDerivative(&src[0].fX, t), |
| 544 | eval_cubic_2ndDerivative(&src[0].fY, t)); |
| 545 | } |
| 546 | |
| 547 | /** Cubic'(t) = At^2 + Bt + C, where |
| 548 | A = 3(-a + 3(b - c) + d) |
| 549 | B = 6(a - 2b + c) |
| 550 | C = 3(b - a) |
| 551 | Solve for t, keeping only those that fit betwee 0 < t < 1 |
| 552 | */ |
| 553 | int SkFindCubicExtrema(SkScalar a, SkScalar b, SkScalar c, SkScalar d, SkScalar tValues[2]) |
| 554 | { |
| 555 | #ifdef SK_SCALAR_IS_FIXED |
| 556 | if (!is_not_monotonic(a, b, c, d)) |
| 557 | return 0; |
| 558 | #endif |
| 559 | |
| 560 | // we divide A,B,C by 3 to simplify |
| 561 | SkScalar A = d - a + 3*(b - c); |
| 562 | SkScalar B = 2*(a - b - b + c); |
| 563 | SkScalar C = b - a; |
| 564 | |
| 565 | return SkFindUnitQuadRoots(A, B, C, tValues); |
| 566 | } |
| 567 | |
| 568 | static void interp_cubic_coords(const SkScalar* src, SkScalar* dst, SkScalar t) |
| 569 | { |
| 570 | SkScalar ab = SkScalarInterp(src[0], src[2], t); |
| 571 | SkScalar bc = SkScalarInterp(src[2], src[4], t); |
| 572 | SkScalar cd = SkScalarInterp(src[4], src[6], t); |
| 573 | SkScalar abc = SkScalarInterp(ab, bc, t); |
| 574 | SkScalar bcd = SkScalarInterp(bc, cd, t); |
| 575 | SkScalar abcd = SkScalarInterp(abc, bcd, t); |
| 576 | |
| 577 | dst[0] = src[0]; |
| 578 | dst[2] = ab; |
| 579 | dst[4] = abc; |
| 580 | dst[6] = abcd; |
| 581 | dst[8] = bcd; |
| 582 | dst[10] = cd; |
| 583 | dst[12] = src[6]; |
| 584 | } |
| 585 | |
| 586 | void SkChopCubicAt(const SkPoint src[4], SkPoint dst[7], SkScalar t) |
| 587 | { |
| 588 | SkASSERT(t > 0 && t < SK_Scalar1); |
| 589 | |
| 590 | interp_cubic_coords(&src[0].fX, &dst[0].fX, t); |
| 591 | interp_cubic_coords(&src[0].fY, &dst[0].fY, t); |
| 592 | } |
| 593 | |
reed@android.com | a964028 | 2009-08-28 20:06:54 +0000 | [diff] [blame] | 594 | /* http://code.google.com/p/skia/issues/detail?id=32 |
| 595 | |
| 596 | This test code would fail when we didn't check the return result of |
| 597 | valid_unit_divide in SkChopCubicAt(... tValues[], int roots). The reason is |
| 598 | that after the first chop, the parameters to valid_unit_divide are equal |
| 599 | (thanks to finite float precision and rounding in the subtracts). Thus |
| 600 | even though the 2nd tValue looks < 1.0, after we renormalize it, we end |
| 601 | up with 1.0, hence the need to check and just return the last cubic as |
| 602 | a degenerate clump of 4 points in the sampe place. |
| 603 | |
| 604 | static void test_cubic() { |
| 605 | SkPoint src[4] = { |
| 606 | { 556.25000, 523.03003 }, |
| 607 | { 556.23999, 522.96002 }, |
| 608 | { 556.21997, 522.89001 }, |
| 609 | { 556.21997, 522.82001 } |
| 610 | }; |
| 611 | SkPoint dst[10]; |
| 612 | SkScalar tval[] = { 0.33333334f, 0.99999994f }; |
| 613 | SkChopCubicAt(src, dst, tval, 2); |
| 614 | } |
| 615 | */ |
| 616 | |
reed@android.com | 8a1c16f | 2008-12-17 15:59:43 +0000 | [diff] [blame] | 617 | void SkChopCubicAt(const SkPoint src[4], SkPoint dst[], const SkScalar tValues[], int roots) |
| 618 | { |
| 619 | #ifdef SK_DEBUG |
| 620 | { |
| 621 | for (int i = 0; i < roots - 1; i++) |
| 622 | { |
| 623 | SkASSERT(is_unit_interval(tValues[i])); |
| 624 | SkASSERT(is_unit_interval(tValues[i+1])); |
| 625 | SkASSERT(tValues[i] < tValues[i+1]); |
| 626 | } |
| 627 | } |
| 628 | #endif |
| 629 | |
| 630 | if (dst) |
| 631 | { |
| 632 | if (roots == 0) // nothing to chop |
| 633 | memcpy(dst, src, 4*sizeof(SkPoint)); |
| 634 | else |
| 635 | { |
| 636 | SkScalar t = tValues[0]; |
| 637 | SkPoint tmp[4]; |
| 638 | |
| 639 | for (int i = 0; i < roots; i++) |
| 640 | { |
| 641 | SkChopCubicAt(src, dst, t); |
| 642 | if (i == roots - 1) |
| 643 | break; |
| 644 | |
reed@android.com | 8a1c16f | 2008-12-17 15:59:43 +0000 | [diff] [blame] | 645 | dst += 3; |
reed@android.com | a964028 | 2009-08-28 20:06:54 +0000 | [diff] [blame] | 646 | // have src point to the remaining cubic (after the chop) |
reed@android.com | 8a1c16f | 2008-12-17 15:59:43 +0000 | [diff] [blame] | 647 | memcpy(tmp, dst, 4 * sizeof(SkPoint)); |
| 648 | src = tmp; |
reed@android.com | a964028 | 2009-08-28 20:06:54 +0000 | [diff] [blame] | 649 | |
| 650 | // watch out in case the renormalized t isn't in range |
| 651 | if (!valid_unit_divide(tValues[i+1] - tValues[i], |
| 652 | SK_Scalar1 - tValues[i], &t)) { |
| 653 | // if we can't, just create a degenerate cubic |
| 654 | dst[4] = dst[5] = dst[6] = src[3]; |
| 655 | break; |
| 656 | } |
reed@android.com | 8a1c16f | 2008-12-17 15:59:43 +0000 | [diff] [blame] | 657 | } |
| 658 | } |
| 659 | } |
| 660 | } |
| 661 | |
| 662 | void SkChopCubicAtHalf(const SkPoint src[4], SkPoint dst[7]) |
| 663 | { |
| 664 | SkScalar x01 = SkScalarAve(src[0].fX, src[1].fX); |
| 665 | SkScalar y01 = SkScalarAve(src[0].fY, src[1].fY); |
| 666 | SkScalar x12 = SkScalarAve(src[1].fX, src[2].fX); |
| 667 | SkScalar y12 = SkScalarAve(src[1].fY, src[2].fY); |
| 668 | SkScalar x23 = SkScalarAve(src[2].fX, src[3].fX); |
| 669 | SkScalar y23 = SkScalarAve(src[2].fY, src[3].fY); |
| 670 | |
| 671 | SkScalar x012 = SkScalarAve(x01, x12); |
| 672 | SkScalar y012 = SkScalarAve(y01, y12); |
| 673 | SkScalar x123 = SkScalarAve(x12, x23); |
| 674 | SkScalar y123 = SkScalarAve(y12, y23); |
| 675 | |
| 676 | dst[0] = src[0]; |
| 677 | dst[1].set(x01, y01); |
| 678 | dst[2].set(x012, y012); |
| 679 | dst[3].set(SkScalarAve(x012, x123), SkScalarAve(y012, y123)); |
| 680 | dst[4].set(x123, y123); |
| 681 | dst[5].set(x23, y23); |
| 682 | dst[6] = src[3]; |
| 683 | } |
| 684 | |
| 685 | static void flatten_double_cubic_extrema(SkScalar coords[14]) |
| 686 | { |
| 687 | coords[4] = coords[8] = coords[6]; |
| 688 | } |
| 689 | |
| 690 | /** Given 4 points on a cubic bezier, chop it into 1, 2, 3 beziers such that |
| 691 | the resulting beziers are monotonic in Y. This is called by the scan converter. |
| 692 | Depending on what is returned, dst[] is treated as follows |
| 693 | 0 dst[0..3] is the original cubic |
| 694 | 1 dst[0..3] and dst[3..6] are the two new cubics |
| 695 | 2 dst[0..3], dst[3..6], dst[6..9] are the three new cubics |
| 696 | If dst == null, it is ignored and only the count is returned. |
| 697 | */ |
reed@android.com | bb13586 | 2009-11-18 13:47:40 +0000 | [diff] [blame] | 698 | int SkChopCubicAtYExtrema(const SkPoint src[4], SkPoint dst[10]) { |
reed@android.com | 8a1c16f | 2008-12-17 15:59:43 +0000 | [diff] [blame] | 699 | SkScalar tValues[2]; |
reed@android.com | bb13586 | 2009-11-18 13:47:40 +0000 | [diff] [blame] | 700 | int roots = SkFindCubicExtrema(src[0].fY, src[1].fY, src[2].fY, |
| 701 | src[3].fY, tValues); |
| 702 | |
reed@android.com | 8a1c16f | 2008-12-17 15:59:43 +0000 | [diff] [blame] | 703 | SkChopCubicAt(src, dst, tValues, roots); |
reed@android.com | bb13586 | 2009-11-18 13:47:40 +0000 | [diff] [blame] | 704 | if (dst && roots > 0) { |
reed@android.com | 8a1c16f | 2008-12-17 15:59:43 +0000 | [diff] [blame] | 705 | // we do some cleanup to ensure our Y extrema are flat |
| 706 | flatten_double_cubic_extrema(&dst[0].fY); |
reed@android.com | bb13586 | 2009-11-18 13:47:40 +0000 | [diff] [blame] | 707 | if (roots == 2) { |
reed@android.com | 8a1c16f | 2008-12-17 15:59:43 +0000 | [diff] [blame] | 708 | flatten_double_cubic_extrema(&dst[3].fY); |
reed@android.com | bb13586 | 2009-11-18 13:47:40 +0000 | [diff] [blame] | 709 | } |
| 710 | } |
| 711 | return roots; |
| 712 | } |
| 713 | |
| 714 | int SkChopCubicAtXExtrema(const SkPoint src[4], SkPoint dst[10]) { |
| 715 | SkScalar tValues[2]; |
| 716 | int roots = SkFindCubicExtrema(src[0].fX, src[1].fX, src[2].fX, |
| 717 | src[3].fX, tValues); |
| 718 | |
| 719 | SkChopCubicAt(src, dst, tValues, roots); |
| 720 | if (dst && roots > 0) { |
| 721 | // we do some cleanup to ensure our Y extrema are flat |
| 722 | flatten_double_cubic_extrema(&dst[0].fX); |
| 723 | if (roots == 2) { |
| 724 | flatten_double_cubic_extrema(&dst[3].fX); |
| 725 | } |
reed@android.com | 8a1c16f | 2008-12-17 15:59:43 +0000 | [diff] [blame] | 726 | } |
| 727 | return roots; |
| 728 | } |
| 729 | |
| 730 | /** http://www.faculty.idc.ac.il/arik/quality/appendixA.html |
| 731 | |
| 732 | Inflection means that curvature is zero. |
| 733 | Curvature is [F' x F''] / [F'^3] |
| 734 | So we solve F'x X F''y - F'y X F''y == 0 |
| 735 | After some canceling of the cubic term, we get |
| 736 | A = b - a |
| 737 | B = c - 2b + a |
| 738 | C = d - 3c + 3b - a |
| 739 | (BxCy - ByCx)t^2 + (AxCy - AyCx)t + AxBy - AyBx == 0 |
| 740 | */ |
| 741 | int SkFindCubicInflections(const SkPoint src[4], SkScalar tValues[]) |
| 742 | { |
| 743 | SkScalar Ax = src[1].fX - src[0].fX; |
| 744 | SkScalar Ay = src[1].fY - src[0].fY; |
| 745 | SkScalar Bx = src[2].fX - 2 * src[1].fX + src[0].fX; |
| 746 | SkScalar By = src[2].fY - 2 * src[1].fY + src[0].fY; |
| 747 | SkScalar Cx = src[3].fX + 3 * (src[1].fX - src[2].fX) - src[0].fX; |
| 748 | SkScalar Cy = src[3].fY + 3 * (src[1].fY - src[2].fY) - src[0].fY; |
| 749 | int count; |
| 750 | |
| 751 | #ifdef SK_SCALAR_IS_FLOAT |
| 752 | count = SkFindUnitQuadRoots(Bx*Cy - By*Cx, Ax*Cy - Ay*Cx, Ax*By - Ay*Bx, tValues); |
| 753 | #else |
| 754 | Sk64 A, B, C, tmp; |
| 755 | |
| 756 | A.setMul(Bx, Cy); |
| 757 | tmp.setMul(By, Cx); |
| 758 | A.sub(tmp); |
| 759 | |
| 760 | B.setMul(Ax, Cy); |
| 761 | tmp.setMul(Ay, Cx); |
| 762 | B.sub(tmp); |
| 763 | |
| 764 | C.setMul(Ax, By); |
| 765 | tmp.setMul(Ay, Bx); |
| 766 | C.sub(tmp); |
| 767 | |
| 768 | count = Sk64FindFixedQuadRoots(A, B, C, tValues); |
| 769 | #endif |
| 770 | |
| 771 | return count; |
| 772 | } |
| 773 | |
| 774 | int SkChopCubicAtInflections(const SkPoint src[], SkPoint dst[10]) |
| 775 | { |
| 776 | SkScalar tValues[2]; |
| 777 | int count = SkFindCubicInflections(src, tValues); |
| 778 | |
| 779 | if (dst) |
| 780 | { |
| 781 | if (count == 0) |
| 782 | memcpy(dst, src, 4 * sizeof(SkPoint)); |
| 783 | else |
| 784 | SkChopCubicAt(src, dst, tValues, count); |
| 785 | } |
| 786 | return count + 1; |
| 787 | } |
| 788 | |
| 789 | template <typename T> void bubble_sort(T array[], int count) |
| 790 | { |
| 791 | for (int i = count - 1; i > 0; --i) |
| 792 | for (int j = i; j > 0; --j) |
| 793 | if (array[j] < array[j-1]) |
| 794 | { |
| 795 | T tmp(array[j]); |
| 796 | array[j] = array[j-1]; |
| 797 | array[j-1] = tmp; |
| 798 | } |
| 799 | } |
| 800 | |
| 801 | #include "SkFP.h" |
| 802 | |
| 803 | // newton refinement |
| 804 | #if 0 |
| 805 | static SkScalar refine_cubic_root(const SkFP coeff[4], SkScalar root) |
| 806 | { |
| 807 | // x1 = x0 - f(t) / f'(t) |
| 808 | |
| 809 | SkFP T = SkScalarToFloat(root); |
| 810 | SkFP N, D; |
| 811 | |
| 812 | // f' = 3*coeff[0]*T^2 + 2*coeff[1]*T + coeff[2] |
| 813 | D = SkFPMul(SkFPMul(coeff[0], SkFPMul(T,T)), 3); |
| 814 | D = SkFPAdd(D, SkFPMulInt(SkFPMul(coeff[1], T), 2)); |
| 815 | D = SkFPAdd(D, coeff[2]); |
| 816 | |
| 817 | if (D == 0) |
| 818 | return root; |
| 819 | |
| 820 | // f = coeff[0]*T^3 + coeff[1]*T^2 + coeff[2]*T + coeff[3] |
| 821 | N = SkFPMul(SkFPMul(SkFPMul(T, T), T), coeff[0]); |
| 822 | N = SkFPAdd(N, SkFPMul(SkFPMul(T, T), coeff[1])); |
| 823 | N = SkFPAdd(N, SkFPMul(T, coeff[2])); |
| 824 | N = SkFPAdd(N, coeff[3]); |
| 825 | |
| 826 | if (N) |
| 827 | { |
| 828 | SkScalar delta = SkFPToScalar(SkFPDiv(N, D)); |
| 829 | |
| 830 | if (delta) |
| 831 | root -= delta; |
| 832 | } |
| 833 | return root; |
| 834 | } |
| 835 | #endif |
| 836 | |
| 837 | #if defined _WIN32 && _MSC_VER >= 1300 && defined SK_SCALAR_IS_FIXED // disable warning : unreachable code if building fixed point for windows desktop |
| 838 | #pragma warning ( disable : 4702 ) |
| 839 | #endif |
| 840 | |
| 841 | /* Solve coeff(t) == 0, returning the number of roots that |
| 842 | lie withing 0 < t < 1. |
| 843 | coeff[0]t^3 + coeff[1]t^2 + coeff[2]t + coeff[3] |
| 844 | */ |
| 845 | static int solve_cubic_polynomial(const SkFP coeff[4], SkScalar tValues[3]) |
| 846 | { |
| 847 | #ifndef SK_SCALAR_IS_FLOAT |
| 848 | return 0; // this is not yet implemented for software float |
| 849 | #endif |
| 850 | |
| 851 | if (SkScalarNearlyZero(coeff[0])) // we're just a quadratic |
| 852 | { |
| 853 | return SkFindUnitQuadRoots(coeff[1], coeff[2], coeff[3], tValues); |
| 854 | } |
| 855 | |
| 856 | SkFP a, b, c, Q, R; |
| 857 | |
| 858 | { |
| 859 | SkASSERT(coeff[0] != 0); |
| 860 | |
| 861 | SkFP inva = SkFPInvert(coeff[0]); |
| 862 | a = SkFPMul(coeff[1], inva); |
| 863 | b = SkFPMul(coeff[2], inva); |
| 864 | c = SkFPMul(coeff[3], inva); |
| 865 | } |
| 866 | Q = SkFPDivInt(SkFPSub(SkFPMul(a,a), SkFPMulInt(b, 3)), 9); |
| 867 | // R = (2*a*a*a - 9*a*b + 27*c) / 54; |
| 868 | R = SkFPMulInt(SkFPMul(SkFPMul(a, a), a), 2); |
| 869 | R = SkFPSub(R, SkFPMulInt(SkFPMul(a, b), 9)); |
| 870 | R = SkFPAdd(R, SkFPMulInt(c, 27)); |
| 871 | R = SkFPDivInt(R, 54); |
| 872 | |
| 873 | SkFP Q3 = SkFPMul(SkFPMul(Q, Q), Q); |
| 874 | SkFP R2MinusQ3 = SkFPSub(SkFPMul(R,R), Q3); |
| 875 | SkFP adiv3 = SkFPDivInt(a, 3); |
| 876 | |
| 877 | SkScalar* roots = tValues; |
| 878 | SkScalar r; |
| 879 | |
| 880 | if (SkFPLT(R2MinusQ3, 0)) // we have 3 real roots |
| 881 | { |
| 882 | #ifdef SK_SCALAR_IS_FLOAT |
| 883 | float theta = sk_float_acos(R / sk_float_sqrt(Q3)); |
| 884 | float neg2RootQ = -2 * sk_float_sqrt(Q); |
| 885 | |
| 886 | r = neg2RootQ * sk_float_cos(theta/3) - adiv3; |
| 887 | if (is_unit_interval(r)) |
| 888 | *roots++ = r; |
| 889 | |
| 890 | r = neg2RootQ * sk_float_cos((theta + 2*SK_ScalarPI)/3) - adiv3; |
| 891 | if (is_unit_interval(r)) |
| 892 | *roots++ = r; |
| 893 | |
| 894 | r = neg2RootQ * sk_float_cos((theta - 2*SK_ScalarPI)/3) - adiv3; |
| 895 | if (is_unit_interval(r)) |
| 896 | *roots++ = r; |
| 897 | |
| 898 | // now sort the roots |
| 899 | bubble_sort(tValues, (int)(roots - tValues)); |
| 900 | #endif |
| 901 | } |
| 902 | else // we have 1 real root |
| 903 | { |
| 904 | SkFP A = SkFPAdd(SkFPAbs(R), SkFPSqrt(R2MinusQ3)); |
| 905 | A = SkFPCubeRoot(A); |
| 906 | if (SkFPGT(R, 0)) |
| 907 | A = SkFPNeg(A); |
| 908 | |
| 909 | if (A != 0) |
| 910 | A = SkFPAdd(A, SkFPDiv(Q, A)); |
| 911 | r = SkFPToScalar(SkFPSub(A, adiv3)); |
| 912 | if (is_unit_interval(r)) |
| 913 | *roots++ = r; |
| 914 | } |
| 915 | |
| 916 | return (int)(roots - tValues); |
| 917 | } |
| 918 | |
| 919 | /* Looking for F' dot F'' == 0 |
| 920 | |
| 921 | A = b - a |
| 922 | B = c - 2b + a |
| 923 | C = d - 3c + 3b - a |
| 924 | |
| 925 | F' = 3Ct^2 + 6Bt + 3A |
| 926 | F'' = 6Ct + 6B |
| 927 | |
| 928 | F' dot F'' -> CCt^3 + 3BCt^2 + (2BB + CA)t + AB |
| 929 | */ |
| 930 | static void formulate_F1DotF2(const SkScalar src[], SkFP coeff[4]) |
| 931 | { |
| 932 | SkScalar a = src[2] - src[0]; |
| 933 | SkScalar b = src[4] - 2 * src[2] + src[0]; |
| 934 | SkScalar c = src[6] + 3 * (src[2] - src[4]) - src[0]; |
| 935 | |
| 936 | SkFP A = SkScalarToFP(a); |
| 937 | SkFP B = SkScalarToFP(b); |
| 938 | SkFP C = SkScalarToFP(c); |
| 939 | |
| 940 | coeff[0] = SkFPMul(C, C); |
| 941 | coeff[1] = SkFPMulInt(SkFPMul(B, C), 3); |
| 942 | coeff[2] = SkFPMulInt(SkFPMul(B, B), 2); |
| 943 | coeff[2] = SkFPAdd(coeff[2], SkFPMul(C, A)); |
| 944 | coeff[3] = SkFPMul(A, B); |
| 945 | } |
| 946 | |
| 947 | // EXPERIMENTAL: can set this to zero to accept all t-values 0 < t < 1 |
| 948 | //#define kMinTValueForChopping (SK_Scalar1 / 256) |
| 949 | #define kMinTValueForChopping 0 |
| 950 | |
| 951 | /* Looking for F' dot F'' == 0 |
| 952 | |
| 953 | A = b - a |
| 954 | B = c - 2b + a |
| 955 | C = d - 3c + 3b - a |
| 956 | |
| 957 | F' = 3Ct^2 + 6Bt + 3A |
| 958 | F'' = 6Ct + 6B |
| 959 | |
| 960 | F' dot F'' -> CCt^3 + 3BCt^2 + (2BB + CA)t + AB |
| 961 | */ |
| 962 | int SkFindCubicMaxCurvature(const SkPoint src[4], SkScalar tValues[3]) |
| 963 | { |
| 964 | SkFP coeffX[4], coeffY[4]; |
| 965 | int i; |
| 966 | |
| 967 | formulate_F1DotF2(&src[0].fX, coeffX); |
| 968 | formulate_F1DotF2(&src[0].fY, coeffY); |
| 969 | |
| 970 | for (i = 0; i < 4; i++) |
| 971 | coeffX[i] = SkFPAdd(coeffX[i],coeffY[i]); |
| 972 | |
| 973 | SkScalar t[3]; |
| 974 | int count = solve_cubic_polynomial(coeffX, t); |
| 975 | int maxCount = 0; |
| 976 | |
| 977 | // now remove extrema where the curvature is zero (mins) |
| 978 | // !!!! need a test for this !!!! |
| 979 | for (i = 0; i < count; i++) |
| 980 | { |
| 981 | // if (not_min_curvature()) |
| 982 | if (t[i] > kMinTValueForChopping && t[i] < SK_Scalar1 - kMinTValueForChopping) |
| 983 | tValues[maxCount++] = t[i]; |
| 984 | } |
| 985 | return maxCount; |
| 986 | } |
| 987 | |
| 988 | int SkChopCubicAtMaxCurvature(const SkPoint src[4], SkPoint dst[13], SkScalar tValues[3]) |
| 989 | { |
| 990 | SkScalar t_storage[3]; |
| 991 | |
| 992 | if (tValues == NULL) |
| 993 | tValues = t_storage; |
| 994 | |
| 995 | int count = SkFindCubicMaxCurvature(src, tValues); |
| 996 | |
| 997 | if (dst) |
| 998 | { |
| 999 | if (count == 0) |
| 1000 | memcpy(dst, src, 4 * sizeof(SkPoint)); |
| 1001 | else |
| 1002 | SkChopCubicAt(src, dst, tValues, count); |
| 1003 | } |
| 1004 | return count + 1; |
| 1005 | } |
| 1006 | |
kbr@chromium.org | 2e08619 | 2010-07-07 22:20:35 +0000 | [diff] [blame] | 1007 | bool SkXRayCrossesMonotonicCubic(const SkXRay& pt, const SkPoint cubic[4], bool* ambiguous) { |
| 1008 | if (ambiguous) { |
| 1009 | *ambiguous = false; |
| 1010 | } |
| 1011 | |
reed@android.com | 945a139 | 2010-02-05 20:41:02 +0000 | [diff] [blame] | 1012 | // Find the minimum and maximum y of the extrema, which are the |
| 1013 | // first and last points since this cubic is monotonic |
| 1014 | SkScalar min_y = SkMinScalar(cubic[0].fY, cubic[3].fY); |
| 1015 | SkScalar max_y = SkMaxScalar(cubic[0].fY, cubic[3].fY); |
| 1016 | |
| 1017 | if (pt.fY == cubic[0].fY |
| 1018 | || pt.fY < min_y |
| 1019 | || pt.fY > max_y) { |
| 1020 | // The query line definitely does not cross the curve |
kbr@chromium.org | 2e08619 | 2010-07-07 22:20:35 +0000 | [diff] [blame] | 1021 | if (ambiguous) { |
| 1022 | *ambiguous = (pt.fY == cubic[0].fY); |
| 1023 | } |
reed@android.com | 945a139 | 2010-02-05 20:41:02 +0000 | [diff] [blame] | 1024 | return false; |
| 1025 | } |
| 1026 | |
kbr@chromium.org | 2e08619 | 2010-07-07 22:20:35 +0000 | [diff] [blame] | 1027 | bool pt_at_extremum = (pt.fY == cubic[3].fY); |
| 1028 | |
reed@android.com | 945a139 | 2010-02-05 20:41:02 +0000 | [diff] [blame] | 1029 | SkScalar min_x = |
| 1030 | SkMinScalar( |
| 1031 | SkMinScalar( |
| 1032 | SkMinScalar(cubic[0].fX, cubic[1].fX), |
| 1033 | cubic[2].fX), |
| 1034 | cubic[3].fX); |
| 1035 | if (pt.fX < min_x) { |
| 1036 | // The query line definitely crosses the curve |
kbr@chromium.org | 2e08619 | 2010-07-07 22:20:35 +0000 | [diff] [blame] | 1037 | if (ambiguous) { |
| 1038 | *ambiguous = pt_at_extremum; |
| 1039 | } |
reed@android.com | 945a139 | 2010-02-05 20:41:02 +0000 | [diff] [blame] | 1040 | return true; |
| 1041 | } |
| 1042 | |
| 1043 | SkScalar max_x = |
| 1044 | SkMaxScalar( |
| 1045 | SkMaxScalar( |
| 1046 | SkMaxScalar(cubic[0].fX, cubic[1].fX), |
| 1047 | cubic[2].fX), |
| 1048 | cubic[3].fX); |
| 1049 | if (pt.fX > max_x) { |
| 1050 | // The query line definitely does not cross the curve |
| 1051 | return false; |
| 1052 | } |
| 1053 | |
| 1054 | // Do a binary search to find the parameter value which makes y as |
| 1055 | // close as possible to the query point. See whether the query |
| 1056 | // line's origin is to the left of the associated x coordinate. |
| 1057 | |
| 1058 | // kMaxIter is chosen as the number of mantissa bits for a float, |
| 1059 | // since there's no way we are going to get more precision by |
| 1060 | // iterating more times than that. |
| 1061 | const int kMaxIter = 23; |
| 1062 | SkPoint eval; |
| 1063 | int iter = 0; |
| 1064 | SkScalar upper_t; |
| 1065 | SkScalar lower_t; |
| 1066 | // Need to invert direction of t parameter if cubic goes up |
| 1067 | // instead of down |
| 1068 | if (cubic[3].fY > cubic[0].fY) { |
| 1069 | upper_t = SK_Scalar1; |
| 1070 | lower_t = SkFloatToScalar(0); |
| 1071 | } else { |
| 1072 | upper_t = SkFloatToScalar(0); |
| 1073 | lower_t = SK_Scalar1; |
| 1074 | } |
| 1075 | do { |
| 1076 | SkScalar t = SkScalarAve(upper_t, lower_t); |
| 1077 | SkEvalCubicAt(cubic, t, &eval, NULL, NULL); |
| 1078 | if (pt.fY > eval.fY) { |
| 1079 | lower_t = t; |
| 1080 | } else { |
| 1081 | upper_t = t; |
| 1082 | } |
| 1083 | } while (++iter < kMaxIter |
| 1084 | && !SkScalarNearlyZero(eval.fY - pt.fY)); |
| 1085 | if (pt.fX <= eval.fX) { |
kbr@chromium.org | 2e08619 | 2010-07-07 22:20:35 +0000 | [diff] [blame] | 1086 | if (ambiguous) { |
| 1087 | *ambiguous = pt_at_extremum; |
| 1088 | } |
reed@android.com | 945a139 | 2010-02-05 20:41:02 +0000 | [diff] [blame] | 1089 | return true; |
| 1090 | } |
| 1091 | return false; |
| 1092 | } |
| 1093 | |
kbr@chromium.org | 2e08619 | 2010-07-07 22:20:35 +0000 | [diff] [blame] | 1094 | int SkNumXRayCrossingsForCubic(const SkXRay& pt, const SkPoint cubic[4], bool* ambiguous) { |
reed@android.com | 945a139 | 2010-02-05 20:41:02 +0000 | [diff] [blame] | 1095 | int num_crossings = 0; |
| 1096 | SkPoint monotonic_cubics[10]; |
| 1097 | int num_monotonic_cubics = SkChopCubicAtYExtrema(cubic, monotonic_cubics); |
kbr@chromium.org | 2e08619 | 2010-07-07 22:20:35 +0000 | [diff] [blame] | 1098 | if (ambiguous) { |
| 1099 | *ambiguous = false; |
| 1100 | } |
| 1101 | bool locally_ambiguous; |
| 1102 | if (SkXRayCrossesMonotonicCubic(pt, &monotonic_cubics[0], &locally_ambiguous)) |
reed@android.com | 945a139 | 2010-02-05 20:41:02 +0000 | [diff] [blame] | 1103 | ++num_crossings; |
kbr@chromium.org | 2e08619 | 2010-07-07 22:20:35 +0000 | [diff] [blame] | 1104 | if (ambiguous) { |
| 1105 | *ambiguous |= locally_ambiguous; |
| 1106 | } |
reed@android.com | 945a139 | 2010-02-05 20:41:02 +0000 | [diff] [blame] | 1107 | if (num_monotonic_cubics > 0) |
kbr@chromium.org | 2e08619 | 2010-07-07 22:20:35 +0000 | [diff] [blame] | 1108 | if (SkXRayCrossesMonotonicCubic(pt, &monotonic_cubics[3], &locally_ambiguous)) |
reed@android.com | 945a139 | 2010-02-05 20:41:02 +0000 | [diff] [blame] | 1109 | ++num_crossings; |
kbr@chromium.org | 2e08619 | 2010-07-07 22:20:35 +0000 | [diff] [blame] | 1110 | if (ambiguous) { |
| 1111 | *ambiguous |= locally_ambiguous; |
| 1112 | } |
reed@android.com | 945a139 | 2010-02-05 20:41:02 +0000 | [diff] [blame] | 1113 | if (num_monotonic_cubics > 1) |
kbr@chromium.org | 2e08619 | 2010-07-07 22:20:35 +0000 | [diff] [blame] | 1114 | if (SkXRayCrossesMonotonicCubic(pt, &monotonic_cubics[6], &locally_ambiguous)) |
reed@android.com | 945a139 | 2010-02-05 20:41:02 +0000 | [diff] [blame] | 1115 | ++num_crossings; |
kbr@chromium.org | 2e08619 | 2010-07-07 22:20:35 +0000 | [diff] [blame] | 1116 | if (ambiguous) { |
| 1117 | *ambiguous |= locally_ambiguous; |
| 1118 | } |
reed@android.com | 945a139 | 2010-02-05 20:41:02 +0000 | [diff] [blame] | 1119 | return num_crossings; |
| 1120 | } |
| 1121 | |
reed@android.com | 8a1c16f | 2008-12-17 15:59:43 +0000 | [diff] [blame] | 1122 | //////////////////////////////////////////////////////////////////////////////// |
| 1123 | |
| 1124 | /* Find t value for quadratic [a, b, c] = d. |
| 1125 | Return 0 if there is no solution within [0, 1) |
| 1126 | */ |
| 1127 | static SkScalar quad_solve(SkScalar a, SkScalar b, SkScalar c, SkScalar d) |
| 1128 | { |
| 1129 | // At^2 + Bt + C = d |
| 1130 | SkScalar A = a - 2 * b + c; |
| 1131 | SkScalar B = 2 * (b - a); |
| 1132 | SkScalar C = a - d; |
| 1133 | |
| 1134 | SkScalar roots[2]; |
| 1135 | int count = SkFindUnitQuadRoots(A, B, C, roots); |
| 1136 | |
| 1137 | SkASSERT(count <= 1); |
| 1138 | return count == 1 ? roots[0] : 0; |
| 1139 | } |
| 1140 | |
| 1141 | /* given a quad-curve and a point (x,y), chop the quad at that point and return |
| 1142 | the new quad's offCurve point. Should only return false if the computed pos |
| 1143 | is the start of the curve (i.e. root == 0) |
| 1144 | */ |
| 1145 | static bool quad_pt2OffCurve(const SkPoint quad[3], SkScalar x, SkScalar y, SkPoint* offCurve) |
| 1146 | { |
| 1147 | const SkScalar* base; |
| 1148 | SkScalar value; |
| 1149 | |
| 1150 | if (SkScalarAbs(x) < SkScalarAbs(y)) { |
| 1151 | base = &quad[0].fX; |
| 1152 | value = x; |
| 1153 | } else { |
| 1154 | base = &quad[0].fY; |
| 1155 | value = y; |
| 1156 | } |
| 1157 | |
| 1158 | // note: this returns 0 if it thinks value is out of range, meaning the |
| 1159 | // root might return something outside of [0, 1) |
| 1160 | SkScalar t = quad_solve(base[0], base[2], base[4], value); |
| 1161 | |
| 1162 | if (t > 0) |
| 1163 | { |
| 1164 | SkPoint tmp[5]; |
| 1165 | SkChopQuadAt(quad, tmp, t); |
| 1166 | *offCurve = tmp[1]; |
| 1167 | return true; |
| 1168 | } else { |
| 1169 | /* t == 0 means either the value triggered a root outside of [0, 1) |
| 1170 | For our purposes, we can ignore the <= 0 roots, but we want to |
| 1171 | catch the >= 1 roots (which given our caller, will basically mean |
| 1172 | a root of 1, give-or-take numerical instability). If we are in the |
| 1173 | >= 1 case, return the existing offCurve point. |
| 1174 | |
| 1175 | The test below checks to see if we are close to the "end" of the |
| 1176 | curve (near base[4]). Rather than specifying a tolerance, I just |
| 1177 | check to see if value is on to the right/left of the middle point |
| 1178 | (depending on the direction/sign of the end points). |
| 1179 | */ |
| 1180 | if ((base[0] < base[4] && value > base[2]) || |
| 1181 | (base[0] > base[4] && value < base[2])) // should root have been 1 |
| 1182 | { |
| 1183 | *offCurve = quad[1]; |
| 1184 | return true; |
| 1185 | } |
| 1186 | } |
| 1187 | return false; |
| 1188 | } |
| 1189 | |
| 1190 | static const SkPoint gQuadCirclePts[kSkBuildQuadArcStorage] = { |
| 1191 | { SK_Scalar1, 0 }, |
| 1192 | { SK_Scalar1, SK_ScalarTanPIOver8 }, |
| 1193 | { SK_ScalarRoot2Over2, SK_ScalarRoot2Over2 }, |
| 1194 | { SK_ScalarTanPIOver8, SK_Scalar1 }, |
| 1195 | |
| 1196 | { 0, SK_Scalar1 }, |
| 1197 | { -SK_ScalarTanPIOver8, SK_Scalar1 }, |
| 1198 | { -SK_ScalarRoot2Over2, SK_ScalarRoot2Over2 }, |
| 1199 | { -SK_Scalar1, SK_ScalarTanPIOver8 }, |
| 1200 | |
| 1201 | { -SK_Scalar1, 0 }, |
| 1202 | { -SK_Scalar1, -SK_ScalarTanPIOver8 }, |
| 1203 | { -SK_ScalarRoot2Over2, -SK_ScalarRoot2Over2 }, |
| 1204 | { -SK_ScalarTanPIOver8, -SK_Scalar1 }, |
| 1205 | |
| 1206 | { 0, -SK_Scalar1 }, |
| 1207 | { SK_ScalarTanPIOver8, -SK_Scalar1 }, |
| 1208 | { SK_ScalarRoot2Over2, -SK_ScalarRoot2Over2 }, |
| 1209 | { SK_Scalar1, -SK_ScalarTanPIOver8 }, |
| 1210 | |
| 1211 | { SK_Scalar1, 0 } |
| 1212 | }; |
| 1213 | |
| 1214 | int SkBuildQuadArc(const SkVector& uStart, const SkVector& uStop, |
| 1215 | SkRotationDirection dir, const SkMatrix* userMatrix, |
| 1216 | SkPoint quadPoints[]) |
| 1217 | { |
| 1218 | // rotate by x,y so that uStart is (1.0) |
| 1219 | SkScalar x = SkPoint::DotProduct(uStart, uStop); |
| 1220 | SkScalar y = SkPoint::CrossProduct(uStart, uStop); |
| 1221 | |
| 1222 | SkScalar absX = SkScalarAbs(x); |
| 1223 | SkScalar absY = SkScalarAbs(y); |
| 1224 | |
| 1225 | int pointCount; |
| 1226 | |
| 1227 | // check for (effectively) coincident vectors |
| 1228 | // this can happen if our angle is nearly 0 or nearly 180 (y == 0) |
| 1229 | // ... we use the dot-prod to distinguish between 0 and 180 (x > 0) |
| 1230 | if (absY <= SK_ScalarNearlyZero && x > 0 && |
| 1231 | ((y >= 0 && kCW_SkRotationDirection == dir) || |
| 1232 | (y <= 0 && kCCW_SkRotationDirection == dir))) { |
| 1233 | |
| 1234 | // just return the start-point |
| 1235 | quadPoints[0].set(SK_Scalar1, 0); |
| 1236 | pointCount = 1; |
| 1237 | } else { |
| 1238 | if (dir == kCCW_SkRotationDirection) |
| 1239 | y = -y; |
| 1240 | |
| 1241 | // what octant (quadratic curve) is [xy] in? |
| 1242 | int oct = 0; |
| 1243 | bool sameSign = true; |
| 1244 | |
| 1245 | if (0 == y) |
| 1246 | { |
| 1247 | oct = 4; // 180 |
| 1248 | SkASSERT(SkScalarAbs(x + SK_Scalar1) <= SK_ScalarNearlyZero); |
| 1249 | } |
| 1250 | else if (0 == x) |
| 1251 | { |
| 1252 | SkASSERT(absY - SK_Scalar1 <= SK_ScalarNearlyZero); |
| 1253 | if (y > 0) |
| 1254 | oct = 2; // 90 |
| 1255 | else |
| 1256 | oct = 6; // 270 |
| 1257 | } |
| 1258 | else |
| 1259 | { |
| 1260 | if (y < 0) |
| 1261 | oct += 4; |
| 1262 | if ((x < 0) != (y < 0)) |
| 1263 | { |
| 1264 | oct += 2; |
| 1265 | sameSign = false; |
| 1266 | } |
| 1267 | if ((absX < absY) == sameSign) |
| 1268 | oct += 1; |
| 1269 | } |
| 1270 | |
| 1271 | int wholeCount = oct << 1; |
| 1272 | memcpy(quadPoints, gQuadCirclePts, (wholeCount + 1) * sizeof(SkPoint)); |
| 1273 | |
| 1274 | const SkPoint* arc = &gQuadCirclePts[wholeCount]; |
| 1275 | if (quad_pt2OffCurve(arc, x, y, &quadPoints[wholeCount + 1])) |
| 1276 | { |
| 1277 | quadPoints[wholeCount + 2].set(x, y); |
| 1278 | wholeCount += 2; |
| 1279 | } |
| 1280 | pointCount = wholeCount + 1; |
| 1281 | } |
| 1282 | |
| 1283 | // now handle counter-clockwise and the initial unitStart rotation |
| 1284 | SkMatrix matrix; |
| 1285 | matrix.setSinCos(uStart.fY, uStart.fX); |
| 1286 | if (dir == kCCW_SkRotationDirection) { |
| 1287 | matrix.preScale(SK_Scalar1, -SK_Scalar1); |
| 1288 | } |
| 1289 | if (userMatrix) { |
| 1290 | matrix.postConcat(*userMatrix); |
| 1291 | } |
| 1292 | matrix.mapPoints(quadPoints, pointCount); |
| 1293 | return pointCount; |
| 1294 | } |
| 1295 | |