> 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].
I wonder which method Apple is using for their recently introduced Bézier curve primitives for real-time 3D rendering in Metal. From their WWDC 2023 presentation [1]:
> Geometry such as hair, fur, and vegetation can have thousands or even millions of primitives. These are typically modeled as fine, smooth curves. Instead of using triangles to approximate these curves, you can use Metal's new curve primitives. These curves will remain smooth even as the camera zooms in. And compared to triangles, curves have a more compact memory footprint and allow faster acceleration structure builds.
> A full curve is made of a series of connected curve segments. Every segment on a curve is its own primitive, and Metal assigns each segment a unique primitive ID. Each of these segments is defined by a series of control points, which control the shape of the curve. These control points are interpolated using a set of basis functions. Depending on the basis function, each curve segment can have 2, 3, or 4 control points. Metal offers four different curve basis functions: Bezier, Catmull-Rom, B-Spline, and Linear. (...)
Finding the sign of the distance has been extremely challenging to me in many ways, so I'm very curious about the approach you're presenting. The snippet you shared has a "a³-bcd ≤ 0" formula which is all I get without more context. Can you elaborate on it or provide resources?
The winding number logic is usually super involved, especially when multiple sub-shapes start overlap and subtracting each other. Is this covered or orthogonal to what you are talking about?
> The winding number logic is usually super involved, especially when multiple sub-shapes start overlap and subtracting each other. Is this covered or orthogonal to what you are talking about?
Orthogonal: The implicit curve only tells you if you are inside or outside (the sign of the SDF), so that is technically sufficient, but usually you want more things: Some kind of anti-aliasing, composite shapes of more than one bezier curve and boolean operators for masking / clipping. Using the stencil buffer for counting the winding number allows to do all of that very easily without tessellation or decomposition at path intersections.
> Can you elaborate on it or provide resources?
If you are interested in the theory behind implicit curve rendering and how to handle the edge cases of cubic bezier curves checkout these papers:
> The distance is only irrelevant for plain 2D text rendering, right?
Yes, as I said it is relevant for text rendering, but not necessarily 2D. It can also be embedded in a 3D perspective as long as the text itself is planar. Meaning you can directly render text in a 3D scene this way without rendering to a texture first.
> But real shadows and lighting would require the distance aspect, no?
I think the difference is in stroke vs fill, not the illumination (as you could still use shadow mapping / projection). In stroking you need to calculate an offset curve either explicitly or implicitly sample it from a signed distance field. Thus the exact distance matters for stroking, for filling it does not.
Couldn’t you do stroking by doing a second fill operation on a slightly scaled down version of the first with the negative space color as the interior?
Yep, stroking is just filling of the space between offset curves (aka. parallel curves), and that "slightly scaled down version of the first" is the "calculate an offset curve explicitly" approach I mentioned.
Though it is very unpractical because the offset curve of a cubic bezier curve is not a cubic bezier curve anymore, instead it is an analytic curve of degree 10. Thus, in practice the offset curves for stroking are either approximated by polygons or implicitly sampled from signed distance fields.
One more thing: Offset curves are different form classical scaling from a center point in all but the most trivial cases where there exists such a center; namely regular polygons. And cubic bezier curves can be concave, even have a self intersecting loop or form a cusp.
Pretty cool. I haven't dived into the details of this post, but it's a problem I've been wrangling with from time to time in my game.
Below is a series of writings on this topic that I enjoyed, by James Blinn, of "Blinn-Phong reflection" fame (and more). Not state-of-the-art, but an interesting read. It's just for cubics, which is what you need to solve the distance formula for quadratic bezier curves (my particular case), rather than the harder cubic curves of the linked article.
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.
Thanks, that the same approach I was suggesting in my other comment in this thread https://news.ycombinator.com/item?id=45628245 . But couldn't find literature specific to the Bezier curves to help break the communication gap, and my specific knowledge of Bezier Curve isn't deep enough.
I am happy to see other people use the approach I consider more natural.
It's a generic global optimization approach and geometry is always full of pathological edge cases, so it's hard to tell if you miss any. Getting to work in the average case is usually easy, but to be sure it always work is much harder.
Your comment motivated me to also comment :) I agree: Bézier curves and b-splines have a lot of rich geometry baked it, so it makes sense to use it, especially if you can avoid second derivatives. I think I misunderstood your comment about the extra control point work to find an initial guess. It looks like we’re saying something pretty close though: you can split the curve smartly to remove the single crossing in this case, which skips all the recursion that I’m suggesting, which is probably better.
Yes the real trouble is true optimality guarantees. I remember that there were edge cases of the above approach needed that a lot of subdivision steps to succeed for general degree curves, so it might end up worse than a more rigorously justified approach in these cases.
Thank you! Do we have a guarantee that these subcurves are solvable with Newton's method? The approach with derivatives has this because we know there is one crossing, and also clipping them to the zero-derivatives makes sure there won't be multiple curve "pits".
The trick is that you only solve Newton on a single subcurve: you recursively find the closest control point on all subcurves, split that subcurve and repeat until the closest control point doesn’t move much, or just however many subdivision steps work in practice. So the last curve that you apply Newton to should be smooth enough to succeed. I think there are edge cases with cusps, but I can’t exactly remember the theoretical guarantees anymore.
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.
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.
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.
This is quite a complicated approach that misses out on most of the basic numerical and geometric methods that will make the problem simpler. I would recommend looking outside of recent SIGGRAPH papers. Brushing up on the basics of Bernstein polynomials, B-splines, and rootfinding methods will lead you to develop simpler algorithms which are likely faster.
I'm extremely curious what those basic methods are. We're in the process of replacing the higher order rootfinding in kurbo with a new solver based on Yuksel's method[1]. If you know of simpler, faster techniques that would be quite interesting.
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
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.
I'm interested in the approach you're describing but it's hard to follow a comment in the margin. Is there a paper or an implementation example somewhere?
The general technique is not recent I was taught it in school in global optimisation class more than 15 years ago.
Here there is a small number of local minimum, the idea is to iterate over them in increasing order.
Can't remember the exact name but here is a more recent paper proposing "Sequential Gradient Descent" https://arxiv.org/abs/2011.04866 which features a similar idea.
There is not really something special to it, it just standard local non linear minimization techniques with constraints Sequential Least Squares Quadratic Programming (SLSQP).
It's just about framing it as an optimization problem looking for "Points" with constraints and applying standard optimization toolbox, and recognizing which type of problem your specific problem is. You can write it as basic gradient descent if you don't care about performance.
The problem of finding a minimum of a quadratic function inside a disk is commonly known as the "Trust Region SubProblem" https://cran.r-project.org/web/packages/trust/vignettes/trus... but in this specific case of distances to curve we are on the easy case of Positive Definite.
What you described in your first message seemed similar to the approach used in the degree N root solving algorithm by Cem Yuksel; splitting the curve in simpler segments, then bisect into them. I'd be happy to explore what you suggested, but I'm not mathematically literate, so I'll be honest with you; what you're saying here is complete gibberish to me, and it's very hard to follow your point. It will take me weeks to figure out your suggestion and make a call as to whether it's actually simpler or more performant than what is proposed in the article.
I have written some gist to illustrate the approach I suggest. The code run but there may be bugs, and it don't use the appropriate optimizer. The purpose is to illustrate the optimisation approach.
The approach I suggest will need more work, and there are probably problematic edge cases to consider and numerical stability issues. Proper proofs have not been done. It's typically some high work-low reward situation.
It's mostly interesting because it highlight the link between roots and local mimimum. And because it respect the structure of the problem more.
To find roots we can find a first root then divide the polynomial by (x-root). And find a root again.
If you are not mathematically literate, it'll probably be hard to do the details necessary to make it performant. But if you use a standard black-box optimizer with constraints it should be able to do it in few iterations.
You can simplify the problem by considering piece-wise segments instead of splines. The extension to chains of segment is roughly the same, and the spatial acceleration structure based on branch-and-bound are easier.
> 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...
I wonder which method Apple is using for their recently introduced Bézier curve primitives for real-time 3D rendering in Metal. From their WWDC 2023 presentation [1]:
> Geometry such as hair, fur, and vegetation can have thousands or even millions of primitives. These are typically modeled as fine, smooth curves. Instead of using triangles to approximate these curves, you can use Metal's new curve primitives. These curves will remain smooth even as the camera zooms in. And compared to triangles, curves have a more compact memory footprint and allow faster acceleration structure builds.
> A full curve is made of a series of connected curve segments. Every segment on a curve is its own primitive, and Metal assigns each segment a unique primitive ID. Each of these segments is defined by a series of control points, which control the shape of the curve. These control points are interpolated using a set of basis functions. Depending on the basis function, each curve segment can have 2, 3, or 4 control points. Metal offers four different curve basis functions: Bezier, Catmull-Rom, B-Spline, and Linear. (...)
1: https://developer.apple.com/videos/play/wwdc2023/10128/?time...
Finding the sign of the distance has been extremely challenging to me in many ways, so I'm very curious about the approach you're presenting. The snippet you shared has a "a³-bcd ≤ 0" formula which is all I get without more context. Can you elaborate on it or provide resources?
The winding number logic is usually super involved, especially when multiple sub-shapes start overlap and subtracting each other. Is this covered or orthogonal to what you are talking about?
> "a³-bcd ≤ 0" formula
These are the coefficients of the implicit curve, finding them can be done once upfront.
For integral quadratic bezier curves that is trivial as they are constant, see: https://www.shadertoy.com/view/fsXcDj
For rational cubic bezier curves it is more involved, see: https://www.shadertoy.com/view/ttVfzh
And for the full complexity of dealing with self intersecting loops and cusps see: https://github.com/Lichtso/contrast_renderer/blob/main/src/f...
> The winding number logic is usually super involved, especially when multiple sub-shapes start overlap and subtracting each other. Is this covered or orthogonal to what you are talking about?
Orthogonal: The implicit curve only tells you if you are inside or outside (the sign of the SDF), so that is technically sufficient, but usually you want more things: Some kind of anti-aliasing, composite shapes of more than one bezier curve and boolean operators for masking / clipping. Using the stencil buffer for counting the winding number allows to do all of that very easily without tessellation or decomposition at path intersections.
> Can you elaborate on it or provide resources?
If you are interested in the theory behind implicit curve rendering and how to handle the edge cases of cubic bezier curves checkout these papers:
Loop, Charles, and Jim Blinn. "Resolution independent curve rendering using programmable graphics hardware." https://www.microsoft.com/en-us/research/wp-content/uploads/...
BARROWCLOUGH, Oliver JD. "A basis for the implicit representation of planar rational cubic Bézier curves." https://arxiv.org/abs/1605.08669
Thanks, I'll look into this. BTW, your 2nd shadertoy link is off (maybe it's private? Edit: seems you fixed it, thanks)
But real shadows and lighting would require the distance aspect, no? The distance is only irrelevant for plain 2D text rendering, right?
> The distance is only irrelevant for plain 2D text rendering, right?
Yes, as I said it is relevant for text rendering, but not necessarily 2D. It can also be embedded in a 3D perspective as long as the text itself is planar. Meaning you can directly render text in a 3D scene this way without rendering to a texture first.
> But real shadows and lighting would require the distance aspect, no?
I think the difference is in stroke vs fill, not the illumination (as you could still use shadow mapping / projection). In stroking you need to calculate an offset curve either explicitly or implicitly sample it from a signed distance field. Thus the exact distance matters for stroking, for filling it does not.
Couldn’t you do stroking by doing a second fill operation on a slightly scaled down version of the first with the negative space color as the interior?
Yep, stroking is just filling of the space between offset curves (aka. parallel curves), and that "slightly scaled down version of the first" is the "calculate an offset curve explicitly" approach I mentioned.
Though it is very unpractical because the offset curve of a cubic bezier curve is not a cubic bezier curve anymore, instead it is an analytic curve of degree 10. Thus, in practice the offset curves for stroking are either approximated by polygons or implicitly sampled from signed distance fields.
Raph Levien has a good blog post about it:
https://raphlinus.github.io/curves/2022/09/09/parallel-bezie...
One more thing: Offset curves are different form classical scaling from a center point in all but the most trivial cases where there exists such a center; namely regular polygons. And cubic bezier curves can be concave, even have a self intersecting loop or form a cusp.
Pretty cool. I haven't dived into the details of this post, but it's a problem I've been wrangling with from time to time in my game.
Below is a series of writings on this topic that I enjoyed, by James Blinn, of "Blinn-Phong reflection" fame (and more). Not state-of-the-art, but an interesting read. It's just for cubics, which is what you need to solve the distance formula for quadratic bezier curves (my particular case), rather than the harder cubic curves of the linked article.
* https://courses.cs.washington.edu/courses/cse590b/13au/lectu...
* https://courses.cs.washington.edu/courses/cse590b/13au/lectu...
* https://courses.cs.washington.edu/courses/cse590b/13au/lectu...
* https://courses.cs.washington.edu/courses/cse590b/13au/lectu...
* https://courses.cs.washington.edu/courses/cse590b/13au/lectu...
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...
Thanks, that the same approach I was suggesting in my other comment in this thread https://news.ycombinator.com/item?id=45628245 . But couldn't find literature specific to the Bezier curves to help break the communication gap, and my specific knowledge of Bezier Curve isn't deep enough.
I am happy to see other people use the approach I consider more natural.
It's a generic global optimization approach and geometry is always full of pathological edge cases, so it's hard to tell if you miss any. Getting to work in the average case is usually easy, but to be sure it always work is much harder.
Your comment motivated me to also comment :) I agree: Bézier curves and b-splines have a lot of rich geometry baked it, so it makes sense to use it, especially if you can avoid second derivatives. I think I misunderstood your comment about the extra control point work to find an initial guess. It looks like we’re saying something pretty close though: you can split the curve smartly to remove the single crossing in this case, which skips all the recursion that I’m suggesting, which is probably better.
Yes the real trouble is true optimality guarantees. I remember that there were edge cases of the above approach needed that a lot of subdivision steps to succeed for general degree curves, so it might end up worse than a more rigorously justified approach in these cases.
Thank you! Do we have a guarantee that these subcurves are solvable with Newton's method? The approach with derivatives has this because we know there is one crossing, and also clipping them to the zero-derivatives makes sure there won't be multiple curve "pits".
The trick is that you only solve Newton on a single subcurve: you recursively find the closest control point on all subcurves, split that subcurve and repeat until the closest control point doesn’t move much, or just however many subdivision steps work in practice. So the last curve that you apply Newton to should be smooth enough to succeed. I think there are edge cases with cusps, but I can’t exactly remember the theoretical guarantees anymore.
I think this is the main reference for this algorithm (should be able to search for the pdf): https://www.sciencedirect.com/science/article/abs/pii/S01678...
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
Your link 404s- private repo?
Oops, updated.
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.
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.
This is quite a complicated approach that misses out on most of the basic numerical and geometric methods that will make the problem simpler. I would recommend looking outside of recent SIGGRAPH papers. Brushing up on the basics of Bernstein polynomials, B-splines, and rootfinding methods will lead you to develop simpler algorithms which are likely faster.
I'm extremely curious what those basic methods are. We're in the process of replacing the higher order rootfinding in kurbo with a new solver based on Yuksel's method[1]. If you know of simpler, faster techniques that would be quite interesting.
[1]: https://crates.io/crates/polycool
This is a pretty good foundational reference for Bézier curves and b-splines https://www.sciencedirect.com/book/9780122490521/curves-and-...
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
this is a good question since for font rendering the length of each curve will usually be only a couple of pixels long
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.
I'm interested in the approach you're describing but it's hard to follow a comment in the margin. Is there a paper or an implementation example somewhere?
The general technique is not recent I was taught it in school in global optimisation class more than 15 years ago.
Here there is a small number of local minimum, the idea is to iterate over them in increasing order.
Can't remember the exact name but here is a more recent paper proposing "Sequential Gradient Descent" https://arxiv.org/abs/2011.04866 which features a similar idea.
Sequential convex programming : http://web.stanford.edu/class/ee364b/lectures/seq_notes.pdf
There is not really something special to it, it just standard local non linear minimization techniques with constraints Sequential Least Squares Quadratic Programming (SLSQP).
It's just about framing it as an optimization problem looking for "Points" with constraints and applying standard optimization toolbox, and recognizing which type of problem your specific problem is. You can write it as basic gradient descent if you don't care about performance.
The problem of finding a minimum of a quadratic function inside a disk is commonly known as the "Trust Region SubProblem" https://cran.r-project.org/web/packages/trust/vignettes/trus... but in this specific case of distances to curve we are on the easy case of Positive Definite.
What you described in your first message seemed similar to the approach used in the degree N root solving algorithm by Cem Yuksel; splitting the curve in simpler segments, then bisect into them. I'd be happy to explore what you suggested, but I'm not mathematically literate, so I'll be honest with you; what you're saying here is complete gibberish to me, and it's very hard to follow your point. It will take me weeks to figure out your suggestion and make a call as to whether it's actually simpler or more performant than what is proposed in the article.
I have written some gist to illustrate the approach I suggest. The code run but there may be bugs, and it don't use the appropriate optimizer. The purpose is to illustrate the optimisation approach.
https://gist.github.com/unrealwill/1ad0e50e8505fd191b617903b...
Point 33 "intersection between bezier curve with a circle" may be useful to find the feasible regions of the subproblems https://pomax.github.io/bezierinfo/#introduction
The approach I suggest will need more work, and there are probably problematic edge cases to consider and numerical stability issues. Proper proofs have not been done. It's typically some high work-low reward situation.
It's mostly interesting because it highlight the link between roots and local mimimum. And because it respect the structure of the problem more.
To find roots we can find a first root then divide the polynomial by (x-root). And find a root again.
If you are not mathematically literate, it'll probably be hard to do the details necessary to make it performant. But if you use a standard black-box optimizer with constraints it should be able to do it in few iterations.
You can simplify the problem by considering piece-wise segments instead of splines. The extension to chains of segment is roughly the same, and the spatial acceleration structure based on branch-and-bound are easier.
Would an undergrad-level algorithm be faster here?
1. Find the derivative, and use an analytic formula for the roots of a quartic equation.
2. Then evaluate the value at the zeros to look for the intervals that cross the [0,1] segment.
[dead]
Ok now all we need is an AI transformer model which can handle bezier curve embeddings as primitives instead of vector embeddings.