Planetary Atmosphere Multi-Scattering in GLSL

Simulating realistic planetary atmospheres requires accounting for light that bounces multiple times between air molecules and aerosol particles before reaching the camera. This article outlines the primary algorithm used to compute multi-scattering diffuse approximations in real-time GLSL shaders: Sébastien Hillaire’s energy-conserving multiple scattering approximation. It details the underlying mathematical formulation, the generation of the precomputed multi-scattering Look-Up Table (LUT), and how GLSL shaders sample this approximation during standard single-scattering raymarching.

The Algorithm: Hillaire’s Multi-Scattering Approximation

Single-scattering models (such as pure Rayleigh and Mie implementations) fail to reproduce the ambient brightness of twilight, the diffuse glow on the planet’s unlit side, and true energy conservation. While simulating infinite scattering bounces traditionally requires offline Monte Carlo integration or complex 4D LUTs (such as Eric Bruneton’s classic precomputed atmospheric scattering), modern real-time rendering relies on Sébastien Hillaire’s 2020 algorithm: "A Scalable and Production Ready Sky and Atmosphere Rendering Technique."

Hillaire's algorithm simplifies multiple scattering by assuming that light after the first bounce becomes increasingly isotropic (diffuse). Instead of tracking directional phase functions across infinite bounces, it treats secondary and higher-order scattering as a diffuse energy transfer represented by an infinite geometric series.

Mathematical Formulation

When light scatters once, a fraction of the incoming radiance is absorbed or scattered away. The fraction of scattered light that remains in the atmosphere across an entire sphere of directions is denoted as \(f_{\text{ms}}\) (the multi-scattering factor).

Assuming each subsequent bounce behaves isotropically, the total radiance \(L_{\text{total}}\) accumulated along a path is the sum of the primary single scattering \(L_{\text{single}}\) and the multi-scattered radiance \(L_{\text{multi}}\):

\[L_{\text{multi}} = L_{2\text{nd}} \cdot \frac{1}{1 - f_{\text{ms}}}\]

Where:

  • \(L_{2\text{nd}}\) is the radiance produced by second-order scattering across a spherical integration.
  • \(f_{\text{ms}}\) is the average single-scattering albedo integrated over the unit sphere.
  • \(\frac{1}{1 - f_{\text{ms}}}\) is the sum of the infinite geometric series representing orders 3 through infinity.

Generating the Multi-Scattering LUT in GLSL

Rather than evaluating the spherical integral dynamically every frame, the algorithm precomputes this diffuse contribution into a low-resolution 2D texture (typically \(32 \times 32\) texels).

The 2D LUT is parameterized by:

  1. Altitude (\(r\)): The height of the sample point from the planet surface to the atmosphere boundary.
  2. Sun Zenith Angle (\(\cos \theta_s\)): The cosine of the angle between the local surface normal and the directional light vector.

A compute or fragment shader calculates the 2D LUT by performing the following steps:

  1. Map UV coordinates to altitude \(r\) and sun zenith angle \(\mu_s = \cos \theta_s\).
  2. Cast a uniform distribution of sample directions over the entire sphere (typically 64 to 128 directions via Fibonacci sphere sampling).
  3. For each direction, raymarch through the atmosphere to integrate single scattering and transmittance along that sub-ray.
  4. Accumulate the incoming radiance into two buffers: the total scattered radiance (\(L_{2\text{nd}}\)) and the transfer factor (\(f_{\text{ms}}\)).
  5. Apply the infinite series factor: \(L_{\text{multi\_LUT}} = \frac{L_{2\text{nd}}}{1.0 - f_{\text{ms}}}\).

Sampling in the Primary Sky Shader

During real-time rendering, the primary camera raymarches through the atmosphere as usual. At each step along the ray:

// Sample single scattering transmittance and densities
vec3 samplePos = cameraPos + rayDir * t;
float altitude = length(samplePos) - groundRadius;
vec3 upVector = normalize(samplePos);
float sunCosTheta = dot(upVector, sunDirection);

// Compute 2D LUT coordinates
vec2 uv = mapAltitudeAndSunZenithToUV(altitude, sunCosTheta);
vec3 multiScatteringLuminance = texture(u_MultiScatteringLUT, uv).rgb;

// Combine phase-dependent single scattering with isotropic multi-scattering
vec3 rayleighPhase = vec3(rayleighPhaseFunction(dot(rayDir, sunDirection)));
vec3 miePhase = vec3(miePhaseFunction(dot(rayDir, sunDirection), mieG));

vec3 singleScattering = sunTransmittance * (rayleighScattering * rayleighPhase + mieScattering * miePhase);
vec3 multiScattering = multiScatteringLuminance * (rayleighScattering + mieScattering);

// Total in-scattered radiance for the step
vec3 stepScattering = (singleScattering + multiScattering) * stepTransmittance;

Because the multi-scattering component uses the isotropic approximation, it is multiplied directly by the scattering coefficients rather than the directional Rayleigh or Mie phase functions. This yields physically consistent ambient lighting, deep atmospheric hues, and realistic horizons in GLSL with minimal runtime overhead.