aboutsummaryrefslogtreecommitdiffhomepage
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/core/SkGeometry.cpp201
-rw-r--r--src/gpu/ccpr/GrCCCubicShader.cpp44
-rw-r--r--src/pathops/SkPathOpsCubic.cpp2
3 files changed, 101 insertions, 146 deletions
diff --git a/src/core/SkGeometry.cpp b/src/core/SkGeometry.cpp
index d42a3b1961..6ed8911320 100644
--- a/src/core/SkGeometry.cpp
+++ b/src/core/SkGeometry.cpp
@@ -542,63 +542,28 @@ static double calc_dot_cross_cubic(const SkPoint& p0, const SkPoint& p1, const S
return (xComp + yComp + wComp);
}
-// Calc coefficients of I(s,t) where roots of I are inflection points of curve
-// I(s,t) = t*(3*d0*s^2 - 3*d1*s*t + d2*t^2)
-// d0 = a1 - 2*a2+3*a3
-// d1 = -a2 + 3*a3
-// d2 = 3*a3
-// a1 = p0 . (p3 x p2)
-// 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], 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 * a3;
- d[2] = d[3] - a2;
- d[1] = d[2] - a2 + a1;
- d[0] = 0;
-}
-
-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]); // classify_cubic should not call this method when s[i] is 0 or NaN.
-
- uint64_t bitsT, bitsS;
- memcpy(&bitsT, &t[i], sizeof(double));
- memcpy(&bitsS, &s[i], sizeof(double));
-
- uint64_t maxExponent = SkTMax(bitsT & 0x7ff0000000000000, bitsS & 0x7ff0000000000000);
+// Returns a positive power of 2 that, when multiplied by n, and excepting the two edge cases listed
+// below, shifts the exponent of n to yield a magnitude somewhere inside [1..2).
+// Returns 2^1023 if abs(n) < 2^-1022 (including 0).
+// Returns NaN if n is Inf or NaN.
+inline static double previous_inverse_pow2(double n) {
+ uint64_t bits;
+ memcpy(&bits, &n, sizeof(double));
+ bits = ((1023llu*2 << 52) + ((1llu << 52) - 1)) - bits; // exp=-exp
+ bits &= (0x7ffllu) << 52; // mantissa=1.0, sign=0
+ memcpy(&n, &bits, sizeof(double));
+ return n;
+}
+
+inline static void write_cubic_inflection_roots(double t0, double s0, double t1, double s1,
+ double* t, double* s) {
+ t[0] = t0;
+ s[0] = s0;
-#ifdef SK_DEBUG
- uint64_t maxExponentValue = maxExponent >> 52;
- // Ensure max(absT,absS) is NOT in denormalized form. SkClassifyCubic is given fp32 points,
- // and does not call this method when s==0, so this should never happen.
- SkASSERT(0 != maxExponentValue);
- // Ensure 1/max(absT,absS) will NOT be in denormalized form. SkClassifyCubic is given fp32
- // points, so this should never happen.
- SkASSERT(2046 != maxExponentValue);
-#endif
-
- // Pick a normalizer that scales the larger exponent to 1 (aka 1023 in biased form), but
- // does NOT change the mantissa (thus preserving accuracy).
- double normalizer;
- uint64_t normalizerExponent = (uint64_t(1023 * 2) << 52) - maxExponent;
- memcpy(&normalizer, &normalizerExponent, sizeof(double));
-
- t[i] *= normalizer;
- s[i] *= normalizer;
- }
-}
-
-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] = -copysign(t[1], t[1] * s[1]);
- s[1] = -fabs(s[1]);
+ t[1] = -copysign(t1, t1 * s1);
+ s[1] = -fabs(s1);
// Ensure t[0]/s[0] <= t[1]/s[1] (s[1] is negative from above).
if (copysign(s[1], s[0]) * t[0] > -fabs(s[0]) * t[1]) {
@@ -607,78 +572,78 @@ static void sort_and_orient_t_s(double t[2], double s[2]) {
}
}
-// See "Resolution Independent Curve Rendering using Programmable Graphics Hardware"
-// https://www.microsoft.com/en-us/research/wp-content/uploads/2005/01/p1000-loop.pdf
-// discr(I) = 3*d2^2 - 4*d1*d3
-// Classification:
-// d1 != 0, discr(I) > 0 Serpentine
-// d1 != 0, discr(I) < 0 Loop
-// d1 != 0, discr(I) = 0 Cusp (with inflection at infinity)
-// 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 double d[4], double t[2], double s[2]) {
- if (0 == d[1]) {
- if (0 == d[2]) {
- if (t && s) {
- t[0] = t[1] = 1;
- s[0] = s[1] = 0; // infinity
- }
- return 0 == d[3] ? SkCubicType::kLineOrPoint : SkCubicType::kQuadratic;
- }
- if (t && s) {
- t[0] = d[3];
- s[0] = 3 * d[2];
- normalize_t_s(t, s, 1);
- t[1] = 1;
- s[1] = 0; // infinity
- }
- return SkCubicType::kCuspAtInfinity;
+SkCubicType SkClassifyCubic(const SkPoint P[4], double t[2], double s[2], double d[4]) {
+ // Find the cubic's inflection function, I = [T^3 -3T^2 3T -1] dot D. (D0 will always be 0
+ // for integral cubics.)
+ //
+ // See "Resolution Independent Curve Rendering using Programmable Graphics Hardware",
+ // 4.2 Curve Categorization:
+ //
+ // https://www.microsoft.com/en-us/research/wp-content/uploads/2005/01/p1000-loop.pdf
+ double A1 = calc_dot_cross_cubic(P[0], P[3], P[2]);
+ double A2 = calc_dot_cross_cubic(P[1], P[0], P[3]);
+ double A3 = calc_dot_cross_cubic(P[2], P[1], P[0]);
+
+ double D3 = 3 * A3;
+ double D2 = D3 - A2;
+ double D1 = D2 - A2 + A1;
+
+ // Shift the exponents in D so the largest magnitude falls somewhere in 1..2. This protects us
+ // from overflow down the road while solving for roots and KLM functionals.
+ double Dmax = std::max(std::max(fabs(D1), fabs(D2)), fabs(D3));
+ double norm = previous_inverse_pow2(Dmax);
+ D1 *= norm;
+ D2 *= norm;
+ D3 *= norm;
+
+ if (d) {
+ d[3] = D3;
+ d[2] = D2;
+ d[1] = D1;
+ d[0] = 0;
}
- const double discr = 3 * d[2] * d[2] - 4 * d[1] * d[3];
- if (discr > 0) {
- if (t && s) {
- 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];
- s[1] = q;
- normalize_t_s(t, s, 2);
- sort_and_orient_t_s(t, s);
- }
- return SkCubicType::kSerpentine;
- } else if (discr < 0) {
- if (t && s) {
- 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]);
- s[1] = d[1] * q;
- normalize_t_s(t, s, 2);
- sort_and_orient_t_s(t, s);
+ // Now use the inflection function to classify the cubic.
+ //
+ // See "Resolution Independent Curve Rendering using Programmable Graphics Hardware",
+ // 4.4 Integral Cubics:
+ //
+ // https://www.microsoft.com/en-us/research/wp-content/uploads/2005/01/p1000-loop.pdf
+ if (0 != D1) {
+ double discr = 3*D2*D2 - 4*D1*D3;
+ if (discr > 0) { // Serpentine.
+ if (t && s) {
+ double q = 3*D2 + copysign(sqrt(3*discr), D2);
+ write_cubic_inflection_roots(q, 6*D1, 2*D3, q, t, s);
+ }
+ return SkCubicType::kSerpentine;
+ } else if (discr < 0) { // Loop.
+ if (t && s) {
+ double q = D2 + copysign(sqrt(-discr), D2);
+ write_cubic_inflection_roots(q, 2*D1, 2*(D2*D2 - D3*D1), D1*q, t, s);
+ }
+ return SkCubicType::kLoop;
+ } else { // Cusp.
+ if (t && s) {
+ write_cubic_inflection_roots(D2, 2*D1, D2, 2*D1, t, s);
+ }
+ return SkCubicType::kLocalCusp;
}
- return SkCubicType::kLoop;
} else {
- if (t && s) {
- t[0] = d[2];
- s[0] = 2 * d[1];
- normalize_t_s(t, s, 1);
- t[1] = t[0];
- s[1] = s[0];
- sort_and_orient_t_s(t, s);
+ if (0 != D2) { // Cusp at T=infinity.
+ if (t && s) {
+ write_cubic_inflection_roots(D3, 3*D2, 1, 0, t, s); // T1=infinity.
+ }
+ return SkCubicType::kCuspAtInfinity;
+ } else { // Degenerate.
+ if (t && s) {
+ write_cubic_inflection_roots(1, 0, 1, 0, t, s); // T0=T1=infinity.
+ }
+ return 0 != D3 ? SkCubicType::kQuadratic : SkCubicType::kLineOrPoint;
}
- return SkCubicType::kLocalCusp;
}
}
-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);
-}
-
template <typename T> void bubble_sort(T array[], int count) {
for (int i = count - 1; i > 0; --i)
for (int j = i; j > 0; --j)
diff --git a/src/gpu/ccpr/GrCCCubicShader.cpp b/src/gpu/ccpr/GrCCCubicShader.cpp
index f8a8352383..9a57c1cd44 100644
--- a/src/gpu/ccpr/GrCCCubicShader.cpp
+++ b/src/gpu/ccpr/GrCCCubicShader.cpp
@@ -15,28 +15,6 @@ using Shader = GrCCCoverageProcessor::Shader;
void GrCCCubicShader::emitSetupCode(GrGLSLVertexGeoBuilder* s, const char* pts,
const char* wind, const char** /*outHull4*/) const {
- // Define a function that normalizes the homogeneous coordinates T=t/s in order to avoid
- // exponent overflow.
- SkString normalizeHomogCoordFn;
- GrShaderVar coord("coord", kFloat2_GrSLType);
- s->emitFunction(kFloat2_GrSLType, "normalize_homogeneous_coord", 1, &coord,
- s->getProgramBuilder()->shaderCaps()->fpManipulationSupport()
- // Exponent manipulation version: Scale the exponents so the larger
- // component has a magnitude in 1..2.
- // (Neither component should be infinity because ccpr crops big paths.)
- ? "int exp;"
- "frexp(max(abs(coord.t), abs(coord.s)), exp);"
- "return coord * ldexp(1, 1 - exp);"
-
- // Division version: Divide by the component with the larger magnitude.
- // (Both should not be 0 because ccpr catches degenerate cubics.)
- : "bool swap = abs(coord.t) > abs(coord.s);"
- "coord = swap ? coord.ts : coord;"
- "coord = float2(1, coord.t/coord.s);"
- "return swap ? coord.ts : coord;",
-
- &normalizeHomogCoordFn);
-
// Find the cubic's power basis coefficients.
s->codeAppendf("float2x4 C = float4x4(-1, 3, -3, 1, "
" 3, -6, 3, 0, "
@@ -48,6 +26,21 @@ void GrCCCubicShader::emitSetupCode(GrGLSLVertexGeoBuilder* s, const char* pts,
s->codeAppend ("float D2 = -determinant(float2x2(C[0].xz, C[1].xz));");
s->codeAppend ("float D1 = +determinant(float2x2(C));");
+ // Shift the exponents in D so the largest magnitude falls somewhere in 1..2. This protects us
+ // from overflow while solving for roots and KLM functionals.
+ s->codeAppend ("float Dmax = max(max(abs(D1), abs(D2)), abs(D3));");
+ s->codeAppend ("float norm;");
+ if (s->getProgramBuilder()->shaderCaps()->fpManipulationSupport()) {
+ s->codeAppend ("int exp;");
+ s->codeAppend ("frexp(Dmax, exp);");
+ s->codeAppend ("norm = ldexp(1, 1 - exp);");
+ } else {
+ s->codeAppend ("norm = 1/Dmax;"); // Dmax will not be 0 because we cull line cubics on CPU.
+ }
+ s->codeAppend ("D3 *= norm;");
+ s->codeAppend ("D2 *= norm;");
+ s->codeAppend ("D1 *= norm;");
+
// Calculate the KLM matrix.
s->declareGlobal(fKLMMatrix);
s->codeAppend ("float discr = 3*D2*D2 - 4*D1*D3;");
@@ -56,10 +49,9 @@ void GrCCCubicShader::emitSetupCode(GrGLSLVertexGeoBuilder* s, const char* pts,
s->codeAppend ("q = x*D2 + (D2 >= 0 ? q : -q);");
s->codeAppend ("float2 l, m;");
- s->codeAppendf("l.ts = %s(float2(q, 2*x * D1));", normalizeHomogCoordFn.c_str());
- s->codeAppendf("m.ts = %s(float2(2, q) * (discr >= 0 ? float2(D3, 1) "
- ": float2(D2*D2 - D3*D1, D1)));",
- normalizeHomogCoordFn.c_str());
+ s->codeAppend ("l.ts = float2(q, 2*x * D1);");
+ s->codeAppend ("m.ts = float2(2, q) * (discr >= 0 ? float2(D3, 1) "
+ ": float2(D2*D2 - D3*D1, D1));");
s->codeAppend ("float4 K;");
s->codeAppend ("float4 lm = l.sstt * m.stst;");
diff --git a/src/pathops/SkPathOpsCubic.cpp b/src/pathops/SkPathOpsCubic.cpp
index eb2ddb7c7d..8d3baa038b 100644
--- a/src/pathops/SkPathOpsCubic.cpp
+++ b/src/pathops/SkPathOpsCubic.cpp
@@ -253,9 +253,7 @@ int SkDCubic::ComplexBreak(const SkPoint pointsPtr[4], SkScalar* t) {
case SkCubicType::kLoop: {
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] = static_cast<SkScalar>((td * se + te * sd) / (2 * sd * se));
- SkASSERT(roughly_between(0, *t, 1));
return (int) (t[0] > 0 && t[0] < 1);
}
}
generated by cgit v1.2.3 (git 2.39.1) at 2024-07-01 19:51:28 +0000