Lommel-Seeliger Scattering in GLSL

This article provides a concise guide to implementing the Lommel-Seeliger scattering law in GLSL for rendering low-albedo, particulate celestial bodies such as asteroids, the Moon, and Mercury. You will learn the core mathematical formulation, understand the physical justification behind its use over standard Lambertian reflectance, and see a practical GLSL code snippet ready for integration into a fragment shader.

The Lommel-Seeliger Formula

Unlike Lambertian reflectance, which models an ideal diffuse surface and exhibits pronounced limb darkening, the Lommel-Seeliger model accounts for single-scattering in porous, particulate regoliths where low albedo makes multiple scattering negligible. It simulates the uniform disc brightness observed on the Moon and rocky asteroids.

The basic intensity \(I\) defined by the Lommel-Seeliger law is:

\[I = \frac{\mu_0}{\mu_0 + \mu}\]

Where:

• \(\mu_0 = \cos(i) = \mathbf{N} \cdot \mathbf{L}\) (the cosine of the angle of incidence between the surface normal \(\mathbf{N}\) and light direction \(\mathbf{L}\))
• \(\mu = \cos(e) = \mathbf{N} \cdot \mathbf{V}\) (the cosine of the emission angle between the surface normal \(\mathbf{N}\) and view direction \(\mathbf{V}\))

GLSL Implementation

In a GLSL fragment shader, evaluating this formula requires taking the clamped dot products of the normalized surface normal with both the light vector and view vector, while avoiding division by zero at grazing angles.

// Evaluates the Lommel-Seeliger scattering factor
float evalLommelSeeliger(vec3 normal, vec3 lightDir, vec3 viewDir) {
    float NdotL = max(dot(normal, lightDir), 0.0);
    float NdotV = max(dot(normal, viewDir), 0.0);
    
    float denom = NdotL + NdotV;
    
    // Guard against division by zero on unlit or back-facing areas
    if (denom <= 0.00001) {
        return 0.0;
    }
    
    return NdotL / denom;
}

Applying Albedo and Phase Functions

To accurately render a planetary surface, the Lommel-Seeliger term is typically combined with the body's geometric albedo (\(\varpi\)) and an empirical phase function \(P(\alpha)\), where \(\alpha\) is the phase angle (\(\cos(\alpha) = \mathbf{L} \cdot \mathbf{V}\)):

vec3 computePlanetaryReflectance(vec3 N, vec3 L, vec3 V, vec3 albedo, float phaseFactor) {
    float scattering = evalLommelSeeliger(N, L, V);
    return albedo * scattering * phaseFactor;
}

By computing \(N \cdot L / (N \cdot L + N \cdot V)\), the shader prevents artificial darkening toward the edges of spherical planetary meshes, properly reproducing the flat, disc-like appearance of low-albedo regolith under full illumination.