3D Fractals: Mandelbulb, Menger Sponge & Distance-Field Rendering

1. Why 3D Fractals Are Hard

The Mandelbrot set lives in the complex plane and is defined by the iteration z → z² + c. In 2D the squaring operation is simply multiplication of complex numbers: (a + bi)² = a² − b² + 2abi. The key property that makes this work fractally is analyticity — the derivative exists and obeys the Cauchy-Riemann equations, which means the map is conformal (angle-preserving) and produces infinitely intricate self-similar boundaries.

The natural extension to three dimensions would require a normed division algebra over ℝ³. By Hurwitz's theorem, such an algebra does not exist — the only valid normed division algebras have dimension 1, 2, 4 (quaternions), or 8 (octonions). Quaternions give a 4D fractal (the quaternion Julia / Mandelbrot set), but projecting it back to 3D produces relatively smooth shapes without the intricate bulb structure we expect.

Two successful approaches exist:

• Algebraic hacks — define a "power" operation on 3D vectors via spherical coordinates (the Mandelbulb approach), sacrificing true algebraicity but achieving visually spectacular results.
• Geometric iteration — apply affine or folding transformations repeatedly in 3D (Menger sponge, Mandelbox), which requires no algebraic structure at all.

Real-time rendering of these implicit 3D objects is made practical by signed distance field ray-marching, which uses distance estimator functions derived from the fractal's escape-time iteration.

2. The Mandelbulb: Spherical Power Formula

The Mandelbulb was introduced by Daniel White and Paul Nylander in 2009. The idea is to represent a 3D point v = (x, y, z) in spherical coordinates (r, θ, φ), raise it to the n-th power by multiplying the radius n times and multiplying the angles, then convert back:

The Mandelbulb iteration is then z → z^n + c (in the above 3D sense) starting from z₀ = 0. A point c is in the set if |z| remains bounded. The canonical choice is n = 8, which yields intricate bulbous structures with elephant-valley and seahorse-valley analogues.

Why n = 8?

For n = 2 the result looks like a smooth ball — not interesting. As n increases the surface develops more self-similar detail. At n = 8, the fractal displays a convincing 3D analogue of the Mandelbrot set with deeply nested spiralling features and a rich boundary Hausdorff dimension estimated at around 2.97 — nearly space-filling.

Limitations

