Creating a Bezier Curve SDF

An example of this SDF being used can be found here.

I should start by making it clear this absolutely is not the best way to implement this SDF. I would strongly recommend using Inigo Quilez's version, which you can find here. The math in this article is correct, to the best of my knowledge, but my implementation has some issues with floating point precision errors.

A second point before we begin: this SDF is for a quadratic Bezier curve. I think creating an SDF of a cubic Bezier is probably not feasible, since it will require finding the roots of a quartic function.

I'm going to assume you are somewhat familiar with Bezier curves. If you aren't, then I highly recommend Freya Holmér's video.

Starting with the usual lerp function

Given three points $\overrightarrow{A}\text{,}\overrightarrow{B}\text{,}\overrightarrow{C}$, a quadratic Bezier curve is defined as

If we expand out the lerps, we get

For convenience, and to save me from typing so much, let $\overrightarrow{M}\text{=}\overrightarrow{A}\text{- 2}\overrightarrow{B}\text{+}\overrightarrow{C}$, and $\overrightarrow{N}\text{= -2}\overrightarrow{A}\text{+ 2}\overrightarrow{B}$. So our function is now simply

Ok, we have our function. Now what? Our goal is that, given a point $\overrightarrow{p}$, we want to find the distance between $\overrightarrow{p}$ and the curve. To do this, we will first find the closest point on the curve to $\overrightarrow{p}$.

If we imagine a circle centred on our point $\overrightarrow{p}$, we can expand the circle outward until it intersects our curve at exactly one point. The curve will be tangent to the circle that that point since, if it wasn't, the curve would have to cross into the circle, and thus cross out of it again, meaning the circle would be intersecting the curve at two points.

Now consider a vector $\overrightarrow{d}$ that points from the point $\overrightarrow{p}$ to the point on the curve given by $\text{f(t)}$.

It follows that the tangent of the curve at its closest point to $\overrightarrow{p}$ will be orthogonal to $\overrightarrow{d}$.

The tangent of our curve is given by its derivative with respect to $\text{t}$.

