Why 3D Fractals Need Different Rendering

Traditional 3D rendering uses triangle meshes - surfaces broken into thousands of small triangles that a GPU can rasterize quickly. This works for characters, buildings, and terrain because those surfaces can be described by a finite number of polygons.

Fractal surfaces are infinitely detailed. No matter how many triangles you use, you lose the fractal's self-similar detail at smaller scales. More fundamentally, fractal surfaces have no closed-form equation for where a ray intersects them - you cannot solve "where does this ray hit the Mandelbulb?" algebraically.

Instead, 3D fractals are rendered using ray marching - a technique where rays are cast from a virtual camera through each pixel, then advanced step by step through space until they find the fractal surface. The key insight that makes this practical is distance estimation: at each step, a function tells you how far away the nearest surface is, so you can take the largest safe step without passing through anything.

Sphere Tracing

Sphere tracing is the specific ray marching algorithm used for rendering distance-field geometry. It was formalized by John C. Hart in his 1996 paper "Sphere tracing: a geometric method for the antialiased ray tracing of implicit surfaces", published in The Visual Computer.

The algorithm is elegant: at each point along a ray, you evaluate a distance estimator (DE) - a function that returns a conservative lower bound on the distance to the nearest surface. You then advance the ray by exactly that amount. The name "sphere tracing" comes from the geometric interpretation: the DE defines a sphere around the current point that is guaranteed to be empty, and the ray advances to the sphere's boundary.

Ray marching loop: p = ray_origin + t × ray_direction d = DE(p) // distance to nearest surface t = t + d // advance by that distance // repeat until d < ε (hit) or t > max_distance (miss)

This approach is fundamentally more efficient than fixed-step ray marching (which uses uniform step sizes). In open space far from the fractal, steps are large and traversal is fast. Near the surface, steps shrink automatically for precision. And critically, the algorithm never overshoots the surface - if the distance estimate is truly a lower bound, you are guaranteed not to miss features.

Hart's foundational work actually began earlier: his 1989 SIGGRAPH paper "Ray Tracing Deterministic 3-D Fractals" (with Sandin and Kauffman) was the first to demonstrate distance estimation for rendering quaternion Julia sets - the first practical method for visualizing 3D fractals.

Sources:

Distance Estimation

The distance estimator is the heart of fractal ray marching. For simple geometry - a sphere, a box, a torus - signed distance functions (SDFs) have exact closed-form expressions. For fractals, the SDF must be approximated from the escape-time iteration itself.

The key insight comes from the Hubbard-Douady potential. For escape-time fractals (Mandelbrot, Mandelbulb, quaternion Julia sets), you iterate a function z → f(z) and track two values: the orbit zn and its running derivative z'n. When the orbit escapes (|zn| exceeds a bailout radius), the distance to the fractal surface can be estimated as:

DE(z) = 0.5 × |z_n| × log(|z_n|) / |z'_n|

This formula derives from the Green's function of the complement of the fractal set. The factor of 0.5 ensures the estimate is a lower bound - safe for sphere tracing, because it guarantees you never step through the surface.

The Running Derivative

The running derivative z'n is computed alongside the main fractal iteration. For a Mandelbrot-type power-n iteration z → zn + c (differentiating with respect to c):

z'_0 = 1 z'_{k+1} = n × z_k^(n-1) × z'_k + 1

For Julia sets, where c is fixed and you differentiate with respect to the initial point z0, the constant term is absent: z'k+1 = n × zkn-1 × z'k.

For 3D fractals, this generalizes to tracking the Jacobian matrix of the transformation. In practice, a scalar running derivative (tracking only the magnitude) is sufficient and much cheaper to compute. This scalar approach was first proposed by community member Buddhi on the Fractal Forums for the Mandelbox.

Fractal-Specific Distance Functions

Each type of 3D fractal uses a different mathematical operation for its iteration, which means each requires its own distance estimation approach. Here are the techniques used by the fractals you can explore on Synaptic Spiral.

Mandelbulb - Spherical Power Mapping

The Mandelbulb extends the Mandelbrot set into 3D using spherical coordinates. Instead of complex multiplication (which doubles angles and squares radii in 2D), the Mandelbulb raises a 3D point to a power n using a "triplex algebra":

1. Convert (x, y, z) → (r, θ, φ) [Cartesian to spherical] 2. Compute v^n = r^n × ( sin(nθ)cos(nφ), sin(nθ)sin(nφ), cos(nθ) ) 3. Iterate: z_{k+1} = z_k^n + c

