Hey, thanks for the nice post. I really enjoyed reading it; it’s good see this kind of thing on the front page.

Since you’re interested in doing this on GPU, an approach that might be interesting to you (although not necessarily more efficient) would be to leverage the intrinsic properties of Bezier curves to feed a near-optimal initial guess to Newton. Some useful facts about Bezier curves: i) Bezier control points form a convex hull of the curve they define ii) Bezier curves defined on [0,1] can be split into two bezier curves, each defined on [0,t] and [t,1] that define the same curve, with a tighter control polygon. iii) This Bezier curve splitting can be done using repeated linear combinations of Bezier control points, so you can skip evaluating Bernstein polynomials directly. iv) there is a mapping from Bezier control points to their corresponding value in the [0,1] parameter space (the term for this for B-Splines is greville abcissae, I’m not sure that there is an explicit name for the equivalent for Bezier curves, but basically the preimage of control point b_i of a degree d curve is i/d, i=0,…,d+1).

These things together sort of imply an algorithm: 1. Subdivide the Bezier curve c into 2 or 3 curves c_1, c_2, c_3 2. Find the closest control point b_j to the target point x 3. Choose the curve c_i corresponding to b_j: this subcurve contains the closest point to x 4. Go to step 1 and repeat this loop several times with c = c_i 5. Then, compute the preimage of the closest control point b_j to x on c (j/d plus some shift and rescaling). This value, t’, will be the initial guess to Newton’s method. 6. Solve for the closest point on the selected subcurve c to x with Newton’s method; this should converge in very few steps because your initial guess is so good, quadratic convergence, blah, blah blah.

The break-even point for this kind of algorithm vs. a derivative based algorithm is very unclear on CPU. But, for GPU, I think the computation can be structured in an architecture friendly way; since computing the euclidean distance between x with all control points and the bezier curve splitting can written in a vectorizable manner, you will probably see a decent speed up. I’ve only really worked with CUDA though, so I’m not sure if this idea maps very cleanly to GLSL.

Here’s an example of the algorithm above for CPU if you are interested: https://github.com/qnzhou/nanospline/commit/5ac97722414dbc75...

> The next step is to work with chains of Bézier curves to make up complex shapes (such as font glyphs). It will lead us to build a signed distance field. This is not trivial at all and mandates one or several dedicated articles. We will hopefully study these subjects in the not-so-distant future.

If you only want to fill a path of bezier curves (e.g. for text rendering) you can do without the "distance" part from "signed distance field" [0], leaving you with a "signed field" aka. an implicit curve [1].

Meaning not having to calculate the exact distance but only the sign (inside or outside) can be done without all the crazy iterative root finding in an actually cheap manner with only four multiplications and one addition per pixel / fragment / sample for a rational cubic curve [3].

[0]: https://en.wikipedia.org/wiki/Signed_distance_function

[1]: https://en.wikipedia.org/wiki/Implicit_curve

[2]: https://github.com/Lichtso/contrast_renderer/blob/a189d64a13...

Last time I found a paper in Graphics Gems titled Solving the Nearest-Point-on-Curve Problem, which transforms the problem into a Bernstein polynomial form. Then an exact solution can be obtained using A Bézier Curve-Based Root-Finder. This is my implementation [1], but it's not very robust for high-degree cases.

[1] https://github.com/Long0x0/distance-to-bezier

Cool. Maybe the next step can be to compute a set of bezier curves that (by a good approximation) closely cover the point-set that is exactly a given distance away from a given set of bezier curves.

Possibly naive question, but at least in the context of using distance fields to store font glyphs, what's the cost of the analytical solution (distance field of N combined bezier curves) vs rasterize at "high enough" resolution and then perform jump flood

I had to deal with the same problem on a GPU once. What worked for me was subdividing the cubic Bezier curve into smaller quadratic ones and then finding the roots of a cubic polynomial for each.

This looks like new developments. Has any of it been applied to rational beziers of degree 3?

This is fundamentally a geometric problem and the author completely missed the geometry aspect by transforming everything into polynomials and root finding.

The naive generic way of finding distances from point P to curve C([0,1]) is a procedure quite standard for global minimization : repeat "find a local minimum on the space constrained to be better than any previous minimum"

- Find a point P0 local minimum of d(P,P0) subject to the constraint P0=C(t0) (aka P0 is on C)

- define d0 = d(P,P0)

- Bisect the curve at t1 into two curves C1 = C([0,t0]) and C2([t0,1])

- Find a point P1 local minimum of d(P,P1) subject to the constraint P1=C(t1) with 0< t1 < t0 (aka P1 on c1) and the additional constraint d(P,P1) < d0 (note here the inequality is strict so that we won't find P0 again) if it exist.

- define d1 = d(P,P1)

- Find a point P2 local minimum of d(P,P2) subject to the constraint P2=C(t2) with t0 < t2 < 1 (aka P2 on c2) and the additional constraint d(P,P2) < d0 (note here the inequality is strict so that we won't find P0 again) if it exist.

- define d2 = d(P,P2)

Here in the case of cubic Bezier the curve have only one loop so you don't need to bissect anymore. If curves where higher order like spirals, you would need to cascade the problems with a shrinking distance constraint (aka looking only for points in R^2 inside the circle

So the distance is min over all local minimum found = min( d0,d1,d2)

Here the local minimization problems are in R^n (and inside disks centered around P) with n the dimension of the space, and because this is numerical approximation, you can either use slack (dual) variables to find the tx which express the on Curve constraint or barrier methods to express the disk constraints once a t parametrization has been chosen.

Multivariate Newton is guaranteed to work because distances are positive so the Sequential Quadratic Programming problems are convex (scipy minimize "SLSQP"). (Whereas the author needed 5 polynomials root, you can only need to solve for 3 points because each solve solves for two coordinates).

A local minimum is a point which satisfy the KKT conditions.

This procedure is quite standard : it's for example use to find the eigenvalues of matrices iteratively. Or finding all solutions to a minimization problem.

Where this procedure shines, is when you have multiple splines, and you want to find the minimal distance to them : you can partition the space efficiently and not compute distances to part of curves which are to far away to have a chance to be a minimum. This will scale independently of the resolution of your chain spline. (imagine a spiral the number of local minimum you need to encounter are proportional to the complexity of the scene and not the resolution of the spiral)

But when you are speaking about computing the whole signed distance field, you should often take the step of solving the Eikonal equation over the space instead of computing individual distances.

Ok now all we need is an AI transformer model which can handle bezier curve embeddings as primitives instead of vector embeddings.