Implementing Dipole Diffusion SSS in GLSL
This article provides an overview and technical breakdown of implementing the dipole diffusion approximation for subsurface light transport in GLSL. Developed by Jensen et al., this model simulates translucent materials—such as skin, wax, and marble—by replacing complex volumetric scattering with a pair of virtual point sources. Below, you will find the fundamental mathematical formulations, the translation of these equations into shader code, and techniques for integrating the evaluation into real-time rendering pipelines.
Mathematical Foundations
The dipole diffusion model approximates the diffuse reflectance profile \(R(r)\) as a function of the radial distance \(r\) from the point of light incidence. It accounts for multiple scattering by placing two point sources: one positive source inside the medium at distance \(z_r\) beneath the surface, and one negative virtual source above the surface at distance \(z_v\) to satisfy boundary conditions.
The key optical properties required are:
- Absorption coefficient (\(\sigma_a\))
- Scattering coefficient (\(\sigma_s\))
- Mean cosine of the scattering angle (\(g\))
From these, derive the reduced coefficients and parameters:
- Reduced scattering coefficient: \(\sigma'_s = \sigma_s (1 - g)\)
- Reduced extinction coefficient: \(\sigma'_t = \sigma_a + \sigma'_s\)
- Effective transport coefficient: \(\sigma_{tr} = \sqrt{3 \sigma_a \sigma'_t}\)
- Reduced albedo: \(\alpha' = \frac{\sigma'_s}{\sigma'_t}\)
- Diffusion constant: \(D = \frac{1}{3 \sigma'_t}\)
The source distances relative to the surface interface are defined by:
- \(z_r = \frac{1}{\sigma'_t}\)
- \(z_v = z_r + 4 A D\)
Here, \(A\) represents the boundary condition parameter derived from the index of refraction (\(\eta\)), determining internal reflection: \[A = \frac{1 + F_{dr}}{1 - F_{dr}}\] Where \(F_{dr} \approx -1.440 / \eta^2 + 0.710 / \eta + 0.668 + 0.0636 \eta\).
The distance from the evaluation point to the positive and negative sources are:
- \(d_r = \sqrt{r^2 + z_r^2}\)
- \(d_v = \sqrt{r^2 + z_v^2}\)
The diffuse reflectance profile \(R(r)\) is computed as: \[R(r) = \frac{\alpha'}{4\pi} \left[ z_r \left(\sigma_{tr} + \frac{1}{d_r}\right) \frac{e^{-\sigma_{tr} d_r}}{d_r^2} + z_v \left(\sigma_{tr} + \frac{1}{d_v}\right) \frac{e^{-\sigma_{tr} d_v}}{d_v^2} \right]\]
GLSL Implementation of the Reflectance Profile
The core mathematical function can be implemented in a reusable GLSL function that calculates the spectral reflectance given a distance \(r\):
// Computes internal reflection parameter A from index of refraction (eta)
float computeA(float eta) {
float Fdr = -1.440 / (eta * eta) + 0.710 / eta + 0.668 + 0.0636 * eta;
return (1.0 + Fdr) / (1.0 - Fdr);
}
// Evaluates the dipole diffusion profile R(r) for a single channel
float dipoleReflectance(float r, float sigma_a, float sigma_s_prime, float eta) {
float sigma_t_prime = sigma_a + sigma_s_prime;
float alpha_prime = sigma_s_prime / sigma_t_prime;
float sigma_tr = sqrt(3.0 * sigma_a * sigma_t_prime);
float D = 1.0 / (3.0 * sigma_t_prime);
float A = computeA(eta);
float z_r = 1.0 / sigma_t_prime;
float z_v = z_r + 4.0 * A * D;
float d_r = sqrt(r * r + z_r * z_r);
float d_v = sqrt(r * r + z_v * z_v);
float factor = alpha_prime / (4.0 * 3.14159265359);
float term_r = z_r * (sigma_tr + 1.0 / d_r) * (exp(-sigma_tr * d_r) / (d_r * d_r));
float term_v = z_v * (sigma_tr + 1.0 / d_v) * (exp(-sigma_tr * d_v) / (d_v * d_v));
return factor * (term_r + term_v);
}
// Vectorized wrapper for RGB optical properties
vec3 evaluateDipoleRGB(float r, vec3 sigma_a, vec3 sigma_s_prime, float eta) {
return vec3(
dipoleReflectance(r, sigma_a.r, sigma_s_prime.r, eta),
dipoleReflectance(r, sigma_a.g, sigma_s_prime.g, eta),
dipoleReflectance(r, sigma_a.b, sigma_s_prime.b, eta)
);
}Applying the Dipole Model in Shaders
In real-time rendering, running a brute-force spatial integration over neighboring pixels for every fragment is computationally expensive. Dipole diffusion is typically integrated in one of two ways:
1. Texture-Space or Screen-Space Importance Sampling
Sample adjacent points around the current surface position using a disk sampling pattern, weighting the incoming irradiance \(E(x_i)\) by the dipole profile value \(R(\|x - x_i\|)\) and the sample area:
uniform sampler2D irradianceTexture;
uniform sampler2D worldPosTexture;
uniform vec3 uSigmaA;
uniform vec3 uSigmaSPrime;
uniform float uEta;
in vec2 vTexCoord;
out vec4 fragColor;
const int NUM_SAMPLES = 16;
const vec2 POISSON_DISK[16] = vec2[](
vec2(-0.326, -0.406), vec2(-0.840, -0.074),
vec2(-0.696, 0.457), vec2(-0.203, 0.621),
vec2( 0.962, -0.195), vec2( 0.473, -0.480),
vec2( 0.519, 0.767), vec2( 0.185, -0.893),
vec2( 0.507, 0.144), vec2( 0.896, 0.412),
vec2(-0.322, -0.933), vec2(-0.792, -0.598),
vec2(-0.165, 0.039), vec2( 0.239, 0.379),
vec2( 0.284, -0.187), vec2(-0.028, 0.923)
);
void main() {
vec3 currentPos = texture(worldPosTexture, vTexCoord).xyz;
vec3 accumulatedLight = vec3(0.0);
vec3 totalWeight = vec3(0.0);
float diskRadius = 0.05; // World-space sampling radius
for (int i = 0; i < NUM_SAMPLES; ++i) {
vec2 sampleUV = vTexCoord + POISSON_DISK[i] * 0.01; // Screen/Texture UV offset
vec3 samplePos = texture(worldPosTexture, sampleUV).xyz;
vec3 sampleIrradiance = texture(irradianceTexture, sampleUV).rgb;
float r = distance(currentPos, samplePos);
vec3 weight = evaluateDipoleRGB(r, uSigmaA, uSigmaSPrime, uEta);
accumulatedLight += sampleIrradiance * weight;
totalWeight += weight;
}
vec3 finalRadiance = accumulatedLight / max(totalWeight, vec3(0.0001));
fragColor = vec4(finalRadiance, 1.0);
}2. Precomputed Kernels for Separable SSS
For production pipelines requiring high performance, the non-linear dipole profile \(R(r)\) is frequently pre-evaluated into a 1D lookup table or approximated as a sum of Gaussians. This permits execution via two separable blur passes (horizontal and vertical) across the screen buffer, depth-weighted to prevent bleeding across geometric discontinuities.