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>
This commit is contained in:
Chris Dalton 2017-06-12 11:22:54 -06:00 committed by Skia Commit-Bot
parent e3a0be73a6
commit 390f6cd6f9
5 changed files with 59 additions and 56 deletions

View File

@ -533,10 +533,10 @@ int SkChopCubicAtInflections(const SkPoint src[], SkPoint dst[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 @@ static SkScalar calc_dot_cross_cubic(const SkPoint& p0, const SkPoint& p1, const
// 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 @@ static void sort_and_orient_t_s(SkScalar t[2], SkScalar s[2]) {
// 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 @@ static SkCubicType classify_cubic(const SkScalar d[4], SkScalar t[2], SkScalar s
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 @@ static SkCubicType classify_cubic(const SkScalar d[4], SkScalar t[2], SkScalar s
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 @@ static SkCubicType classify_cubic(const SkScalar d[4], SkScalar t[2], SkScalar s
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);
}

View File

@ -182,8 +182,8 @@ enum class SkCubicType {
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);
///////////////////////////////////////////////////////////////////////////////

View File

@ -763,7 +763,7 @@ static void calc_inf_cusp_klm(const SkPoint pts[4], SkScalar tn, SkScalar sn, Sk
// | ..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 @@ static void calc_line_klm(const SkPoint pts[4], SkMatrix* klm) {
-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 @@ int GrPathUtils::chopCubicAtLoopIntersection(const SkPoint src[4], SkPoint dst[1
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;
}
}

View File

@ -146,7 +146,7 @@ namespace GrPathUtils {
// 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

View File

@ -248,14 +248,14 @@ int SkDCubic::ComplexBreak(const SkPoint pointsPtr[4], SkScalar* t) {
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);
}