Contact-Hardening Shadows for Point Lights in GLSL

This article details how to implement contact-hardening soft shadows for omnidirectional point lights in GLSL using Percentage-Closer Soft Shadows (PCSS) adapted for cubemaps. It covers the mathematical principles of penumbra estimation using radial distance, the generation of 3D sample directions around a cubemap ray, the blocker search phase, and the adaptive percentage-closer filtering (PCF) pass to produce realistic, distance-dependent shadow softening.

Fundamentals of Point Light PCSS

Traditional PCSS assumes a 2D planar projection from a directional or spot light. Point lights emit light radially in all directions, requiring omnidirectional shadow mapping typically backed by an OpenGL or Vulkan samplerCube.

Standard depth buffers store non-linear projected \(z\)-values, but point light shadow algorithms rely on linear Euclidean distance:

\[\text{distance} = \| \mathbf{p}_{\text{fragment}} - \mathbf{p}_{\text{light}} \|\]

The cubemap stores either this distance normalized by the light's maximum range \([0, 1]\) or raw floating-point distances in a texture format such as GL_R32F.

Contact hardening scales the penumbra width (\(w_{\text{penumbra}}\)) using similar triangles:

\[w_{\text{penumbra}} = \frac{d_{\text{receiver}} - d_{\text{blocker}}}{d_{\text{blocker}}} \cdot w_{\text{light}}\]

Where:

\(d_{\text{receiver}}\) is the linear distance from the point light to the current fragment.
\(d_{\text{blocker}}\) is the average linear distance of occluding surfaces between the light and the fragment.
\(w_{\text{light}}\) is the physical or simulated radius of the spherical light source.

Constructing the Sampling Basis

Sampling a cubemap requires a 3D direction vector. Offsetting this direction across a virtual disk requires building an orthonormal tangent basis around the normalized direction vector \(\mathbf{L} = \frac{\mathbf{p}_{\text{fragment}} - \mathbf{p}_{\text{light}}}{\| \mathbf{p}_{\text{fragment}} - \mathbf{p}_{\text{light}} \|}\).

void createOrthonormalBasis(in vec3 dir, out vec3 tangent, out vec3 bitangent) {
    vec3 up = abs(dir.y) < 0.999 ? vec3(0.0, 1.0, 0.0) : vec3(0.0, 0.0, 1.0);
    tangent = normalize(cross(up, dir));
    bitangent = cross(dir, tangent);
}

Predefined 2D sampling patterns, such as Poisson disks or Vogel disks, map onto this local frame to yield perturbed 3D sample directions:

vec3 getSampleDirection(vec3 dir, vec3 tangent, vec3 bitangent, vec2 diskSample, float radius) {
    return normalize(dir + (tangent * diskSample.x + bitangent * diskSample.y) * radius);
}

Step 1: Cubemap Blocker Search

The blocker search samples the cubemap within a search cone proportional to the light size and receiver distance. The search radius is determined by:

\[R_{\text{search}} = \frac{w_{\text{light}}}{d_{\text{receiver}}}\]

During this pass, samples where the stored depth is less than \(d_{\text{receiver}} - \text{bias}\) represent blockers.

const int BLOCKER_SAMPLES = 16;
const vec2 POISSON_DISK[16] = vec2[](
    vec2(-0.94201624, -0.39906216), vec2(0.94558609, -0.76890725),
    vec2(-0.094184101, -0.92938870), vec2(0.34495938, 0.29387760),
    vec2(-0.91588581, 0.45771432), vec2(-0.81544232, -0.87912464),
    vec2(-0.38277543, 0.27676845), vec2(0.97484398, 0.75648379),
    vec2(0.44323325, -0.97511554), vec2(0.53742981, -0.47373420),
    vec2(-0.26496911, -0.41893023), vec2(0.79197514, 0.19090188),
    vec2(-0.24188840, 0.99706507), vec2(-0.81409955, 0.91437590),
    vec2(0.19984126, 0.78641367), vec2(0.14383161, -0.14100790)
);

