1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/gfx/skia/trunk/src/core/SkGeometry.cpp Wed Dec 31 06:09:35 2014 +0100
1.3 @@ -0,0 +1,1468 @@
1.4 +/*
1.5 + * Copyright 2006 The Android Open Source Project
1.6 + *
1.7 + * Use of this source code is governed by a BSD-style license that can be
1.8 + * found in the LICENSE file.
1.9 + */
1.10 +
1.11 +#include "SkGeometry.h"
1.12 +#include "SkMatrix.h"
1.13 +
1.14 +bool SkXRayCrossesLine(const SkXRay& pt,
1.15 + const SkPoint pts[2],
1.16 + bool* ambiguous) {
1.17 + if (ambiguous) {
1.18 + *ambiguous = false;
1.19 + }
1.20 + // Determine quick discards.
1.21 + // Consider query line going exactly through point 0 to not
1.22 + // intersect, for symmetry with SkXRayCrossesMonotonicCubic.
1.23 + if (pt.fY == pts[0].fY) {
1.24 + if (ambiguous) {
1.25 + *ambiguous = true;
1.26 + }
1.27 + return false;
1.28 + }
1.29 + if (pt.fY < pts[0].fY && pt.fY < pts[1].fY)
1.30 + return false;
1.31 + if (pt.fY > pts[0].fY && pt.fY > pts[1].fY)
1.32 + return false;
1.33 + if (pt.fX > pts[0].fX && pt.fX > pts[1].fX)
1.34 + return false;
1.35 + // Determine degenerate cases
1.36 + if (SkScalarNearlyZero(pts[0].fY - pts[1].fY))
1.37 + return false;
1.38 + if (SkScalarNearlyZero(pts[0].fX - pts[1].fX)) {
1.39 + // We've already determined the query point lies within the
1.40 + // vertical range of the line segment.
1.41 + if (pt.fX <= pts[0].fX) {
1.42 + if (ambiguous) {
1.43 + *ambiguous = (pt.fY == pts[1].fY);
1.44 + }
1.45 + return true;
1.46 + }
1.47 + return false;
1.48 + }
1.49 + // Ambiguity check
1.50 + if (pt.fY == pts[1].fY) {
1.51 + if (pt.fX <= pts[1].fX) {
1.52 + if (ambiguous) {
1.53 + *ambiguous = true;
1.54 + }
1.55 + return true;
1.56 + }
1.57 + return false;
1.58 + }
1.59 + // Full line segment evaluation
1.60 + SkScalar delta_y = pts[1].fY - pts[0].fY;
1.61 + SkScalar delta_x = pts[1].fX - pts[0].fX;
1.62 + SkScalar slope = SkScalarDiv(delta_y, delta_x);
1.63 + SkScalar b = pts[0].fY - SkScalarMul(slope, pts[0].fX);
1.64 + // Solve for x coordinate at y = pt.fY
1.65 + SkScalar x = SkScalarDiv(pt.fY - b, slope);
1.66 + return pt.fX <= x;
1.67 +}
1.68 +
1.69 +/** If defined, this makes eval_quad and eval_cubic do more setup (sometimes
1.70 + involving integer multiplies by 2 or 3, but fewer calls to SkScalarMul.
1.71 + May also introduce overflow of fixed when we compute our setup.
1.72 +*/
1.73 +// #define DIRECT_EVAL_OF_POLYNOMIALS
1.74 +
1.75 +////////////////////////////////////////////////////////////////////////
1.76 +
1.77 +static int is_not_monotonic(SkScalar a, SkScalar b, SkScalar c) {
1.78 + SkScalar ab = a - b;
1.79 + SkScalar bc = b - c;
1.80 + if (ab < 0) {
1.81 + bc = -bc;
1.82 + }
1.83 + return ab == 0 || bc < 0;
1.84 +}
1.85 +
1.86 +////////////////////////////////////////////////////////////////////////
1.87 +
1.88 +static bool is_unit_interval(SkScalar x) {
1.89 + return x > 0 && x < SK_Scalar1;
1.90 +}
1.91 +
1.92 +static int valid_unit_divide(SkScalar numer, SkScalar denom, SkScalar* ratio) {
1.93 + SkASSERT(ratio);
1.94 +
1.95 + if (numer < 0) {
1.96 + numer = -numer;
1.97 + denom = -denom;
1.98 + }
1.99 +
1.100 + if (denom == 0 || numer == 0 || numer >= denom) {
1.101 + return 0;
1.102 + }
1.103 +
1.104 + SkScalar r = SkScalarDiv(numer, denom);
1.105 + if (SkScalarIsNaN(r)) {
1.106 + return 0;
1.107 + }
1.108 + SkASSERT(r >= 0 && r < SK_Scalar1);
1.109 + if (r == 0) { // catch underflow if numer <<<< denom
1.110 + return 0;
1.111 + }
1.112 + *ratio = r;
1.113 + return 1;
1.114 +}
1.115 +
1.116 +/** From Numerical Recipes in C.
1.117 +
1.118 + Q = -1/2 (B + sign(B) sqrt[B*B - 4*A*C])
1.119 + x1 = Q / A
1.120 + x2 = C / Q
1.121 +*/
1.122 +int SkFindUnitQuadRoots(SkScalar A, SkScalar B, SkScalar C, SkScalar roots[2]) {
1.123 + SkASSERT(roots);
1.124 +
1.125 + if (A == 0) {
1.126 + return valid_unit_divide(-C, B, roots);
1.127 + }
1.128 +
1.129 + SkScalar* r = roots;
1.130 +
1.131 + SkScalar R = B*B - 4*A*C;
1.132 + if (R < 0 || SkScalarIsNaN(R)) { // complex roots
1.133 + return 0;
1.134 + }
1.135 + R = SkScalarSqrt(R);
1.136 +
1.137 + SkScalar Q = (B < 0) ? -(B-R)/2 : -(B+R)/2;
1.138 + r += valid_unit_divide(Q, A, r);
1.139 + r += valid_unit_divide(C, Q, r);
1.140 + if (r - roots == 2) {
1.141 + if (roots[0] > roots[1])
1.142 + SkTSwap<SkScalar>(roots[0], roots[1]);
1.143 + else if (roots[0] == roots[1]) // nearly-equal?
1.144 + r -= 1; // skip the double root
1.145 + }
1.146 + return (int)(r - roots);
1.147 +}
1.148 +
1.149 +///////////////////////////////////////////////////////////////////////////////
1.150 +///////////////////////////////////////////////////////////////////////////////
1.151 +
1.152 +static SkScalar eval_quad(const SkScalar src[], SkScalar t) {
1.153 + SkASSERT(src);
1.154 + SkASSERT(t >= 0 && t <= SK_Scalar1);
1.155 +
1.156 +#ifdef DIRECT_EVAL_OF_POLYNOMIALS
1.157 + SkScalar C = src[0];
1.158 + SkScalar A = src[4] - 2 * src[2] + C;
1.159 + SkScalar B = 2 * (src[2] - C);
1.160 + return SkScalarMulAdd(SkScalarMulAdd(A, t, B), t, C);
1.161 +#else
1.162 + SkScalar ab = SkScalarInterp(src[0], src[2], t);
1.163 + SkScalar bc = SkScalarInterp(src[2], src[4], t);
1.164 + return SkScalarInterp(ab, bc, t);
1.165 +#endif
1.166 +}
1.167 +
1.168 +static SkScalar eval_quad_derivative(const SkScalar src[], SkScalar t) {
1.169 + SkScalar A = src[4] - 2 * src[2] + src[0];
1.170 + SkScalar B = src[2] - src[0];
1.171 +
1.172 + return 2 * SkScalarMulAdd(A, t, B);
1.173 +}
1.174 +
1.175 +static SkScalar eval_quad_derivative_at_half(const SkScalar src[]) {
1.176 + SkScalar A = src[4] - 2 * src[2] + src[0];
1.177 + SkScalar B = src[2] - src[0];
1.178 + return A + 2 * B;
1.179 +}
1.180 +
1.181 +void SkEvalQuadAt(const SkPoint src[3], SkScalar t, SkPoint* pt,
1.182 + SkVector* tangent) {
1.183 + SkASSERT(src);
1.184 + SkASSERT(t >= 0 && t <= SK_Scalar1);
1.185 +
1.186 + if (pt) {
1.187 + pt->set(eval_quad(&src[0].fX, t), eval_quad(&src[0].fY, t));
1.188 + }
1.189 + if (tangent) {
1.190 + tangent->set(eval_quad_derivative(&src[0].fX, t),
1.191 + eval_quad_derivative(&src[0].fY, t));
1.192 + }
1.193 +}
1.194 +
1.195 +void SkEvalQuadAtHalf(const SkPoint src[3], SkPoint* pt, SkVector* tangent) {
1.196 + SkASSERT(src);
1.197 +
1.198 + if (pt) {
1.199 + SkScalar x01 = SkScalarAve(src[0].fX, src[1].fX);
1.200 + SkScalar y01 = SkScalarAve(src[0].fY, src[1].fY);
1.201 + SkScalar x12 = SkScalarAve(src[1].fX, src[2].fX);
1.202 + SkScalar y12 = SkScalarAve(src[1].fY, src[2].fY);
1.203 + pt->set(SkScalarAve(x01, x12), SkScalarAve(y01, y12));
1.204 + }
1.205 + if (tangent) {
1.206 + tangent->set(eval_quad_derivative_at_half(&src[0].fX),
1.207 + eval_quad_derivative_at_half(&src[0].fY));
1.208 + }
1.209 +}
1.210 +
1.211 +static void interp_quad_coords(const SkScalar* src, SkScalar* dst, SkScalar t) {
1.212 + SkScalar ab = SkScalarInterp(src[0], src[2], t);
1.213 + SkScalar bc = SkScalarInterp(src[2], src[4], t);
1.214 +
1.215 + dst[0] = src[0];
1.216 + dst[2] = ab;
1.217 + dst[4] = SkScalarInterp(ab, bc, t);
1.218 + dst[6] = bc;
1.219 + dst[8] = src[4];
1.220 +}
1.221 +
1.222 +void SkChopQuadAt(const SkPoint src[3], SkPoint dst[5], SkScalar t) {
1.223 + SkASSERT(t > 0 && t < SK_Scalar1);
1.224 +
1.225 + interp_quad_coords(&src[0].fX, &dst[0].fX, t);
1.226 + interp_quad_coords(&src[0].fY, &dst[0].fY, t);
1.227 +}
1.228 +
1.229 +void SkChopQuadAtHalf(const SkPoint src[3], SkPoint dst[5]) {
1.230 + SkScalar x01 = SkScalarAve(src[0].fX, src[1].fX);
1.231 + SkScalar y01 = SkScalarAve(src[0].fY, src[1].fY);
1.232 + SkScalar x12 = SkScalarAve(src[1].fX, src[2].fX);
1.233 + SkScalar y12 = SkScalarAve(src[1].fY, src[2].fY);
1.234 +
1.235 + dst[0] = src[0];
1.236 + dst[1].set(x01, y01);
1.237 + dst[2].set(SkScalarAve(x01, x12), SkScalarAve(y01, y12));
1.238 + dst[3].set(x12, y12);
1.239 + dst[4] = src[2];
1.240 +}
1.241 +
1.242 +/** Quad'(t) = At + B, where
1.243 + A = 2(a - 2b + c)
1.244 + B = 2(b - a)
1.245 + Solve for t, only if it fits between 0 < t < 1
1.246 +*/
1.247 +int SkFindQuadExtrema(SkScalar a, SkScalar b, SkScalar c, SkScalar tValue[1]) {
1.248 + /* At + B == 0
1.249 + t = -B / A
1.250 + */
1.251 + return valid_unit_divide(a - b, a - b - b + c, tValue);
1.252 +}
1.253 +
1.254 +static inline void flatten_double_quad_extrema(SkScalar coords[14]) {
1.255 + coords[2] = coords[6] = coords[4];
1.256 +}
1.257 +
1.258 +/* Returns 0 for 1 quad, and 1 for two quads, either way the answer is
1.259 + stored in dst[]. Guarantees that the 1/2 quads will be monotonic.
1.260 + */
1.261 +int SkChopQuadAtYExtrema(const SkPoint src[3], SkPoint dst[5]) {
1.262 + SkASSERT(src);
1.263 + SkASSERT(dst);
1.264 +
1.265 + SkScalar a = src[0].fY;
1.266 + SkScalar b = src[1].fY;
1.267 + SkScalar c = src[2].fY;
1.268 +
1.269 + if (is_not_monotonic(a, b, c)) {
1.270 + SkScalar tValue;
1.271 + if (valid_unit_divide(a - b, a - b - b + c, &tValue)) {
1.272 + SkChopQuadAt(src, dst, tValue);
1.273 + flatten_double_quad_extrema(&dst[0].fY);
1.274 + return 1;
1.275 + }
1.276 + // if we get here, we need to force dst to be monotonic, even though
1.277 + // we couldn't compute a unit_divide value (probably underflow).
1.278 + b = SkScalarAbs(a - b) < SkScalarAbs(b - c) ? a : c;
1.279 + }
1.280 + dst[0].set(src[0].fX, a);
1.281 + dst[1].set(src[1].fX, b);
1.282 + dst[2].set(src[2].fX, c);
1.283 + return 0;
1.284 +}
1.285 +
1.286 +/* Returns 0 for 1 quad, and 1 for two quads, either way the answer is
1.287 + stored in dst[]. Guarantees that the 1/2 quads will be monotonic.
1.288 + */
1.289 +int SkChopQuadAtXExtrema(const SkPoint src[3], SkPoint dst[5]) {
1.290 + SkASSERT(src);
1.291 + SkASSERT(dst);
1.292 +
1.293 + SkScalar a = src[0].fX;
1.294 + SkScalar b = src[1].fX;
1.295 + SkScalar c = src[2].fX;
1.296 +
1.297 + if (is_not_monotonic(a, b, c)) {
1.298 + SkScalar tValue;
1.299 + if (valid_unit_divide(a - b, a - b - b + c, &tValue)) {
1.300 + SkChopQuadAt(src, dst, tValue);
1.301 + flatten_double_quad_extrema(&dst[0].fX);
1.302 + return 1;
1.303 + }
1.304 + // if we get here, we need to force dst to be monotonic, even though
1.305 + // we couldn't compute a unit_divide value (probably underflow).
1.306 + b = SkScalarAbs(a - b) < SkScalarAbs(b - c) ? a : c;
1.307 + }
1.308 + dst[0].set(a, src[0].fY);
1.309 + dst[1].set(b, src[1].fY);
1.310 + dst[2].set(c, src[2].fY);
1.311 + return 0;
1.312 +}
1.313 +
1.314 +// F(t) = a (1 - t) ^ 2 + 2 b t (1 - t) + c t ^ 2
1.315 +// F'(t) = 2 (b - a) + 2 (a - 2b + c) t
1.316 +// F''(t) = 2 (a - 2b + c)
1.317 +//
1.318 +// A = 2 (b - a)
1.319 +// B = 2 (a - 2b + c)
1.320 +//
1.321 +// Maximum curvature for a quadratic means solving
1.322 +// Fx' Fx'' + Fy' Fy'' = 0
1.323 +//
1.324 +// t = - (Ax Bx + Ay By) / (Bx ^ 2 + By ^ 2)
1.325 +//
1.326 +SkScalar SkFindQuadMaxCurvature(const SkPoint src[3]) {
1.327 + SkScalar Ax = src[1].fX - src[0].fX;
1.328 + SkScalar Ay = src[1].fY - src[0].fY;
1.329 + SkScalar Bx = src[0].fX - src[1].fX - src[1].fX + src[2].fX;
1.330 + SkScalar By = src[0].fY - src[1].fY - src[1].fY + src[2].fY;
1.331 + SkScalar t = 0; // 0 means don't chop
1.332 +
1.333 + (void)valid_unit_divide(-(Ax * Bx + Ay * By), Bx * Bx + By * By, &t);
1.334 + return t;
1.335 +}
1.336 +
1.337 +int SkChopQuadAtMaxCurvature(const SkPoint src[3], SkPoint dst[5]) {
1.338 + SkScalar t = SkFindQuadMaxCurvature(src);
1.339 + if (t == 0) {
1.340 + memcpy(dst, src, 3 * sizeof(SkPoint));
1.341 + return 1;
1.342 + } else {
1.343 + SkChopQuadAt(src, dst, t);
1.344 + return 2;
1.345 + }
1.346 +}
1.347 +
1.348 +#define SK_ScalarTwoThirds (0.666666666f)
1.349 +
1.350 +void SkConvertQuadToCubic(const SkPoint src[3], SkPoint dst[4]) {
1.351 + const SkScalar scale = SK_ScalarTwoThirds;
1.352 + dst[0] = src[0];
1.353 + dst[1].set(src[0].fX + SkScalarMul(src[1].fX - src[0].fX, scale),
1.354 + src[0].fY + SkScalarMul(src[1].fY - src[0].fY, scale));
1.355 + dst[2].set(src[2].fX + SkScalarMul(src[1].fX - src[2].fX, scale),
1.356 + src[2].fY + SkScalarMul(src[1].fY - src[2].fY, scale));
1.357 + dst[3] = src[2];
1.358 +}
1.359 +
1.360 +//////////////////////////////////////////////////////////////////////////////
1.361 +///// CUBICS // CUBICS // CUBICS // CUBICS // CUBICS // CUBICS // CUBICS /////
1.362 +//////////////////////////////////////////////////////////////////////////////
1.363 +
1.364 +static void get_cubic_coeff(const SkScalar pt[], SkScalar coeff[4]) {
1.365 + coeff[0] = pt[6] + 3*(pt[2] - pt[4]) - pt[0];
1.366 + coeff[1] = 3*(pt[4] - pt[2] - pt[2] + pt[0]);
1.367 + coeff[2] = 3*(pt[2] - pt[0]);
1.368 + coeff[3] = pt[0];
1.369 +}
1.370 +
1.371 +void SkGetCubicCoeff(const SkPoint pts[4], SkScalar cx[4], SkScalar cy[4]) {
1.372 + SkASSERT(pts);
1.373 +
1.374 + if (cx) {
1.375 + get_cubic_coeff(&pts[0].fX, cx);
1.376 + }
1.377 + if (cy) {
1.378 + get_cubic_coeff(&pts[0].fY, cy);
1.379 + }
1.380 +}
1.381 +
1.382 +static SkScalar eval_cubic(const SkScalar src[], SkScalar t) {
1.383 + SkASSERT(src);
1.384 + SkASSERT(t >= 0 && t <= SK_Scalar1);
1.385 +
1.386 + if (t == 0) {
1.387 + return src[0];
1.388 + }
1.389 +
1.390 +#ifdef DIRECT_EVAL_OF_POLYNOMIALS
1.391 + SkScalar D = src[0];
1.392 + SkScalar A = src[6] + 3*(src[2] - src[4]) - D;
1.393 + SkScalar B = 3*(src[4] - src[2] - src[2] + D);
1.394 + SkScalar C = 3*(src[2] - D);
1.395 +
1.396 + return SkScalarMulAdd(SkScalarMulAdd(SkScalarMulAdd(A, t, B), t, C), t, D);
1.397 +#else
1.398 + SkScalar ab = SkScalarInterp(src[0], src[2], t);
1.399 + SkScalar bc = SkScalarInterp(src[2], src[4], t);
1.400 + SkScalar cd = SkScalarInterp(src[4], src[6], t);
1.401 + SkScalar abc = SkScalarInterp(ab, bc, t);
1.402 + SkScalar bcd = SkScalarInterp(bc, cd, t);
1.403 + return SkScalarInterp(abc, bcd, t);
1.404 +#endif
1.405 +}
1.406 +
1.407 +/** return At^2 + Bt + C
1.408 +*/
1.409 +static SkScalar eval_quadratic(SkScalar A, SkScalar B, SkScalar C, SkScalar t) {
1.410 + SkASSERT(t >= 0 && t <= SK_Scalar1);
1.411 +
1.412 + return SkScalarMulAdd(SkScalarMulAdd(A, t, B), t, C);
1.413 +}
1.414 +
1.415 +static SkScalar eval_cubic_derivative(const SkScalar src[], SkScalar t) {
1.416 + SkScalar A = src[6] + 3*(src[2] - src[4]) - src[0];
1.417 + SkScalar B = 2*(src[4] - 2 * src[2] + src[0]);
1.418 + SkScalar C = src[2] - src[0];
1.419 +
1.420 + return eval_quadratic(A, B, C, t);
1.421 +}
1.422 +
1.423 +static SkScalar eval_cubic_2ndDerivative(const SkScalar src[], SkScalar t) {
1.424 + SkScalar A = src[6] + 3*(src[2] - src[4]) - src[0];
1.425 + SkScalar B = src[4] - 2 * src[2] + src[0];
1.426 +
1.427 + return SkScalarMulAdd(A, t, B);
1.428 +}
1.429 +
1.430 +void SkEvalCubicAt(const SkPoint src[4], SkScalar t, SkPoint* loc,
1.431 + SkVector* tangent, SkVector* curvature) {
1.432 + SkASSERT(src);
1.433 + SkASSERT(t >= 0 && t <= SK_Scalar1);
1.434 +
1.435 + if (loc) {
1.436 + loc->set(eval_cubic(&src[0].fX, t), eval_cubic(&src[0].fY, t));
1.437 + }
1.438 + if (tangent) {
1.439 + tangent->set(eval_cubic_derivative(&src[0].fX, t),
1.440 + eval_cubic_derivative(&src[0].fY, t));
1.441 + }
1.442 + if (curvature) {
1.443 + curvature->set(eval_cubic_2ndDerivative(&src[0].fX, t),
1.444 + eval_cubic_2ndDerivative(&src[0].fY, t));
1.445 + }
1.446 +}
1.447 +
1.448 +/** Cubic'(t) = At^2 + Bt + C, where
1.449 + A = 3(-a + 3(b - c) + d)
1.450 + B = 6(a - 2b + c)
1.451 + C = 3(b - a)
1.452 + Solve for t, keeping only those that fit betwee 0 < t < 1
1.453 +*/
1.454 +int SkFindCubicExtrema(SkScalar a, SkScalar b, SkScalar c, SkScalar d,
1.455 + SkScalar tValues[2]) {
1.456 + // we divide A,B,C by 3 to simplify
1.457 + SkScalar A = d - a + 3*(b - c);
1.458 + SkScalar B = 2*(a - b - b + c);
1.459 + SkScalar C = b - a;
1.460 +
1.461 + return SkFindUnitQuadRoots(A, B, C, tValues);
1.462 +}
1.463 +
1.464 +static void interp_cubic_coords(const SkScalar* src, SkScalar* dst,
1.465 + SkScalar t) {
1.466 + SkScalar ab = SkScalarInterp(src[0], src[2], t);
1.467 + SkScalar bc = SkScalarInterp(src[2], src[4], t);
1.468 + SkScalar cd = SkScalarInterp(src[4], src[6], t);
1.469 + SkScalar abc = SkScalarInterp(ab, bc, t);
1.470 + SkScalar bcd = SkScalarInterp(bc, cd, t);
1.471 + SkScalar abcd = SkScalarInterp(abc, bcd, t);
1.472 +
1.473 + dst[0] = src[0];
1.474 + dst[2] = ab;
1.475 + dst[4] = abc;
1.476 + dst[6] = abcd;
1.477 + dst[8] = bcd;
1.478 + dst[10] = cd;
1.479 + dst[12] = src[6];
1.480 +}
1.481 +
1.482 +void SkChopCubicAt(const SkPoint src[4], SkPoint dst[7], SkScalar t) {
1.483 + SkASSERT(t > 0 && t < SK_Scalar1);
1.484 +
1.485 + interp_cubic_coords(&src[0].fX, &dst[0].fX, t);
1.486 + interp_cubic_coords(&src[0].fY, &dst[0].fY, t);
1.487 +}
1.488 +
1.489 +/* http://code.google.com/p/skia/issues/detail?id=32
1.490 +
1.491 + This test code would fail when we didn't check the return result of
1.492 + valid_unit_divide in SkChopCubicAt(... tValues[], int roots). The reason is
1.493 + that after the first chop, the parameters to valid_unit_divide are equal
1.494 + (thanks to finite float precision and rounding in the subtracts). Thus
1.495 + even though the 2nd tValue looks < 1.0, after we renormalize it, we end
1.496 + up with 1.0, hence the need to check and just return the last cubic as
1.497 + a degenerate clump of 4 points in the sampe place.
1.498 +
1.499 + static void test_cubic() {
1.500 + SkPoint src[4] = {
1.501 + { 556.25000, 523.03003 },
1.502 + { 556.23999, 522.96002 },
1.503 + { 556.21997, 522.89001 },
1.504 + { 556.21997, 522.82001 }
1.505 + };
1.506 + SkPoint dst[10];
1.507 + SkScalar tval[] = { 0.33333334f, 0.99999994f };
1.508 + SkChopCubicAt(src, dst, tval, 2);
1.509 + }
1.510 + */
1.511 +
1.512 +void SkChopCubicAt(const SkPoint src[4], SkPoint dst[],
1.513 + const SkScalar tValues[], int roots) {
1.514 +#ifdef SK_DEBUG
1.515 + {
1.516 + for (int i = 0; i < roots - 1; i++)
1.517 + {
1.518 + SkASSERT(is_unit_interval(tValues[i]));
1.519 + SkASSERT(is_unit_interval(tValues[i+1]));
1.520 + SkASSERT(tValues[i] < tValues[i+1]);
1.521 + }
1.522 + }
1.523 +#endif
1.524 +
1.525 + if (dst) {
1.526 + if (roots == 0) { // nothing to chop
1.527 + memcpy(dst, src, 4*sizeof(SkPoint));
1.528 + } else {
1.529 + SkScalar t = tValues[0];
1.530 + SkPoint tmp[4];
1.531 +
1.532 + for (int i = 0; i < roots; i++) {
1.533 + SkChopCubicAt(src, dst, t);
1.534 + if (i == roots - 1) {
1.535 + break;
1.536 + }
1.537 +
1.538 + dst += 3;
1.539 + // have src point to the remaining cubic (after the chop)
1.540 + memcpy(tmp, dst, 4 * sizeof(SkPoint));
1.541 + src = tmp;
1.542 +
1.543 + // watch out in case the renormalized t isn't in range
1.544 + if (!valid_unit_divide(tValues[i+1] - tValues[i],
1.545 + SK_Scalar1 - tValues[i], &t)) {
1.546 + // if we can't, just create a degenerate cubic
1.547 + dst[4] = dst[5] = dst[6] = src[3];
1.548 + break;
1.549 + }
1.550 + }
1.551 + }
1.552 + }
1.553 +}
1.554 +
1.555 +void SkChopCubicAtHalf(const SkPoint src[4], SkPoint dst[7]) {
1.556 + SkScalar x01 = SkScalarAve(src[0].fX, src[1].fX);
1.557 + SkScalar y01 = SkScalarAve(src[0].fY, src[1].fY);
1.558 + SkScalar x12 = SkScalarAve(src[1].fX, src[2].fX);
1.559 + SkScalar y12 = SkScalarAve(src[1].fY, src[2].fY);
1.560 + SkScalar x23 = SkScalarAve(src[2].fX, src[3].fX);
1.561 + SkScalar y23 = SkScalarAve(src[2].fY, src[3].fY);
1.562 +
1.563 + SkScalar x012 = SkScalarAve(x01, x12);
1.564 + SkScalar y012 = SkScalarAve(y01, y12);
1.565 + SkScalar x123 = SkScalarAve(x12, x23);
1.566 + SkScalar y123 = SkScalarAve(y12, y23);
1.567 +
1.568 + dst[0] = src[0];
1.569 + dst[1].set(x01, y01);
1.570 + dst[2].set(x012, y012);
1.571 + dst[3].set(SkScalarAve(x012, x123), SkScalarAve(y012, y123));
1.572 + dst[4].set(x123, y123);
1.573 + dst[5].set(x23, y23);
1.574 + dst[6] = src[3];
1.575 +}
1.576 +
1.577 +static void flatten_double_cubic_extrema(SkScalar coords[14]) {
1.578 + coords[4] = coords[8] = coords[6];
1.579 +}
1.580 +
1.581 +/** Given 4 points on a cubic bezier, chop it into 1, 2, 3 beziers such that
1.582 + the resulting beziers are monotonic in Y. This is called by the scan
1.583 + converter. Depending on what is returned, dst[] is treated as follows:
1.584 + 0 dst[0..3] is the original cubic
1.585 + 1 dst[0..3] and dst[3..6] are the two new cubics
1.586 + 2 dst[0..3], dst[3..6], dst[6..9] are the three new cubics
1.587 + If dst == null, it is ignored and only the count is returned.
1.588 +*/
1.589 +int SkChopCubicAtYExtrema(const SkPoint src[4], SkPoint dst[10]) {
1.590 + SkScalar tValues[2];
1.591 + int roots = SkFindCubicExtrema(src[0].fY, src[1].fY, src[2].fY,
1.592 + src[3].fY, tValues);
1.593 +
1.594 + SkChopCubicAt(src, dst, tValues, roots);
1.595 + if (dst && roots > 0) {
1.596 + // we do some cleanup to ensure our Y extrema are flat
1.597 + flatten_double_cubic_extrema(&dst[0].fY);
1.598 + if (roots == 2) {
1.599 + flatten_double_cubic_extrema(&dst[3].fY);
1.600 + }
1.601 + }
1.602 + return roots;
1.603 +}
1.604 +
1.605 +int SkChopCubicAtXExtrema(const SkPoint src[4], SkPoint dst[10]) {
1.606 + SkScalar tValues[2];
1.607 + int roots = SkFindCubicExtrema(src[0].fX, src[1].fX, src[2].fX,
1.608 + src[3].fX, tValues);
1.609 +
1.610 + SkChopCubicAt(src, dst, tValues, roots);
1.611 + if (dst && roots > 0) {
1.612 + // we do some cleanup to ensure our Y extrema are flat
1.613 + flatten_double_cubic_extrema(&dst[0].fX);
1.614 + if (roots == 2) {
1.615 + flatten_double_cubic_extrema(&dst[3].fX);
1.616 + }
1.617 + }
1.618 + return roots;
1.619 +}
1.620 +
1.621 +/** http://www.faculty.idc.ac.il/arik/quality/appendixA.html
1.622 +
1.623 + Inflection means that curvature is zero.
1.624 + Curvature is [F' x F''] / [F'^3]
1.625 + So we solve F'x X F''y - F'y X F''y == 0
1.626 + After some canceling of the cubic term, we get
1.627 + A = b - a
1.628 + B = c - 2b + a
1.629 + C = d - 3c + 3b - a
1.630 + (BxCy - ByCx)t^2 + (AxCy - AyCx)t + AxBy - AyBx == 0
1.631 +*/
1.632 +int SkFindCubicInflections(const SkPoint src[4], SkScalar tValues[]) {
1.633 + SkScalar Ax = src[1].fX - src[0].fX;
1.634 + SkScalar Ay = src[1].fY - src[0].fY;
1.635 + SkScalar Bx = src[2].fX - 2 * src[1].fX + src[0].fX;
1.636 + SkScalar By = src[2].fY - 2 * src[1].fY + src[0].fY;
1.637 + SkScalar Cx = src[3].fX + 3 * (src[1].fX - src[2].fX) - src[0].fX;
1.638 + SkScalar Cy = src[3].fY + 3 * (src[1].fY - src[2].fY) - src[0].fY;
1.639 +
1.640 + return SkFindUnitQuadRoots(Bx*Cy - By*Cx,
1.641 + Ax*Cy - Ay*Cx,
1.642 + Ax*By - Ay*Bx,
1.643 + tValues);
1.644 +}
1.645 +
1.646 +int SkChopCubicAtInflections(const SkPoint src[], SkPoint dst[10]) {
1.647 + SkScalar tValues[2];
1.648 + int count = SkFindCubicInflections(src, tValues);
1.649 +
1.650 + if (dst) {
1.651 + if (count == 0) {
1.652 + memcpy(dst, src, 4 * sizeof(SkPoint));
1.653 + } else {
1.654 + SkChopCubicAt(src, dst, tValues, count);
1.655 + }
1.656 + }
1.657 + return count + 1;
1.658 +}
1.659 +
1.660 +template <typename T> void bubble_sort(T array[], int count) {
1.661 + for (int i = count - 1; i > 0; --i)
1.662 + for (int j = i; j > 0; --j)
1.663 + if (array[j] < array[j-1])
1.664 + {
1.665 + T tmp(array[j]);
1.666 + array[j] = array[j-1];
1.667 + array[j-1] = tmp;
1.668 + }
1.669 +}
1.670 +
1.671 +/**
1.672 + * Given an array and count, remove all pair-wise duplicates from the array,
1.673 + * keeping the existing sorting, and return the new count
1.674 + */
1.675 +static int collaps_duplicates(SkScalar array[], int count) {
1.676 + for (int n = count; n > 1; --n) {
1.677 + if (array[0] == array[1]) {
1.678 + for (int i = 1; i < n; ++i) {
1.679 + array[i - 1] = array[i];
1.680 + }
1.681 + count -= 1;
1.682 + } else {
1.683 + array += 1;
1.684 + }
1.685 + }
1.686 + return count;
1.687 +}
1.688 +
1.689 +#ifdef SK_DEBUG
1.690 +
1.691 +#define TEST_COLLAPS_ENTRY(array) array, SK_ARRAY_COUNT(array)
1.692 +
1.693 +static void test_collaps_duplicates() {
1.694 + static bool gOnce;
1.695 + if (gOnce) { return; }
1.696 + gOnce = true;
1.697 + const SkScalar src0[] = { 0 };
1.698 + const SkScalar src1[] = { 0, 0 };
1.699 + const SkScalar src2[] = { 0, 1 };
1.700 + const SkScalar src3[] = { 0, 0, 0 };
1.701 + const SkScalar src4[] = { 0, 0, 1 };
1.702 + const SkScalar src5[] = { 0, 1, 1 };
1.703 + const SkScalar src6[] = { 0, 1, 2 };
1.704 + const struct {
1.705 + const SkScalar* fData;
1.706 + int fCount;
1.707 + int fCollapsedCount;
1.708 + } data[] = {
1.709 + { TEST_COLLAPS_ENTRY(src0), 1 },
1.710 + { TEST_COLLAPS_ENTRY(src1), 1 },
1.711 + { TEST_COLLAPS_ENTRY(src2), 2 },
1.712 + { TEST_COLLAPS_ENTRY(src3), 1 },
1.713 + { TEST_COLLAPS_ENTRY(src4), 2 },
1.714 + { TEST_COLLAPS_ENTRY(src5), 2 },
1.715 + { TEST_COLLAPS_ENTRY(src6), 3 },
1.716 + };
1.717 + for (size_t i = 0; i < SK_ARRAY_COUNT(data); ++i) {
1.718 + SkScalar dst[3];
1.719 + memcpy(dst, data[i].fData, data[i].fCount * sizeof(dst[0]));
1.720 + int count = collaps_duplicates(dst, data[i].fCount);
1.721 + SkASSERT(data[i].fCollapsedCount == count);
1.722 + for (int j = 1; j < count; ++j) {
1.723 + SkASSERT(dst[j-1] < dst[j]);
1.724 + }
1.725 + }
1.726 +}
1.727 +#endif
1.728 +
1.729 +static SkScalar SkScalarCubeRoot(SkScalar x) {
1.730 + return SkScalarPow(x, 0.3333333f);
1.731 +}
1.732 +
1.733 +/* Solve coeff(t) == 0, returning the number of roots that
1.734 + lie withing 0 < t < 1.
1.735 + coeff[0]t^3 + coeff[1]t^2 + coeff[2]t + coeff[3]
1.736 +
1.737 + Eliminates repeated roots (so that all tValues are distinct, and are always
1.738 + in increasing order.
1.739 +*/
1.740 +static int solve_cubic_poly(const SkScalar coeff[4], SkScalar tValues[3]) {
1.741 + if (SkScalarNearlyZero(coeff[0])) { // we're just a quadratic
1.742 + return SkFindUnitQuadRoots(coeff[1], coeff[2], coeff[3], tValues);
1.743 + }
1.744 +
1.745 + SkScalar a, b, c, Q, R;
1.746 +
1.747 + {
1.748 + SkASSERT(coeff[0] != 0);
1.749 +
1.750 + SkScalar inva = SkScalarInvert(coeff[0]);
1.751 + a = coeff[1] * inva;
1.752 + b = coeff[2] * inva;
1.753 + c = coeff[3] * inva;
1.754 + }
1.755 + Q = (a*a - b*3) / 9;
1.756 + R = (2*a*a*a - 9*a*b + 27*c) / 54;
1.757 +
1.758 + SkScalar Q3 = Q * Q * Q;
1.759 + SkScalar R2MinusQ3 = R * R - Q3;
1.760 + SkScalar adiv3 = a / 3;
1.761 +
1.762 + SkScalar* roots = tValues;
1.763 + SkScalar r;
1.764 +
1.765 + if (R2MinusQ3 < 0) { // we have 3 real roots
1.766 + SkScalar theta = SkScalarACos(R / SkScalarSqrt(Q3));
1.767 + SkScalar neg2RootQ = -2 * SkScalarSqrt(Q);
1.768 +
1.769 + r = neg2RootQ * SkScalarCos(theta/3) - adiv3;
1.770 + if (is_unit_interval(r)) {
1.771 + *roots++ = r;
1.772 + }
1.773 + r = neg2RootQ * SkScalarCos((theta + 2*SK_ScalarPI)/3) - adiv3;
1.774 + if (is_unit_interval(r)) {
1.775 + *roots++ = r;
1.776 + }
1.777 + r = neg2RootQ * SkScalarCos((theta - 2*SK_ScalarPI)/3) - adiv3;
1.778 + if (is_unit_interval(r)) {
1.779 + *roots++ = r;
1.780 + }
1.781 + SkDEBUGCODE(test_collaps_duplicates();)
1.782 +
1.783 + // now sort the roots
1.784 + int count = (int)(roots - tValues);
1.785 + SkASSERT((unsigned)count <= 3);
1.786 + bubble_sort(tValues, count);
1.787 + count = collaps_duplicates(tValues, count);
1.788 + roots = tValues + count; // so we compute the proper count below
1.789 + } else { // we have 1 real root
1.790 + SkScalar A = SkScalarAbs(R) + SkScalarSqrt(R2MinusQ3);
1.791 + A = SkScalarCubeRoot(A);
1.792 + if (R > 0) {
1.793 + A = -A;
1.794 + }
1.795 + if (A != 0) {
1.796 + A += Q / A;
1.797 + }
1.798 + r = A - adiv3;
1.799 + if (is_unit_interval(r)) {
1.800 + *roots++ = r;
1.801 + }
1.802 + }
1.803 +
1.804 + return (int)(roots - tValues);
1.805 +}
1.806 +
1.807 +/* Looking for F' dot F'' == 0
1.808 +
1.809 + A = b - a
1.810 + B = c - 2b + a
1.811 + C = d - 3c + 3b - a
1.812 +
1.813 + F' = 3Ct^2 + 6Bt + 3A
1.814 + F'' = 6Ct + 6B
1.815 +
1.816 + F' dot F'' -> CCt^3 + 3BCt^2 + (2BB + CA)t + AB
1.817 +*/
1.818 +static void formulate_F1DotF2(const SkScalar src[], SkScalar coeff[4]) {
1.819 + SkScalar a = src[2] - src[0];
1.820 + SkScalar b = src[4] - 2 * src[2] + src[0];
1.821 + SkScalar c = src[6] + 3 * (src[2] - src[4]) - src[0];
1.822 +
1.823 + coeff[0] = c * c;
1.824 + coeff[1] = 3 * b * c;
1.825 + coeff[2] = 2 * b * b + c * a;
1.826 + coeff[3] = a * b;
1.827 +}
1.828 +
1.829 +/* Looking for F' dot F'' == 0
1.830 +
1.831 + A = b - a
1.832 + B = c - 2b + a
1.833 + C = d - 3c + 3b - a
1.834 +
1.835 + F' = 3Ct^2 + 6Bt + 3A
1.836 + F'' = 6Ct + 6B
1.837 +
1.838 + F' dot F'' -> CCt^3 + 3BCt^2 + (2BB + CA)t + AB
1.839 +*/
1.840 +int SkFindCubicMaxCurvature(const SkPoint src[4], SkScalar tValues[3]) {
1.841 + SkScalar coeffX[4], coeffY[4];
1.842 + int i;
1.843 +
1.844 + formulate_F1DotF2(&src[0].fX, coeffX);
1.845 + formulate_F1DotF2(&src[0].fY, coeffY);
1.846 +
1.847 + for (i = 0; i < 4; i++) {
1.848 + coeffX[i] += coeffY[i];
1.849 + }
1.850 +
1.851 + SkScalar t[3];
1.852 + int count = solve_cubic_poly(coeffX, t);
1.853 + int maxCount = 0;
1.854 +
1.855 + // now remove extrema where the curvature is zero (mins)
1.856 + // !!!! need a test for this !!!!
1.857 + for (i = 0; i < count; i++) {
1.858 + // if (not_min_curvature())
1.859 + if (t[i] > 0 && t[i] < SK_Scalar1) {
1.860 + tValues[maxCount++] = t[i];
1.861 + }
1.862 + }
1.863 + return maxCount;
1.864 +}
1.865 +
1.866 +int SkChopCubicAtMaxCurvature(const SkPoint src[4], SkPoint dst[13],
1.867 + SkScalar tValues[3]) {
1.868 + SkScalar t_storage[3];
1.869 +
1.870 + if (tValues == NULL) {
1.871 + tValues = t_storage;
1.872 + }
1.873 +
1.874 + int count = SkFindCubicMaxCurvature(src, tValues);
1.875 +
1.876 + if (dst) {
1.877 + if (count == 0) {
1.878 + memcpy(dst, src, 4 * sizeof(SkPoint));
1.879 + } else {
1.880 + SkChopCubicAt(src, dst, tValues, count);
1.881 + }
1.882 + }
1.883 + return count + 1;
1.884 +}
1.885 +
1.886 +bool SkXRayCrossesMonotonicCubic(const SkXRay& pt, const SkPoint cubic[4],
1.887 + bool* ambiguous) {
1.888 + if (ambiguous) {
1.889 + *ambiguous = false;
1.890 + }
1.891 +
1.892 + // Find the minimum and maximum y of the extrema, which are the
1.893 + // first and last points since this cubic is monotonic
1.894 + SkScalar min_y = SkMinScalar(cubic[0].fY, cubic[3].fY);
1.895 + SkScalar max_y = SkMaxScalar(cubic[0].fY, cubic[3].fY);
1.896 +
1.897 + if (pt.fY == cubic[0].fY
1.898 + || pt.fY < min_y
1.899 + || pt.fY > max_y) {
1.900 + // The query line definitely does not cross the curve
1.901 + if (ambiguous) {
1.902 + *ambiguous = (pt.fY == cubic[0].fY);
1.903 + }
1.904 + return false;
1.905 + }
1.906 +
1.907 + bool pt_at_extremum = (pt.fY == cubic[3].fY);
1.908 +
1.909 + SkScalar min_x =
1.910 + SkMinScalar(
1.911 + SkMinScalar(
1.912 + SkMinScalar(cubic[0].fX, cubic[1].fX),
1.913 + cubic[2].fX),
1.914 + cubic[3].fX);
1.915 + if (pt.fX < min_x) {
1.916 + // The query line definitely crosses the curve
1.917 + if (ambiguous) {
1.918 + *ambiguous = pt_at_extremum;
1.919 + }
1.920 + return true;
1.921 + }
1.922 +
1.923 + SkScalar max_x =
1.924 + SkMaxScalar(
1.925 + SkMaxScalar(
1.926 + SkMaxScalar(cubic[0].fX, cubic[1].fX),
1.927 + cubic[2].fX),
1.928 + cubic[3].fX);
1.929 + if (pt.fX > max_x) {
1.930 + // The query line definitely does not cross the curve
1.931 + return false;
1.932 + }
1.933 +
1.934 + // Do a binary search to find the parameter value which makes y as
1.935 + // close as possible to the query point. See whether the query
1.936 + // line's origin is to the left of the associated x coordinate.
1.937 +
1.938 + // kMaxIter is chosen as the number of mantissa bits for a float,
1.939 + // since there's no way we are going to get more precision by
1.940 + // iterating more times than that.
1.941 + const int kMaxIter = 23;
1.942 + SkPoint eval;
1.943 + int iter = 0;
1.944 + SkScalar upper_t;
1.945 + SkScalar lower_t;
1.946 + // Need to invert direction of t parameter if cubic goes up
1.947 + // instead of down
1.948 + if (cubic[3].fY > cubic[0].fY) {
1.949 + upper_t = SK_Scalar1;
1.950 + lower_t = 0;
1.951 + } else {
1.952 + upper_t = 0;
1.953 + lower_t = SK_Scalar1;
1.954 + }
1.955 + do {
1.956 + SkScalar t = SkScalarAve(upper_t, lower_t);
1.957 + SkEvalCubicAt(cubic, t, &eval, NULL, NULL);
1.958 + if (pt.fY > eval.fY) {
1.959 + lower_t = t;
1.960 + } else {
1.961 + upper_t = t;
1.962 + }
1.963 + } while (++iter < kMaxIter
1.964 + && !SkScalarNearlyZero(eval.fY - pt.fY));
1.965 + if (pt.fX <= eval.fX) {
1.966 + if (ambiguous) {
1.967 + *ambiguous = pt_at_extremum;
1.968 + }
1.969 + return true;
1.970 + }
1.971 + return false;
1.972 +}
1.973 +
1.974 +int SkNumXRayCrossingsForCubic(const SkXRay& pt,
1.975 + const SkPoint cubic[4],
1.976 + bool* ambiguous) {
1.977 + int num_crossings = 0;
1.978 + SkPoint monotonic_cubics[10];
1.979 + int num_monotonic_cubics = SkChopCubicAtYExtrema(cubic, monotonic_cubics);
1.980 + if (ambiguous) {
1.981 + *ambiguous = false;
1.982 + }
1.983 + bool locally_ambiguous;
1.984 + if (SkXRayCrossesMonotonicCubic(pt,
1.985 + &monotonic_cubics[0],
1.986 + &locally_ambiguous))
1.987 + ++num_crossings;
1.988 + if (ambiguous) {
1.989 + *ambiguous |= locally_ambiguous;
1.990 + }
1.991 + if (num_monotonic_cubics > 0)
1.992 + if (SkXRayCrossesMonotonicCubic(pt,
1.993 + &monotonic_cubics[3],
1.994 + &locally_ambiguous))
1.995 + ++num_crossings;
1.996 + if (ambiguous) {
1.997 + *ambiguous |= locally_ambiguous;
1.998 + }
1.999 + if (num_monotonic_cubics > 1)
1.1000 + if (SkXRayCrossesMonotonicCubic(pt,
1.1001 + &monotonic_cubics[6],
1.1002 + &locally_ambiguous))
1.1003 + ++num_crossings;
1.1004 + if (ambiguous) {
1.1005 + *ambiguous |= locally_ambiguous;
1.1006 + }
1.1007 + return num_crossings;
1.1008 +}
1.1009 +
1.1010 +///////////////////////////////////////////////////////////////////////////////
1.1011 +
1.1012 +/* Find t value for quadratic [a, b, c] = d.
1.1013 + Return 0 if there is no solution within [0, 1)
1.1014 +*/
1.1015 +static SkScalar quad_solve(SkScalar a, SkScalar b, SkScalar c, SkScalar d) {
1.1016 + // At^2 + Bt + C = d
1.1017 + SkScalar A = a - 2 * b + c;
1.1018 + SkScalar B = 2 * (b - a);
1.1019 + SkScalar C = a - d;
1.1020 +
1.1021 + SkScalar roots[2];
1.1022 + int count = SkFindUnitQuadRoots(A, B, C, roots);
1.1023 +
1.1024 + SkASSERT(count <= 1);
1.1025 + return count == 1 ? roots[0] : 0;
1.1026 +}
1.1027 +
1.1028 +/* given a quad-curve and a point (x,y), chop the quad at that point and place
1.1029 + the new off-curve point and endpoint into 'dest'.
1.1030 + Should only return false if the computed pos is the start of the curve
1.1031 + (i.e. root == 0)
1.1032 +*/
1.1033 +static bool truncate_last_curve(const SkPoint quad[3], SkScalar x, SkScalar y,
1.1034 + SkPoint* dest) {
1.1035 + const SkScalar* base;
1.1036 + SkScalar value;
1.1037 +
1.1038 + if (SkScalarAbs(x) < SkScalarAbs(y)) {
1.1039 + base = &quad[0].fX;
1.1040 + value = x;
1.1041 + } else {
1.1042 + base = &quad[0].fY;
1.1043 + value = y;
1.1044 + }
1.1045 +
1.1046 + // note: this returns 0 if it thinks value is out of range, meaning the
1.1047 + // root might return something outside of [0, 1)
1.1048 + SkScalar t = quad_solve(base[0], base[2], base[4], value);
1.1049 +
1.1050 + if (t > 0) {
1.1051 + SkPoint tmp[5];
1.1052 + SkChopQuadAt(quad, tmp, t);
1.1053 + dest[0] = tmp[1];
1.1054 + dest[1].set(x, y);
1.1055 + return true;
1.1056 + } else {
1.1057 + /* t == 0 means either the value triggered a root outside of [0, 1)
1.1058 + For our purposes, we can ignore the <= 0 roots, but we want to
1.1059 + catch the >= 1 roots (which given our caller, will basically mean
1.1060 + a root of 1, give-or-take numerical instability). If we are in the
1.1061 + >= 1 case, return the existing offCurve point.
1.1062 +
1.1063 + The test below checks to see if we are close to the "end" of the
1.1064 + curve (near base[4]). Rather than specifying a tolerance, I just
1.1065 + check to see if value is on to the right/left of the middle point
1.1066 + (depending on the direction/sign of the end points).
1.1067 + */
1.1068 + if ((base[0] < base[4] && value > base[2]) ||
1.1069 + (base[0] > base[4] && value < base[2])) // should root have been 1
1.1070 + {
1.1071 + dest[0] = quad[1];
1.1072 + dest[1].set(x, y);
1.1073 + return true;
1.1074 + }
1.1075 + }
1.1076 + return false;
1.1077 +}
1.1078 +
1.1079 +static const SkPoint gQuadCirclePts[kSkBuildQuadArcStorage] = {
1.1080 +// The mid point of the quadratic arc approximation is half way between the two
1.1081 +// control points. The float epsilon adjustment moves the on curve point out by
1.1082 +// two bits, distributing the convex test error between the round rect
1.1083 +// approximation and the convex cross product sign equality test.
1.1084 +#define SK_MID_RRECT_OFFSET \
1.1085 + (SK_Scalar1 + SK_ScalarTanPIOver8 + FLT_EPSILON * 4) / 2
1.1086 + { SK_Scalar1, 0 },
1.1087 + { SK_Scalar1, SK_ScalarTanPIOver8 },
1.1088 + { SK_MID_RRECT_OFFSET, SK_MID_RRECT_OFFSET },
1.1089 + { SK_ScalarTanPIOver8, SK_Scalar1 },
1.1090 +
1.1091 + { 0, SK_Scalar1 },
1.1092 + { -SK_ScalarTanPIOver8, SK_Scalar1 },
1.1093 + { -SK_MID_RRECT_OFFSET, SK_MID_RRECT_OFFSET },
1.1094 + { -SK_Scalar1, SK_ScalarTanPIOver8 },
1.1095 +
1.1096 + { -SK_Scalar1, 0 },
1.1097 + { -SK_Scalar1, -SK_ScalarTanPIOver8 },
1.1098 + { -SK_MID_RRECT_OFFSET, -SK_MID_RRECT_OFFSET },
1.1099 + { -SK_ScalarTanPIOver8, -SK_Scalar1 },
1.1100 +
1.1101 + { 0, -SK_Scalar1 },
1.1102 + { SK_ScalarTanPIOver8, -SK_Scalar1 },
1.1103 + { SK_MID_RRECT_OFFSET, -SK_MID_RRECT_OFFSET },
1.1104 + { SK_Scalar1, -SK_ScalarTanPIOver8 },
1.1105 +
1.1106 + { SK_Scalar1, 0 }
1.1107 +#undef SK_MID_RRECT_OFFSET
1.1108 +};
1.1109 +
1.1110 +int SkBuildQuadArc(const SkVector& uStart, const SkVector& uStop,
1.1111 + SkRotationDirection dir, const SkMatrix* userMatrix,
1.1112 + SkPoint quadPoints[]) {
1.1113 + // rotate by x,y so that uStart is (1.0)
1.1114 + SkScalar x = SkPoint::DotProduct(uStart, uStop);
1.1115 + SkScalar y = SkPoint::CrossProduct(uStart, uStop);
1.1116 +
1.1117 + SkScalar absX = SkScalarAbs(x);
1.1118 + SkScalar absY = SkScalarAbs(y);
1.1119 +
1.1120 + int pointCount;
1.1121 +
1.1122 + // check for (effectively) coincident vectors
1.1123 + // this can happen if our angle is nearly 0 or nearly 180 (y == 0)
1.1124 + // ... we use the dot-prod to distinguish between 0 and 180 (x > 0)
1.1125 + if (absY <= SK_ScalarNearlyZero && x > 0 &&
1.1126 + ((y >= 0 && kCW_SkRotationDirection == dir) ||
1.1127 + (y <= 0 && kCCW_SkRotationDirection == dir))) {
1.1128 +
1.1129 + // just return the start-point
1.1130 + quadPoints[0].set(SK_Scalar1, 0);
1.1131 + pointCount = 1;
1.1132 + } else {
1.1133 + if (dir == kCCW_SkRotationDirection) {
1.1134 + y = -y;
1.1135 + }
1.1136 + // what octant (quadratic curve) is [xy] in?
1.1137 + int oct = 0;
1.1138 + bool sameSign = true;
1.1139 +
1.1140 + if (0 == y) {
1.1141 + oct = 4; // 180
1.1142 + SkASSERT(SkScalarAbs(x + SK_Scalar1) <= SK_ScalarNearlyZero);
1.1143 + } else if (0 == x) {
1.1144 + SkASSERT(absY - SK_Scalar1 <= SK_ScalarNearlyZero);
1.1145 + oct = y > 0 ? 2 : 6; // 90 : 270
1.1146 + } else {
1.1147 + if (y < 0) {
1.1148 + oct += 4;
1.1149 + }
1.1150 + if ((x < 0) != (y < 0)) {
1.1151 + oct += 2;
1.1152 + sameSign = false;
1.1153 + }
1.1154 + if ((absX < absY) == sameSign) {
1.1155 + oct += 1;
1.1156 + }
1.1157 + }
1.1158 +
1.1159 + int wholeCount = oct << 1;
1.1160 + memcpy(quadPoints, gQuadCirclePts, (wholeCount + 1) * sizeof(SkPoint));
1.1161 +
1.1162 + const SkPoint* arc = &gQuadCirclePts[wholeCount];
1.1163 + if (truncate_last_curve(arc, x, y, &quadPoints[wholeCount + 1])) {
1.1164 + wholeCount += 2;
1.1165 + }
1.1166 + pointCount = wholeCount + 1;
1.1167 + }
1.1168 +
1.1169 + // now handle counter-clockwise and the initial unitStart rotation
1.1170 + SkMatrix matrix;
1.1171 + matrix.setSinCos(uStart.fY, uStart.fX);
1.1172 + if (dir == kCCW_SkRotationDirection) {
1.1173 + matrix.preScale(SK_Scalar1, -SK_Scalar1);
1.1174 + }
1.1175 + if (userMatrix) {
1.1176 + matrix.postConcat(*userMatrix);
1.1177 + }
1.1178 + matrix.mapPoints(quadPoints, pointCount);
1.1179 + return pointCount;
1.1180 +}
1.1181 +
1.1182 +
1.1183 +///////////////////////////////////////////////////////////////////////////////
1.1184 +//
1.1185 +// NURB representation for conics. Helpful explanations at:
1.1186 +//
1.1187 +// http://citeseerx.ist.psu.edu/viewdoc/
1.1188 +// download?doi=10.1.1.44.5740&rep=rep1&type=ps
1.1189 +// and
1.1190 +// http://www.cs.mtu.edu/~shene/COURSES/cs3621/NOTES/spline/NURBS/RB-conics.html
1.1191 +//
1.1192 +// F = (A (1 - t)^2 + C t^2 + 2 B (1 - t) t w)
1.1193 +// ------------------------------------------
1.1194 +// ((1 - t)^2 + t^2 + 2 (1 - t) t w)
1.1195 +//
1.1196 +// = {t^2 (P0 + P2 - 2 P1 w), t (-2 P0 + 2 P1 w), P0}
1.1197 +// ------------------------------------------------
1.1198 +// {t^2 (2 - 2 w), t (-2 + 2 w), 1}
1.1199 +//
1.1200 +
1.1201 +static SkScalar conic_eval_pos(const SkScalar src[], SkScalar w, SkScalar t) {
1.1202 + SkASSERT(src);
1.1203 + SkASSERT(t >= 0 && t <= SK_Scalar1);
1.1204 +
1.1205 + SkScalar src2w = SkScalarMul(src[2], w);
1.1206 + SkScalar C = src[0];
1.1207 + SkScalar A = src[4] - 2 * src2w + C;
1.1208 + SkScalar B = 2 * (src2w - C);
1.1209 + SkScalar numer = SkScalarMulAdd(SkScalarMulAdd(A, t, B), t, C);
1.1210 +
1.1211 + B = 2 * (w - SK_Scalar1);
1.1212 + C = SK_Scalar1;
1.1213 + A = -B;
1.1214 + SkScalar denom = SkScalarMulAdd(SkScalarMulAdd(A, t, B), t, C);
1.1215 +
1.1216 + return SkScalarDiv(numer, denom);
1.1217 +}
1.1218 +
1.1219 +// F' = 2 (C t (1 + t (-1 + w)) - A (-1 + t) (t (-1 + w) - w) + B (1 - 2 t) w)
1.1220 +//
1.1221 +// t^2 : (2 P0 - 2 P2 - 2 P0 w + 2 P2 w)
1.1222 +// t^1 : (-2 P0 + 2 P2 + 4 P0 w - 4 P1 w)
1.1223 +// t^0 : -2 P0 w + 2 P1 w
1.1224 +//
1.1225 +// We disregard magnitude, so we can freely ignore the denominator of F', and
1.1226 +// divide the numerator by 2
1.1227 +//
1.1228 +// coeff[0] for t^2
1.1229 +// coeff[1] for t^1
1.1230 +// coeff[2] for t^0
1.1231 +//
1.1232 +static void conic_deriv_coeff(const SkScalar src[],
1.1233 + SkScalar w,
1.1234 + SkScalar coeff[3]) {
1.1235 + const SkScalar P20 = src[4] - src[0];
1.1236 + const SkScalar P10 = src[2] - src[0];
1.1237 + const SkScalar wP10 = w * P10;
1.1238 + coeff[0] = w * P20 - P20;
1.1239 + coeff[1] = P20 - 2 * wP10;
1.1240 + coeff[2] = wP10;
1.1241 +}
1.1242 +
1.1243 +static SkScalar conic_eval_tan(const SkScalar coord[], SkScalar w, SkScalar t) {
1.1244 + SkScalar coeff[3];
1.1245 + conic_deriv_coeff(coord, w, coeff);
1.1246 + return t * (t * coeff[0] + coeff[1]) + coeff[2];
1.1247 +}
1.1248 +
1.1249 +static bool conic_find_extrema(const SkScalar src[], SkScalar w, SkScalar* t) {
1.1250 + SkScalar coeff[3];
1.1251 + conic_deriv_coeff(src, w, coeff);
1.1252 +
1.1253 + SkScalar tValues[2];
1.1254 + int roots = SkFindUnitQuadRoots(coeff[0], coeff[1], coeff[2], tValues);
1.1255 + SkASSERT(0 == roots || 1 == roots);
1.1256 +
1.1257 + if (1 == roots) {
1.1258 + *t = tValues[0];
1.1259 + return true;
1.1260 + }
1.1261 + return false;
1.1262 +}
1.1263 +
1.1264 +struct SkP3D {
1.1265 + SkScalar fX, fY, fZ;
1.1266 +
1.1267 + void set(SkScalar x, SkScalar y, SkScalar z) {
1.1268 + fX = x; fY = y; fZ = z;
1.1269 + }
1.1270 +
1.1271 + void projectDown(SkPoint* dst) const {
1.1272 + dst->set(fX / fZ, fY / fZ);
1.1273 + }
1.1274 +};
1.1275 +
1.1276 +// We only interpolate one dimension at a time (the first, at +0, +3, +6).
1.1277 +static void p3d_interp(const SkScalar src[7], SkScalar dst[7], SkScalar t) {
1.1278 + SkScalar ab = SkScalarInterp(src[0], src[3], t);
1.1279 + SkScalar bc = SkScalarInterp(src[3], src[6], t);
1.1280 + dst[0] = ab;
1.1281 + dst[3] = SkScalarInterp(ab, bc, t);
1.1282 + dst[6] = bc;
1.1283 +}
1.1284 +
1.1285 +static void ratquad_mapTo3D(const SkPoint src[3], SkScalar w, SkP3D dst[]) {
1.1286 + dst[0].set(src[0].fX * 1, src[0].fY * 1, 1);
1.1287 + dst[1].set(src[1].fX * w, src[1].fY * w, w);
1.1288 + dst[2].set(src[2].fX * 1, src[2].fY * 1, 1);
1.1289 +}
1.1290 +
1.1291 +void SkConic::evalAt(SkScalar t, SkPoint* pt, SkVector* tangent) const {
1.1292 + SkASSERT(t >= 0 && t <= SK_Scalar1);
1.1293 +
1.1294 + if (pt) {
1.1295 + pt->set(conic_eval_pos(&fPts[0].fX, fW, t),
1.1296 + conic_eval_pos(&fPts[0].fY, fW, t));
1.1297 + }
1.1298 + if (tangent) {
1.1299 + tangent->set(conic_eval_tan(&fPts[0].fX, fW, t),
1.1300 + conic_eval_tan(&fPts[0].fY, fW, t));
1.1301 + }
1.1302 +}
1.1303 +
1.1304 +void SkConic::chopAt(SkScalar t, SkConic dst[2]) const {
1.1305 + SkP3D tmp[3], tmp2[3];
1.1306 +
1.1307 + ratquad_mapTo3D(fPts, fW, tmp);
1.1308 +
1.1309 + p3d_interp(&tmp[0].fX, &tmp2[0].fX, t);
1.1310 + p3d_interp(&tmp[0].fY, &tmp2[0].fY, t);
1.1311 + p3d_interp(&tmp[0].fZ, &tmp2[0].fZ, t);
1.1312 +
1.1313 + dst[0].fPts[0] = fPts[0];
1.1314 + tmp2[0].projectDown(&dst[0].fPts[1]);
1.1315 + tmp2[1].projectDown(&dst[0].fPts[2]); dst[1].fPts[0] = dst[0].fPts[2];
1.1316 + tmp2[2].projectDown(&dst[1].fPts[1]);
1.1317 + dst[1].fPts[2] = fPts[2];
1.1318 +
1.1319 + // to put in "standard form", where w0 and w2 are both 1, we compute the
1.1320 + // new w1 as sqrt(w1*w1/w0*w2)
1.1321 + // or
1.1322 + // w1 /= sqrt(w0*w2)
1.1323 + //
1.1324 + // However, in our case, we know that for dst[0]:
1.1325 + // w0 == 1, and for dst[1], w2 == 1
1.1326 + //
1.1327 + SkScalar root = SkScalarSqrt(tmp2[1].fZ);
1.1328 + dst[0].fW = tmp2[0].fZ / root;
1.1329 + dst[1].fW = tmp2[2].fZ / root;
1.1330 +}
1.1331 +
1.1332 +static SkScalar subdivide_w_value(SkScalar w) {
1.1333 + return SkScalarSqrt(SK_ScalarHalf + w * SK_ScalarHalf);
1.1334 +}
1.1335 +
1.1336 +void SkConic::chop(SkConic dst[2]) const {
1.1337 + SkScalar scale = SkScalarInvert(SK_Scalar1 + fW);
1.1338 + SkScalar p1x = fW * fPts[1].fX;
1.1339 + SkScalar p1y = fW * fPts[1].fY;
1.1340 + SkScalar mx = (fPts[0].fX + 2 * p1x + fPts[2].fX) * scale * SK_ScalarHalf;
1.1341 + SkScalar my = (fPts[0].fY + 2 * p1y + fPts[2].fY) * scale * SK_ScalarHalf;
1.1342 +
1.1343 + dst[0].fPts[0] = fPts[0];
1.1344 + dst[0].fPts[1].set((fPts[0].fX + p1x) * scale,
1.1345 + (fPts[0].fY + p1y) * scale);
1.1346 + dst[0].fPts[2].set(mx, my);
1.1347 +
1.1348 + dst[1].fPts[0].set(mx, my);
1.1349 + dst[1].fPts[1].set((p1x + fPts[2].fX) * scale,
1.1350 + (p1y + fPts[2].fY) * scale);
1.1351 + dst[1].fPts[2] = fPts[2];
1.1352 +
1.1353 + dst[0].fW = dst[1].fW = subdivide_w_value(fW);
1.1354 +}
1.1355 +
1.1356 +/*
1.1357 + * "High order approximation of conic sections by quadratic splines"
1.1358 + * by Michael Floater, 1993
1.1359 + */
1.1360 +#define AS_QUAD_ERROR_SETUP \
1.1361 + SkScalar a = fW - 1; \
1.1362 + SkScalar k = a / (4 * (2 + a)); \
1.1363 + SkScalar x = k * (fPts[0].fX - 2 * fPts[1].fX + fPts[2].fX); \
1.1364 + SkScalar y = k * (fPts[0].fY - 2 * fPts[1].fY + fPts[2].fY);
1.1365 +
1.1366 +void SkConic::computeAsQuadError(SkVector* err) const {
1.1367 + AS_QUAD_ERROR_SETUP
1.1368 + err->set(x, y);
1.1369 +}
1.1370 +
1.1371 +bool SkConic::asQuadTol(SkScalar tol) const {
1.1372 + AS_QUAD_ERROR_SETUP
1.1373 + return (x * x + y * y) <= tol * tol;
1.1374 +}
1.1375 +
1.1376 +int SkConic::computeQuadPOW2(SkScalar tol) const {
1.1377 + AS_QUAD_ERROR_SETUP
1.1378 + SkScalar error = SkScalarSqrt(x * x + y * y) - tol;
1.1379 +
1.1380 + if (error <= 0) {
1.1381 + return 0;
1.1382 + }
1.1383 + uint32_t ierr = (uint32_t)error;
1.1384 + return (34 - SkCLZ(ierr)) >> 1;
1.1385 +}
1.1386 +
1.1387 +static SkPoint* subdivide(const SkConic& src, SkPoint pts[], int level) {
1.1388 + SkASSERT(level >= 0);
1.1389 +
1.1390 + if (0 == level) {
1.1391 + memcpy(pts, &src.fPts[1], 2 * sizeof(SkPoint));
1.1392 + return pts + 2;
1.1393 + } else {
1.1394 + SkConic dst[2];
1.1395 + src.chop(dst);
1.1396 + --level;
1.1397 + pts = subdivide(dst[0], pts, level);
1.1398 + return subdivide(dst[1], pts, level);
1.1399 + }
1.1400 +}
1.1401 +
1.1402 +int SkConic::chopIntoQuadsPOW2(SkPoint pts[], int pow2) const {
1.1403 + SkASSERT(pow2 >= 0);
1.1404 + *pts = fPts[0];
1.1405 + SkDEBUGCODE(SkPoint* endPts =) subdivide(*this, pts + 1, pow2);
1.1406 + SkASSERT(endPts - pts == (2 * (1 << pow2) + 1));
1.1407 + return 1 << pow2;
1.1408 +}
1.1409 +
1.1410 +bool SkConic::findXExtrema(SkScalar* t) const {
1.1411 + return conic_find_extrema(&fPts[0].fX, fW, t);
1.1412 +}
1.1413 +
1.1414 +bool SkConic::findYExtrema(SkScalar* t) const {
1.1415 + return conic_find_extrema(&fPts[0].fY, fW, t);
1.1416 +}
1.1417 +
1.1418 +bool SkConic::chopAtXExtrema(SkConic dst[2]) const {
1.1419 + SkScalar t;
1.1420 + if (this->findXExtrema(&t)) {
1.1421 + this->chopAt(t, dst);
1.1422 + // now clean-up the middle, since we know t was meant to be at
1.1423 + // an X-extrema
1.1424 + SkScalar value = dst[0].fPts[2].fX;
1.1425 + dst[0].fPts[1].fX = value;
1.1426 + dst[1].fPts[0].fX = value;
1.1427 + dst[1].fPts[1].fX = value;
1.1428 + return true;
1.1429 + }
1.1430 + return false;
1.1431 +}
1.1432 +
1.1433 +bool SkConic::chopAtYExtrema(SkConic dst[2]) const {
1.1434 + SkScalar t;
1.1435 + if (this->findYExtrema(&t)) {
1.1436 + this->chopAt(t, dst);
1.1437 + // now clean-up the middle, since we know t was meant to be at
1.1438 + // an Y-extrema
1.1439 + SkScalar value = dst[0].fPts[2].fY;
1.1440 + dst[0].fPts[1].fY = value;
1.1441 + dst[1].fPts[0].fY = value;
1.1442 + dst[1].fPts[1].fY = value;
1.1443 + return true;
1.1444 + }
1.1445 + return false;
1.1446 +}
1.1447 +
1.1448 +void SkConic::computeTightBounds(SkRect* bounds) const {
1.1449 + SkPoint pts[4];
1.1450 + pts[0] = fPts[0];
1.1451 + pts[1] = fPts[2];
1.1452 + int count = 2;
1.1453 +
1.1454 + SkScalar t;
1.1455 + if (this->findXExtrema(&t)) {
1.1456 + this->evalAt(t, &pts[count++]);
1.1457 + }
1.1458 + if (this->findYExtrema(&t)) {
1.1459 + this->evalAt(t, &pts[count++]);
1.1460 + }
1.1461 + bounds->set(pts, count);
1.1462 +}
1.1463 +
1.1464 +void SkConic::computeFastBounds(SkRect* bounds) const {
1.1465 + bounds->set(fPts, 3);
1.1466 +}
1.1467 +
1.1468 +bool SkConic::findMaxCurvature(SkScalar* t) const {
1.1469 + // TODO: Implement me
1.1470 + return false;
1.1471 +}