Because the spherical power does not obey the product rule for derivatives in the true algebraic sense (it's not a holomorphic map), the Mandelbulb is not a true 3D Mandelbrot set. Interior regions and the boundary do not have the same mathematical guarantees as the 2D case — but visually and computationally it is one of the most rewarding 3D fractals known.

Distance Estimator

To ray-march the Mandelbulb we need a distance estimator (DE) — a lower bound on the distance from a point to the fractal surface. Derived from the derivative magnitue along the iteration:

This is Íñigo Quílez's formula. It gives a conservative (undershooting) distance bound — always safe to step by a fraction (typically 0.8–0.99) of the DE to prevent tunnelling through thin features.

3. Menger Sponge: Iterated Function System in 3D

The Menger sponge is constructed by starting with a unit cube and recursively dividing it into a 3×3×3 grid of 27 subcubes, then removing the central cube and the 6 face-centre cubes (leaving 20 subcubes), and repeating the process on each remaining subcube.

After k iterations the sponge has 20^k cubes each of side length (1/3)^k. Its Hausdorff dimension is:

The sponge has zero volume (in the limit) but infinite surface area — a characteristic property of fractals with dimension between 2 and 3.

SDF for the Menger Sponge (Folding Approach)

The Menger sponge can be computed with a folding SDF that mirrors the recursive construction without branching:

The sortings (xy, xz, yz swaps) fold the point into a canonical octant at each level, and the scale factor accumulates through the chain rule just as with the Mandelbulb DE.

Comparison of Approaches

• Escape-time (Mandelbulb): suitable for fractals defined by complex-like iteration; natural smooth colouring via orbit traps.
• Iterated folding (Menger, Mandelbox): suitable for fractals defined by affine maps; typically cheaper to evaluate per step.
• Both yield a valid DE for sphere-tracing; the choice is purely definitional.

4. Signed Distance Fields & Ray Marching

A signed distance field (SDF) is a scalar function f(p) that returns the signed distance from point p to the nearest surface: negative inside, positive outside, zero on the surface. Ray marching (sphere tracing) exploits the SDF to advance a ray safely:

1. Cast a ray from camera origin o in direction d (normalised).
2. At the current position p, evaluate the SDF: t = f(p).
3. Jump forward by t along the ray: p += t · d.
4. Repeat until |t| < ε (hit) or accumulated distance > max_dist (miss).

Normal Estimation

Surface normals are computed via the gradient of the SDF using finite differences:

This requires 6 SDF evaluations per hit. Tetrahedron-based finite differences (Quílez trick) reduce this to 4 evaluations while preserving accuracy.

Performance Considerations

• MAX_STEPS is the dominant cost. Thin fractal features require many steps; typical range 64–256.
• Step size can be multiplied by a relaxation factor ω > 1.0 (overrelaxation / "enhanced sphere tracing") to reduce step count in open regions, checking for backtrack when the ray overshoots.
• Bounding volumes (e.g., sphere of radius 1.2 around the fractal) eliminate background pixels cheaply.

5. SDF Distance Estimators for Fractals

For escape-time fractals the DE is derived from the modulus of the derivative along the orbit. The general formula for a map z → f(z) starting at c is:

A key subtlety: the DE must be divided by the accumulated scale factor for folding-based fractals. For the Menger sponge, each iteration multiplies coordinates by 3, so the DE is divided by 3^ITERATIONS. Forgetting this produces a shape that appears geometrically correct but with normals that point in random directions — a common debugging headache.

Ambient Occlusion from Step Count

A cheap but effective ambient occlusion estimate is the number of ray-marching steps taken divided by MAX_STEPS. Regions requiring many small steps (near deep crevices) are darker; open flat regions (few steps) are brighter. This correctly darkens fractal fjords and deep creases without extra rays.

6. Orbit-Trap Colouring & Smooth Iteration Count

Smooth Iteration Count

The raw escape-time iteration count is an integer, producing visible bands on the rendered image. Smooth colouring replaces the band index n with a real value:

This formula comes from the potential theory of the Mandelbrot set: the Green's function G(c) = lim log|z_n| / 2^n is continuous across band boundaries. The correction logâ‚‚(logâ‚‚(r)) removes the integer staircasing.

Orbit Traps

Orbit traps colour a point based on how close the orbit passes to a geometric shape:

• Point trap: record min(|z − trap_point|) over the orbit. Map this minimum to a colour.
• Plane trap: min absolute value of one coordinate over all iterations.
• XY-cross trap: min(|x|, |y|) — produces vibrant cross-shaped colour patterns.
• Circle trap: min | |z| − r_trap | — produces ring patterns around interior features.

Multiple traps can be blended by weighted minimum or by tracking which trap was triggered at the closest approach. The result is far richer than simple gradient colouring.

Normal-Based Lighting

Once a hit point and its normal are known, standard Phong or Blinn-Phong shading applies directly. Adding a specular highlight with exponent 64–128 dramatically enhances the perception of the complex surface. Soft shadows can be approximated by marching a secondary ray toward the light and tracking how close it comes to the surface (the "shadow softness" penumbra technique by Quílez).

7. GLSL Ray-Marcher: Complete Shader

Below is a self-contained GLSL fragment shader for the power-8 Mandelbulb. It runs in WebGL 1.0 with appropriate precision hints.

// ── Mandelbulb Ray-Marcher (GLSL ES 3.0) ──────────────────────────
precision highp float;
uniform vec2 u_resolution;
uniform float u_time;

const int ITER = 8;
const int STEPS = 120;
const float BAILOUT = 4.0;
const int N = 8;

// Distance estimator for Mandelbulb power-8
float mandelbulbDE(vec3 pos) {
  vec3 z = pos;
  float dr = 1.0, r = 0.0;
  for (int i = 0; i < ITER; i++) {
    r = length(z);
    if (r > BAILOUT) break;
    float theta = acos(clamp(z.z / r, -1.0, 1.0));
    float phi   = atan(z.y, z.x);
    dr = pow(r, float(N - 1)) * float(N) * dr + 1.0;
    float zr = pow(r, float(N));
    z = zr * vec3(
      sin(float(N)*theta) * cos(float(N)*phi),
      sin(float(N)*theta) * sin(float(N)*phi),
      cos(float(N)*theta)
    ) + pos;
  }
  return 0.5 * log(r) * r / dr;
}

// Finite-difference normal (tetrahedron method)
vec3 calcNormal(vec3 p) {
  const vec2 k = vec2(1.0, -1.0) * 5e-4;
  return normalize(
    k.xyy * mandelbulbDE(p + k.xyy) +
    k.yyx * mandelbulbDE(p + k.yyx) +
    k.yxy * mandelbulbDE(p + k.yxy) +
    k.xxx * mandelbulbDE(p + k.xxx)
  );
}

// Sphere-trace
float trace(vec3 ro, vec3 rd) {
  float t = 0.0;
  for (int i = 0; i < STEPS; i++) {
    float d = mandelbulbDE(ro + t * rd);
    if (d < 1e-4) return t;
    t += d * 0.8;
    if (t > 8.0) break;
  }
  return -1.0;
}

void main() {
  vec2 uv = (gl_FragCoord.xy - 0.5 * u_resolution) / u_resolution.y;

  // Orbiting camera
  float angle = u_time * 0.3;
  vec3 ro = vec3(2.5 * sin(angle), 0.5, 2.5 * cos(angle));
  vec3 ta = vec3(0.0);
  vec3 ww = normalize(ta - ro);
  vec3 uu = normalize(cross(ww, vec3(0, 1, 0)));
  vec3 vv = cross(uu, ww);
  vec3 rd = normalize(uv.x * uu + uv.y * vv + 1.5 * ww);

  float t = trace(ro, rd);
  vec3 col = vec3(0.05, 0.05, 0.08);  // background

  if (t > 0.0) {
    vec3 p    = ro + t * rd;
    vec3 n    = calcNormal(p);
    vec3 lig  = normalize(vec3(0.8, 0.6, 0.4));
    float diff = clamp(dot(n, lig), 0.0, 1.0);
    float spec = pow(clamp(dot(reflect(-lig, n), -rd), 0.0, 1.0), 64.0);
    float trap = clamp(0.1 / abs(p.x * p.y), 0.0, 1.0);
    vec3 baseCol = mix(vec3(0.1, 0.2, 0.8), vec3(0.9, 0.4, 0.1), trap);
    col = baseCol * diff + vec3(0.5) * spec;
    col *= 0.8 + 0.2 * (1.0 - t / 8.0);
  }

  gl_FragColor = vec4(pow(col, vec3(0.4545)), 1.0);
}

8. Extensions: Mandelbox, Juliabulb & Folding Space

Mandelbox

The Mandelbox (Tom Lowe, 2010) replaces the spherical power with a box fold followed by a sphere fold:

The Mandelbox produces infinitely nested box-like structures — recursive hallways, city-like skylines, fractal coral. The scale parameter (typically −2.0 to −1.5) governs overall character; values near −1.5 produce wispy fern-like structures while −2.0 yields massive angular cliffs.

Juliabulb

Just as the 2D Julia set freezes c and varies zâ‚€, the Juliabulb fixes a 3D parameter c and uses varying start points. This slices a different region of parameter space, producing distinctly symmetric shapes. Animating c through a path creates morphing fractal animations.

Kleinian Groups & Sphere-Inversion Fractals

More general folding fractals apply Möbius transformations — inversions through spheres — recursively. Related to 3D Kleinian limit sets and hyperbolic geometry, these produce translucent bubble-cluster fractals with an accurate DE derivable from the chain rule.

Hybrid Fractals

Alternating iteration steps from different formulas — one Mandelbulb step then one box-fold step — creates "hybrid" fractals with entirely new morphologies. Most celebrated fractal renders online are precisely such hybrids, exploiting the limitless composability of distance estimators.