I stumbled upon this interesting paper that demonstrates a new way to evaluate Bézier curves using additional points which are called Seiler points. Then using this points and just three linear interpolations the curve can be evaluated. By storing the Seiler points instead of the middle two Bézier control points it becomes an efficient way to evaluate the curves.

Given a cubic Bézier curve with control points \(b_0\), \(b_1\), \(b_2\) and $b_3$ the Seiler points \(s_1\) and \(s_2\) can be computed in the following way: \(s_1 = 3 b_1 - b_0 - b_3,\) and \(s_2 = 3 b_2 - b_3 - b_0.\) Then a curve evaluation \(C(t)\) boils down to:

\[b_{03} = \texttt{lerp}(b_0, b_3, t),\] \[s_{12} = \texttt{lerp}(s_1, s_2, t),\] \[C(t) = \texttt{lerp}(b_{03}, s_{12}, (1-t)t)\]

Geometrically this can be interpreted in the following way:

Taken from https://www.cemyuksel.com/research/seilers_interpolation/

\(b_{03}\) is a point on \(b_0b_3\) and \(s_{12}\) is point on \(s_1s_2\). Then the final evaluated point is another linear combination between those points. From this it seems that evaluating with Seiler points is akin to evaluating a bilinear surface. With this it seems possible to efficiently evaluate a Bezier curve in shaders. We can use the quadrilateral \(b_0b_1s_2s_1\) as your bilinear surface. Then given your pixel position \(p\) invert it to obtain bilinear coordinates \(u\) and \(v\). With this we can easily check if we are evaluating a point on the curve because only on the curve it holds that:

\(t = u\) and \((1-t)t = v.\)

Then we can substitute and check if

\((1-u)u = v\).

However, chances are not very likely for that to happen for any given pixel so we need to add some \(\epsilon\) to this comparison. We can make it so:

\[d = \texttt{abs}(v - (1-u)u)\]

then we can check whether \(d < \epsilon\) to determine if we are close enough to the curve. The effect of this is shown in the following ShaderToy:

Unfortunately, this does not give the result that one wants. The thickness is not uniform along the curve and for inflected curves it seems to dissapear altogether. This seems to be because we are evaluating based on the \(uv\) coordinates of the Seiler quad. Still, it seems to be useful in drawing the area under (or above) the curve as checking whether \(v-(1-u)u < 0\) will give this quite nicely.