Double Henyey-Greenstein Phase Function in GLSL

The double Henyey-Greenstein phase function simulates realistic light scattering in volumetric media such as clouds, fog, and biological tissue by balancing forward and backward scattering lobes. While a standard Henyey-Greenstein function models scattering biased in a single direction, combining two instances via a linear blend allows shaders to simultaneously capture sharp forward glare and subtle backscattering effects. In OpenGL Shading Language (GLSL), this is achieved by evaluating two separate asymmetry lobes and interpolating them using a blending factor.

The Mathematical Formulation

The single Henyey-Greenstein phase function defines the angular distribution of scattered light based on the cosine of the scattering angle (\(\cos\theta\)) and an asymmetry parameter \(g \in (-1, 1)\):

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

When \(g > 0\), light predominantly scatters forward; when \(g < 0\), it predominantly scatters backward.

The double Henyey-Greenstein phase function combines two parameters—a forward lobe parameter \(g_1 > 0\) and a backward lobe parameter \(g_2 < 0\)—weighted by an interpolation factor \(w \in [0, 1]\):

\[P_{2HG}(\cos\theta, g_1, g_2, w) = w \cdot P(\cos\theta, g_1) + (1 - w) \cdot P(\cos\theta, g_2)\]

GLSL Implementation

In a GLSL fragment or compute shader, the angle cosine \(\cos\theta\) is calculated using the dot product between the incoming light ray direction and the outgoing viewing ray direction.

#define PI 3.14159265359

// Evaluates a single Henyey-Greenstein phase function lobe
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));
}

// Combines forward and backward lobes
float doubleHenyeyGreenstein(float cosTheta, float g1, float g2, float w) {
    float forwardScattering = henyeyGreenstein(cosTheta, g1);
    float backwardScattering = henyeyGreenstein(cosTheta, g2);
    return mix(backwardScattering, forwardScattering, w);
}

How the Combination Works in Shaders

  1. Directional Input: The cosine of the angle between the normalized incoming light vector L and the normalized view direction V is computed using float cosTheta = dot(L, V).
  2. Forward Lobe Evaluation: The forward scattering component uses a positive asymmetry parameter (e.g., g1 = 0.7). This creates high-intensity forward scattering around \(\cos\theta \approx 1\), mimicking the bright halo observed around light sources in foggy environments.
  3. Backward Lobe Evaluation: The backward scattering component uses a negative asymmetry parameter (e.g., g2 = -0.3). This boosts reflectance near \(\cos\theta \approx -1\), mimicking retroreflection or diffuse back-lighting.
  4. Energy Blending: The GLSL mix() function blends the two components according to weight w. Because both sub-functions integrate to 1 over the unit sphere, their convex combination also integrates to 1, preserving energy conservation across the medium.