The power-8 formula produces the iconic Mandelbulb shape discovered by Daniel White and Paul Nylander in 2009. The distance estimator uses a scalar running derivative:

dr_{k+1} = n × r^(n-1) × dr_k + 1 DE = 0.5 × log(r) × r / dr

Explore the Mandelbulb live

Sources:

Mandelbox - Box Fold + Sphere Fold

The Mandelbox, discovered by Tom Lowe (Tglad) in 2010, takes a completely different approach. Instead of raising points to a power, it applies two geometric folding operations:

Box Fold: if component > 1, reflect to (2 - component) if component < -1, reflect to (-2 - component) Sphere Fold: if |z| < minRadius², scale by (fixedRadius/minRadius)² if |z| < fixedRadius², scale by (fixedRadius/|z|)² Iteration: z_{k+1} = scale × sphereFold(boxFold(z_k)) + c

The box fold is an isometry (it preserves distances), while the sphere fold scales space non-uniformly. The distance estimator tracks how these operations scale the running derivative at each step. Values of the scale parameter near −2.5 to −1.5 typically produce the most complex structures.

Explore the Mandelbox live

Sources:

Quaternion Julia Sets - 4D Algebra

Quaternion Julia sets extend the classic 2D Julia set into 4D using quaternion algebra (a + bi + cj + dk). A 3D cross-section is taken for visualization. The iteration is qk+1 = qk² + c using quaternion multiplication, and the distance estimator follows the same Hubbard-Douady formula:

q'_0 = 1 q'_{k+1} = 2 × q_k × q'_k [quaternion product] DE = 0.5 × |q_n| × log(|q_n|) / |q'_n|

Keenan Crane's 2005 paper "Ray Tracing Quaternion Julia Sets on the GPU" demonstrated that this could be done at interactive frame rates using GPU fragment shaders - an early milestone in real-time fractal rendering.

Explore Quaternion Julia Sets live

Benesi Pine Tree - Trigonometric Variants

The Benesi formulas are alternative 3D fractal constructions that use different trigonometric mappings than the Mandelbulb. When no analytical distance estimator is available, a numerical approach called the Makin/Buddhi 4-point Delta-DE is used - it probes the escape function at four nearby points to numerically approximate the gradient, similar to computing normals by central differences.

Explore the Benesi Pine Tree live

Lighting & Shading in Fractal Space

Once a ray hits the fractal surface, the scene needs lighting to look three-dimensional. Since fractal surfaces are defined implicitly by the distance estimator rather than by explicit geometry, all lighting calculations must be derived from the DE function itself.

Normal Estimation via Central Differences

The surface normal (which determines how light reflects) is the gradient of the distance field. Since we have no analytical expression for this gradient, it is computed numerically by sampling the DE at six points surrounding the hit point:

normal(p) = normalize( DE(p + (ε,0,0)) - DE(p - (ε,0,0)), DE(p + (0,ε,0)) - DE(p - (0,ε,0)), DE(p + (0,0,ε)) - DE(p - (0,0,ε)) )

This requires six additional DE evaluations per surface hit - expensive for complex fractals. Inigo Quilez documented an optimized tetrahedron technique that achieves the same quality with only four evaluations using cleverly chosen sample offsets.

Ambient Occlusion from the Distance Field

Ambient occlusion (AO) - the darkening of crevices and enclosed areas - can be approximated cheaply by sampling the DE at increasing distances along the surface normal. If the DE returns a value smaller than the sample distance, nearby geometry is blocking light from that direction:

ao = 0 for i = 1 to 5: expected = step_size × i actual = DE(p + normal × expected) ao += (expected - actual) / 2^i

An even simpler approach (and one commonly seen in real-time fractal renderers) colors surfaces darker based on how many ray march steps were needed to reach them. More steps means the ray was navigating narrow crevices - a natural proxy for occlusion.

Soft Shadows via Ray Marching

Inigo Quilez developed a widely-used soft shadow technique that works naturally with distance fields. For each surface point, a shadow ray is marched toward the light. At each step, the ratio of the DE value to the distance traveled gives a penumbra estimate:

shadow = 1.0 for each step along shadow ray: h = DE(current_position) shadow = min(shadow, k × h / distance_traveled) // k controls shadow softness (higher = harder shadows)

