Atmospheric Rayleigh and Mie Scattering in GLSL

Atmospheric rendering in real-time computer graphics relies on simulating how sunlight interacts with particles suspended in a planetary atmosphere. Rayleigh scattering describes light interacting with air molecules smaller than the wavelength of light, producing the blue sky and reddish sunsets, while Mie scattering describes interaction with larger aerosol particles, yielding the haze and white glare around the sun. In GLSL, this phenomenon is simulated through spherical numerical integration—commonly known as ray marching—where shaders evaluate the radiative transfer equation across view rays passing through concentric planetary and atmospheric spheres to compute optical depth, transmittance, and in-scattered radiance.

Planetary Geometry and Atmospheric Density

The physical model treats the planet as a solid sphere of radius \(R_p\) concentric with an outer atmospheric boundary sphere of radius \(R_a\). The altitude \(h\) of any point \(\mathbf{x}\) above the planetary surface is calculated as:

\[h = \|\mathbf{x}\| - R_p\]

Atmospheric density decreases exponentially with altitude according to the barometric formula. Two distinct particle densities are computed using their respective scale heights (\(H_R \approx 8000\,\text{m}\) for air molecules, and \(H_M \approx 1200\,\text{m}\) for aerosols):

\[\rho_R(h) = \exp\left(-\frac{h}{H_R}\right), \quad \rho_M(h) = \exp\left(-\frac{h}{H_M}\right)\]

In GLSL, ray-sphere intersection tests determine the entry point \(\mathbf{p}_{\text{start}}\) and exit point \(\mathbf{p}_{\text{end}}\) of the camera ray within the atmospheric shell:

bool raySphereIntersect(vec3 ro, vec3 rd, float rad, out float t0, out float t1) {
    float b = dot(ro, rd);
    float c = dot(ro, ro) - rad * rad;
    float d = b * b - c;
    if (d < 0.0) return false;
    float sqrtD = sqrt(d);
    t0 = -b - sqrtD;
    t1 = -b + sqrtD;
    return true;
}

Phase Functions and Scattering Coefficients

The directional distribution of scattered light is governed by phase functions dependent on the scattering angle \(\theta\) (the angle between the view direction \(\mathbf{v}\) and the sun direction \(\mathbf{l}\), where \(\mu = \cos\theta = \mathbf{v} \cdot \mathbf{l}\)).

Rayleigh Phase Function

Rayleigh scattering radiates almost symmetrically forward and backward:

\[P_R(\mu) = \frac{3}{16\pi} (1 + \mu^2)\]

The wavelength-dependent scattering coefficient \(\beta_R\) dictates the color of the atmosphere:

\[\beta_R = (5.8 \times 10^{-6}, 13.5 \times 10^{-6}, 33.1 \times 10^{-6})\,\text{m}^{-1}\]

Mie Phase Function

Mie scattering produces strong forward scattering, typically approximated using the Cornette-Shanks or Henyey-Greenstein phase function with an asymmetry factor \(g \in [0.75, 0.99]\):

\[P_M(\mu) = \frac{3}{8\pi} \frac{(1 - g^2)(1 + \mu^2)}{(2 + g^2)(1 + g^2 - 2g\mu)^{3/2}}\]

The Mie scattering coefficient \(\beta_M\) is generally constant across visible wavelengths:

\[\beta_M \approx 21.0 \times 10^{-6}\,\text{m}^{-1}\]

Spherical Integration via Ray Marching

The total in-scattered light reaching the camera along a ray of length \(L\) requires computing a double integral. The outer integral samples points along the camera view ray, while the inner integral computes the optical depth along the light ray from each sample point to the top of the atmosphere.

For a point along the view ray, optical depth represents the integral of density:

\[\tau(A, B) = \beta_R \int_A^B \rho_R(h(s))\,ds + \beta_M \cdot 1.11 \int_A^B \rho_M(h(s))\,ds\]

Transmittance (attenuation) is derived using the Beer-Lambert law:

\[T(A, B) = \exp(-\tau(A, B))\]

The accumulated in-scattering integral is formulated as:

