blob: d8470cc4e4b01021f3d64fb83d6187593a9d29d1 [file] [log] [blame]
// Graphite-specific vertex shader code
const float $PI = 3.141592653589793238;
// Wang's formula gives the minimum number of evenly spaced (in the parametric sense) line segments
// that a bezier curve must be chopped into in order to guarantee all lines stay within a distance
// of "1/precision" pixels from the true curve. Its definition for a bezier curve of degree "n" is
// as follows:
//
// maxLength = max([length(p[i+2] - 2p[i+1] + p[i]) for (0 <= i <= n-2)])
// numParametricSegments = sqrt(maxLength * precision * n*(n - 1)/8)
//
// (Goldman, Ron. (2003). 5.6.3 Wang's Formula. "Pyramid Algorithms: A Dynamic Programming Approach
// to Curves and Surfaces for Geometric Modeling". Morgan Kaufmann Publishers.)
const float $kDegree = 3;
const float $kPrecision = 4; // Must match skgpu::tess::kPrecision
const float $kLengthTerm = ($kDegree * ($kDegree - 1) / 8.0) * $kPrecision;
const float $kLengthTermPow2 = (($kDegree * $kDegree) * (($kDegree - 1) * ($kDegree - 1)) / 64.0) *
($kPrecision * $kPrecision);
// Returns the length squared of the largest forward difference from Wang's cubic formula.
float wangs_formula_max_fdiff_p2(float2 p0, float2 p1, float2 p2, float2 p3,
float2x2 matrix) {
float2 d0 = matrix * (fma(float2(-2), p1, p2) + p0);
float2 d1 = matrix * (fma(float2(-2), p2, p3) + p1);
return max(dot(d0,d0), dot(d1,d1));
}
float wangs_formula_cubic(float2 p0, float2 p1, float2 p2, float2 p3,
float2x2 matrix) {
float m = wangs_formula_max_fdiff_p2(p0, p1, p2, p3, matrix);
return max(ceil(sqrt($kLengthTerm * sqrt(m))), 1.0);
}
float wangs_formula_cubic_log2(float2 p0, float2 p1, float2 p2, float2 p3,
float2x2 matrix) {
float m = wangs_formula_max_fdiff_p2(p0, p1, p2, p3, matrix);
return ceil(log2(max($kLengthTermPow2 * m, 1.0)) * .25);
}
float wangs_formula_conic_p2(float2 p0, float2 p1, float2 p2, float w) {
// Translate the bounding box center to the origin.
float2 C = (min(min(p0, p1), p2) + max(max(p0, p1), p2)) * 0.5;
p0 -= C;
p1 -= C;
p2 -= C;
// Compute max length.
float m = sqrt(max(max(dot(p0,p0), dot(p1,p1)), dot(p2,p2)));
// Compute forward differences.
float2 dp = fma(float2(-2.0 * w), p1, p0) + p2;
float dw = abs(fma(-2.0, w, 2.0));
// Compute numerator and denominator for parametric step size of linearization. Here, the
// epsilon referenced from the cited paper is 1/precision.
float rp_minus_1 = max(0.0, fma(m, $kPrecision, -1.0));
float numer = length(dp) * $kPrecision + rp_minus_1 * dw;
float denom = 4 * min(w, 1.0);
return numer/denom;
}
float wangs_formula_conic(float2 p0, float2 p1, float2 p2, float w) {
float n2 = wangs_formula_conic_p2(p0, p1, p2, w);
return max(ceil(sqrt(n2)), 1.0);
}
float wangs_formula_conic_log2(float2 p0, float2 p1, float2 p2, float w) {
float n2 = wangs_formula_conic_p2(p0, p1, p2, w);
return ceil(log2(max(n2, 1.0)) * .5);
}
// Returns the normalized difference between a and b, i.e. normalize(a - b), with care taken for
// if 'a' and/or 'b' have large coordinates.
float2 robust_normalize_diff(float2 a, float2 b) {
float2 diff = a - b;
if (diff == float2(0.0)) {
return float2(0.0);
} else {
float invMag = 1.0 / max(abs(diff.x), abs(diff.y));
return normalize(invMag * diff);
}
}
// Returns the cosine of the angle between a and b, assuming a and b are unit vectors already.
// Guaranteed to be between [-1, 1].
float cosine_between_unit_vectors(float2 a, float2 b) {
// Since a and b are assumed to be normalized, the cosine is equal to the dot product, although
// we clamp that to ensure it falls within the expected range of [-1, 1].
return clamp(dot(a, b), -1.0, 1.0);
}
// Extends the middle radius to either the miter point, or the bevel edge if we surpassed the
// miter limit and need to revert to a bevel join.
float miter_extent(float cosTheta, float miterLimit) {
float x = fma(cosTheta, .5, .5);
return (x * miterLimit * miterLimit >= 1.0) ? inversesqrt(x) : sqrt(x);
}
// Returns the number of radial segments required for each radian of rotation, in order for the
// curve to appear "smooth" as defined by the approximate device-space stroke radius.
float num_radial_segments_per_radian(float approxDevStrokeRadius) {
return .5 / acos(max(1.0 - (1.0 / $kPrecision) / approxDevStrokeRadius, -1.0));
}
// Unlike mix(), this does not return b when t==1. But it otherwise seems to get better
// precision than "a*(1 - t) + b*t" for things like chopping cubics on exact cusp points.
// We override this result anyway when t==1 so it shouldn't be a problem.
float unchecked_mix(float a, float b, float T) {
return fma(b - a, T, a);
}
float2 unchecked_mix(float2 a, float2 b, float T) {
return fma(b - a, float2(T), a);
}
float4 unchecked_mix(float4 a, float4 b, float4 T) {
return fma(b - a, T, a);
}
// Compute a vertex position for the curve described by p01 and p23 packed control points,
// tessellated to the given resolve level, and assuming it will be drawn as a filled curve.
float2 tessellate_filled_curve(float2x2 vectorXform,
float resolveLevel, float idxInResolveLevel,
float4 p01, float4 p23) {
float2 localcoord;
if (isinf(p23.z)) {
// This patch is an exact triangle.
localcoord = (resolveLevel != 0) ? p01.zw
: (idxInResolveLevel != 0) ? p23.xy
: p01.xy;
} else {
float2 p0=p01.xy, p1=p01.zw, p2=p23.xy, p3=p23.zw;
float w = -1; // w < 0 tells us to treat the instance as an integral cubic.
float maxResolveLevel;
if (isinf(p23.w)) {
// Conics are 3 points, with the weight in p3.
w = p3.x;
maxResolveLevel = wangs_formula_conic_log2(vectorXform*p0,
vectorXform*p1,
vectorXform*p2, w);
p1 *= w; // Unproject p1.
p3 = p2; // Duplicate the endpoint for shared code that also runs on cubics.
} else {
// The patch is an integral cubic.
maxResolveLevel = wangs_formula_cubic_log2(p0, p1, p2, p3, vectorXform);
}
if (resolveLevel > maxResolveLevel) {
// This vertex is at a higher resolve level than we need. Demote to a lower
// resolveLevel, which will produce a degenerate triangle.
idxInResolveLevel = floor(ldexp(idxInResolveLevel,
int(maxResolveLevel - resolveLevel)));
resolveLevel = maxResolveLevel;
}
// Promote our location to a discrete position in the maximum fixed resolve level.
// This is extra paranoia to ensure we get the exact same fp32 coordinates for
// colocated points from different resolve levels (e.g., the vertices T=3/4 and
// T=6/8 should be exactly colocated).
float fixedVertexID = floor(.5 + ldexp(idxInResolveLevel, int(5 - resolveLevel)));
if (0 < fixedVertexID && fixedVertexID < 32) {
float T = fixedVertexID * (1 / 32.0);
// Evaluate at T. Use De Casteljau's for its accuracy and stability.
float2 ab = mix(p0, p1, T);
float2 bc = mix(p1, p2, T);
float2 cd = mix(p2, p3, T);
float2 abc = mix(ab, bc, T);
float2 bcd = mix(bc, cd, T);
float2 abcd = mix(abc, bcd, T);
// Evaluate the conic weight at T.
float u = mix(1.0, w, T);
float v = w + 1 - u; // == mix(w, 1, T)
float uv = mix(u, v, T);
localcoord = (w < 0) ? /*cubic*/ abcd : /*conic*/ abc/uv;
} else {
localcoord = (fixedVertexID == 0) ? p0.xy : p3.xy;
}
}
return localcoord;
}
// Device coords are in xy, local coords are in zw, since for now perspective isn't supported.
float4 tessellate_stroked_curve(float edgeID, float maxEdges,
float2x2 affineMatrix,
float2 translate,
float maxScale /* derived from affineMatrix */,
float4 p01, float4 p23,
float2 lastControlPoint,
float2 strokeParams) {
float2 p0=p01.xy, p1=p01.zw, p2=p23.xy, p3=p23.zw;
float w = -1; // w<0 means the curve is an integral cubic.
if (isinf(p23.w)) {
// Conics are 3 points, with the weight in p3.
w = p3.x;
p3 = p2; // Setting p3 equal to p2 works for the remaining rotational logic.
}
// Call Wang's formula to determine parametric segments before transform points for hairlines
// so that it is consistent with how the CPU tested the control points for chopping.
float numParametricSegments;
if (w < 0) {
if (p0 == p1 && p2 == p3) {
numParametricSegments = 1; // a line
} else {
numParametricSegments = wangs_formula_cubic(p0, p1, p2, p3, affineMatrix);
}
} else {
numParametricSegments = wangs_formula_conic(affineMatrix * p0,
affineMatrix * p1,
affineMatrix * p2, w);
}
// Matches skgpu::tess::StrokeParams
float strokeRadius = strokeParams.x;
float joinType = strokeParams.y; // <0 = round join, ==0 = bevel join, >0 encodes miter limit
bool isHairline = strokeParams.x == 0.0;
float numRadialSegmentsPerRadian;
if (isHairline) {
numRadialSegmentsPerRadian = num_radial_segments_per_radian(1.0);
strokeRadius = 0.5;
} else {
numRadialSegmentsPerRadian = num_radial_segments_per_radian(maxScale * strokeParams.x);
}
if (isHairline) {
// Hairline case. Transform the points before tessellation. We can still hold off on the
// translate until the end; we just need to perform the scale and skew right now.
p0 = affineMatrix * p0;
p1 = affineMatrix * p1;
p2 = affineMatrix * p2;
p3 = affineMatrix * p3;
lastControlPoint = affineMatrix * lastControlPoint;
}
// Find the starting and ending tangents.
float2 tan0 = robust_normalize_diff((p0 == p1) ? ((p1 == p2) ? p3 : p2) : p1, p0);
float2 tan1 = robust_normalize_diff(p3, (p3 == p2) ? ((p2 == p1) ? p0 : p1) : p2);
if (tan0 == float2(0)) {
// The stroke is a point. This special case tells us to draw a stroke-width circle as a
// 180 degree point stroke instead.
tan0 = float2(1,0);
tan1 = float2(-1,0);
}
// Determine how many edges to give to the join. We emit the first and final edges
// of the join twice: once full width and once restricted to half width. This guarantees
// perfect seaming by matching the vertices from the join as well as from the strokes on
// either side.
float numEdgesInJoin;
if (joinType >= 0 /*Is the join not a round type?*/) {
// Bevel(0) and miter(+) joins get 1 and 2 segments respectively.
// +2 because we emit the beginning and ending edges twice (see above comments).
numEdgesInJoin = sign(joinType) + 1 + 2;
} else {
float2 prevTan = robust_normalize_diff(p0, lastControlPoint);
float joinRads = acos(cosine_between_unit_vectors(prevTan, tan0));
float numRadialSegmentsInJoin = max(ceil(joinRads * numRadialSegmentsPerRadian), 1);
// +2 because we emit the beginning and ending edges twice (see above comment).
numEdgesInJoin = numRadialSegmentsInJoin + 2;
// The stroke section needs at least two edges. Don't assign more to the join than
// "maxEdges - 2". (This is only relevant when the ideal max edge count calculated
// on the CPU had to be limited to maxEdges in the draw call).
numEdgesInJoin = min(numEdgesInJoin, maxEdges - 2);
}
// Find which direction the curve turns.
// NOTE: Since the curve is not allowed to inflect, we can just check F'(.5) x F''(.5).
// NOTE: F'(.5) x F''(.5) has the same sign as (P2 - P0) x (P3 - P1)
float turn = cross_length_2d(p2 - p0, p3 - p1);
float combinedEdgeID = abs(edgeID) - numEdgesInJoin;
if (combinedEdgeID < 0) {
tan1 = tan0;
// Don't let tan0 become zero. The code as-is isn't built to handle that case. tan0=0
// means the join is disabled, and to disable it with the existing code we can leave
// tan0 equal to tan1.
if (lastControlPoint != p0) {
tan0 = robust_normalize_diff(p0, lastControlPoint);
}
turn = cross_length_2d(tan0, tan1);
}
// Calculate the curve's starting angle and rotation.
float cosTheta = cosine_between_unit_vectors(tan0, tan1);
float rotation = acos(cosTheta);
if (turn < 0) {
// Adjust sign of rotation to match the direction the curve turns.
rotation = -rotation;
}
float numRadialSegments;
float strokeOutset = sign(edgeID);
if (combinedEdgeID < 0) {
// We belong to the preceding join. The first and final edges get duplicated, so we only
// have "numEdgesInJoin - 2" segments.
numRadialSegments = numEdgesInJoin - 2;
numParametricSegments = 1; // Joins don't have parametric segments.
p3 = p2 = p1 = p0; // Colocate all points on the junction point.
// Shift combinedEdgeID to the range [-1, numRadialSegments]. This duplicates the first
// edge and lands one edge at the very end of the join. (The duplicated final edge will
// actually come from the section of our strip that belongs to the stroke.)
combinedEdgeID += numRadialSegments + 1;
// We normally restrict the join on one side of the junction, but if the tangents are
// nearly equivalent this could theoretically result in bad seaming and/or cracks on the
// side we don't put it on. If the tangents are nearly equivalent then we leave the join
// double-sided.
float sinEpsilon = 1e-2; // ~= sin(180deg / 3000)
bool tangentsNearlyParallel =
(abs(turn) * inversesqrt(dot(tan0, tan0) * dot(tan1, tan1))) < sinEpsilon;
if (!tangentsNearlyParallel || dot(tan0, tan1) < 0) {
// There are two edges colocated at the beginning. Leave the first one double sided
// for seaming with the previous stroke. (The double sided edge at the end will
// actually come from the section of our strip that belongs to the stroke.)
if (combinedEdgeID >= 0) {
strokeOutset = (turn < 0) ? min(strokeOutset, 0) : max(strokeOutset, 0);
}
}
combinedEdgeID = max(combinedEdgeID, 0);
} else {
// We belong to the stroke. Unless numRadialSegmentsPerRadian is incredibly high,
// clamping to maxCombinedSegments will be a no-op because the draw call was invoked with
// sufficient vertices to cover the worst case scenario of 180 degree rotation.
float maxCombinedSegments = maxEdges - numEdgesInJoin - 1;
numRadialSegments = max(ceil(abs(rotation) * numRadialSegmentsPerRadian), 1);
numRadialSegments = min(numRadialSegments, maxCombinedSegments);
numParametricSegments = min(numParametricSegments,
maxCombinedSegments - numRadialSegments + 1);
}
// Additional parameters for final tessellation evaluation.
float radsPerSegment = rotation / numRadialSegments;
float numCombinedSegments = numParametricSegments + numRadialSegments - 1;
bool isFinalEdge = (combinedEdgeID >= numCombinedSegments);
if (combinedEdgeID > numCombinedSegments) {
strokeOutset = 0; // The strip has more edges than we need. Drop this one.
}
// Edge #2 extends to the miter point.
if (abs(edgeID) == 2 && joinType > 0/*Is the join a miter type?*/) {
strokeOutset *= miter_extent(cosTheta, joinType/*miterLimit*/);
}
float2 tangent, strokeCoord;
if (combinedEdgeID != 0 && !isFinalEdge) {
// Compute the location and tangent direction of the stroke edge with the integral id
// "combinedEdgeID", where combinedEdgeID is the sorted-order index of parametric and radial
// edges. Start by finding the tangent function's power basis coefficients. These define a
// tangent direction (scaled by some uniform value) as:
// |T^2|
// Tangent_Direction(T) = dx,dy = |A 2B C| * |T |
// |. . .| |1 |
float2 A, B, C = p1 - p0;
float2 D = p3 - p0;
if (w >= 0.0) {
// P0..P2 represent a conic and P3==P2. The derivative of a conic has a cumbersome
// order-4 denominator. However, this isn't necessary if we are only interested in a
// vector in the same *direction* as a given tangent line. Since the denominator scales
// dx and dy uniformly, we can throw it out completely after evaluating the derivative
// with the standard quotient rule. This leaves us with a simpler quadratic function
// that we use to find a tangent.
C *= w;
B = .5*D - C;
A = (w - 1.0) * D;
p1 *= w;
} else {
float2 E = p2 - p1;
B = E - C;
A = fma(float2(-3), E, D);
}
// FIXME(crbug.com/800804,skbug.com/11268): Consider normalizing the exponents in A,B,C at
// this point in order to prevent fp32 overflow.
// Now find the coefficients that give a tangent direction from a parametric edge ID:
//
// |parametricEdgeID^2|
// Tangent_Direction(parametricEdgeID) = dx,dy = |A B_ C_| * |parametricEdgeID |
// |. . .| |1 |
//
float2 B_ = B * (numParametricSegments * 2.0);
float2 C_ = C * (numParametricSegments * numParametricSegments);
// Run a binary search to determine the highest parametric edge that is located on or before
// the combinedEdgeID. A combined ID is determined by the sum of complete parametric and
// radial segments behind it. i.e., find the highest parametric edge where:
//
// parametricEdgeID + floor(numRadialSegmentsAtParametricT) <= combinedEdgeID
//
float lastParametricEdgeID = 0.0;
float maxParametricEdgeID = min(numParametricSegments - 1.0, combinedEdgeID);
float negAbsRadsPerSegment = -abs(radsPerSegment);
float maxRotation0 = (1.0 + combinedEdgeID) * abs(radsPerSegment);
for (int exp = 5 /*max resolve level*/ - 1; exp >= 0; --exp) {
// Test the parametric edge at lastParametricEdgeID + 2^exp.
float testParametricID = lastParametricEdgeID + exp2(float(exp));
if (testParametricID <= maxParametricEdgeID) {
float2 testTan = fma(float2(testParametricID), A, B_);
testTan = fma(float2(testParametricID), testTan, C_);
float cosRotation = dot(normalize(testTan), tan0);
float maxRotation = fma(testParametricID, negAbsRadsPerSegment, maxRotation0);
maxRotation = min(maxRotation, $PI);
// Is rotation <= maxRotation? (i.e., is the number of complete radial segments
// behind testT, + testParametricID <= combinedEdgeID?)
if (cosRotation >= cos(maxRotation)) {
// testParametricID is on or before the combinedEdgeID. Keep it!
lastParametricEdgeID = testParametricID;
}
}
}
// Find the T value of the parametric edge at lastParametricEdgeID.
float parametricT = lastParametricEdgeID / numParametricSegments;
// Now that we've identified the highest parametric edge on or before the
// combinedEdgeID, the highest radial edge is easy:
float lastRadialEdgeID = combinedEdgeID - lastParametricEdgeID;
// Find the angle of tan0, i.e. the angle between tan0 and the positive x axis.
float angle0 = acos(clamp(tan0.x, -1.0, 1.0));
angle0 = tan0.y >= 0.0 ? angle0 : -angle0;
// Find the tangent vector on the edge at lastRadialEdgeID. By construction it is already
// normalized.
float radialAngle = fma(lastRadialEdgeID, radsPerSegment, angle0);
tangent = float2(cos(radialAngle), sin(radialAngle));
float2 norm = float2(-tangent.y, tangent.x);
// Find the T value where the tangent is orthogonal to norm. This is a quadratic:
//
// dot(norm, Tangent_Direction(T)) == 0
//
// |T^2|
// norm * |A 2B C| * |T | == 0
// |. . .| |1 |
//
float a=dot(norm,A), b_over_2=dot(norm,B), c=dot(norm,C);
float discr_over_4 = max(b_over_2*b_over_2 - a*c, 0.0);
float q = sqrt(discr_over_4);
if (b_over_2 > 0.0) {
q = -q;
}
q -= b_over_2;
// Roots are q/a and c/q. Since each curve section does not inflect or rotate more than 180
// degrees, there can only be one tangent orthogonal to "norm" inside 0..1. Pick the root
// nearest .5.
float _5qa = -.5*q*a;
float2 root = (abs(fma(q,q,_5qa)) < abs(fma(a,c,_5qa))) ? float2(q,a) : float2(c,q);
float radialT = (root.t != 0.0) ? root.s / root.t : 0.0;
radialT = clamp(radialT, 0.0, 1.0);
if (lastRadialEdgeID == 0.0) {
// The root finder above can become unstable when lastRadialEdgeID == 0 (e.g., if
// there are roots at exatly 0 and 1 both). radialT should always == 0 in this case.
radialT = 0.0;
}
// Now that we've identified the T values of the last parametric and radial edges, our final
// T value for combinedEdgeID is whichever is larger.
float T = max(parametricT, radialT);
// Evaluate the cubic at T. Use De Casteljau's for its accuracy and stability.
float2 ab = unchecked_mix(p0, p1, T);
float2 bc = unchecked_mix(p1, p2, T);
float2 cd = unchecked_mix(p2, p3, T);
float2 abc = unchecked_mix(ab, bc, T);
float2 bcd = unchecked_mix(bc, cd, T);
float2 abcd = unchecked_mix(abc, bcd, T);
// Evaluate the conic weight at T.
float u = unchecked_mix(1.0, w, T);
float v = w + 1 - u; // == mix(w, 1, T)
float uv = unchecked_mix(u, v, T);
// If we went with T=parametricT, then update the tangent. Otherwise leave it at the radial
// tangent found previously. (In the event that parametricT == radialT, we keep the radial
// tangent.)
if (T != radialT) {
// We must re-normalize here because the tangent is determined by the curve coefficients
tangent = w >= 0.0 ? robust_normalize_diff(bc*u, ab*v)
: robust_normalize_diff(bcd, abc);
}
strokeCoord = (w >= 0.0) ? abc/uv : abcd;
} else {
// Edges at the beginning and end of the strip use exact endpoints and tangents. This
// ensures crack-free seaming between instances.
tangent = (combinedEdgeID == 0) ? tan0 : tan1;
strokeCoord = (combinedEdgeID == 0) ? p0 : p3;
}
// At this point 'tangent' is normalized, so the orthogonal vector is also normalized.
float2 ortho = float2(tangent.y, -tangent.x);
strokeCoord += ortho * (strokeRadius * strokeOutset);
if (isHairline) {
// Hairline case. The scale and skew already happened before tessellation.
// TODO: There's probably a more efficient way to tessellate the hairline that lets us
// avoid inverting the affine matrix to get back to local coords, but it's just a 2x2 so
// this works for now.
return float4(strokeCoord + translate, inverse(affineMatrix) * strokeCoord);
} else {
// Normal case. Do the transform after tessellation.
return float4(affineMatrix * strokeCoord + translate, strokeCoord);
}
}