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

Related Articles

• Mandelbulb - Practical distance estimation example
• Hybrid Fractals - Combining distance estimates
• GLSL 3D Fractals Series