[HN Gopher] Fast calculation of the distance to cubic Bezier cur...
___________________________________________________________________
Fast calculation of the distance to cubic Bezier curves on the GPU
Author : ux
Score : 111 points
Date : 2025-10-18 09:25 UTC (13 hours ago)
(HTM) web link (blog.pkh.me)
(TXT) w3m dump (blog.pkh.me)
| amelius wrote:
| 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.
| jongjong wrote:
| Ok now all we need is an AI transformer model which can handle
| bezier curve embeddings as primitives instead of vector
| embeddings.
| Lichtso wrote:
| > The next step is to work with chains of Bezier 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...
| ux wrote:
| 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 "a3-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?
| Lichtso wrote:
| > "a3-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_rend
| erer/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 Bezier curves."
| https://arxiv.org/abs/1605.08669
| ux wrote:
| Thanks, I'll look into this. BTW, your 2nd shadertoy link
| is off (maybe it's private? Edit: seems you fixed it,
| thanks)
| vlovich123 wrote:
| But real shadows and lighting would require the distance
| aspect, no? The distance is only irrelevant for plain 2D text
| rendering, right?
| Lichtso wrote:
| > 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.
| vlovich123 wrote:
| 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?
| Lichtso wrote:
| 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.
| cubefox wrote:
| I wonder which method Apple is using for their recently
| introduced Bezier 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...
| 0xml wrote:
| 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 Bezier 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
| jasonjmcghee wrote:
| Your link 404s- private repo?
| 0xml wrote:
| Oops, updated.
| jasonjmcghee wrote:
| 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
| WithinReason wrote:
| this is a good question since for font rendering the length of
| each curve will usually be only a couple of pixels long
| GistNoesis wrote:
| 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.
| ux wrote:
| 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?
| GistNoesis wrote:
| 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/vig
| nettes/trus... but in this specific case of distances to
| curve we are on the easy case of Positive Definite.
| ux wrote:
| 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.
| GistNoesis wrote:
| 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/1ad0e50e8505fd191b6179
| 03b...
|
| 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.
| atilimcetin wrote:
| 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.
| phkahler wrote:
| This looks like new developments. Has any of it been applied to
| rational beziers of degree 3?
| mmorse1217 wrote:
| 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...
| ux wrote:
| 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".
___________________________________________________________________
(page generated 2025-10-18 23:01 UTC)