GLSL 3D Fractals - Distance Estimation Methods

Distance estimation is the cornerstone of efficient 3D fractal rendering. Unlike traditional ray tracing, which requires expensive intersection tests, distance estimation allows us to step along rays efficiently using guaranteed lower bounds on distances. This article explores the mathematical foundations and practical implementations of distance estimation for fractals.

Why Distance Estimation?

Traditional ray-fractal intersection would require:

Testing every iteration of the fractal formula
Complex geometric intersection calculations
Potentially infinite computation for fractals

Distance estimation provides:

Guaranteed lower bounds: We know we're at least distance $d$ away
Efficient stepping: Step by $d$ without missing the surface
Early termination: Stop when close enough or too far
Real-time performance: Enables 60 FPS rendering

Mathematical Foundation

Distance Estimation Formula

For fractals defined by iteration $z_{n+1} = f(z_n)$ , the distance estimate is:

$$d \approx \frac{|z_n| \cdot \log(|z_n|)}{|z'_n|}$$

Where:

$z_n$ is the iterated value after $n$ iterations
$z'_n$ is the derivative of the iteration
The formula assumes $|z_n| > 1$ (outside the fractal)

Derivative Accumulation

The derivative $z'_n$ is computed alongside the iteration:

$$z'_{n+1} = f'(z_n) \cdot z'_n$$

Starting with $z'_0 = 1$ .

Proof Sketch

The distance estimate comes from the relationship between the iteration and the distance to the fractal set. For a point $\mathbf{p}$ outside the fractal:

$$d(\mathbf{p}) \propto \frac{|\mathbf{p} - \mathbf{p}_0|}{|\nabla f(\mathbf{p}_0)|}$$

Where $\mathbf{p}_0$ is the closest point on the fractal. The iteration $z_n$ approximates this relationship.

Implementation Patterns

Basic Distance Estimation

float fractalDE(vec3 pos) {
    vec3 z = pos;
    float dr = 1.0;  // Derivative
    float r = 0.0;
    
    for (int i = 0; i < maxIterations; i++) {
        r = length(z);
        if (r > bailout) break;
        
        // Update derivative
        dr = fDerivative(z, dr);
        
        // Iterate
        z = f(z);
    }
    
    // Distance estimate
    return 0.5 * log(r) * r / dr;
}

Mandelbulb Distance Estimation

For the Mandelbulb with power $n$ :

float mandelbulbDE(vec3 pos, float power) {
    vec3 z = pos;
    float dr = 1.0;
    float r = 0.0;
    
    for (int i = 0; i < iterations; i++) {
        r = length(z);
        if (r > 4.0) break;
        
        // Derivative: dr = r^(n-1) * n * dr + 1
        dr = pow(r, power - 1.0) * power * dr + 1.0;
        
        // Iteration: z = z^n + c
        z = powN(z, power) + pos;
    }
    
    return 0.5 * log(r) * r / dr;
}

Quaternion Julia Distance Estimation

For quaternion Julia sets:

float quaternionJuliaDE(vec3 pos, vec4 c) {
    vec4 z = vec4(pos, 0.0);
    float dr = 1.0;
    float r = 0.0;
    
    for (int i = 0; i < iterations; i++) {
        r = length(z);
        if (r > 4.0) break;
        
        // Derivative: dr = 2 * r * dr
        dr = 2.0 * r * dr;
        
        // Iteration: z = z^2 + c
        z = qsqr(z) + c;
    }
    
    return 0.5 * log(r) * r / dr;
}

Optimization Techniques

Bounding Spheres

Use bounding spheres to skip empty space:

float boundingSphere(vec3 pos, vec3 center, float radius) {
    return length(pos - center) - radius;
}

float sceneDE(vec3 pos) {
    // Quick rejection
    float bound = boundingSphere(pos, vec3(0.0), 2.0);
    if (bound > 0.1) return bound;
    
    // Full distance estimation
    return fractalDE(pos);
}

Adaptive Step Sizes

Adjust step size based on distance:

float t = 0.0;
for (int i = 0; i < maxSteps; i++) {
    vec3 p = ro + rd * t;
    float d = sceneDE(p);
    
    if (d < epsilon) {
        // Hit!
        return t;
    }
    
    // Adaptive step: use distance but clamp
    t += max(d * 0.8, minStep);
    
    if (t > maxDist) break;
}

Early Termination

Stop iterating when we know we're outside:

for (int i = 0; i < maxIterations; i++) {
    r = length(z);
    
    // Early termination
    if (r > bailout) break;
    if (r < minRadius) break;  // Inside, but diverged
    
    // Continue iteration...
}

Level of Detail

Reduce iterations for distant objects:

int iterations = maxIterations;
float dist = length(cameraPos - objectPos);
if (dist > farDistance) {
    iterations = maxIterations / 2;
}

Accuracy vs Performance

Trade-offs

More iterations: More accurate, slower
Higher bailout: More accurate, slower
Smaller epsilon: More accurate, slower
Adaptive steps: Better performance, similar accuracy

Practical Guidelines

Mandelbulb: 8-12 iterations, bailout 4.0
Julia sets: 12-16 iterations, bailout 4.0
Geometric fractals: 4-8 iterations (exact distance)
IFS fractals: 6-10 iterations

Common Pitfalls

Division by Zero

Always check for small values:

if (r < 0.0001) return 0.0;
return 0.5 * log(r) * r / dr;

Derivative Overflow

Clamp derivative to prevent overflow:

dr = min(dr, 1e10);

Precision Issues

Use appropriate precision qualifiers:

precision highp float;  // For detailed fractals

Advanced Techniques

Hybrid Distance Estimation

Combine multiple estimation methods:

float hybridDE(vec3 pos) {
    float d1 = mandelbulbDE(pos);
    float d2 = juliaDE(pos);
    return min(d1, d2);  // Union
}

Distance Field Smoothing

Smooth transitions between regions:

float smoothMin(float a, float b, float k) {
    float h = clamp(0.5 + 0.5 * (b - a) / k, 0.0, 1.0);
    return mix(b, a, h) - k * h * (1.0 - h);
}

References

"Distance Estimation Method for Rendering 3D Fractals" by Inigo Quilez
"The Science of Fractal Images" by Peitgen and Saupe
Inigo Quilez - Distance Estimation

Glass