Add cubic clipping, similar to that of quad clipping. Both Newton-Raphson and Bisection are implemented, because it is not clear which one will yield the highest performance on a given platform. Bisection is turned on as the default.
git-svn-id: http://skia.googlecode.com/svn/trunk@101 2bbb7eff-a529-9590-31e7-b0007b416f81
diff --git a/src/core/SkCubicClipper.cpp b/src/core/SkCubicClipper.cpp
new file mode 100644
index 0000000..9022ff3
--- /dev/null
+++ b/src/core/SkCubicClipper.cpp
@@ -0,0 +1,166 @@
+/*
+ * Copyright (C) 2009 The Android Open Source Project
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#include "SkCubicClipper.h"
+#include "SkGeometry.h"
+
+SkCubicClipper::SkCubicClipper() {}
+
+void SkCubicClipper::setClip(const SkIRect& clip) {
+ // conver to scalars, since that's where we'll see the points
+ fClip.set(clip);
+}
+
+
+static bool chopMonoCubicAtY(SkPoint pts[4], SkScalar y, SkScalar* t) {
+ SkScalar ycrv[4];
+ ycrv[0] = pts[0].fY - y;
+ ycrv[1] = pts[1].fY - y;
+ ycrv[2] = pts[2].fY - y;
+ ycrv[3] = pts[3].fY - y;
+
+#ifdef NEWTON_RAPHSON // Quadratic convergence, typically <= 3 iterations.
+ // Initial guess.
+ // TODO(turk): Check for zero denominator? Shouldn't happen unless the curve
+ // is not only monotonic but degenerate.
+#ifdef SK_SCALAR_IS_FLOAT
+ SkScalar t1 = ycrv[0] / (ycrv[0] - ycrv[3]);
+#else // !SK_SCALAR_IS_FLOAT
+ SkScalar t1 = SkDivBits(ycrv[0], ycrv[0] - ycrv[3], 16);
+#endif // !SK_SCALAR_IS_FLOAT
+
+ // Newton's iterations.
+ const SkScalar tol = SK_Scalar1 / 16384; // This leaves 2 fixed noise bits.
+ SkScalar t0;
+ const int maxiters = 5;
+ int iters = 0;
+ bool converged;
+ do {
+ t0 = t1;
+ SkScalar y01 = SkScalarInterp(ycrv[0], ycrv[1], t0);
+ SkScalar y12 = SkScalarInterp(ycrv[1], ycrv[2], t0);
+ SkScalar y23 = SkScalarInterp(ycrv[2], ycrv[3], t0);
+ SkScalar y012 = SkScalarInterp(y01, y12, t0);
+ SkScalar y123 = SkScalarInterp(y12, y23, t0);
+ SkScalar y0123 = SkScalarInterp(y012, y123, t0);
+ SkScalar yder = (y123 - y012) * 3;
+ // TODO(turk): check for yder==0: horizontal.
+#ifdef SK_SCALAR_IS_FLOAT
+ t1 -= y0123 / yder;
+#else // !SK_SCALAR_IS_FLOAT
+ t1 -= SkDivBits(y0123, yder, 16);
+#endif // !SK_SCALAR_IS_FLOAT
+ converged = SkScalarAbs(t1 - t0) <= tol; // NaN-safe
+ ++iters;
+ } while (!converged && (iters < maxiters));
+ *t = t1; // Return the result.
+
+ // The result might be valid, even if outside of the range [0, 1], but
+ // we never evaluate a Bezier outside this interval, so we return false.
+ if (t1 < 0 || t1 > SK_Scalar1)
+ return false; // This shouldn't happen, but check anyway.
+ return converged;
+
+#else // BISECTION // Linear convergence, typically 16 iterations.
+
+ // Check that the endpoints straddle zero.
+ SkScalar tNeg, tPos; // Parameter where the function is negative and positive, respectively.
+ if (ycrv[0] < 0) {
+ if (ycrv[3] < 0)
+ return false;
+ tNeg = 0;
+ tPos = SK_Scalar1;
+ } else if (ycrv[0] > 0) {
+ if (ycrv[3] > 0)
+ return false;
+ tNeg = SK_Scalar1;
+ tPos = 0;
+ } else {
+ *t = 0;
+ return true;
+ }
+
+ const SkScalar tol = SK_Scalar1 / 65536; // 1 for fixed, 1e-5 for float.
+ int iters = 0;
+ do {
+ SkScalar tMid = (tPos + tNeg) / 2;
+ SkScalar y01 = SkScalarInterp(ycrv[0], ycrv[1], tMid);
+ SkScalar y12 = SkScalarInterp(ycrv[1], ycrv[2], tMid);
+ SkScalar y23 = SkScalarInterp(ycrv[2], ycrv[3], tMid);
+ SkScalar y012 = SkScalarInterp(y01, y12, tMid);
+ SkScalar y123 = SkScalarInterp(y12, y23, tMid);
+ SkScalar y0123 = SkScalarInterp(y012, y123, tMid);
+ if (y0123 == 0) {
+ *t = tMid;
+ return true;
+ }
+ if (y0123 < 0) tNeg = tMid;
+ else tPos = tMid;
+ ++iters;
+ } while (!(SkScalarAbs(tPos - tNeg) <= tol)); // Nan-safe
+
+ *t = (tNeg + tPos) / 2;
+ return true;
+#endif // BISECTION
+}
+
+
+bool SkCubicClipper::clipCubic(const SkPoint srcPts[4], SkPoint dst[4]) {
+ bool reverse;
+
+ // we need the data to be monotonically descending in Y
+ if (srcPts[0].fY > srcPts[2].fY) {
+ dst[0] = srcPts[3];
+ dst[1] = srcPts[2];
+ dst[2] = srcPts[1];
+ dst[3] = srcPts[0];
+ reverse = true;
+ } else {
+ memcpy(dst, srcPts, 3 * sizeof(SkPoint));
+ reverse = false;
+ }
+
+ // are we completely above or below
+ const SkScalar ctop = fClip.fTop;
+ const SkScalar cbot = fClip.fBottom;
+ if (dst[2].fY <= ctop || dst[0].fY >= cbot) {
+ return false;
+ }
+
+ SkScalar t;
+ SkPoint tmp[5]; // for SkChopCubicAt
+
+ // are we partially above
+ if (dst[0].fY < ctop && chopMonoCubicAtY(dst, ctop, &t)) {
+ SkChopCubicAt(dst, tmp, t);
+ dst[0] = tmp[2];
+ dst[1] = tmp[3];
+ }
+
+ // are we partially below
+ if (dst[2].fY > cbot && chopMonoCubicAtY(dst, cbot, &t)) {
+ SkChopCubicAt(dst, tmp, t);
+ dst[1] = tmp[1];
+ dst[2] = tmp[2];
+ }
+
+ if (reverse) {
+ SkTSwap<SkPoint>(dst[0], dst[3]);
+ SkTSwap<SkPoint>(dst[1], dst[2]);
+ }
+ return true;
+}
+