And so, what we want is the value of $\text{t}$ where $\overrightarrow{d}\text{(t) · f '(t) = 0}$.

If you expand out the dot product, you get a cubic equation with the following coefficients

The coefficients are all scalar values, so our goal is now to find the roots of this cubic equation. I relied entirely on this Wikipedia article to learn how to do so. Deriving a way to find the roots of a cubic equation is well beyond me.

We need to manipulate our cubic equation into a "depressed" cubic of the form

We can ignore $\text{x}$ for now, since it's just t with an offset, which is what we're trying to find. $\text{m}$ and $\text{n}$ are defined as

There are two situations we need to consider. The first is when there is only a single root to our cubic equation, like in the image above. But consider a situation like this

In this case, there are three values for $\mathrm{t}$ where our cubic will be zero.

Fortunately, there is an easy way to know how many roots our cubic equation will have.

If the discriminant, $\text{dis}$, is greater than zero, then there's exactly one real root. Otherwise, there are three real roots, one of which might be a double root.

As stated earlier, to convert the root of our depressed cubic back to the root of our starting cubic, we just need to adjust the result by a constant value, $\text{adj}$.

Let's start with the case where there's only one root, since it's simpler. This is just copied from the previously linked wikipedia article, but split up to more closely match the GLSL code.

The $\text{clamp(v, min, max)}$ function constrains the input, $\mathrm{v}$, to be between the min and max values.

Now that we have our root, we just have to calculate the distance between $\text{f(t)}$ and $\overrightarrow{p}$

The case where there are three roots is more complex.

It is worth noting that, while this is a correct way to calculate θ, it tends to produce significant floating point errors. I'm not sure why that is, but my best guess is that it is a combination of $\mathrm{m}$ and $\mathrm{n}$ potentially being significantly different orders of magnitude, arccos only being approximated in GLSL, and WebGL using 16 bit floats, all together resulting in significant compounding errors.

Google Gemini gave me the following, which seems to produce better results by keeping the values of each step of the calculation closer in magnitude. I'll admit it's unlikely I would have been able to figure this out on my own, so chalk up a win for AI, much to my annoyance.

My typesetting is really falling apart here. Now we can find our roots.

And finally, we get the distance between our point and the point given by the closest root

Note that there is a case we didn't consider. If $\overrightarrow{A}\text{,}\overrightarrow{B}\text{,}\overrightarrow{C}$, are co-linear, we can end up with a quadratic equation rather than a cubic equation. In that case, we just use the SDF for a straight line with $\overrightarrow{A}$ and $\overrightarrow{C}$ as the end points.

Here is the final SDF written out in GLSL.

float lineSDF(vec2 p, vec2 A, vec2 B)
{
    p = p - A;
    B = B - A;

    float s = dot(p,B)/dot(B,B);

    return length(p - B * clamp(s, 0., 1.));
}

float bezierSDF(vec2 p, vec2 A, vec2 B, vec2 C)
{
    vec2 M = A - (2.0 * B) + C;
    vec2 N = (-2.0 * A) + (2.0 * B);

    float a = dot(-2.0 * M, M);

    // If 'a' is zero (or close enough to it), then we have a straight line. We have to handle that
    // differently, otherwise we'll get a division by zero
    if(abs(a) < 0.0001)
    {
        return lineSDF(p, A, C);
    }

    float b = -dot(N, M) - dot(2.0 * M, N);

    float c = dot(2.0 * M, p) - dot(N, N) - dot(2.0 * M, A);

    float d = dot(N, p) - dot(N, A);

    float m = (3.0 * a * c - b * b) / (3.0 * a * a);
    float n = (2.0 * b * b * b - 9.0 * a * b * c + 27.0 * a * a * d) / (27.0 * a * a * a);

    float dis = (n*n)/4. + (m*m*m)/27.;
    float adj = b/(3.*a);

    if(dis > 0.)
    {
        float s = sqrt(dis);
        float u1 = -n/2. + s;
        float u2 = -n/2. - s;
        float t = sign(u1)*pow(abs(u1), 1./3.) + sign(u2)*pow(abs(u2), 1./3.) - adj;
        t = clamp(t, 0., 1.);
        return distance(p, M*t*t + N*t + A);
    }
    else
    {
        //float ac = acos(((3.*n)/(2.*m)) * sqrt(-m/3.));
        float ac = acos(-n/(2.*sqrt(-m*m*m/27.)));

        float s = 2. * sqrt(-m/3.);
        float r1 = clamp(s * cos(ac/3.) - adj, 0., 1.);
        float r2 = clamp(s * cos((ac - 2.*PI)/3.) - adj, 0., 1.);
        float r3 = clamp(s * cos((ac - 4.*PI)/3.) - adj, 0., 1.);

        vec2 l1 = p - (M*r1*r1 + N*r1 + A);
        vec2 l2 = p - (M*r2*r2 + N*r2 + A);
        vec2 l3 = p - (M*r3*r3 + N*r3 + A);

        return sqrt(min(min(dot(l1,l1), dot(l2,l2)), dot(l3,l3)));
    }
}

Addendum

I mentioned at the beginning of the article that an SDF of a cubic Bezier curve is likely unfeasible. However, since then, I've started to think otherwise. We had to find the roots of a cubic function for the quadratic Bezier SDF. Similarly, we would need to find the roots of a quartic function for the SDF of a cubic Bezier. Wikipedia has the full equation written out. While it looks extremely long and complicated, it only uses square and cube roots, and large portions of it are reused several times. It also doesn't use the arccos function, which is where I think the issues with imprecision crop up.

Speaking of imprecision. While the equation given to me by Gemini produced better results, there are still issues when the control points of the curve in certain posistions. Parts of the curve will disappear, or it will be cut short.