Pre-Integrated Skin Shading with Curvature in GLSL

Pre-integrated skin shading is an efficient real-time rendering technique that approximates subsurface scattering by decoupling light transport into pre-computed lookups. Instead of performing expensive diffusion simulations at runtime, this method calculates the local surface curvature dynamically and uses it alongside the light incidence angle to sample a pre-integrated two-dimensional lookup table (LUT). This article outlines the mathematical basis of the algorithm, explains how to estimate surface curvature dynamically using screen-space derivatives, and provides a direct, production-ready GLSL implementation.

The Algorithm Concept

Subsurface scattering causes light to bleed through translucent surfaces such as human skin, softening harsh shadow boundaries and shifting shadow penumbras toward red due to blood absorption. Eric Penner’s pre-integrated skin shading algorithm models this behavior by integrating the diffusion profile over spheres of varying radii.

The scattering response is stored in a 2D texture (LUT):

• Horizontal Axis (U): Light incidence angle, mapped from the dot product of the surface normal and the light vector (\(N \cdot L\)).
• Vertical Axis (V): Local surface curvature (\(\kappa = 1 / r\)), representing how sharply the surface curves.

At runtime, the shader determines the local curvature, computes the lighting angle, and samples this LUT to yield the scattered diffuse color.

Step 1: Estimating Curvature in Screen Space

Curvature can be authored offline in texture maps or computed dynamically in the fragment shader. Dynamic estimation uses screen-space partial derivatives (dFdx and dFdy) to measure how quickly the geometric surface normal changes relative to world-space position:

\[\kappa \approx \frac{\|\nabla N\|}{\|\nabla P\|}\]

In GLSL, this is derived by comparing the length of the normal delta to the length of the position delta:

float estimateCurvature(vec3 worldNormal, vec3 worldPosition) {
    vec3 dNdx = dFdx(worldNormal);
    vec3 dNdy = dFdy(worldNormal);
    vec3 dPdx = dFdx(worldPosition);
    vec3 dPdy = dFdy(worldPosition);

    float deltaNormal = length(dNdx) + length(dNdy);
    float deltaPos = length(dPdx) + length(dPdy);

    // Prevent division by zero and scale to adjust scattering responsiveness
    return deltaNormal / max(deltaPos, 0.0001);
}

Step 2: Preparing the LUT Coordinates

The horizontal coordinate is parameterized by mapping \(N \cdot L\) from the range \([-1.0, 1.0]\) to \([0.0, 1.0]\):

\[u = 0.5 \times (N \cdot L) + 0.5\]

The vertical coordinate is parameterized by scaling the estimated curvature into the normalized \([0.0, 1.0]\) range based on the model's unit scale:

\[v = \text{clamp}(\kappa \times \text{curvatureScale}, 0.0, 1.0)\]

Step 3: Complete GLSL Shader Implementation

Below is the complete fragment shader logic for calculating the pre-integrated diffuse scattering term:

#version 330 core

in vec3 vWorldPos;
in vec3 vWorldNormal;

out vec4 fragColor;

uniform sampler2D uPreIntegratedLUT;
uniform vec3 uLightDir;       // Direction towards the light
uniform vec3 uLightColor;     // Light intensity and color
uniform float uCurvatureScale; // Global tuning multiplier (e.g., 0.01 - 0.1)

float computeCurvature(vec3 normal, vec3 position) {
    vec3 dNdx = dFdx(normal);
    vec3 dNdy = dFdy(normal);
    vec3 dPdx = dFdx(position);
    vec3 dPdy = dFdy(position);

    float deltaNormal = length(dNdx) + length(dNdy);
    float deltaPosition = length(dPdx) + length(dPdy);

    return deltaNormal / max(deltaPosition, 0.00001);
}

void main() {
    vec3 N = normalize(vWorldNormal);
    vec3 L = normalize(uLightDir);

    // 1. Calculate N dot L
    float NdotL = dot(N, L);

    // 2. Map NdotL to [0.0, 1.0] for the LUT's U coordinate
    float u = NdotL * 0.5 + 0.5;

    // 3. Compute screen-space curvature for the LUT's V coordinate
    float rawCurvature = computeCurvature(N, vWorldPos);
    float v = clamp(rawCurvature * uCurvatureScale, 0.0, 1.0);

    // 4. Sample the pre-integrated scattering profile (RGB represents wavelength diffusion)
    vec3 sssDiffuse = texture(uPreIntegratedLUT, vec2(u, v)).rgb;

    // 5. Combine with light properties
    vec3 finalDiffuse = sssDiffuse * uLightColor;

    fragColor = vec4(finalDiffuse, 1.0);
}

Integration Notes

• Shadow Handling: Standard shadow mapping produces hard edges that break subsurface illusions. To soften shadows, blend a blurred shadow map or apply Penner's separate shadow penumbra scattering pass.
• Derivative Artifacts: Screen-space derivatives can produce high-frequency noise across polygon seams or silhouette edges. Filtering the calculated curvature or combining it with an offline-baked curvature map eliminates visual noise while preserving dynamic deformation scattering.