Spotlight #62 – WebGL 3D — GLSL Raymarching, Signed Distance Fields & GPGPU Simulation

I. Raymarching and Signed Distance Fields — The Core Technique

Before looking at individual simulations, it is worth understanding the shared foundation they all build on. Traditional 3D rendering works by transforming triangles through a rasterisation pipeline: every surface must be described as a mesh of polygons. This works beautifully for engineering models but is ill-suited to objects defined implicitly — fractals, fluid surfaces, algebraic shapes — where no triangle mesh is practical.

A signed distance field (SDF) is a function d(p) that returns the signed distance from any point p in 3D space to the nearest surface: positive outside, negative inside, zero on the surface. For a sphere centred at the origin with radius r:

d(p)  =  length(p)  -  r

SDFs can be combined with simple arithmetic: the union of two shapes is min(d1, d2); the intersection is max(d1, d2); subtraction (hollowing one shape from another) is max(d1, -d2). This algebra of shapes, called constructive solid geometry (CSG), lets you build arbitrarily complex objects from primitives without ever constructing a mesh.

The Raymarching Loop

Given an SDF, raymarching renders it by shooting a ray from the camera through each screen pixel and marching along the ray in steps sized by the SDF value:

float t = 0.0;
for (int i = 0; i < MAX_STEPS; i++) {
    vec3 p = ray_origin + t * ray_dir;
    float d = sdf(p);
    if (d < EPSILON) { /* hit surface */ break; }
    if (t > MAX_DIST) { /* missed, sky colour */ break; }
    t += d;  // safe step: won't overshoot surface
}

The key insight is that the SDF value is a safe step size: the ray cannot overshoot the surface because d(p) is exactly the minimum distance to any surface from point p. This sphere-tracing algorithm (Hart, 1996) is both safe and efficient: in open space, large steps are taken; near a surface, steps shrink to sub-pixel precision. A typical complex scene converges in 50–100 iterations per pixel, each iteration performing one SDF evaluation — a workload that maps perfectly to massively parallel GPU execution, with one thread per pixel.

Surface Normals and Lighting

Once a hit is found, the surface normal is computed by finite differences — sampling the SDF at six nearby points and taking the gradient:

vec3 normal(vec3 p) {
    float e = 0.001;
    return normalize(vec3(
        sdf(p + vec3(e,0,0)) - sdf(p - vec3(e,0,0)),
        sdf(p + vec3(0,e,0)) - sdf(p - vec3(0,e,0)),
        sdf(p + vec3(0,0,e)) - sdf(p - vec3(0,0,e))
    ));
}

This requires only six additional SDF calls per pixel and works for any implicit surface without explicitly computing a mesh normal. Standard Blinn-Phong or physically-based lighting then uses the normal, viewing direction, and light direction in the usual way. Soft shadows and ambient occlusion are computed by additional short ray marches toward the light source: if the shadow ray's minimum SDF value is very small, the point is partially occluded, producing smooth penumbras impossible with hard shadow maps.

II. Mandelbulb — Distance-Estimated 3D Fractal

The Mandelbrot set lives in the complex plane: for each point c, iterate z → z² + c and test whether the orbit escapes. Extending this to three dimensions requires a notion of "multiplication" for 3D vectors. The Mandelbulb (White and Nylander, 2009) uses a triplex power in spherical coordinates (r, θ, φ):

z^n:  r_new = r^n
      θ_new = n · θ
      φ_new = n · φ

This is not an algebra in the strict sense (it lacks associativity for n ≠ 2), but it produces visually extraordinary results for n = 8: a bulbous structure covered in intricate spherical bumps, each of which, on closer inspection, is covered in smaller copies of the same bumps — a self-similar fractal at all scales.

The Distance Estimator

The Mandelbulb has no closed-form SDF, but it has a distance estimator (DE) derived from the derivative of the iteration map. Track both the current iterate z and its derivative dz/dc during the iteration loop:

// Each iteration:
dz = n * pow(r, n-1) * dz + 1.0;
z  = triplex_pow(z, n) + c;
if (length(z) > BAILOUT) break;

After the iteration escapes, the distance estimate is:

de  =  0.5 · log(|z|) · |z| / |dz/dc|

This estimate is accurate to within a small constant factor near the surface, which is sufficient for raymarching: the sphere-tracing loop can safely use it as a step size. The cost is that each raymarching step requires running the full fractal iteration to determine the DE, making the Mandelbulb one of the most GPU-intensive simulations on the site — yet still interactive at 60 fps on modern hardware.

III. Menger Sponge — IFS via SDF Folding

The Menger sponge is constructed by repeatedly removing the centre cube and the centres of each face from a cube — an iterated function system (IFS) that produces a fractal with Hausdorff dimension log(20)/log(3) ≈ 2.727. Constructing a mesh for even 4 iterations requires millions of polygons. The SDF approach, using domain folding, evaluates any iteration depth with a constant-cost function.

The key trick is modular folding: instead of recursively subdividing the domain, fold space so that a single SDF evaluation of the base cube represents the entire fractal. At each iteration, the coordinates are folded into the fundamental domain of the symmetry group, scaled up by a factor of 3, and the cross-shaped cavity is subtracted:

