Cornette-Shanks Phase Function in GLSL

This article provides an overview and practical implementation of the Cornette-Shanks phase function used for simulating atmospheric aerosol scattering in real-time computer graphics. It details the underlying mathematical equation, contrasts its behavior with other scattering models, and presents a production-ready, optimized GLSL implementation suitable for volumetric rendering and sky shaders.

Mathematical Formulation

The Cornette-Shanks phase function (introduced in 1992) is an analytical approximation of Mie scattering designed to address shortcomings in the widely used Henyey-Greenstein phase function. While Henyey-Greenstein models forward and backward scattering effectively, it does not reduce to the Rayleigh phase function when the asymmetry parameter approaches zero. Cornette-Shanks incorporates a \((1 + \cos^2\theta)\) term to properly mirror Rayleigh-like dipolar behavior for small particles while retaining directional forward-scattering lobes for larger aerosols.

The normalized formulation of the phase function \(P(\theta)\) over a sphere of \(4\pi\) steradians is:

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

Where:

  • \(\theta\) is the scattering angle between the incoming light direction and the viewing direction.
  • \(\cos\theta\) is the dot product of the light ray direction vector and the observer ray direction vector.
  • \(g\) is the asymmetry parameter, defined on the interval \((-1, 1)\). Positive values represent forward scattering (typical for aerosols and mist), zero represents isotropic/Rayleigh-like symmetric scattering, and negative values represent backward scattering.

GLSL Implementation

Below is an optimized implementation of the Cornette-Shanks phase function written in GLSL:

#define PI 3.14159265359

/**
 * Computes the Cornette-Shanks phase function for aerosol scattering.
 * 
 * @param cosTheta The cosine of the angle between the light direction and the view direction.
 * @param g The asymmetry parameter in the range (-1.0, 1.0).
 * @return The scattering intensity value normalized over 4*PI.
 */
float cornetteShanksPhase(float cosTheta, float g) {
    float g2 = g * g;
    float cosTheta2 = cosTheta * cosTheta;
    
    // Normalization and shape coefficient
    float factor = (3.0 / (8.0 * PI)) * ((1.0 - g2) / (2.0 + g2));
    
    // Angular numerator and denominator
    float num = 1.0 + cosTheta2;
    float denom = pow(1.0 + g2 - 2.0 * g * cosTheta, 1.5);
    
    return factor * (num / max(denom, 1e-6));
}

Precomputed Variant

If \(g\) is constant across the entire screen or scene (such as in a global atmosphere setting), the terms independent of \(\cos\theta\) can be precalculated on the CPU and passed as uniforms:

// Uniforms calculated once on the CPU:
// uniform float u_CSFactor; // (3.0 / (8.0 * PI)) * ((1.0 - g2) / (2.0 + g2))
// uniform float u_G;        // g
// uniform float u_G2Plus1;  // 1.0 + g * g

float cornetteShanksPhaseOptimized(float cosTheta, float csFactor, float g, float g2Plus1) {
    float num = 1.0 + cosTheta * cosTheta;
    float denom = pow(g2Plus1 - 2.0 * g * cosTheta, 1.5);
    return csFactor * (num / max(denom, 1e-6));
}

Practical Considerations

  1. Direction Vectors: Ensure the vectors used to calculate \(\cos\theta\) are normalized. When computing the dot product, cosTheta = dot(rayDir, lightDir). Depending on your light direction convention (pointing toward the light source vs. pointing along the ray of light), you may need to invert the sign of cosTheta.
  2. Singularities: When \(|g| \to 1\) and \(\cos\theta \to 1\), the denominator approaches zero. Clamping the denominator using max(denom, 1e-6) prevents division by zero artifacts when looking directly into a concentrated light source.
  3. Parameter Tuning: Aerosols and haze typically use \(g\) values between 0.7 and 0.85, which create a pronounced forward-scattering halo around light sources without severe numerical instability.