\[I = \int_0^L T(\mathbf{x}_0, \mathbf{x}(s)) \cdot T(\mathbf{x}(s), \mathbf{x}_{\text{sun}}) \cdot \left[ \beta_R \rho_R(h(s)) P_R(\mu) + \beta_M \rho_M(h(s)) P_M(\mu) \right] ds\]

GLSL Implementation

The following fragment shader function executes the spherical integration using numerical Riemann sums:

#define PRIMARY_STEPS 32
#define LIGHT_STEPS 8

const float PI = 3.14159265359;
const vec3 BETA_R = vec3(5.8e-6, 13.5e-6, 33.1e-6);
const vec3 BETA_M = vec3(21.0e-6);
const float H_R = 8000.0;
const float H_M = 1200.0;
const float G = 0.76;

vec3 integrateAtmosphere(vec3 rayOrigin, vec3 rayDir, vec3 sunDir, float tMin, float tMax, float rPlanet, float rAtmo) {
    float stepSize = (tMax - tMin) / float(PRIMARY_STEPS);
    vec3 currentPos = rayOrigin + rayDir * (tMin + stepSize * 0.5);
    
    vec3 totalRayleigh = vec3(0.0);
    vec3 totalMie = vec3(0.0);
    
    float optDepthR = 0.0;
    float optDepthM = 0.0;
    
    float cosTheta = dot(rayDir, sunDir);
    float cosTheta2 = cosTheta * cosTheta;
    
    // Phase functions
    float phaseR = (3.0 / (16.0 * PI)) * (1.0 + cosTheta2);
    float phaseM = (3.0 / (8.0 * PI)) * ((1.0 - G * G) * (1.0 + cosTheta2)) / 
                   ((2.0 + G * G) * pow(1.0 + G * G - 2.0 * G * cosTheta, 1.5));
    
    for (int i = 0; i < PRIMARY_STEPS; ++i) {
        float h = length(currentPos) - rPlanet;
        if (h < 0.0) break; // Intersected planet
        
        float hr = exp(-h / H_R) * stepSize;
        float hm = exp(-h / H_M) * stepSize;
        
        optDepthR += hr;
        optDepthM += hm;
        
        // Inner ray march toward the sun
        float tSun0, tSun1;
        raySphereIntersect(currentPos, sunDir, rAtmo, tSun0, tSun1);
        float sunStepSize = tSun1 / float(LIGHT_STEPS);
        vec3 sunPos = currentPos + sunDir * (sunStepSize * 0.5);
        
        float sunOptDepthR = 0.0;
        float sunOptDepthM = 0.0;
        
        for (int j = 0; j < LIGHT_STEPS; ++j) {
            float sunH = length(sunPos) - rPlanet;
            sunOptDepthR += exp(-sunH / H_R) * sunStepSize;
            sunOptDepthM += exp(-sunH / H_M) * sunStepSize;
            sunPos += sunDir * sunStepSize;
        }
        
        vec3 tau = BETA_R * (optDepthR + sunOptDepthR) + BETA_M * 1.11 * (optDepthM + sunOptDepthM);
        vec3 attenuation = exp(-tau);
        
        totalRayleigh += attenuation * hr;
        totalMie += attenuation * hm;
        
        currentPos += rayDir * stepSize;
    }
    
    return (totalRayleigh * BETA_R * phaseR + totalMie * BETA_M * phaseM);
}

Performance Optimizations

Direct execution of \(N \times M\) nested integration loops per pixel is computationally intensive for real-time applications. Common real-time optimizations include:

  1. Look-Up Tables (LUTs): Precomputing the optical depth \(\tau(\mathbf{x}, \mathbf{l})\) as a 2D texture parameterized by altitude and sun zenith angle eliminates the inner ray march entirely, reducing the runtime complexity from \(\mathcal{O}(N \times M)\) to \(\mathcal{O}(N)\).
  2. Multiple Scattering Approximations: Single scattering ignores secondary illumination from the atmosphere and ground reflections. Real-time engines account for higher-order scattering by compounding the precomputed transmittance using dual-scattering terms or parameterizing an ambient energy-preservation factor.
  3. Temporal Sampling: Offsetting the ray march sample positions with interleaved gradient noise or blue noise, combined with temporal anti-aliasing (TAA), allows equivalent visual quality using as few as 8 to 16 primary steps.