Henyey-Greenstein Phase Function in GLSL Fog

This article explains the mathematical formula for the Henyey-Greenstein phase function and details its implementation in GLSL for real-time volumetric fog rendering. The Henyey-Greenstein phase function is widely used in computer graphics to approximate anisotropic light scattering—known as Mie scattering—in participating media like clouds, haze, and dust. Below, you will find the exact mathematical formulation, parameter explanations, and ready-to-use GLSL shader code optimized for volumetric raymarching.

Mathematical Formula

The standard Henyey-Greenstein phase function defines the probability distribution of light scattering at an angle \(\theta\) relative to the incoming light ray:

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

Where:

• \(\theta\) is the scattering angle between the direction of the light ray and the view direction (the direction light travels toward the camera).
• \(\cos\theta\) is computed using the dot product of the normalized light direction vector \(\vec{L}\) and the view ray direction \(\vec{V}\): \(\cos\theta = \vec{V} \cdot \vec{L}\).
• \(g\) is the asymmetry parameter, defined in the range \([-1, 1]\):
  • \(g > 0\): Produces forward scattering (light scatters forward, common in fog and water droplets).
  • \(g < 0\): Produces backward scattering (light bounces back toward the light source).
  • \(g = 0\): Produces isotropic scattering, reducing the formula to \(P(\theta) = \frac{1}{4\pi}\).

GLSL Implementation

In GLSL shaders, the constant factor \(\frac{1}{4\pi} \approx 0.07957747\) can be precalculated to reduce computational overhead during volumetric raymarching.

#define PI 3.14159265359
#define INV_4PI 0.07957747154594767

float henyeyGreenstein(float cosTheta, float g) {
    float g2 = g * g;
    float denom = 1.0 + g2 - 2.0 * g * cosTheta;
    return INV_4PI * ((1.0 - g2) / max(pow(denom, 1.5), 0.0001));
}

Integration in Volumetric Raymarching

To evaluate scattering at a specific sample point along a ray:

// rayDir: normalized camera ray direction (pointing away from camera)
// lightDir: normalized direction towards the light source
// Since scattering points toward the camera, we compare rayDir with lightDir
float cosTheta = dot(rayDir, lightDir);

// Evaluate phase function with an asymmetry factor (e.g., g = 0.6 for forward scattering fog)
float phase = henyeyGreenstein(cosTheta, 0.6);

// Multiply the accumulated in-scattered light by the phase factor
vec3 inScattering = lightColor * lightTransmittance * fogDensity * phase;

Dual Henyey-Greenstein Phase Function

Atmospheric phenomena often exhibit both strong forward scattering around light sources and a subtle backward scattering rim. A common industry technique is the two-lobe (dual) Henyey-Greenstein function, which combines two phase evaluations using a blend factor \(k\):

\[P_{dual}(\theta) = k \cdot P(\theta, g_1) + (1 - k) \cdot P(\theta, g_2)\]

In GLSL:

float dualHenyeyGreenstein(float cosTheta, float g1, float g2, float k) {
    return mix(henyeyGreenstein(cosTheta, g2), henyeyGreenstein(cosTheta, g1), k);
}

Using \(g_1 > 0\) (such as \(0.7\)), \(g_2 < 0\) (such as \(-0.2\)), and \(k \approx 0.8\) produces realistic atmospheric fog with both a pronounced forward silver lining and a slight backward glow.