float menger_sdf(vec3 p, int iterations) {
    float d = box_sdf(p, vec3(1.0));
    float scale = 1.0;
    for (int i = 0; i < iterations; i++) {
        // Fold into unit cell, reflect across axes
        vec3 q = mod(p * scale, 2.0) - 1.0;
        q = abs(q);
        // Cross cavity: remove centre tube in each axis pair
        float cross = min(max(q.x, q.y),
                     min(max(q.y, q.z),
                         max(q.z, q.x)));
        float cavity = cross - 1.0/3.0;
        d = max(d, -cavity / scale);
        scale *= 3.0;
    }
    return d;
}

Each loop iteration corresponds to one level of the fractal hierarchy. The division by scale corrects the SDF to account for the spatial compression introduced by the fold. The beauty of this approach is that iterations 1 through 6+ are available at the cost of a handful of additional loop iterations, with no geometry change — just a slider controlling the iteration count.

IV. Wormhole 3D — Morris-Thorne Geometry in GLSL

A Morris-Thorne wormhole is a solution to Einstein's field equations describing a traversable tunnel between two regions of spacetime. The metric in spherical coordinates is:

ds²  =  -e^(2Φ) dt²  +  dl²  +  (b₀² + l²)(dθ² + sin²θ dφ²)

where l is the proper distance along the tunnel axis, b₀ is the throat radius (minimum circumference radius / 2π), and Φ(l) is the redshift function. The embedding surface — a 2D slice that captures the spatial geometry — is the surface of revolution of z(r) = ±b₀ · ln(r/b₀ + √((r/b₀)²-1)), which flares outward from the throat like two connected trumpets.

The GLSL simulation does not solve the full general-relativistic ray equations (which would require integrating coupled ODEs per pixel per frame). Instead, it uses a geometric approximation: rays are bent as they approach the throat by a space-warp shader that smoothly maps incoming directions to outgoing directions according to the embedding's intrinsic curvature. Stars and nebulae appear to wrap around the throat due to gravitational lensing. This is the same technique used in the film Interstellar (2014), for which Kip Thorne and the Double Negative visual effects team developed a physically accurate GLSL wormhole renderer.

V. Supershape — Gielis Superformula on the GPU

Johan Gielis introduced the superformula in 2003 as a generalisation of the circle that can describe an extraordinary range of natural and mathematical shapes with just six parameters (m, n1, n2, n3, a, b):

r(θ)  =  [ |cos(mθ/4)/a|^n2  +  |sin(mθ/4)/b|^n3 ]^(-1/n1)

In polar coordinates, r(θ) gives the radius at angle θ. For a 3D supershape, two independent superformula evaluations are applied as latitude and longitude radii, producing a sphere-like surface deformed by the product of the two radial functions. The GPU maps each pixel directly to a (θ, φ) coordinate, evaluates the two formulas, computes the 3D position on the surface, and renders via a standard lighting model — no mesh, no geometry shader, just a pure parametric surface evaluated per pixel.

VI. GPGPU Particles — 250,000 Particles via Render-to-Texture

The central challenge of GPU particle simulation is a WebGL constraint: a shader cannot read from and write to the same texture simultaneously. The ping-pong technique solves this by maintaining two pairs of render targets (A and B). On odd frames, the compute pass reads from A and writes to B; on even frames it reads from B and writes to A. The CPU never touches the particle data after initialisation.

// Pseudo-code for one frame:
gl.bindFramebuffer(writeFBO);     // bind target B
gl.useProgram(computeShader);
gl.bindTexture(readTexture);      // bind target A
gl.drawArrays(TRIANGLES, 0, 6);  // full-screen quad: runs physics on GPU
swap(A, B);                       // next frame reads from B

// Render pass:
gl.bindFramebuffer(null);         // back buffer
gl.useProgram(renderShader);
gl.bindTexture(positionTexture);  // current positions (just written)
gl.drawArrays(POINTS, 0, N);     // each point vertex samples its own texel

Each particle's position is stored as an RGBA32F texel (four 32-bit floats: xyz position, w for age or mass). A 512×512 texture stores 262,144 particles. The compute shader is a full-screen quad whose fragment shader looks up the particle position at gl_FragCoord, applies a physics update (gravity, curl noise, or attractor), and writes the new position as the output colour. The render shader is a points draw call: each point vertex index maps to a texel coordinate, reads position via a texture fetch in the vertex shader, and outputs a point sprite.

Curl Noise and Divergence-Free Flow

The most visually compelling particle behaviour uses curl noise for the velocity field. Given any 3D vector potential A(p), the curl F = ∇ × A is mathematically guaranteed to be divergence-free:

∇ · (∇ × A)  =  0   (identity in vector calculus)

A divergence-free field has no sources or sinks — particles neither clump into dense blobs nor fly apart. Instead they form incompressible swirling streams that look strikingly like smoke, clouds, and ocean currents. The vector potential A is taken as a 3D simplex noise field, and the curl is computed analytically (the noise gradient is differentiable). Advecting 250,000 particles through curl noise at 60 fps requires no CPU involvement after the shader is compiled — the GPU handles the entire simulation and rendering pipeline.