Disc Area Light Irradiance with GLSL Polynomials

Evaluating the diffuse irradiance from a disc area light in real-time rendering requires integrating the clamped cosine of incoming radiance across the disc's subtended solid angle. While the analytical solution involves computationally heavy operations like elliptic integrals and trigonometric horizon-clipping, real-time engines replace these evaluations with low-order polynomial curve fits. This article explains how disc light geometry is parameterized and how polynomial approximations are implemented in GLSL to compute diffuse irradiance efficiently.

The Geometric Formulation

To approximate disc irradiance, the disc's geometry is reduced to a set of scalar parameters relative to the shaded surface point \(P\) with surface normal \(N\):

  1. Light Vector (\(L\)) and Distance (\(d\)): The vector from \(P\) to the disc center \(C\), where \(d = \|C - P\|\) and \(L = (C - P) / d\).
  2. Disc Orientation Factor: The alignment of the disc normal \(N_{\text{disc}}\) relative to \(-L\), given by \(\cos \theta_{\text{disc}} = \max(0.0, -N_{\text{disc}} \cdot L)\).
  3. Subtended Angular Radius (\(\sin \alpha\)): For a disc of radius \(R\), the sine of the half-angle subtended by the light source is: \[\sin \alpha = \frac{R}{\sqrt{d^2 + R^2}}\]
  4. Surface Alignment (\(N \cdot L\)): The primary cosine factor between the surface normal and the direction to the disc center.

Polynomial Fitting Strategy

The total irradiance received by a differential surface element from a disc light is governed by:

\[E = \int_{\Omega} L_i \max(0, N \cdot \omega) \, d\omega\]

Direct integration is piecewise due to the horizon boundary (\(N \cdot \omega = 0\)). To avoid runtime branching and expensive trigonometric operations, the integral is precomputed offline over normalized ranges of \(x = N \cdot L\) and \(y = \sin \alpha\), then fitted to a polynomial surface.

A widely used formulation fits the diffuse response using a low-degree polynomial that approximates the smoothing effect of the area light near the horizon (\(N \cdot L \approx 0\)). Instead of an abrupt cut-off, the light wraps smoothly around the surface based on the angular size \(\sin \alpha\):

\[E \approx \frac{\max(0.0, x + y)}{1.0 + y} \cdot \left( a x^2 + b x y + c y^2 + d x + e y \right)\]

Alternatively, the simplified form popularized by real-time physically based shading pipelines uses a clamped quadratic or cubic fit to model the attenuation and horizon factor:

\[f(x, y) = \text{clamp}\left(\frac{x + y^2}{1.0 + y^2}, 0.0, 1.0\right) \cdot \left(1.0 - y^2 \cdot (1.0 - x)\right)\]

This scales the irradiance smoothly from a point-source response (\(y \to 0\)) to a hemispherical wrap (\(y \to 1\)).

GLSL Implementation

Below is a complete GLSL function that computes the diffuse irradiance from a disc area light using a polynomial approximation:

// Computes diffuse irradiance from a disc area light
// P: Shading point world position
// N: Shading point normal (normalized)
// discPos: Center of disc light
// discNorm: Normal vector of the disc (normalized, pointing forward)
// discRadius: Radius of the disc light in world units
// discLuminance: Radiometric power / emission of the disc
vec3 evaluateDiscLightIrradiance(
    vec3 P, 
    vec3 N, 
    vec3 discPos, 
    vec3 discNorm, 
    float discRadius, 
    vec3 discLuminance
) {
    vec3 dir = discPos - P;
    float distSq = dot(dir, dir);
    float dist = sqrt(distSq);
    vec3 L = dir / dist;

    // Cosine of angle between disc normal and incoming direction to surface
    float cosDisc = max(dot(-discNorm, L), 0.0);
    if (cosDisc <= 0.0) {
        return vec3(0.0);
    }

    // Subtended half-angle sine approximation: sin(alpha)
    float sinAlpha = discRadius / sqrt(distSq + discRadius * discRadius);
    float sinAlphaSq = sinAlpha * sinAlpha;

    // Angle between surface normal and vector to disc center
    float NoL = dot(N, L);

    // Polynomial approximation of the clamped cosine integral over the disc
    // Horizon wrap factor parameterized by angular radius
    float x = NoL;
    float y = sinAlpha;

    // Polynomial fit: smoothly extends irradiance over the geometric horizon
    float horizonTerm = clamp((x + y) / (1.0 + y), 0.0, 1.0);
    
    // Quadratic correction factor to account for solid angle expansion
    float polyWeight = horizonTerm * (0.5 * x * y + (1.0 - 0.5 * y));
    float diffuseFactor = max(0.0, x) + polyWeight * (1.0 - max(0.0, x));

    // Distance falloff with radius offset to prevent singularity at zero distance
    float attenuation = 1.0 / (distSq + discRadius * discRadius);

    // Projected solid angle scale factor
    float solidAngleFactor = sinAlphaSq * cosDisc;

    return discLuminance * (diffuseFactor * solidAngleFactor * attenuation);
}

Performance Benefits in Shaders

Using polynomial fits for disc lights offers three main performance advantages:

  • Elimination of Transcendental Functions: Eliminates expensive acos, atan, and square roots within integration loops, reducing ALU cycle counts.
  • Vectorization and Execution Uniformity: Polynomial evaluations contain no branch divergence, maintaining uniform execution across shader warps/wavefronts.
  • Continuous Horizon Handling: The smooth mathematical boundary of polynomials eliminates hard clipping artifacts at the geometric horizon without requiring manual super-sampling.