Linearly Transformed Cosines for Polygonal Lights in GLSL
This article provides a mathematical breakdown and GLSL shader implementation of Linearly Transformed Cosines (LTC) for real-time polygonal area lights. It explains how spherical BRDF distributions are mapped to clamped cosine distributions through linear transformations, covers the application of Stokes' theorem to reduce surface integrals to analytical edge evaluations, and presents the corresponding GLSL code for horizon clipping and lighting evaluation.
Core Mathematical Formulation
Linearly Transformed Cosines (LTC) approximate complex BRDF lobes by applying a linear transformation matrix \(M\) to a standard clamped cosine distribution \(D_0\):
\[D(\omega) = \frac{D_0(M^{-1} \omega)}{\det(M) \| M^{-1} \omega \|^3}\]
Where:
- \(D_0(\omega_0) = \frac{1}{\pi} \max(0, \cos\theta_0)\) is the clamped cosine distribution.
- \(M\) is a \(3 \times 3\) matrix representing roughness, anisotropy, and view angle.
- \(\omega\) is the incoming light direction.
Evaluating the reflected radiance \(L_o\) over an arbitrary polygonal domain \(P\) involves computing:
\[E = \int_P D(\omega) (\mathbf{n} \cdot \omega) \, d\omega\]
By substituting the transformed integration domain, the integral simplifies to evaluating the clamped cosine distribution \(D_0\) over a transformed polygon \(P' = M^{-1} P\):
\[E = \int_{P'} D_0(\omega_0) (\mathbf{n}_0 \cdot \omega_0) \, d\omega_0\]
This operation shifts the complexity from evaluating an intricate BRDF profile to computing the geometric form factor of a transformed polygon under a simple Lambertian distribution.
Analytical Evaluation via Stokes' Theorem
By applying Stokes' theorem, the area integral over the spherical polygon \(P'\) projects down to a boundary line integral along the polygon’s directed edges:
\[\mathbf{F} = \frac{1}{2\pi} \sum_{i=1}^{N} \mathbf{v}_i \times \mathbf{v}_{i+1} \frac{\theta_i}{\|\mathbf{v}_i \times \mathbf{v}_{i+1}\|}\]
Where:
- \(\mathbf{v}_i\) are the normalized vertices of the transformed polygon projected onto the unit sphere: \(\mathbf{v}_i = \frac{M^{-1}(P_i - X)}{\|M^{-1}(P_i - X)\|}\) for shading point \(X\).
- \(\theta_i = \arccos(\mathbf{v}_i \cdot \mathbf{v}_{i+1})\) is the angle subtended by edge \((\mathbf{v}_i, \mathbf{v}_{i+1})\) at the origin.
- The scalar irradiance is the component aligned with the local normal \(\mathbf{n}_0 = (0, 0, 1)\):
\[E = \max(0, \mathbf{F}_z)\]
Horizon Clipping
Vertices of the transformed polygon located below the local tangent horizon (\(z < 0\)) cause erroneous back-facing illumination. To avoid integration over the lower hemisphere, edges crossing the \(z = 0\) plane must be clipped.
If an edge transitions between \(z_i \ge 0\) and \(z_{i+1} < 0\), a replacement boundary vertex is interpolated:
\[\mathbf{v}_{clip} = \mathbf{v}_i + \frac{-z_i}{z_{i+1} - z_i} (\mathbf{v}_{i+1} - \mathbf{v}_i)\]
Clipping a quad can yield a polygon with three, four, or five vertices.
GLSL Implementation
1. Edge Integration Function
This function evaluates the line integral for a single directed edge between two normalized points on the unit sphere:
vec3 IntegrateEdge(vec3 v1, vec3 v2) {
float cosTheta = dot(v1, v2);
cosTheta = clamp(cosTheta, -0.9999, 0.9999);
float theta = acos(cosTheta);
vec3 crossVec = cross(v1, v2);
float len = length(crossVec);
if (len < 1e-6) {
return vec3(0.0);
}
return crossVec * (theta / len);
}2. Horizon Clipping Implementation
This routine clips an arbitrary quad against the horizon plane \(z = 0\):
int ClipPolygonToHorizon(vec3 L[4], out vec3 L_clipped[5]) {
int n = 0;
for (int i = 0; i < 4; ++i) {
vec3 p1 = L[i];
vec3 p2 = L[(i + 1) % 4];
if (p1.z >= 0.0) {
L_clipped[n++] = p1;
}
if ((p1.z >= 0.0 && p2.z < 0.0) || (p1.z < 0.0 && p2.z >= 0.0)) {
float t = -p1.z / (p2.z - p1.z);
L_clipped[n++] = mix(p1, p2, t);
}
}
return n;
}3. Evaluating Irradiance
The complete shading function accepts the transformed polygon, executes clipping, projects points to the unit sphere, and calculates the total vector irradiance:
vec3 EvaluateLTC(mat3 Minv, vec3 points[4], vec3 P, vec3 N, vec3 V) {
// Construct local tangent space (T, B, N)
vec3 T = normalize(V - N * dot(V, N));
vec3 B = cross(N, T);
mat3 tangentToWorld = mat3(T, B, N);
mat3 worldToTangent = transpose(tangentToWorld);
// Transform polygon vertices into local LTC space
mat3 M = Minv * worldToTangent;
vec3 L[4];
for (int i = 0; i < 4; ++i) {
L[i] = M * (points[i] - P);
}
// Clip against horizon
vec3 L_clipped[5];
int n = ClipPolygonToHorizon(L, L_clipped);
if (n < 3) {
return vec3(0.0);
}
// Normalize clipped vertices onto unit sphere
for (int i = 0; i < n; ++i) {
L_clipped[i] = normalize(L_clipped[i]);
}
// Integrate edges using Stokes' theorem
vec3 vectorIrradiance = vec3(0.0);
for (int i = 0; i < n; ++i) {
int next = (i + 1) % n;
vectorIrradiance += IntegrateEdge(L_clipped[i], L_clipped[next]);
}
// Scalar irradiance over the clamped cosine (divided by 2*PI)
float scalarIrradiance = max(0.0, vectorIrradiance.z * (1.0 / (2.0 * 3.14159265)));
return vec3(scalarIrradiance);
}4. Sampling LUTs for the LTC Matrix
The inverse transformation matrix \(M^{-1}\) and the directional Fresnel magnitude are retrieved from precomputed lookup tables (LUTs) parametrized by the roughness \(\alpha\) and the view angle \(\cos\theta = \mathbf{n} \cdot \mathbf{v}\):
uniform sampler2D ltcMatTexture;
uniform sampler2D ltcAmpTexture;
void main() {
float NdotV = clamp(dot(normal, viewDir), 0.0, 1.0);
vec2 uv = vec2(roughness, sqrt(1.0 - NdotV));
uv = uv * (63.0 / 64.0) + (0.5 / 64.0); // Texture coordinate bias
vec4 t1 = texture(ltcMatTexture, uv);
vec4 t2 = texture(ltcAmpTexture, uv);
// Reconstruct M^-1
mat3 Minv = mat3(
vec3(t1.x, 0.0, t1.y),
vec3(0.0, 1.0, 0.0),
vec3(t1.z, 0.0, t1.w)
);
vec3 specularColor = mix(vec3(0.04), albedo, metallic);
vec3 specularFactor = specularColor * t2.x + (1.0 - specularColor) * t2.y;
vec3 spec = EvaluateLTC(Minv, quadVertices, worldPos, normal, viewDir) * specularFactor;
vec3 diff = EvaluateLTC(mat3(1.0), quadVertices, worldPos, normal, viewDir) * albedo;
vec3 finalRadiance = lightRadiance * (spec + diff);
}