bool findBlocker(
    samplerCube shadowMap,
    vec3 lightToFrag,
    float receiverDist,
    float lightSize,
    float bias,
    out float avgBlockerDepth
) {
    vec3 dir = normalize(lightToFrag);
    vec3 tangent, bitangent;
    createOrthonormalBasis(dir, tangent, bitangent);

    float searchRadius = (lightSize / receiverDist);
    float blockerSum = 0.0;
    int numBlockers = 0;

    for (int i = 0; i < BLOCKER_SAMPLES; ++i) {
        vec3 sampleDir = getSampleDirection(dir, tangent, bitangent, POISSON_DISK[i], searchRadius);
        float sampledDist = texture(shadowMap, sampleDir).r;

        if (sampledDist < receiverDist - bias) {
            blockerSum += sampledDist;
            numBlockers++;
        }
    }

    if (numBlockers == 0) {
        avgBlockerDepth = 0.0;
        return false;
    }

    avgBlockerDepth = blockerSum / float(numBlockers);
    return true;
}

Step 2: Penumbra Estimation

If no blockers are detected, the fragment is completely unoccluded (shadow factor equals 1.0). If blockers exist, the ratio between the receiver distance and the average blocker distance determines the dynamic filter kernel size.

float calculatePenumbraRadius(float receiverDist, float avgBlockerDist, float lightSize) {
    return ((receiverDist - avgBlockerDist) / avgBlockerDist) * (lightSize / receiverDist);
}

Dividing by receiverDist adjusts the angular spread of the filter radius to maintain consistent screen-space sizing regardless of the distance between the light source and the geometry.

Step 3: Adaptive Filtering (PCF Pass)

Using the derived penumbra radius, the cubemap is sampled again using a larger sample count. A depth comparison is performed per sample, accumulating the final occlusion factor.

const int PCF_SAMPLES = 32;

float filterShadowPCF(
    samplerCube shadowMap,
    vec3 lightToFrag,
    float receiverDist,
    float filterRadius,
    float bias
) {
    vec3 dir = normalize(lightToFrag);
    vec3 tangent, bitangent;
    createOrthonormalBasis(dir, tangent, bitangent);

    float shadow = 0.0;
    for (int i = 0; i < PCF_SAMPLES; ++i) {
        // Vogel disk or larger Poisson disk sample
        float angle = float(i) * 2.3999632; // Golden angle in radians
        float r = sqrt(float(i) + 0.5) / sqrt(float(PCF_SAMPLES));
        vec2 diskSample = vec2(cos(angle), sin(angle)) * r;

        vec3 sampleDir = getSampleDirection(dir, tangent, bitangent, diskSample, filterRadius);
        float sampledDist = texture(shadowMap, sampleDir).r;

        shadow += (sampledDist < receiverDist - bias) ? 0.0 : 1.0;
    }

    return shadow / float(PCF_SAMPLES);
}

Complete Fragment Integration

Combining these steps produces the complete GLSL execution sequence for contact-hardening point light shadows:

float computePointPCSS(
    samplerCube shadowMap,
    vec3 fragPos,
    vec3 lightPos,
    float lightSize,
    float normalBias
) {
    vec3 lightToFrag = fragPos - lightPos;
    float receiverDist = length(lightToFrag);

    float avgBlockerDist;
    bool hasBlocker = findBlocker(
        shadowMap, 
        lightToFrag, 
        receiverDist, 
        lightSize, 
        normalBias, 
        avgBlockerDist
    );

    if (!hasBlocker) {
        return 1.0; // Fully lit
    }

    float penumbraRadius = calculatePenumbraRadius(receiverDist, avgBlockerDist, lightSize);
    
    // Clamp to prevent sub-pixel undersampling or excessive blur
    penumbraRadius = clamp(penumbraRadius, 0.001, 0.1);

    return filterShadowPCF(
        shadowMap, 
        lightToFrag, 
        receiverDist, 
        penumbraRadius, 
        normalBias
    );
}

Performance and Artifact Mitigation

Self-Shadowing and Dynamic Bias: Point lights suffer from divergence artifacts at the corners of cubemap faces. The shadow bias must scale proportionally with the surface normal and the light vector angle:
```
float bias = max(maxBias * (1.0 - dot(surfaceNormal, -dir)), minBias);
```
Cubemap Seams: Filtering across cubemap edges can cause texture coordinate discontinuity. GLSL handles hardware cubemap seamless filtering via glEnable(GL_TEXTURE_CUBE_MAP_SEAMLESS).
Early Exits: If the average blocker distance is extremely close to the receiver distance or if all samples in the blocker search report occlusion, the PCF pass can early-exit with a value of 0.0 to minimize divergent texture fetches.