/* * Copyright 2006 The Android Open Source Project * * Use of this source code is governed by a BSD-style license that can be * found in the LICENSE file. */ #include "SkGeometry.h" #include "Sk64.h" #include "SkMatrix.h" bool SkXRayCrossesLine(const SkXRay& pt, const SkPoint pts[2], bool* ambiguous) { if (ambiguous) { *ambiguous = false; } // Determine quick discards. // Consider query line going exactly through point 0 to not // intersect, for symmetry with SkXRayCrossesMonotonicCubic. if (pt.fY == pts[0].fY) { if (ambiguous) { *ambiguous = true; } return false; } if (pt.fY < pts[0].fY && pt.fY < pts[1].fY) return false; if (pt.fY > pts[0].fY && pt.fY > pts[1].fY) return false; if (pt.fX > pts[0].fX && pt.fX > pts[1].fX) return false; // Determine degenerate cases if (SkScalarNearlyZero(pts[0].fY - pts[1].fY)) return false; if (SkScalarNearlyZero(pts[0].fX - pts[1].fX)) { // We've already determined the query point lies within the // vertical range of the line segment. if (pt.fX <= pts[0].fX) { if (ambiguous) { *ambiguous = (pt.fY == pts[1].fY); } return true; } return false; } // Ambiguity check if (pt.fY == pts[1].fY) { if (pt.fX <= pts[1].fX) { if (ambiguous) { *ambiguous = true; } return true; } return false; } // Full line segment evaluation SkScalar delta_y = pts[1].fY - pts[0].fY; SkScalar delta_x = pts[1].fX - pts[0].fX; SkScalar slope = SkScalarDiv(delta_y, delta_x); SkScalar b = pts[0].fY - SkScalarMul(slope, pts[0].fX); // Solve for x coordinate at y = pt.fY SkScalar x = SkScalarDiv(pt.fY - b, slope); return pt.fX <= x; } /** If defined, this makes eval_quad and eval_cubic do more setup (sometimes involving integer multiplies by 2 or 3, but fewer calls to SkScalarMul. May also introduce overflow of fixed when we compute our setup. */ #ifdef SK_SCALAR_IS_FIXED #define DIRECT_EVAL_OF_POLYNOMIALS #endif //////////////////////////////////////////////////////////////////////// #ifdef SK_SCALAR_IS_FIXED static int is_not_monotonic(int a, int b, int c, int d) { return (((a - b) | (b - c) | (c - d)) & ((b - a) | (c - b) | (d - c))) >> 31; } static int is_not_monotonic(int a, int b, int c) { return (((a - b) | (b - c)) & ((b - a) | (c - b))) >> 31; } #else static int is_not_monotonic(float a, float b, float c) { float ab = a - b; float bc = b - c; if (ab < 0) bc = -bc; return ab == 0 || bc < 0; } #endif //////////////////////////////////////////////////////////////////////// static bool is_unit_interval(SkScalar x) { return x > 0 && x < SK_Scalar1; } static int valid_unit_divide(SkScalar numer, SkScalar denom, SkScalar* ratio) { SkASSERT(ratio); if (numer < 0) { numer = -numer; denom = -denom; } if (denom == 0 || numer == 0 || numer >= denom) return 0; SkScalar r = SkScalarDiv(numer, denom); if (SkScalarIsNaN(r)) { return 0; } SkASSERT(r >= 0 && r < SK_Scalar1); if (r == 0) // catch underflow if numer <<<< denom return 0; *ratio = r; return 1; } /** From Numerical Recipes in C. Q = -1/2 (B + sign(B) sqrt[B*B - 4*A*C]) x1 = Q / A x2 = C / Q */ int SkFindUnitQuadRoots(SkScalar A, SkScalar B, SkScalar C, SkScalar roots[2]) { SkASSERT(roots); if (A == 0) return valid_unit_divide(-C, B, roots); SkScalar* r = roots; #ifdef SK_SCALAR_IS_FLOAT float R = B*B - 4*A*C; if (R < 0 || SkScalarIsNaN(R)) { // complex roots return 0; } R = sk_float_sqrt(R); #else Sk64 RR, tmp; RR.setMul(B,B); tmp.setMul(A,C); tmp.shiftLeft(2); RR.sub(tmp); if (RR.isNeg()) return 0; SkFixed R = RR.getSqrt(); #endif SkScalar Q = (B < 0) ? -(B-R)/2 : -(B+R)/2; r += valid_unit_divide(Q, A, r); r += valid_unit_divide(C, Q, r); if (r - roots == 2) { if (roots[0] > roots[1]) SkTSwap(roots[0], roots[1]); else if (roots[0] == roots[1]) // nearly-equal? r -= 1; // skip the double root } return (int)(r - roots); } #ifdef SK_SCALAR_IS_FIXED /** Trim A/B/C down so that they are all <= 32bits and then call SkFindUnitQuadRoots() */ static int Sk64FindFixedQuadRoots(const Sk64& A, const Sk64& B, const Sk64& C, SkFixed roots[2]) { int na = A.shiftToMake32(); int nb = B.shiftToMake32(); int nc = C.shiftToMake32(); int shift = SkMax32(na, SkMax32(nb, nc)); SkASSERT(shift >= 0); return SkFindUnitQuadRoots(A.getShiftRight(shift), B.getShiftRight(shift), C.getShiftRight(shift), roots); } #endif ///////////////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////////////// static SkScalar eval_quad(const SkScalar src[], SkScalar t) { SkASSERT(src); SkASSERT(t >= 0 && t <= SK_Scalar1); #ifdef DIRECT_EVAL_OF_POLYNOMIALS SkScalar C = src[0]; SkScalar A = src[4] - 2 * src[2] + C; SkScalar B = 2 * (src[2] - C); return SkScalarMulAdd(SkScalarMulAdd(A, t, B), t, C); #else SkScalar ab = SkScalarInterp(src[0], src[2], t); SkScalar bc = SkScalarInterp(src[2], src[4], t); return SkScalarInterp(ab, bc, t); #endif } static SkScalar eval_quad_derivative(const SkScalar src[], SkScalar t) { SkScalar A = src[4] - 2 * src[2] + src[0]; SkScalar B = src[2] - src[0]; return 2 * SkScalarMulAdd(A, t, B); } static SkScalar eval_quad_derivative_at_half(const SkScalar src[]) { SkScalar A = src[4] - 2 * src[2] + src[0]; SkScalar B = src[2] - src[0]; return A + 2 * B; } void SkEvalQuadAt(const SkPoint src[3], SkScalar t, SkPoint* pt, SkVector* tangent) { SkASSERT(src); SkASSERT(t >= 0 && t <= SK_Scalar1); if (pt) pt->set(eval_quad(&src[0].fX, t), eval_quad(&src[0].fY, t)); if (tangent) tangent->set(eval_quad_derivative(&src[0].fX, t), eval_quad_derivative(&src[0].fY, t)); } void SkEvalQuadAtHalf(const SkPoint src[3], SkPoint* pt, SkVector* tangent) { SkASSERT(src); if (pt) { SkScalar x01 = SkScalarAve(src[0].fX, src[1].fX); SkScalar y01 = SkScalarAve(src[0].fY, src[1].fY); SkScalar x12 = SkScalarAve(src[1].fX, src[2].fX); SkScalar y12 = SkScalarAve(src[1].fY, src[2].fY); pt->set(SkScalarAve(x01, x12), SkScalarAve(y01, y12)); } if (tangent) tangent->set(eval_quad_derivative_at_half(&src[0].fX), eval_quad_derivative_at_half(&src[0].fY)); } static void interp_quad_coords(const SkScalar* src, SkScalar* dst, SkScalar t) { SkScalar ab = SkScalarInterp(src[0], src[2], t); SkScalar bc = SkScalarInterp(src[2], src[4], t); dst[0] = src[0]; dst[2] = ab; dst[4] = SkScalarInterp(ab, bc, t); dst[6] = bc; dst[8] = src[4]; } void SkChopQuadAt(const SkPoint src[3], SkPoint dst[5], SkScalar t) { SkASSERT(t > 0 && t < SK_Scalar1); interp_quad_coords(&src[0].fX, &dst[0].fX, t); interp_quad_coords(&src[0].fY, &dst[0].fY, t); } void SkChopQuadAtHalf(const SkPoint src[3], SkPoint dst[5]) { SkScalar x01 = SkScalarAve(src[0].fX, src[1].fX); SkScalar y01 = SkScalarAve(src[0].fY, src[1].fY); SkScalar x12 = SkScalarAve(src[1].fX, src[2].fX); SkScalar y12 = SkScalarAve(src[1].fY, src[2].fY); dst[0] = src[0]; dst[1].set(x01, y01); dst[2].set(SkScalarAve(x01, x12), SkScalarAve(y01, y12)); dst[3].set(x12, y12); dst[4] = src[2]; } /** Quad'(t) = At + B, where A = 2(a - 2b + c) B = 2(b - a) Solve for t, only if it fits between 0 < t < 1 */ int SkFindQuadExtrema(SkScalar a, SkScalar b, SkScalar c, SkScalar tValue[1]) { /* At + B == 0 t = -B / A */ #ifdef SK_SCALAR_IS_FIXED return is_not_monotonic(a, b, c) && valid_unit_divide(a - b, a - b - b + c, tValue); #else return valid_unit_divide(a - b, a - b - b + c, tValue); #endif } static inline void flatten_double_quad_extrema(SkScalar coords[14]) { coords[2] = coords[6] = coords[4]; } /* Returns 0 for 1 quad, and 1 for two quads, either way the answer is stored in dst[]. Guarantees that the 1/2 quads will be monotonic. */ int SkChopQuadAtYExtrema(const SkPoint src[3], SkPoint dst[5]) { SkASSERT(src); SkASSERT(dst); #if 0 static bool once = true; if (once) { once = false; SkPoint s[3] = { 0, 26398, 0, 26331, 0, 20621428 }; SkPoint d[6]; int n = SkChopQuadAtYExtrema(s, d); 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); } #endif SkScalar a = src[0].fY; SkScalar b = src[1].fY; SkScalar c = src[2].fY; if (is_not_monotonic(a, b, c)) { SkScalar tValue; if (valid_unit_divide(a - b, a - b - b + c, &tValue)) { SkChopQuadAt(src, dst, tValue); flatten_double_quad_extrema(&dst[0].fY); return 1; } // if we get here, we need to force dst to be monotonic, even though // we couldn't compute a unit_divide value (probably underflow). b = SkScalarAbs(a - b) < SkScalarAbs(b - c) ? a : c; } dst[0].set(src[0].fX, a); dst[1].set(src[1].fX, b); dst[2].set(src[2].fX, c); return 0; } /* Returns 0 for 1 quad, and 1 for two quads, either way the answer is stored in dst[]. Guarantees that the 1/2 quads will be monotonic. */ int SkChopQuadAtXExtrema(const SkPoint src[3], SkPoint dst[5]) { SkASSERT(src); SkASSERT(dst); SkScalar a = src[0].fX; SkScalar b = src[1].fX; SkScalar c = src[2].fX; if (is_not_monotonic(a, b, c)) { SkScalar tValue; if (valid_unit_divide(a - b, a - b - b + c, &tValue)) { SkChopQuadAt(src, dst, tValue); flatten_double_quad_extrema(&dst[0].fX); return 1; } // if we get here, we need to force dst to be monotonic, even though // we couldn't compute a unit_divide value (probably underflow). b = SkScalarAbs(a - b) < SkScalarAbs(b - c) ? a : c; } dst[0].set(a, src[0].fY); dst[1].set(b, src[1].fY); dst[2].set(c, src[2].fY); return 0; } // F(t) = a (1 - t) ^ 2 + 2 b t (1 - t) + c t ^ 2 // F'(t) = 2 (b - a) + 2 (a - 2b + c) t // F''(t) = 2 (a - 2b + c) // // A = 2 (b - a) // B = 2 (a - 2b + c) // // Maximum curvature for a quadratic means solving // Fx' Fx'' + Fy' Fy'' = 0 // // t = - (Ax Bx + Ay By) / (Bx ^ 2 + By ^ 2) // int SkChopQuadAtMaxCurvature(const SkPoint src[3], SkPoint dst[5]) { SkScalar Ax = src[1].fX - src[0].fX; SkScalar Ay = src[1].fY - src[0].fY; SkScalar Bx = src[0].fX - src[1].fX - src[1].fX + src[2].fX; SkScalar By = src[0].fY - src[1].fY - src[1].fY + src[2].fY; SkScalar t = 0; // 0 means don't chop #ifdef SK_SCALAR_IS_FLOAT (void)valid_unit_divide(-(Ax * Bx + Ay * By), Bx * Bx + By * By, &t); #else // !!! should I use SkFloat here? seems like it Sk64 numer, denom, tmp; numer.setMul(Ax, -Bx); tmp.setMul(Ay, -By); numer.add(tmp); if (numer.isPos()) // do nothing if numer <= 0 { denom.setMul(Bx, Bx); tmp.setMul(By, By); denom.add(tmp); SkASSERT(!denom.isNeg()); if (numer < denom) { t = numer.getFixedDiv(denom); SkASSERT(t >= 0 && t <= SK_Fixed1); // assert that we're numerically stable (ha!) if ((unsigned)t >= SK_Fixed1) // runtime check for numerical stability t = 0; // ignore the chop } } #endif if (t == 0) { memcpy(dst, src, 3 * sizeof(SkPoint)); return 1; } else { SkChopQuadAt(src, dst, t); return 2; } } #ifdef SK_SCALAR_IS_FLOAT #define SK_ScalarTwoThirds (0.666666666f) #else #define SK_ScalarTwoThirds ((SkFixed)(43691)) #endif void SkConvertQuadToCubic(const SkPoint src[3], SkPoint dst[4]) { const SkScalar scale = SK_ScalarTwoThirds; dst[0] = src[0]; dst[1].set(src[0].fX + SkScalarMul(src[1].fX - src[0].fX, scale), src[0].fY + SkScalarMul(src[1].fY - src[0].fY, scale)); dst[2].set(src[2].fX + SkScalarMul(src[1].fX - src[2].fX, scale), src[2].fY + SkScalarMul(src[1].fY - src[2].fY, scale)); dst[3] = src[2]; } //////////////////////////////////////////////////////////////////////////////////////// ///// CUBICS // CUBICS // CUBICS // CUBICS // CUBICS // CUBICS // CUBICS // CUBICS ///// //////////////////////////////////////////////////////////////////////////////////////// static void get_cubic_coeff(const SkScalar pt[], SkScalar coeff[4]) { coeff[0] = pt[6] + 3*(pt[2] - pt[4]) - pt[0]; coeff[1] = 3*(pt[4] - pt[2] - pt[2] + pt[0]); coeff[2] = 3*(pt[2] - pt[0]); coeff[3] = pt[0]; } void SkGetCubicCoeff(const SkPoint pts[4], SkScalar cx[4], SkScalar cy[4]) { SkASSERT(pts); if (cx) get_cubic_coeff(&pts[0].fX, cx); if (cy) get_cubic_coeff(&pts[0].fY, cy); } static SkScalar eval_cubic(const SkScalar src[], SkScalar t) { SkASSERT(src); SkASSERT(t >= 0 && t <= SK_Scalar1); if (t == 0) return src[0]; #ifdef DIRECT_EVAL_OF_POLYNOMIALS SkScalar D = src[0]; SkScalar A = src[6] + 3*(src[2] - src[4]) - D; SkScalar B = 3*(src[4] - src[2] - src[2] + D); SkScalar C = 3*(src[2] - D); return SkScalarMulAdd(SkScalarMulAdd(SkScalarMulAdd(A, t, B), t, C), t, D); #else SkScalar ab = SkScalarInterp(src[0], src[2], t); SkScalar bc = SkScalarInterp(src[2], src[4], t); SkScalar cd = SkScalarInterp(src[4], src[6], t); SkScalar abc = SkScalarInterp(ab, bc, t); SkScalar bcd = SkScalarInterp(bc, cd, t); return SkScalarInterp(abc, bcd, t); #endif } /** return At^2 + Bt + C */ static SkScalar eval_quadratic(SkScalar A, SkScalar B, SkScalar C, SkScalar t) { SkASSERT(t >= 0 && t <= SK_Scalar1); return SkScalarMulAdd(SkScalarMulAdd(A, t, B), t, C); } static SkScalar eval_cubic_derivative(const SkScalar src[], SkScalar t) { SkScalar A = src[6] + 3*(src[2] - src[4]) - src[0]; SkScalar B = 2*(src[4] - 2 * src[2] + src[0]); SkScalar C = src[2] - src[0]; return eval_quadratic(A, B, C, t); } static SkScalar eval_cubic_2ndDerivative(const SkScalar src[], SkScalar t) { SkScalar A = src[6] + 3*(src[2] - src[4]) - src[0]; SkScalar B = src[4] - 2 * src[2] + src[0]; return SkScalarMulAdd(A, t, B); } void SkEvalCubicAt(const SkPoint src[4], SkScalar t, SkPoint* loc, SkVector* tangent, SkVector* curvature) { SkASSERT(src); SkASSERT(t >= 0 && t <= SK_Scalar1); if (loc) loc->set(eval_cubic(&src[0].fX, t), eval_cubic(&src[0].fY, t)); if (tangent) tangent->set(eval_cubic_derivative(&src[0].fX, t), eval_cubic_derivative(&src[0].fY, t)); if (curvature) curvature->set(eval_cubic_2ndDerivative(&src[0].fX, t), eval_cubic_2ndDerivative(&src[0].fY, t)); } /** Cubic'(t) = At^2 + Bt + C, where A = 3(-a + 3(b - c) + d) B = 6(a - 2b + c) C = 3(b - a) Solve for t, keeping only those that fit betwee 0 < t < 1 */ int SkFindCubicExtrema(SkScalar a, SkScalar b, SkScalar c, SkScalar d, SkScalar tValues[2]) { #ifdef SK_SCALAR_IS_FIXED if (!is_not_monotonic(a, b, c, d)) return 0; #endif // we divide A,B,C by 3 to simplify SkScalar A = d - a + 3*(b - c); SkScalar B = 2*(a - b - b + c); SkScalar C = b - a; return SkFindUnitQuadRoots(A, B, C, tValues); } static void interp_cubic_coords(const SkScalar* src, SkScalar* dst, SkScalar t) { SkScalar ab = SkScalarInterp(src[0], src[2], t); SkScalar bc = SkScalarInterp(src[2], src[4], t); SkScalar cd = SkScalarInterp(src[4], src[6], t); SkScalar abc = SkScalarInterp(ab, bc, t); SkScalar bcd = SkScalarInterp(bc, cd, t); SkScalar abcd = SkScalarInterp(abc, bcd, t); dst[0] = src[0]; dst[2] = ab; dst[4] = abc; dst[6] = abcd; dst[8] = bcd; dst[10] = cd; dst[12] = src[6]; } void SkChopCubicAt(const SkPoint src[4], SkPoint dst[7], SkScalar t) { SkASSERT(t > 0 && t < SK_Scalar1); interp_cubic_coords(&src[0].fX, &dst[0].fX, t); interp_cubic_coords(&src[0].fY, &dst[0].fY, t); } /* http://code.google.com/p/skia/issues/detail?id=32 This test code would fail when we didn't check the return result of valid_unit_divide in SkChopCubicAt(... tValues[], int roots). The reason is that after the first chop, the parameters to valid_unit_divide are equal (thanks to finite float precision and rounding in the subtracts). Thus even though the 2nd tValue looks < 1.0, after we renormalize it, we end up with 1.0, hence the need to check and just return the last cubic as a degenerate clump of 4 points in the sampe place. static void test_cubic() { SkPoint src[4] = { { 556.25000, 523.03003 }, { 556.23999, 522.96002 }, { 556.21997, 522.89001 }, { 556.21997, 522.82001 } }; SkPoint dst[10]; SkScalar tval[] = { 0.33333334f, 0.99999994f }; SkChopCubicAt(src, dst, tval, 2); } */ void SkChopCubicAt(const SkPoint src[4], SkPoint dst[], const SkScalar tValues[], int roots) { #ifdef SK_DEBUG { for (int i = 0; i < roots - 1; i++) { SkASSERT(is_unit_interval(tValues[i])); SkASSERT(is_unit_interval(tValues[i+1])); SkASSERT(tValues[i] < tValues[i+1]); } } #endif if (dst) { if (roots == 0) // nothing to chop memcpy(dst, src, 4*sizeof(SkPoint)); else { SkScalar t = tValues[0]; SkPoint tmp[4]; for (int i = 0; i < roots; i++) { SkChopCubicAt(src, dst, t); if (i == roots - 1) break; dst += 3; // have src point to the remaining cubic (after the chop) memcpy(tmp, dst, 4 * sizeof(SkPoint)); src = tmp; // watch out in case the renormalized t isn't in range if (!valid_unit_divide(tValues[i+1] - tValues[i], SK_Scalar1 - tValues[i], &t)) { // if we can't, just create a degenerate cubic dst[4] = dst[5] = dst[6] = src[3]; break; } } } } } void SkChopCubicAtHalf(const SkPoint src[4], SkPoint dst[7]) { SkScalar x01 = SkScalarAve(src[0].fX, src[1].fX); SkScalar y01 = SkScalarAve(src[0].fY, src[1].fY); SkScalar x12 = SkScalarAve(src[1].fX, src[2].fX); SkScalar y12 = SkScalarAve(src[1].fY, src[2].fY); SkScalar x23 = SkScalarAve(src[2].fX, src[3].fX); SkScalar y23 = SkScalarAve(src[2].fY, src[3].fY); SkScalar x012 = SkScalarAve(x01, x12); SkScalar y012 = SkScalarAve(y01, y12); SkScalar x123 = SkScalarAve(x12, x23); SkScalar y123 = SkScalarAve(y12, y23); dst[0] = src[0]; dst[1].set(x01, y01); dst[2].set(x012, y012); dst[3].set(SkScalarAve(x012, x123), SkScalarAve(y012, y123)); dst[4].set(x123, y123); dst[5].set(x23, y23); dst[6] = src[3]; } static void flatten_double_cubic_extrema(SkScalar coords[14]) { coords[4] = coords[8] = coords[6]; } /** Given 4 points on a cubic bezier, chop it into 1, 2, 3 beziers such that the resulting beziers are monotonic in Y. This is called by the scan converter. Depending on what is returned, dst[] is treated as follows 0 dst[0..3] is the original cubic 1 dst[0..3] and dst[3..6] are the two new cubics 2 dst[0..3], dst[3..6], dst[6..9] are the three new cubics If dst == null, it is ignored and only the count is returned. */ int SkChopCubicAtYExtrema(const SkPoint src[4], SkPoint dst[10]) { SkScalar tValues[2]; int roots = SkFindCubicExtrema(src[0].fY, src[1].fY, src[2].fY, src[3].fY, tValues); SkChopCubicAt(src, dst, tValues, roots); if (dst && roots > 0) { // we do some cleanup to ensure our Y extrema are flat flatten_double_cubic_extrema(&dst[0].fY); if (roots == 2) { flatten_double_cubic_extrema(&dst[3].fY); } } return roots; } int SkChopCubicAtXExtrema(const SkPoint src[4], SkPoint dst[10]) { SkScalar tValues[2]; int roots = SkFindCubicExtrema(src[0].fX, src[1].fX, src[2].fX, src[3].fX, tValues); SkChopCubicAt(src, dst, tValues, roots); if (dst && roots > 0) { // we do some cleanup to ensure our Y extrema are flat flatten_double_cubic_extrema(&dst[0].fX); if (roots == 2) { flatten_double_cubic_extrema(&dst[3].fX); } } return roots; } /** http://www.faculty.idc.ac.il/arik/quality/appendixA.html Inflection means that curvature is zero. Curvature is [F' x F''] / [F'^3] So we solve F'x X F''y - F'y X F''y == 0 After some canceling of the cubic term, we get A = b - a B = c - 2b + a C = d - 3c + 3b - a (BxCy - ByCx)t^2 + (AxCy - AyCx)t + AxBy - AyBx == 0 */ int SkFindCubicInflections(const SkPoint src[4], SkScalar tValues[]) { SkScalar Ax = src[1].fX - src[0].fX; SkScalar Ay = src[1].fY - src[0].fY; SkScalar Bx = src[2].fX - 2 * src[1].fX + src[0].fX; SkScalar By = src[2].fY - 2 * src[1].fY + src[0].fY; SkScalar Cx = src[3].fX + 3 * (src[1].fX - src[2].fX) - src[0].fX; SkScalar Cy = src[3].fY + 3 * (src[1].fY - src[2].fY) - src[0].fY; int count; #ifdef SK_SCALAR_IS_FLOAT count = SkFindUnitQuadRoots(Bx*Cy - By*Cx, Ax*Cy - Ay*Cx, Ax*By - Ay*Bx, tValues); #else Sk64 A, B, C, tmp; A.setMul(Bx, Cy); tmp.setMul(By, Cx); A.sub(tmp); B.setMul(Ax, Cy); tmp.setMul(Ay, Cx); B.sub(tmp); C.setMul(Ax, By); tmp.setMul(Ay, Bx); C.sub(tmp); count = Sk64FindFixedQuadRoots(A, B, C, tValues); #endif return count; } int SkChopCubicAtInflections(const SkPoint src[], SkPoint dst[10]) { SkScalar tValues[2]; int count = SkFindCubicInflections(src, tValues); if (dst) { if (count == 0) memcpy(dst, src, 4 * sizeof(SkPoint)); else SkChopCubicAt(src, dst, tValues, count); } return count + 1; } template void bubble_sort(T array[], int count) { for (int i = count - 1; i > 0; --i) for (int j = i; j > 0; --j) if (array[j] < array[j-1]) { T tmp(array[j]); array[j] = array[j-1]; array[j-1] = tmp; } } #include "SkFP.h" // newton refinement #if 0 static SkScalar refine_cubic_root(const SkFP coeff[4], SkScalar root) { // x1 = x0 - f(t) / f'(t) SkFP T = SkScalarToFloat(root); SkFP N, D; // f' = 3*coeff[0]*T^2 + 2*coeff[1]*T + coeff[2] D = SkFPMul(SkFPMul(coeff[0], SkFPMul(T,T)), 3); D = SkFPAdd(D, SkFPMulInt(SkFPMul(coeff[1], T), 2)); D = SkFPAdd(D, coeff[2]); if (D == 0) return root; // f = coeff[0]*T^3 + coeff[1]*T^2 + coeff[2]*T + coeff[3] N = SkFPMul(SkFPMul(SkFPMul(T, T), T), coeff[0]); N = SkFPAdd(N, SkFPMul(SkFPMul(T, T), coeff[1])); N = SkFPAdd(N, SkFPMul(T, coeff[2])); N = SkFPAdd(N, coeff[3]); if (N) { SkScalar delta = SkFPToScalar(SkFPDiv(N, D)); if (delta) root -= delta; } return root; } #endif /** * Given an array and count, remove all pair-wise duplicates from the array, * keeping the existing sorting, and return the new count */ static int collaps_duplicates(float array[], int count) { for (int n = count; n > 1; --n) { if (array[0] == array[1]) { for (int i = 1; i < n; ++i) { array[i - 1] = array[i]; } count -= 1; } else { array += 1; } } return count; } #ifdef SK_DEBUG #define TEST_COLLAPS_ENTRY(array) array, SK_ARRAY_COUNT(array) static void test_collaps_duplicates() { static bool gOnce; if (gOnce) { return; } gOnce = true; const float src0[] = { 0 }; const float src1[] = { 0, 0 }; const float src2[] = { 0, 1 }; const float src3[] = { 0, 0, 0 }; const float src4[] = { 0, 0, 1 }; const float src5[] = { 0, 1, 1 }; const float src6[] = { 0, 1, 2 }; const struct { const float* fData; int fCount; int fCollapsedCount; } data[] = { { TEST_COLLAPS_ENTRY(src0), 1 }, { TEST_COLLAPS_ENTRY(src1), 1 }, { TEST_COLLAPS_ENTRY(src2), 2 }, { TEST_COLLAPS_ENTRY(src3), 1 }, { TEST_COLLAPS_ENTRY(src4), 2 }, { TEST_COLLAPS_ENTRY(src5), 2 }, { TEST_COLLAPS_ENTRY(src6), 3 }, }; for (size_t i = 0; i < SK_ARRAY_COUNT(data); ++i) { float dst[3]; memcpy(dst, data[i].fData, data[i].fCount * sizeof(dst[0])); int count = collaps_duplicates(dst, data[i].fCount); SkASSERT(data[i].fCollapsedCount == count); for (int j = 1; j < count; ++j) { SkASSERT(dst[j-1] < dst[j]); } } } #endif #if defined _WIN32 && _MSC_VER >= 1300 && defined SK_SCALAR_IS_FIXED // disable warning : unreachable code if building fixed point for windows desktop #pragma warning ( disable : 4702 ) #endif /* Solve coeff(t) == 0, returning the number of roots that lie withing 0 < t < 1. coeff[0]t^3 + coeff[1]t^2 + coeff[2]t + coeff[3] Eliminates repeated roots (so that all tValues are distinct, and are always in increasing order. */ static int solve_cubic_polynomial(const SkFP coeff[4], SkScalar tValues[3]) { #ifndef SK_SCALAR_IS_FLOAT return 0; // this is not yet implemented for software float #endif if (SkScalarNearlyZero(coeff[0])) // we're just a quadratic { return SkFindUnitQuadRoots(coeff[1], coeff[2], coeff[3], tValues); } SkFP a, b, c, Q, R; { SkASSERT(coeff[0] != 0); SkFP inva = SkFPInvert(coeff[0]); a = SkFPMul(coeff[1], inva); b = SkFPMul(coeff[2], inva); c = SkFPMul(coeff[3], inva); } Q = SkFPDivInt(SkFPSub(SkFPMul(a,a), SkFPMulInt(b, 3)), 9); // R = (2*a*a*a - 9*a*b + 27*c) / 54; R = SkFPMulInt(SkFPMul(SkFPMul(a, a), a), 2); R = SkFPSub(R, SkFPMulInt(SkFPMul(a, b), 9)); R = SkFPAdd(R, SkFPMulInt(c, 27)); R = SkFPDivInt(R, 54); SkFP Q3 = SkFPMul(SkFPMul(Q, Q), Q); SkFP R2MinusQ3 = SkFPSub(SkFPMul(R,R), Q3); SkFP adiv3 = SkFPDivInt(a, 3); SkScalar* roots = tValues; SkScalar r; if (SkFPLT(R2MinusQ3, 0)) // we have 3 real roots { #ifdef SK_SCALAR_IS_FLOAT float theta = sk_float_acos(R / sk_float_sqrt(Q3)); float neg2RootQ = -2 * sk_float_sqrt(Q); r = neg2RootQ * sk_float_cos(theta/3) - adiv3; if (is_unit_interval(r)) *roots++ = r; r = neg2RootQ * sk_float_cos((theta + 2*SK_ScalarPI)/3) - adiv3; if (is_unit_interval(r)) *roots++ = r; r = neg2RootQ * sk_float_cos((theta - 2*SK_ScalarPI)/3) - adiv3; if (is_unit_interval(r)) *roots++ = r; SkDEBUGCODE(test_collaps_duplicates();) // now sort the roots int count = (int)(roots - tValues); SkASSERT((unsigned)count <= 3); bubble_sort(tValues, count); count = collaps_duplicates(tValues, count); roots = tValues + count; // so we compute the proper count below #endif } else // we have 1 real root { SkFP A = SkFPAdd(SkFPAbs(R), SkFPSqrt(R2MinusQ3)); A = SkFPCubeRoot(A); if (SkFPGT(R, 0)) A = SkFPNeg(A); if (A != 0) A = SkFPAdd(A, SkFPDiv(Q, A)); r = SkFPToScalar(SkFPSub(A, adiv3)); if (is_unit_interval(r)) *roots++ = r; } return (int)(roots - tValues); } /* Looking for F' dot F'' == 0 A = b - a B = c - 2b + a C = d - 3c + 3b - a F' = 3Ct^2 + 6Bt + 3A F'' = 6Ct + 6B F' dot F'' -> CCt^3 + 3BCt^2 + (2BB + CA)t + AB */ static void formulate_F1DotF2(const SkScalar src[], SkFP coeff[4]) { SkScalar a = src[2] - src[0]; SkScalar b = src[4] - 2 * src[2] + src[0]; SkScalar c = src[6] + 3 * (src[2] - src[4]) - src[0]; SkFP A = SkScalarToFP(a); SkFP B = SkScalarToFP(b); SkFP C = SkScalarToFP(c); coeff[0] = SkFPMul(C, C); coeff[1] = SkFPMulInt(SkFPMul(B, C), 3); coeff[2] = SkFPMulInt(SkFPMul(B, B), 2); coeff[2] = SkFPAdd(coeff[2], SkFPMul(C, A)); coeff[3] = SkFPMul(A, B); } // EXPERIMENTAL: can set this to zero to accept all t-values 0 < t < 1 //#define kMinTValueForChopping (SK_Scalar1 / 256) #define kMinTValueForChopping 0 /* Looking for F' dot F'' == 0 A = b - a B = c - 2b + a C = d - 3c + 3b - a F' = 3Ct^2 + 6Bt + 3A F'' = 6Ct + 6B F' dot F'' -> CCt^3 + 3BCt^2 + (2BB + CA)t + AB */ int SkFindCubicMaxCurvature(const SkPoint src[4], SkScalar tValues[3]) { SkFP coeffX[4], coeffY[4]; int i; formulate_F1DotF2(&src[0].fX, coeffX); formulate_F1DotF2(&src[0].fY, coeffY); for (i = 0; i < 4; i++) coeffX[i] = SkFPAdd(coeffX[i],coeffY[i]); SkScalar t[3]; int count = solve_cubic_polynomial(coeffX, t); int maxCount = 0; // now remove extrema where the curvature is zero (mins) // !!!! need a test for this !!!! for (i = 0; i < count; i++) { // if (not_min_curvature()) if (t[i] > kMinTValueForChopping && t[i] < SK_Scalar1 - kMinTValueForChopping) tValues[maxCount++] = t[i]; } return maxCount; } int SkChopCubicAtMaxCurvature(const SkPoint src[4], SkPoint dst[13], SkScalar tValues[3]) { SkScalar t_storage[3]; if (tValues == NULL) tValues = t_storage; int count = SkFindCubicMaxCurvature(src, tValues); if (dst) { if (count == 0) memcpy(dst, src, 4 * sizeof(SkPoint)); else SkChopCubicAt(src, dst, tValues, count); } return count + 1; } bool SkXRayCrossesMonotonicCubic(const SkXRay& pt, const SkPoint cubic[4], bool* ambiguous) { if (ambiguous) { *ambiguous = false; } // Find the minimum and maximum y of the extrema, which are the // first and last points since this cubic is monotonic SkScalar min_y = SkMinScalar(cubic[0].fY, cubic[3].fY); SkScalar max_y = SkMaxScalar(cubic[0].fY, cubic[3].fY); if (pt.fY == cubic[0].fY || pt.fY < min_y || pt.fY > max_y) { // The query line definitely does not cross the curve if (ambiguous) { *ambiguous = (pt.fY == cubic[0].fY); } return false; } bool pt_at_extremum = (pt.fY == cubic[3].fY); SkScalar min_x = SkMinScalar( SkMinScalar( SkMinScalar(cubic[0].fX, cubic[1].fX), cubic[2].fX), cubic[3].fX); if (pt.fX < min_x) { // The query line definitely crosses the curve if (ambiguous) { *ambiguous = pt_at_extremum; } return true; } SkScalar max_x = SkMaxScalar( SkMaxScalar( SkMaxScalar(cubic[0].fX, cubic[1].fX), cubic[2].fX), cubic[3].fX); if (pt.fX > max_x) { // The query line definitely does not cross the curve return false; } // Do a binary search to find the parameter value which makes y as // close as possible to the query point. See whether the query // line's origin is to the left of the associated x coordinate. // kMaxIter is chosen as the number of mantissa bits for a float, // since there's no way we are going to get more precision by // iterating more times than that. const int kMaxIter = 23; SkPoint eval; int iter = 0; SkScalar upper_t; SkScalar lower_t; // Need to invert direction of t parameter if cubic goes up // instead of down if (cubic[3].fY > cubic[0].fY) { upper_t = SK_Scalar1; lower_t = SkFloatToScalar(0); } else { upper_t = SkFloatToScalar(0); lower_t = SK_Scalar1; } do { SkScalar t = SkScalarAve(upper_t, lower_t); SkEvalCubicAt(cubic, t, &eval, NULL, NULL); if (pt.fY > eval.fY) { lower_t = t; } else { upper_t = t; } } while (++iter < kMaxIter && !SkScalarNearlyZero(eval.fY - pt.fY)); if (pt.fX <= eval.fX) { if (ambiguous) { *ambiguous = pt_at_extremum; } return true; } return false; } int SkNumXRayCrossingsForCubic(const SkXRay& pt, const SkPoint cubic[4], bool* ambiguous) { int num_crossings = 0; SkPoint monotonic_cubics[10]; int num_monotonic_cubics = SkChopCubicAtYExtrema(cubic, monotonic_cubics); if (ambiguous) { *ambiguous = false; } bool locally_ambiguous; if (SkXRayCrossesMonotonicCubic(pt, &monotonic_cubics[0], &locally_ambiguous)) ++num_crossings; if (ambiguous) { *ambiguous |= locally_ambiguous; } if (num_monotonic_cubics > 0) if (SkXRayCrossesMonotonicCubic(pt, &monotonic_cubics[3], &locally_ambiguous)) ++num_crossings; if (ambiguous) { *ambiguous |= locally_ambiguous; } if (num_monotonic_cubics > 1) if (SkXRayCrossesMonotonicCubic(pt, &monotonic_cubics[6], &locally_ambiguous)) ++num_crossings; if (ambiguous) { *ambiguous |= locally_ambiguous; } return num_crossings; } //////////////////////////////////////////////////////////////////////////////// /* Find t value for quadratic [a, b, c] = d. Return 0 if there is no solution within [0, 1) */ static SkScalar quad_solve(SkScalar a, SkScalar b, SkScalar c, SkScalar d) { // At^2 + Bt + C = d SkScalar A = a - 2 * b + c; SkScalar B = 2 * (b - a); SkScalar C = a - d; SkScalar roots[2]; int count = SkFindUnitQuadRoots(A, B, C, roots); SkASSERT(count <= 1); return count == 1 ? roots[0] : 0; } /* given a quad-curve and a point (x,y), chop the quad at that point and return the new quad's offCurve point. Should only return false if the computed pos is the start of the curve (i.e. root == 0) */ static bool quad_pt2OffCurve(const SkPoint quad[3], SkScalar x, SkScalar y, SkPoint* offCurve) { const SkScalar* base; SkScalar value; if (SkScalarAbs(x) < SkScalarAbs(y)) { base = &quad[0].fX; value = x; } else { base = &quad[0].fY; value = y; } // note: this returns 0 if it thinks value is out of range, meaning the // root might return something outside of [0, 1) SkScalar t = quad_solve(base[0], base[2], base[4], value); if (t > 0) { SkPoint tmp[5]; SkChopQuadAt(quad, tmp, t); *offCurve = tmp[1]; return true; } else { /* t == 0 means either the value triggered a root outside of [0, 1) For our purposes, we can ignore the <= 0 roots, but we want to catch the >= 1 roots (which given our caller, will basically mean a root of 1, give-or-take numerical instability). If we are in the >= 1 case, return the existing offCurve point. The test below checks to see if we are close to the "end" of the curve (near base[4]). Rather than specifying a tolerance, I just check to see if value is on to the right/left of the middle point (depending on the direction/sign of the end points). */ if ((base[0] < base[4] && value > base[2]) || (base[0] > base[4] && value < base[2])) // should root have been 1 { *offCurve = quad[1]; return true; } } return false; } #ifdef SK_SCALAR_IS_FLOAT // Due to floating point issues (i.e., 1.0f - SK_ScalarRoot2Over2 != // SK_ScalarRoot2Over2 - SK_ScalarTanPIOver8) a cruder root2over2 // approximation is required to make the quad circle points convex. The // root of the problem is that with the root2over2 value in SkScalar.h // the arcs really are ever so slightly concave. Some alternative fixes // to this problem (besides just arbitrarily pushing out the mid-point as // is done here) are: // Adjust all the points (not just the middle) to both better approximate // the curve and remain convex // Switch over to using cubics rather then quads // Use a different method to create the mid-point (e.g., compute // the two side points, average them, then move it out as needed #ifdef SK_REDEFINE_ROOT2OVER2_TO_MAKE_ARCTOS_CONVEX #define SK_ScalarRoot2Over2_QuadCircle 0.7072f #else #define SK_ScalarRoot2Over2_QuadCircle SK_ScalarRoot2Over2 #endif #else #define SK_ScalarRoot2Over2_QuadCircle SK_FixedRoot2Over2 #endif static const SkPoint gQuadCirclePts[kSkBuildQuadArcStorage] = { { SK_Scalar1, 0 }, { SK_Scalar1, SK_ScalarTanPIOver8 }, { SK_ScalarRoot2Over2_QuadCircle, SK_ScalarRoot2Over2_QuadCircle }, { SK_ScalarTanPIOver8, SK_Scalar1 }, { 0, SK_Scalar1 }, { -SK_ScalarTanPIOver8, SK_Scalar1 }, { -SK_ScalarRoot2Over2_QuadCircle, SK_ScalarRoot2Over2_QuadCircle }, { -SK_Scalar1, SK_ScalarTanPIOver8 }, { -SK_Scalar1, 0 }, { -SK_Scalar1, -SK_ScalarTanPIOver8 }, { -SK_ScalarRoot2Over2_QuadCircle, -SK_ScalarRoot2Over2_QuadCircle }, { -SK_ScalarTanPIOver8, -SK_Scalar1 }, { 0, -SK_Scalar1 }, { SK_ScalarTanPIOver8, -SK_Scalar1 }, { SK_ScalarRoot2Over2_QuadCircle, -SK_ScalarRoot2Over2_QuadCircle }, { SK_Scalar1, -SK_ScalarTanPIOver8 }, { SK_Scalar1, 0 } }; int SkBuildQuadArc(const SkVector& uStart, const SkVector& uStop, SkRotationDirection dir, const SkMatrix* userMatrix, SkPoint quadPoints[]) { // rotate by x,y so that uStart is (1.0) SkScalar x = SkPoint::DotProduct(uStart, uStop); SkScalar y = SkPoint::CrossProduct(uStart, uStop); SkScalar absX = SkScalarAbs(x); SkScalar absY = SkScalarAbs(y); int pointCount; // check for (effectively) coincident vectors // this can happen if our angle is nearly 0 or nearly 180 (y == 0) // ... we use the dot-prod to distinguish between 0 and 180 (x > 0) if (absY <= SK_ScalarNearlyZero && x > 0 && ((y >= 0 && kCW_SkRotationDirection == dir) || (y <= 0 && kCCW_SkRotationDirection == dir))) { // just return the start-point quadPoints[0].set(SK_Scalar1, 0); pointCount = 1; } else { if (dir == kCCW_SkRotationDirection) y = -y; // what octant (quadratic curve) is [xy] in? int oct = 0; bool sameSign = true; if (0 == y) { oct = 4; // 180 SkASSERT(SkScalarAbs(x + SK_Scalar1) <= SK_ScalarNearlyZero); } else if (0 == x) { SkASSERT(absY - SK_Scalar1 <= SK_ScalarNearlyZero); if (y > 0) oct = 2; // 90 else oct = 6; // 270 } else { if (y < 0) oct += 4; if ((x < 0) != (y < 0)) { oct += 2; sameSign = false; } if ((absX < absY) == sameSign) oct += 1; } int wholeCount = oct << 1; memcpy(quadPoints, gQuadCirclePts, (wholeCount + 1) * sizeof(SkPoint)); const SkPoint* arc = &gQuadCirclePts[wholeCount]; if (quad_pt2OffCurve(arc, x, y, &quadPoints[wholeCount + 1])) { quadPoints[wholeCount + 2].set(x, y); wholeCount += 2; } pointCount = wholeCount + 1; } // now handle counter-clockwise and the initial unitStart rotation SkMatrix matrix; matrix.setSinCos(uStart.fY, uStart.fX); if (dir == kCCW_SkRotationDirection) { matrix.preScale(SK_Scalar1, -SK_Scalar1); } if (userMatrix) { matrix.postConcat(*userMatrix); } matrix.mapPoints(quadPoints, pointCount); return pointCount; }