Convert SkClassifyCubic to double precision
Even though it's in homogeneous coordinates, we still get unfortunate
artifacts if this math is not done in double precision. Prioritizing
correctness for now; we can revisit in the future as the need for
performance dictates.
Bug: skia:
Change-Id: If416ef6b70291f1454fcb9f7630d1108644ac2a5
Reviewed-on: https://skia-review.googlesource.com/19501
Commit-Queue: Chris Dalton <csmartdalton@google.com>
Reviewed-by: Greg Daniel <egdaniel@google.com>
diff --git a/src/core/SkGeometry.cpp b/src/core/SkGeometry.cpp
index e140286..2927d18 100644
--- a/src/core/SkGeometry.cpp
+++ b/src/core/SkGeometry.cpp
@@ -533,10 +533,10 @@
// Assumes the third component of points is 1.
// Calcs p0 . (p1 x p2)
-static SkScalar calc_dot_cross_cubic(const SkPoint& p0, const SkPoint& p1, const SkPoint& p2) {
- const SkScalar xComp = p0.fX * (p1.fY - p2.fY);
- const SkScalar yComp = p0.fY * (p2.fX - p1.fX);
- const SkScalar wComp = p1.fX * p2.fY - p1.fY * p2.fX;
+static double calc_dot_cross_cubic(const SkPoint& p0, const SkPoint& p1, const SkPoint& p2) {
+ const double xComp = (double) p0.fX * (double) (p1.fY - p2.fY);
+ const double yComp = (double) p0.fY * (double) (p2.fX - p1.fX);
+ const double wComp = (double) p1.fX * (double) p2.fY - (double) p1.fY * (double) p2.fX;
return (xComp + yComp + wComp);
}
@@ -549,42 +549,42 @@
// a2 = p1 . (p0 x p3)
// a3 = p2 . (p1 x p0)
// Places the values of d1, d2, d3 in array d passed in
-static void calc_cubic_inflection_func(const SkPoint p[4], SkScalar d[4]) {
- SkScalar a1 = calc_dot_cross_cubic(p[0], p[3], p[2]);
- SkScalar a2 = calc_dot_cross_cubic(p[1], p[0], p[3]);
- SkScalar a3 = calc_dot_cross_cubic(p[2], p[1], p[0]);
+static void calc_cubic_inflection_func(const SkPoint p[4], double d[4]) {
+ const double a1 = calc_dot_cross_cubic(p[0], p[3], p[2]);
+ const double a2 = calc_dot_cross_cubic(p[1], p[0], p[3]);
+ const double a3 = calc_dot_cross_cubic(p[2], p[1], p[0]);
- d[3] = 3.f * a3;
+ d[3] = 3 * a3;
d[2] = d[3] - a2;
d[1] = d[2] - a2 + a1;
d[0] = 0;
}
-static void normalize_t_s(float t[], float s[], int count) {
+static void normalize_t_s(double t[], double s[], int count) {
// Keep the exponents at or below zero to avoid overflow down the road.
for (int i = 0; i < count; ++i) {
SkASSERT(0 != s[i]);
- union { float value; int32_t bits; } tt, ss, norm;
+ union { double value; int64_t bits; } tt, ss, norm;
tt.value = t[i];
ss.value = s[i];
- int32_t expT = ((tt.bits >> 23) & 0xff) - 127,
- expS = ((ss.bits >> 23) & 0xff) - 127;
- int32_t expNorm = -SkTMax(expT, expS) + 127;
- SkASSERT(expNorm > 0 && expNorm < 255); // ensure we have a valid non-zero exponent.
- norm.bits = expNorm << 23;
+ int64_t expT = ((tt.bits >> 52) & 0x7ff) - 1023,
+ expS = ((ss.bits >> 52) & 0x7ff) - 1023;
+ int64_t expNorm = -SkTMax(expT, expS) + 1023;
+ SkASSERT(expNorm > 0 && expNorm < 2047); // ensure we have a valid non-zero exponent.
+ norm.bits = expNorm << 52;
t[i] *= norm.value;
s[i] *= norm.value;
}
}
-static void sort_and_orient_t_s(SkScalar t[2], SkScalar s[2]) {
+static void sort_and_orient_t_s(double t[2], double s[2]) {
// This copysign/abs business orients the implicit function so positive values are always on the
// "left" side of the curve.
- t[1] = -SkScalarCopySign(t[1], t[1] * s[1]);
- s[1] = -SkScalarAbs(s[1]);
+ t[1] = -copysign(t[1], t[1] * s[1]);
+ s[1] = -fabs(s[1]);
// Ensure t[0]/s[0] <= t[1]/s[1] (s[1] is negative from above).
- if (SkScalarCopySign(s[1], s[0]) * t[0] > -SkScalarAbs(s[0]) * t[1]) {
+ if (copysign(s[1], s[0]) * t[0] > -fabs(s[0]) * t[1]) {
std::swap(t[0], t[1]);
std::swap(s[0], s[1]);
}
@@ -600,15 +600,15 @@
// d1 = 0, d2 != 0 Cusp (with cusp at infinity)
// d1 = d2 = 0, d3 != 0 Quadratic
// d1 = d2 = d3 = 0 Line or Point
-static SkCubicType classify_cubic(const SkScalar d[4], SkScalar t[2], SkScalar s[2]) {
- SkScalar tolerance = SkTMax(SkScalarAbs(d[1]), SkScalarAbs(d[2]));
- tolerance = SkTMax(tolerance, SkScalarAbs(d[3]));
+static SkCubicType classify_cubic(const double d[4], double t[2], double s[2]) {
+ double tolerance = SkTMax(fabs(d[1]), fabs(d[2]));
+ tolerance = SkTMax(tolerance, fabs(d[3]));
tolerance = tolerance * 1e-5;
- if (!SkScalarNearlyZero(d[1], tolerance)) {
- const SkScalar discr = 3 * d[2] * d[2] - 4 * d[1] * d[3];
+ if (fabs(d[1]) > tolerance) {
+ const double discr = 3 * d[2] * d[2] - 4 * d[1] * d[3];
if (discr > 0) {
if (t && s) {
- const SkScalar q = 3 * d[2] + SkScalarCopySign(SkScalarSqrt(3 * discr), d[2]);
+ const double q = 3 * d[2] + copysign(sqrt(3 * discr), d[2]);
t[0] = q;
s[0] = 6 * d[1];
t[1] = 2 * d[3];
@@ -619,7 +619,7 @@
return SkCubicType::kSerpentine;
} else if (discr < 0) {
if (t && s) {
- const SkScalar q = d[2] + SkScalarCopySign(SkScalarSqrt(-discr), d[2]);
+ const double q = d[2] + copysign(sqrt(-discr), d[2]);
t[0] = q;
s[0] = 2 * d[1];
t[1] = 2 * (d[2] * d[2] - d[3] * d[1]);
@@ -641,7 +641,7 @@
return SkCubicType::kLocalCusp;
}
} else {
- if (!SkScalarNearlyZero(d[2], tolerance)) {
+ if (fabs(d[2]) > tolerance) {
if (t && s) {
t[0] = d[3];
s[0] = 3 * d[2];
@@ -655,15 +655,14 @@
t[0] = t[1] = 1;
s[0] = s[1] = 0; // infinity
}
- return !SkScalarNearlyZero(d[3], tolerance) ? SkCubicType::kQuadratic
- : SkCubicType::kLineOrPoint;
+ return fabs(d[3]) > tolerance ? SkCubicType::kQuadratic : SkCubicType::kLineOrPoint;
}
}
}
-SkCubicType SkClassifyCubic(const SkPoint src[4], SkScalar t[2], SkScalar s[2], SkScalar d[4]) {
- SkScalar localD[4];
- SkScalar* dd = d ? d : localD;
+SkCubicType SkClassifyCubic(const SkPoint src[4], double t[2], double s[2], double d[4]) {
+ double localD[4];
+ double* dd = d ? d : localD;
calc_cubic_inflection_func(src, dd);
return classify_cubic(dd, t, s);
}
diff --git a/src/core/SkGeometry.h b/src/core/SkGeometry.h
index ef789a3..3cb96c3 100644
--- a/src/core/SkGeometry.h
+++ b/src/core/SkGeometry.h
@@ -182,8 +182,8 @@
https://www.microsoft.com/en-us/research/wp-content/uploads/2005/01/p1000-loop.pdf
*/
-SkCubicType SkClassifyCubic(const SkPoint p[4], SkScalar t[2] = nullptr, SkScalar s[2] = nullptr,
- SkScalar d[4] = nullptr);
+SkCubicType SkClassifyCubic(const SkPoint p[4], double t[2] = nullptr, double s[2] = nullptr,
+ double d[4] = nullptr);
///////////////////////////////////////////////////////////////////////////////
diff --git a/src/gpu/GrPathUtils.cpp b/src/gpu/GrPathUtils.cpp
index 022d2c6..6047e60 100644
--- a/src/gpu/GrPathUtils.cpp
+++ b/src/gpu/GrPathUtils.cpp
@@ -763,7 +763,7 @@
// | ..L.. | * | . . . . | == | 0 0 1/3 1 |
// | ..K.. | | 1 1 1 1 | | 0 1/3 2/3 1 |
//
-static void calc_quadratic_klm(const SkPoint pts[4], SkScalar d3, SkMatrix* klm) {
+static void calc_quadratic_klm(const SkPoint pts[4], double d3, SkMatrix* klm) {
SkMatrix klmAtPts;
klmAtPts.setAll(0, 1.f/3, 1,
0, 0, 1,
@@ -798,22 +798,26 @@
-nx, -ny, k);
}
-SkCubicType GrPathUtils::getCubicKLM(const SkPoint src[4], SkMatrix* klm, SkScalar t[2],
- SkScalar s[2]) {
- SkScalar d[4];
+SkCubicType GrPathUtils::getCubicKLM(const SkPoint src[4], SkMatrix* klm, double t[2],
+ double s[2]) {
+ double d[4];
SkCubicType type = SkClassifyCubic(src, t, s, d);
+
+ const SkScalar tt[2] = {static_cast<SkScalar>(t[0]), static_cast<SkScalar>(t[1])};
+ const SkScalar ss[2] = {static_cast<SkScalar>(s[0]), static_cast<SkScalar>(s[1])};
+
switch (type) {
case SkCubicType::kSerpentine:
- calc_serp_klm(src, t[0], s[0], t[1], s[1], klm);
+ calc_serp_klm(src, tt[0], ss[0], tt[1], ss[1], klm);
break;
case SkCubicType::kLoop:
- calc_loop_klm(src, t[0], s[0], t[1], s[1], klm);
+ calc_loop_klm(src, tt[0], ss[0], tt[1], ss[1], klm);
break;
case SkCubicType::kLocalCusp:
- calc_serp_klm(src, t[0], s[0], t[1], s[1], klm);
+ calc_serp_klm(src, tt[0], ss[0], tt[1], ss[1], klm);
break;
case SkCubicType::kCuspAtInfinity:
- calc_inf_cusp_klm(src, t[0], s[0], klm);
+ calc_inf_cusp_klm(src, tt[0], ss[0], klm);
break;
case SkCubicType::kQuadratic:
calc_quadratic_klm(src, d[3], klm);
@@ -831,20 +835,20 @@
SkSTArray<2, SkScalar> chops;
*loopIndex = -1;
- SkScalar t[2], s[2];
+ double t[2], s[2];
if (SkCubicType::kLoop == GrPathUtils::getCubicKLM(src, klm, t, s)) {
- t[0] /= s[0];
- t[1] /= s[1];
- SkASSERT(t[0] <= t[1]); // Technically t0 != t1 in a loop, but there may be FP error.
+ SkScalar t0 = static_cast<SkScalar>(t[0] / s[0]);
+ SkScalar t1 = static_cast<SkScalar>(t[1] / s[1]);
+ SkASSERT(t0 <= t1); // Technically t0 != t1 in a loop, but there may be FP error.
- if (t[0] < 1 && t[1] > 0) {
+ if (t0 < 1 && t1 > 0) {
*loopIndex = 0;
- if (t[0] > 0) {
- chops.push_back(t[0]);
+ if (t0 > 0) {
+ chops.push_back(t0);
*loopIndex = 1;
}
- if (t[1] < 1) {
- chops.push_back(t[1]);
+ if (t1 < 1) {
+ chops.push_back(t1);
*loopIndex = chops.count() - 1;
}
}
diff --git a/src/gpu/GrPathUtils.h b/src/gpu/GrPathUtils.h
index 35f94d0..49c4ee5 100644
--- a/src/gpu/GrPathUtils.h
+++ b/src/gpu/GrPathUtils.h
@@ -146,7 +146,7 @@
// intersect with K (See SkClassifyCubic).
//
// Returns the cubic's classification.
- SkCubicType getCubicKLM(const SkPoint src[4], SkMatrix* klm, SkScalar t[2], SkScalar s[2]);
+ SkCubicType getCubicKLM(const SkPoint src[4], SkMatrix* klm, double t[2], double s[2]);
// Chops the cubic bezier passed in by src, at the double point (intersection point)
// if the curve is a cubic loop. If it is a loop, there will be two parametric values for
diff --git a/src/pathops/SkPathOpsCubic.cpp b/src/pathops/SkPathOpsCubic.cpp
index d7b905c..5c672fa 100644
--- a/src/pathops/SkPathOpsCubic.cpp
+++ b/src/pathops/SkPathOpsCubic.cpp
@@ -248,14 +248,14 @@
if (cubic.monotonicInX() && cubic.monotonicInY()) {
return 0;
}
- SkScalar tt[2], ss[2];
+ double tt[2], ss[2];
SkCubicType cubicType = SkClassifyCubic(pointsPtr, tt, ss);
switch (cubicType) {
case SkCubicType::kLoop: {
- const SkScalar &td = tt[0], &te = tt[1], &sd = ss[0], &se = ss[1];
+ const double &td = tt[0], &te = tt[1], &sd = ss[0], &se = ss[1];
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);
+ t[0] = static_cast<SkScalar>((td * se + te * sd) / (2 * sd * se));
SkASSERT(roughly_between(0, *t, 1));
return (int) (t[0] > 0 && t[0] < 1);
}