Use more stable root finding methods for cubics
Applies the quadratic formula from "Numerical Recipes in C",
Section 5.6, to the homogeneous quadratic equations that find cubic
inflection points and loop intersections. Also addresses KLM
orientation ahead of time, rather than negating K and L after the
fact.
Bug: skia:
Change-Id: Ic7e0818e2fe49b7724f9b583bae52281cfb1aea1
Reviewed-on: https://skia-review.googlesource.com/13481
Commit-Queue: Chris Dalton <csmartdalton@google.com>
Reviewed-by: Cary Clark <caryclark@google.com>
diff --git a/src/gpu/GrPathUtils.cpp b/src/gpu/GrPathUtils.cpp
index 8991500..91a48f7 100644
--- a/src/gpu/GrPathUtils.cpp
+++ b/src/gpu/GrPathUtils.cpp
@@ -648,25 +648,11 @@
return skipRow;
}
-static void negate_kl(SkMatrix* klm) {
- // We could use klm->postScale(-1, -1), but it ends up doing a full matrix multiply.
- for (int i = 0; i < 6; ++i) {
- (*klm)[i] = -(*klm)[i];
- }
-}
-
-static void calc_serp_klm(const SkPoint pts[4], const SkScalar d[4], SkMatrix* klm) {
+static void calc_serp_klm(const SkPoint pts[4], SkScalar tl, SkScalar sl, SkScalar tm, SkScalar sm,
+ SkMatrix* klm) {
SkMatrix CIT;
int skipCol = calc_inverse_transpose_power_basis_matrix(pts, &CIT);
- SkASSERT(d[0] >= 0);
- const SkScalar root = SkScalarSqrt(3 * d[0]);
-
- const SkScalar tl = 3 * d[2] + root;
- const SkScalar sl = 6 * d[1];
- const SkScalar tm = 3 * d[2] - root;
- const SkScalar sm = 6 * d[1];
-
SkMatrix klmCoeffs;
int col = 0;
if (0 != skipCol) {
@@ -694,17 +680,10 @@
klmCoeffs[8] = tm * tm * tm;
klm->setConcat(klmCoeffs, CIT);
-
- // If d1 > 0 we need to flip the orientation of our curve
- // This is done by negating the k and l values
- // We want negative distance values to be on the inside
- if (d[1] > 0) {
- negate_kl(klm);
- }
}
-static void calc_loop_klm(const SkPoint pts[4], SkScalar d1, SkScalar td, SkScalar sd,
- SkScalar te, SkScalar se, SkMatrix* klm) {
+static void calc_loop_klm(const SkPoint pts[4], SkScalar td, SkScalar sd, SkScalar te, SkScalar se,
+ SkMatrix* klm) {
SkMatrix CIT;
int skipCol = calc_inverse_transpose_power_basis_matrix(pts, &CIT);
@@ -738,25 +717,13 @@
klmCoeffs[8] = te * te * td;
klm->setConcat(klmCoeffs, CIT);
-
- // For the general loop curve, we flip the orientation in the same pattern as the serp case
- // above. Thus we only check d1. Technically we should check the value of the hessian as well
- // cause we care about the sign of d1*Hessian. However, the Hessian is always negative outside
- // the loop section and positive inside. We take care of the flipping for the loop sections
- // later on.
- if (d1 > 0) {
- negate_kl(klm);
- }
}
// For the case when we have a cusp at a parameter value of infinity (discr == 0, d1 == 0).
-static void calc_inf_cusp_klm(const SkPoint pts[4], SkScalar d2, SkScalar d3, SkMatrix* klm) {
+static void calc_inf_cusp_klm(const SkPoint pts[4], SkScalar tn, SkScalar sn, SkMatrix* klm) {
SkMatrix CIT;
int skipCol = calc_inverse_transpose_power_basis_matrix(pts, &CIT);
- const SkScalar tn = d3;
- const SkScalar sn = 3 * d2;
-
SkMatrix klmCoeffs;
int col = 0;
if (0 != skipCol) {
@@ -814,7 +781,7 @@
// If d3 > 0 we need to flip the orientation of our curve
// This is done by negating the k and l values
if (d3 > 0) {
- negate_kl(klm);
+ klm->postScale(-1, -1);
}
}
@@ -846,11 +813,11 @@
int chop_count = 0;
if (SkCubicType::kLoop == cType) {
SkASSERT(d[0] < 0);
- const SkScalar tempSqrt = SkScalarSqrt(-d[0]);
- td = d[2] + tempSqrt;
- sd = 2.f * d[1];
- te = d[2] - tempSqrt;
- se = 2.f * d[1];
+ const SkScalar q = d[2] + SkScalarCopySign(SkScalarSqrt(-d[0]), d[2]);
+ td = q;
+ sd = 2 * d[1];
+ te = 2 * (d[2] * d[2] - d[3] * d[1]);
+ se = d[1] * q;
t1 = td / sd;
t2 = te / se;
@@ -895,16 +862,38 @@
if (klm) {
switch (cType) {
- case SkCubicType::kSerpentine:
- case SkCubicType::kLocalCusp:
- calc_serp_klm(src, d, klm);
+ case SkCubicType::kSerpentine: {
+ SkASSERT(d[0] >= 0);
+ const SkScalar q = 3 * d[2] + SkScalarCopySign(SkScalarSqrt(3 * d[0]), d[2]);
+ const SkScalar tl = q;
+ const SkScalar sl = 6 * d[1];
+ const SkScalar tm = 2 * d[3];
+ const SkScalar sm = q;
+ // This copysign/abs business orients the implicit function so positive values are
+ // always on the "left" side of the curve.
+ calc_serp_klm(src, tl, sl, -SkScalarCopySign(tm, tm * sm), -SkScalarAbs(sm), klm);
break;
+ }
+ case SkCubicType::kLocalCusp: {
+ SkASSERT(0 == d[0]);
+ const SkScalar t = d[2];
+ const SkScalar s = 2 * d[1];
+ // This copysign/abs business orients the implicit function so positive values are
+ // always on the "left" side of the curve.
+ calc_serp_klm(src, t, s, -SkScalarCopySign(t, t * s), -SkScalarAbs(s), klm);
+ break;
+ }
case SkCubicType::kLoop:
- calc_loop_klm(src, d[1], td, sd, te, se, klm);
+ // This copysign/abs business orients the implicit function so positive values are
+ // always on the "left" side of the curve.
+ calc_loop_klm(src, td, sd, -SkScalarCopySign(te, te * se), -SkScalarAbs(se), klm);
break;
- case SkCubicType::kInfiniteCusp:
- calc_inf_cusp_klm(src, d[2], d[3], klm);
+ case SkCubicType::kInfiniteCusp: {
+ const SkScalar tn = d[3];
+ const SkScalar sn = 3 * d[2];
+ calc_inf_cusp_klm(src, tn, sn, klm);
break;
+ }
case SkCubicType::kQuadratic:
calc_quadratic_klm(src, d[3], klm);
break;
diff --git a/src/pathops/SkPathOpsCubic.cpp b/src/pathops/SkPathOpsCubic.cpp
index 1e74eb0..794e54f 100644
--- a/src/pathops/SkPathOpsCubic.cpp
+++ b/src/pathops/SkPathOpsCubic.cpp
@@ -255,16 +255,14 @@
// crib code from gpu path utils that finds t values where loop self-intersects
// use it to find mid of t values which should be a friendly place to chop
SkASSERT(d[0] < 0);
- SkScalar tempSqrt = SkScalarSqrt(-d[0]);
- SkScalar ls = d[2] - tempSqrt;
- SkScalar lt = 2.f * d[1];
- SkScalar ms = d[2] + tempSqrt;
- SkScalar mt = 2.f * d[1];
- if (roughly_between(0, ls, lt) && roughly_between(0, ms, mt)) {
- ls = ls / lt;
- ms = ms / mt;
- SkASSERT(roughly_between(0, ls, 1) && roughly_between(0, ms, 1));
- t[0] = (ls + ms) / 2;
+ const SkScalar q = d[2] + SkScalarCopySign(SkScalarSqrt(-d[0]), d[2]);
+ const SkScalar td = q;
+ const SkScalar sd = 2 * d[1];
+ const SkScalar te = 2 * (d[2] * d[2] - d[3] * d[1]);
+ const SkScalar se = d[1] * q;
+ if (roughly_between(0, td, sd) && roughly_between(0, te, se)) {
+ SkASSERT(roughly_between(0, td / sd, 1) && roughly_between(0, te / se, 1));
+ t[0] = (td * se + te * sd) / (2 * sd * se);
SkASSERT(roughly_between(0, *t, 1));
return (int) (t[0] > 0 && t[0] < 1);
}