Cloud Optical Thickness Integration in GLSL
This article explores how to calculate and integrate cloud optical thickness through procedural 3D noise fields using GLSL raymarching. By numerically integrating density values sampled from procedural 3D noise (such as Perlin-Worley noise) along a ray path, real-time shaders compute the optical depth, evaluate transmittance via the Beer-Lambert law, and determine both self-shadowing and light attenuation through participating media.
The Physics of Optical Thickness
Optical thickness (or optical depth, denoted as \(\tau\)) measures how opaque a medium is to light passing through it. Along a ray of length \(d\) parameterized by \(s\), optical thickness is defined by integrating the extinction coefficient \(\sigma_e(s)\) along the path:
\[\tau = \int_{0}^{d} \sigma_e(s) \, ds\]
In homogeneous media, \(\sigma_e\) is constant, but procedural clouds represent heterogeneous media where \(\sigma_e(s)\) varies continuously based on local density \(\rho(s)\). The relationship is typically scaled by an extinction factor: \(\sigma_e(s) = \kappa \cdot \rho(s)\).
Transmittance (\(T\)), representing the fraction of light that penetrates the medium without being scattered or absorbed, is computed using the Beer-Lambert law:
\[T = \exp(-\tau) = \exp\left(-\int_{0}^{d} \kappa \cdot \rho(s) \, ds\right)\]
Sampling the Procedural 3D Noise Field
Clouds derive their shape from procedural noise fields sampled within a defined volume (such as a bounding box or spherical atmosphere layer). Common setups combine low-frequency noise (to define broad cloud billows, often a blend of Perlin and Worley noise) with high-frequency Worley noise (to erode edges and simulate wispy condensation).
In GLSL, density is sampled at a world position p:
float sampleDensity(vec3 p) {
// Apply wind animation and scale
vec3 uvw = p * 0.001 + vec3(0.0, 0.0, u_time * 0.02);
// Low-frequency base shape
float baseNoise = texture(u_perlinWorleyNoise, uvw).r;
// High-frequency detail erosion
float detailNoise = texture(u_worleyNoise, uvw * 4.0).r;
float density = baseNoise - (1.0 - detailNoise) * 0.25;
// Apply altitude attenuation mask and contrast
float heightFraction = clamp((p.y - u_cloudBottom) / (u_cloudTop - u_cloudBottom), 0.0, 1.0);
density *= smoothstep(0.0, 0.2, heightFraction) * smoothstep(1.0, 0.7, heightFraction);
return max(0.0, density - u_coverageThreshold) * u_densityMultiplier;
}Numerical Integration via Raymarching
Because analytical integration over procedural noise is impossible, shaders evaluate the integral via Riemann summation using a uniform step size \(\Delta s\):
\[\tau \approx \sum_{i=0}^{N} \kappa \cdot \rho(\mathbf{x}_i) \cdot \Delta s\]
Where \(\mathbf{x}_i = \mathbf{x}_0 + i \cdot \Delta s \cdot \vec{\omega}\), with \(\vec{\omega}\) being the normalized view ray direction.
Integrating Primary and Shadow Rays
Rendering realistic clouds requires two levels of integration:
- Primary Ray Integration: Accumulates optical depth along the view ray to calculate how much background light reaches the camera.
- Light Ray Integration (In-Scattering): At each step along the primary ray, a secondary ray marches toward the sun/light source to integrate the optical thickness between that sample point and the light. This determines how much direct sunlight illuminates that segment of the cloud.
// Integrate optical thickness toward the light source
float marchLight(vec3 samplePos, vec3 lightDir, float stepSize, int steps) {
float lightOpticalDepth = 0.0;
vec3 currentPos = samplePos;
for (int i = 0; i < steps; ++i) {
currentPos += lightDir * stepSize;
float density = sampleDensity(currentPos);
lightOpticalDepth += density * stepSize;
}
return lightOpticalDepth;
}
// Main integration along the primary view ray
vec4 integrateCloudOpticalThickness(vec3 rayOrigin, vec3 rayDir, float maxDist) {
float primaryStepSize = 25.0;
float lightStepSize = 50.0;
int primarySteps = 64;
int lightSteps = 6;
float totalOpticalDepth = 0.0;
float transmittance = 1.0;
vec3 accumulatedColor = vec3(0.0);
vec3 currentPos = rayOrigin;
for (int i = 0; i < primarySteps; ++i) {
float density = sampleDensity(currentPos);
if (density > 0.001) {
// Integrate step optical thickness
float stepOpticalDepth = density * primaryStepSize * u_extinctionCoefficient;
totalOpticalDepth += stepOpticalDepth;
// Optical thickness to light
float lightOpticalDepth = marchLight(currentPos, u_lightDir, lightStepSize, lightSteps) * u_extinctionCoefficient;
float lightTransmittance = exp(-lightOpticalDepth);
// Phase function approximation (e.g., Henyey-Greenstein)
float phase = 0.25; // Constant placeholder
vec3 inScatteredLight = u_sunColor * lightTransmittance * phase;
// Local transmittance for this step
float stepTransmittance = exp(-stepOpticalDepth);
// Accumulate visible energy
accumulatedColor += transmittance * (1.0 - stepTransmittance) * inScatteredLight;
transmittance *= stepTransmittance;
// Early exit if the medium becomes fully opaque
if (transmittance < 0.01) {
break;
}
}
currentPos += rayDir * primaryStepSize;
}
return vec4(accumulatedColor, 1.0 - transmittance);
}Performance Optimization Strategies
Direct raymarching of 3D noise fields is computationally heavy. Production GLSL implementations rely on three core optimizations:
- Jittering and Temporal Reprojection: Ray origins are jittered using blue noise or interleaved gradients along the ray direction to turn step-aliasing into high-frequency noise, which is then smoothed over multiple frames using temporal anti-aliasing (TAA).
- Coarse-to-Fine Stepping: Shaders use larger step sizes while density is zero. Once a non-zero density value is encountered, the step size decreases to properly resolve edges, reverting to larger steps upon exiting dense regions.
- Low-Resolution Offscreen Buffers: Because volumetric clouds contain predominantly low-frequency lighting variations, the primary integration loop is often rendered at half or quarter resolution and upscaled using bilateral depth-aware filtering to prevent edge bleeding against high-resolution geometry.