Mie Scattering Forward-Lobe Phase Function in GLSL

This article details the mathematical formulation and GLSL implementation used to calculate the forward-lobe intensity of Mie scattering in real-time atmospheric shaders. In atmospheric rendering, large particulate matter such as aerosols, dust, and water droplets causes Mie scattering, which directs light strongly along the incoming path of illumination. You will learn the exact formulas used to approximate this effect—primarily the Henyey-Greenstein phase function and the Cornette-Shanks variation—and how to implement them efficiently in GLSL.

The Henyey-Greenstein Phase Function

Because true Mie scattering involves complex infinite series solutions via Maxwell's equations (Lorentz-Mie theory), real-time graphics rely on empirical phase functions. The standard formula used to evaluate the forward-lobe intensity in GLSL atmospheric simulation is the Henyey-Greenstein phase function:

\[P_{HG}(\theta, g) = \frac{1}{4\pi} \frac{1 - g^2}{(1 + g^2 - 2g \cos\theta)^{3/2}}\]

Variables Defined:

  • \(\theta\) (Theta): The scattering angle between the incoming light direction vector (\(\vec{L}\)) and the view direction vector (\(\vec{V}\)).
  • \(\cos\theta\): The dot product of the normalized light direction and view direction (\(\cos\theta = \vec{L} \cdot \vec{V}\)).
  • \(g\) (Asymmetry Parameter): A coefficient strictly in the range \([-1, 1]\) defining the scattering directionality:
    • \(g > 0\): Forward scattering (characteristic of Mie scattering; typically set between \(0.75\) and \(0.99\) for hazes and sun halos).
    • \(g = 0\): Isotropic scattering (equal in all directions).
    • \(g < 0\): Back-scattering.

The Cornette-Shanks Variation

For higher physical accuracy without significant performance cost, atmospheric renderers often replace the standard Henyey-Greenstein equation with the Cornette-Shanks phase function. It provides a physically plausible approximation of Mie scattering by factoring in the polarization-dependent \((1 + \cos^2\theta)\) term:

\[P_{CS}(\theta, g) = \frac{3}{8\pi} \frac{(1 - g^2)(1 + \cos^2\theta)}{(2 + g^2)(1 + g^2 - 2g \cos\theta)^{3/2}}\]

This function retains the pronounced forward peak of the Henyey-Greenstein equation while correctly blending into Rayleigh-like behaviors at wider angles.

GLSL Implementation

Below are optimized implementations of both phase functions for fragment and compute shaders.

#define PI 3.14159265359

// Standard Henyey-Greenstein phase function
float henyeyGreenstein(float cosTheta, float g) {
    float g2 = g * g;
    float denom = 1.0 + g2 - 2.0 * g * cosTheta;
    return (1.0 / (4.0 * PI)) * ((1.0 - g2) / pow(max(denom, 0.0001), 1.5));
}

// Cornette-Shanks variation for improved atmospheric realism
float cornetteShanks(float cosTheta, float g) {
    float g2 = g * g;
    float denom = 1.0 + g2 - 2.0 * g * cosTheta;
    float numerator = (1.0 - g2) * (1.0 + cosTheta * cosTheta);
    return (3.0 / (8.0 * PI)) * (numerator / ((2.0 + g2) * pow(max(denom, 0.0001), 1.5)));
}

Integration in the Raymarching Loop

In an atmospheric raymarcher, evaluate the scattering angle at every sample point using:

vec3 lightDir = normalize(sunPosition - samplePosition);
vec3 viewDir  = normalize(rayDirection);

// Compute cosine of the scattering angle
float cosTheta = dot(viewDir, lightDir);

// Evaluate Mie phase intensity using an asymmetry factor (e.g., g = 0.8)
float miePhase = cornetteShanks(cosTheta, 0.8);

// Combine with the Mie scattering coefficient and in-scattered light
vec3 inScatteredMie = incomingLight * mieExtinctionCoefficient * miePhase;

A clamping threshold (such as max(denom, 0.0001)) prevents division by zero when looking directly into the sun (\(\cos\theta \to 1\)) when \(g \to 1\).