1 #version 310 es
2 3 // SPDX-License-Identifier: Unlicense OR MIT
4 5 precision mediump float;
6 7 layout(location=0) in vec2 vFrom;
8 layout(location=1) in vec2 vCtrl;
9 layout(location=2) in vec2 vTo;
10 11 layout(location = 0) out vec4 fragCover;
12 13 void main() {
14 float dx = vTo.x - vFrom.x;
15 // Sort from and to in increasing order so the root below
16 // is always the positive square root, if any.
17 // We need the direction of the curve below, so this can't be
18 // done from the vertex shader.
19 bool increasing = vTo.x >= vFrom.x;
20 vec2 left = increasing ? vFrom : vTo;
21 vec2 right = increasing ? vTo : vFrom;
22 23 // The signed horizontal extent of the fragment.
24 vec2 extent = clamp(vec2(vFrom.x, vTo.x), -0.5, 0.5);
25 // Find the t where the curve crosses the middle of the
26 // extent, x₀.
27 // Given the Bézier curve with x coordinates P₀, P₁, P₂
28 // where P₀ is at the origin, its x coordinate in t
29 // is given by:
30 //
31 // x(t) = 2(1-t)tP₁ + t²P₂
32 //
33 // Rearranging:
34 //
35 // x(t) = (P₂ - 2P₁)t² + 2P₁t
36 //
37 // Setting x(t) = x₀ and using Muller's quadratic formula ("Citardauq")
38 // for robustnesss,
39 //
40 // t = 2x₀/(2P₁±√(4P₁²+4(P₂-2P₁)x₀))
41 //
42 // which simplifies to
43 //
44 // t = x₀/(P₁±√(P₁²+(P₂-2P₁)x₀))
45 //
46 // Setting v = P₂-P₁,
47 //
48 // t = x₀/(P₁±√(P₁²+(v-P₁)x₀))
49 //
50 // t lie in [0; 1]; P₂ ≥ P₁ and P₁ ≥ 0 since we split curves where
51 // the control point lies before the start point or after the end point.
52 // It can then be shown that only the positive square root is valid.
53 float midx = mix(extent.x, extent.y, 0.5);
54 float x0 = midx - left.x;
55 vec2 p1 = vCtrl - left;
56 vec2 v = right - vCtrl;
57 float t = x0/(p1.x+sqrt(p1.x*p1.x+(v.x-p1.x)*x0));
58 // Find y(t) on the curve.
59 float y = mix(mix(left.y, vCtrl.y, t), mix(vCtrl.y, right.y, t), t);
60 // And the slope.
61 vec2 d_half = mix(p1, v, t);
62 float dy = d_half.y/d_half.x;
63 // Together, y and dy form a line approximation.
64 65 // Compute the fragment area above the line.
66 // The area is symmetric around dy = 0. Scale slope with extent width.
67 float width = extent.y - extent.x;
68 dy = abs(dy*width);
69 70 vec4 sides = vec4(dy*+0.5 + y, dy*-0.5 + y, (+0.5-y)/dy, (-0.5-y)/dy);
71 sides = clamp(sides+0.5, 0.0, 1.0);
72 73 float area = 0.5*(sides.z - sides.z*sides.y + 1.0 - sides.x+sides.x*sides.w);
74 area *= width;
75 76 // Work around issue #13.
77 if (width == 0.0)
78 area = 0.0;
79 80 fragCover.r = area;
81 }
82