This produces realistic soft shadows with no additional geometry - something that is particularly effective for revealing the intricate three-dimensional structure of fractals.

Color Mapping Techniques

Color is what transforms fractal geometry from mathematical curiosity into visual art. Several techniques exist, often combined for richer results.

Iteration Count Coloring

The simplest method: color each point by the number of iterations before its orbit escapes. Raw integer counts produce visible banding, so a smooth iteration count is used instead:

smooth_iter = n + 1 - log(log(|z_n|)) / log(power)

This fractional value is mapped to a continuous color gradient, eliminating the staircase artifacts of integer-based coloring.

Orbit Traps

Orbit traps color each point based on how close its iterative orbit passes to a predefined geometric shape (the "trap"). During iteration, you record the minimum distance from the orbit to the trap and map it to color. Common trap shapes include points, lines, circles, and crosses. Different trap geometries produce dramatically different visual styles from the same underlying fractal.

Distance-Based Coloring (3D)

For 3D fractals rendered via ray marching, several distance-derived quantities provide color information:

  • Final distance estimate at the surface - reveals fine detail at the boundary
  • Total ray distance traveled - creates depth fog effects
  • Number of ray march steps - reveals surface complexity (more steps = more intricate geometry)
  • Minimum distance the ray passed to any surface during marching - creates glow effects around the fractal
Sources:

GPU Acceleration with WebGL

Ray marching is an embarrassingly parallel algorithm - each pixel's ray is completely independent of every other pixel. This maps perfectly to the GPU execution model, where thousands of shader cores run simultaneously.

In practice, a full-screen quad (two triangles) is drawn, and the fragment shader does all the work: constructing a ray for each pixel, running the entire ray marching loop, evaluating the distance estimator, computing normals, and applying lighting - all per-pixel, all in parallel. A 1920×1080 render dispatches roughly two million parallel ray marching computations simultaneously.

WebGL/GLSL Considerations

  • Loop limits - GLSL requires compile-time-known loop bounds, so maximum iteration counts and ray march steps must be defined as constants
  • Floating-point precision - highp float is essential for fractal detail; medium precision causes visible artifacts, especially at zoom
  • Branching efficiency - GPUs execute fragments in warps/wavefronts; divergent branching (some pixels hitting the surface early, others marching far) reduces throughput
  • Uniforms for interactivity - Camera position, fractal constants, and iteration counts are passed as uniforms, allowing real-time interactive control without recompiling the shader

All five 3D visualizations on Synaptic Spiral - Mandelbulb, Mandelbox, Benesi Pine Tree, Hopf Fibration, and Quaternion Julia - are rendered using these techniques in your browser via WebGL fragment shaders.

Sources:

Key Figures & History

3D fractal visualization is a relatively young field. A handful of individuals - spanning academia, the demoscene, and online fractal communities - built the foundations.

John C. Hart

Foundational figure in fractal rendering. His 1989 SIGGRAPH paper (with Sandin and Kauffman) first demonstrated distance estimation for rendering quaternion Julia sets. His 1996 sphere tracing paper formalized the algorithm now used universally for rendering distance-field geometry.

Daniel White

Amateur fractal enthusiast who began exploring 3D extensions of the Mandelbrot set on the Fractal Forums in 2007. His "triplex" algebra using spherical coordinates, combined with Paul Nylander's suggestion to use power-8, produced the Mandelbulb in 2009. Read his original account at Skytopia.

Tom Lowe (Tglad)

Discovered the Mandelbox in 2010, introducing folding-based fractal iteration as an alternative to the power-based approach of the Mandelbulb. Posted to the Fractal Forums community.

Inigo Quilez

Co-created Shadertoy (with Pol Jeremias) in 2013. His website iquilezles.org is the single most comprehensive resource on SDF techniques, ray marching, soft shadows, ambient occlusion, and fractal rendering. Professional work includes Pixar and Oculus/Meta.

Mikael Hvidtfeldt Christensen

Authored the definitive 8-part blog series "Distance Estimated 3D Fractals" (2011 - 2012) on his Syntopia blog - the most thorough tutorial on the subject. Also created Fragmentarium, a GPU-based IDE for exploring fractals with GLSL.

Keenan Crane

Demonstrated GPU-accelerated ray tracing of quaternion Julia sets in 2005, an early milestone in real-time fractal rendering. Now a professor of computer science at Carnegie Mellon University.

Further Reading