How to raytrace Bezier surfaces?

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.

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.

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.

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.)