How to raytrace Bezier surfaces?

I tried this question on math.SE and surprisingly, the answer was "the equations are too nasty, just feed the function it to a numerical root-finder". But if you consider yourself "a graphics guy" like me, and have played extensively with Bezier curves for design work, I got to believe that better can be done. There is a published algorithm by Kajiya that I don't have the background to understand (Sylvester Matrices), but the related advice on math.SE was that the result is a degree-18 polynomial in t, and you still need to solve that numerically. I had another idea with similar result.

So, is it a total pipe dream to hope to solve the Ray/Bezier-surface intersection algebraically, thus making it possible to code explicitly and have super-fast super-smoothness?

Barring that, what's the fastest method for performing this calculation? Can you "find the wiggles" to get a tight bound (and target) for recursive subdivision? If you have to use a numerical root-finder (sigh), what properties does it need and is there a best choice for speed?

My original thought was about preparing for a specific surface, similar to Laplace expansion as described in the answer to my other math question about triangles. But I'd be interested in general methods, too. I'm just thinking of a fixed set of shapes, like the Utah teapot. But I'd be very interested in ways of optimizing for temporal coherence across animated frames.

• Are you looking for a general method that you can apply to an arbitrary Bezier surface, or a way of preparing a fast method for a specific surface? Will your surface shape be fixed before runtime? Aug 30, 2015 at 22:05
• Note that you can raymarch bezier surfaces a lot easier than raytracing it. You can also raytrace or raymarch univariate surfaces a lot easier than other kinds! blog.demofox.org/2015/07/28/rectangular-bezier-patches Aug 31, 2015 at 4:16

First off, here's the Kajiya method I think you're thinking of: Kajiya, Ray Tracing Parametric Patches, SIGGRAPH 82. The tech report version might be more informative.

What I hope you get from that is that it's not impossible and it's not conceptually difficult if you don't mind getting your hands dirty with some algebraic geometry and complex numbers. However, doing it directly is absurdly expensive.

"Real" ray tracers tend to do some combination of two things:

• Placing a bounding hierarchy (e.g. AABBs) on the patch to get a good "initial value" for a numeric root finder. If you do this well, you can avoid the "wrinkle" problem.
• Tesselating the patch into DDG shells and ray tracing them like polygon meshes.

That last point sounds like it kills the "super-smoothness" requirement, but it isn't nearly as bad as that if you're using ray differentials. Matching the tessellation level to the "size" of the ray bounds the error nicely. Besides, you probably need differentials for texture coordinates anyway, so you might as well use it to control the accuracy of the intersection test too.

Exploiting temporal coherence isn't a bad idea, but exactly how you'd do that depends a lot on your scene graph representation. You might want to look at ray coherence. Ask your favourite search engine about ray packet tracing and ray reordering.

is it a total pipe dream to hope to solve the Ray/Bezier-surface intersection algebraically

Yes, it's a pipe dream. A bicubic Bezier patch is an algebraic surface of degree 18. To intersect a ray with this surface, you have to find the roots of a polynomial of degree 18. There is no formula for these roots -- you have to find them by numerical methods. In fact, there are mathematical results (the Abel-Ruffini theorem) telling us that there can never be formulae for roots of equations beyond degree 4. The math doesn't just say that the formulae have not yet been found; it says that they never will be found, because they can't exist.

If you really want to do analytical (algebraic) ray-tracing of curvaceous shapes, you might try using Steiner patches. These have degree 4, so the ray-patch intersection can be computed by finding the roots of a quartic (i.e. a polynomial of degree 4). There are formulae for finding roots of quartics, but they're pretty nasty, and it's surprisingly hard to write code that implements the formulae reliably.

• I would like to point out that it is not true that there can never be formulae for roots of polynomial equations of degree greater than 4. Although I do agree that it is likely not a good idea to try to find the intersections of a line with a bicubic Bezier patch algebraically, the Abel-Ruffini theorem is at most a partial explanation for this.
– user9485
Jul 19, 2020 at 14:30

Another option, which I used a couple of decades ago (yikes!), is to use Toth's scheme from 1985 that employed interval arithmetic to narrow down the search space. IIRC, eventually it will resort to Newton-Rhapson but, again IIRC, I think it rarely required more than one or two steps to get to a good solution.

Though I haven't looked at it (well, apart from a quick glance) Mitchell has published some more recent work on ray tracing with interval maths.

(I should add that, if you are only doing Bezier surfaces, then the interval method might be a bit "overkill" since you can use tricks like blossoming to get bounds and derivatives. If, however, you combine Bezier curves with other functions, e.g. rotation around an axis, then its generality is more useful.)

https://www.shadertoy.com/results?query=bezier sort by age, in case of compatibility issues:

,...shows many solutions of many spline-subsets, either returning the distance to a 2d spline, or tracing a 3d patch. Splines and patches come in many forms. heavensine beind simplest, bezier being simple, nurbs being overly complex. The more constrains yo add to your spline, the simpler it gets. NURBS is extension overkill; - its Non-Uniform-ness of weights (the "NU") diminishes efficiency in comparison to to more symmetric splines - its Ration-al-ness (the R) also adds some complexity, for segmenting (rationing) and mixing with nearby segments (recursively solved).

bezier-patch-tracing is root-solving and with that comes contextual prioritization on precision; in what order to solve the quadratic equation. this becomes impractical on higher exponents than cubic, due to exponential complexity and precision loss.

ray-marching==sphere-tracking is the simpler heuristic approach to root-solving, which seems to be the simple and most efficient solution to rendering most splines patches.

Lagrange-representation simplifies tracing/marching (as the L-points are ON the spline while the ControlVector-points (of the exact same spline) are rarely on the spline)

The special case of a heavensine-spline, where first derivatives of stat and end are ==0. simplifies continuity and involves less differentials (less subtraction). A heavensine-patch can be traced efficiently in a single pass: https://www.shadertoy.com/view/4djfW3 while other cubic (or higher) splines make the heuristic sphere-tracking/ray-marching approach more efficient (and "precise enough") than daring to analytically calculate all the roots to keep the smallest positive root (with exponentially accumulating precision errors for each root).

In computer graphics, splines and patches have been nearly completely replaced by z-brushing by 2006. z-brushing uses displacement maps with homogeneous coordinates, or even using a "type" that us a union of sphere and linesegments (linesegments have a radius of 0, spheres have a length of 0, a union is simple and useful). For a minor loss in precision for a big gain in performance at relatively low memory cost for a look-up-table, that is easily made dynamic on a gpu.

• Nevermind. all 3d patch solutions are done by sphere tracking.
– ollj
Aug 24, 2017 at 15:53
• it just significantly increases performance and precision when the patch is simpler. to a point where a patch of seamess heavensines gets you very far in a few iterations:
– ollj
Aug 24, 2017 at 15:56
– ollj
Aug 24, 2017 at 15:56

https://www.shadertoy.com/view/MldczM Is ANALYTIC intersection of a ray with a triangular b-spline, 3 quadratic bezier splines, 3 CVs for each spline, 3 corner CVs are shared by 2 splines. All cv heights are orthogonal to the triangle plane. It solves for <=2 roots in barycentric coordinates. There's never 3 intersections, so this is a surprisingly simple tracing of a curved function. The first derivative is super short, too.

ESDF keys for camera movement The parent shader in the comments has less optimized, less precise, more repetitive, naive simpler code.

It still has the error of rendering the negative surface if the camera is inside the prism of the triangle (above or below triangle). This seems to be a simpler special case, with only 1 root, where it fails to ignore the other root. The branch to solve this may be rather simple.

• Interesting, but to be honest, it's only one step up from a triangle, and IIRC (but unfortunately being out of the office I can't access my copy of Farin's "Curves and Surfaces for CAGD" to check) to get C1 (or G1?) continuity I believe you need a quartic Bezier triangle, i.e. 15 control points Jul 21, 2020 at 16:59

simple case == analytic: (4 roots or less)

For the simplest bezier patch; a base-triangle in 3d space (of any shape, just have 3 non-colinear-points in 3d.) (ideally all .z coordinates are == zero ... for triangle-normal-simplicity and faster speed (mote zeroes in matrices and dot-products), with 3 triangle-plane-normals on the triangle corners, and 3 triangle-plane-normals on the triangle-side-midpoints, having one of 6 CV-heights on each of those 6 plane-normals (can be positive or negative or 0), that define a quatratic-bezier patch, that plane-normal-projects exactly on the base-triangle the intersection of a line with such a simple bezier patch is actually analytic and fast and precise (enough within how i use it with ieee 32 bit float .).

https://www.shadertoy.com/view/fsBcWc demonstrates this, but it has a culling-bug, IF the camera is above or below the base-triangle. Bonus: the surface-normal is rather cheap to calculate.

===

general case = likely non analytic, long slow loops, but relatively low O(n) complexity (more than 4 roots, usually around 16 roots):

GENERALLY, solving for intersection of line-or-ray with almost-all-bezier-patches is non-analytic (more than 4 roots);

and usually a slow iteration process (with all the primitive-recursions of a bezier, that can be turned into loops-with-accumulator) to a point where other approaches are done:

• Rasterization into triangle meshes (relatively fast for bezier patches)
• substituting the surface with a height-map or displacement map, with various 3d-estimating rendering approaches (super fast, but textures need more memory)
• substituting the surface with a signed-distance-field and perform sphere-tracking on it (arbitrarily fast/precise, usually less energy inefficient, costs too much energy for mobile-hardware)
• possibly complex numbers, double-integrals and 3rd derivatives, who knows, but this may not be worth the compute cost for insignificant gains in precision.

(shadertoy.com has MANY approaches to "bezier patches that span a square/rectangle") some ray-trace (usually iterating a surface) some sphere-track (ray--march), not many rasterize, because shadertoy is 99% fragment shading and next to null vertex-shading.)

most bezier-tracing is non analytic, 18 roots for a rectangular bezier patch.

while quaretic-root-solving (4 roots (or less)) is analytic, it involves so long fractions for each of the 4 roots (independently, and your first root must be solved this way), that the precision suffers significantly (for finding the 1st root) (compared to calculating 3 cubic roots of a cubic-solver), that 32-bit-floats are likely not precise enough for (naive) quartic-root-solving. naive quartic root solving solves 4 roots independently. non naive solves 1of4 roots this way, and then inserts that as low-precision-constant into a remaining problem of cubic-root-solving, getting you slightly better precision,. but only IF the first root is not-too-far-away from the other 3 roots. if it is far away, the precision is likely even worse, than the naive solver (but the non-naive process is likely still much energy-efficient)

And you ideally know in advance, which of the 4 roots is closest tot he other 3 roots, and solve that root first, and then solve the remaining 3 roots as cubic.

This is why shadertoy.com mostly stops at cubic-root-solving and while it /shows/demoes/knows and has quartc-root solving (opngl opensource), it also knows that it is too slow and too low precision with 32 bit floats.