Raymarching Volumetric Fog in GLSL Compute Shaders

This article explores the technical implementation of real-time volumetric fog using raymarching step accumulation within a OpenGL or Vulkan GLSL compute shader. It breaks down the process of casting rays through participating media, sampling density distributions, evaluating direct lighting and shadow maps, applying phase functions, and numerically integrating transmittance and in-scattering using the Beer-Lambert law to produce cohesive atmospheric lighting.

Mathematical Foundations: The Radiative Transfer Equation

Volumetric fog simulates how light interacts with microscopic particles suspended in the air. The core behavior is governed by the simplified Radiative Transfer Equation (RTE), which balances light extinction (out-scattering and absorption) and in-scattering along a viewing ray.

According to the Beer-Lambert law, the fraction of light that survives traveling a distance \(d\) through a medium with extinction coefficient \(\sigma_t\) is the transmittance \(T\):

\[T = e^{-\int_0^d \sigma_t(x) \, dx}\]

When discretized for raymarching, light accumulation integrates scattered light at discrete samples along the ray while compounding the transmittance decay.

Setting Up the Compute Shader

A compute shader operates over an invocation grid aligned with the render target, processing individual pixels or a low-resolution frustum-aligned voxel grid (often called a "froxel" grid).

#version 450 core
layout(local_size_x = 16, local_size_y = 16, local_size_z = 1) in;

layout(binding = 0, rgba16f) uniform writeonly image2D outVolumetricFog;
layout(binding = 1) uniform sampler2D sceneDepthTex;
layout(binding = 2) uniform sampler2D shadowMap;
layout(binding = 3) uniform sampler3D densityNoiseTex;

uniform mat4 invViewProj;
uniform vec3 cameraPos;
uniform vec3 lightDir;
uniform vec3 lightColor;

Ray Construction and Depth Clamping

For each invocation, the compute shader reconstructs a world-space ray originating at the camera position and traveling toward the scene geometry.

Calculate Screen UVs: Convert gl_GlobalInvocationID.xy into normalized device coordinates (NDC) in the range [-1, 1].
Reconstruct World Ray: Unproject the NDC coordinates using the inverse view-projection matrix to determine the normalized ray direction \(\vec{D}\).
Determine Ray Bounds: Sample the scene depth buffer at the current UV. Linearize this depth to establish the maximum ray distance \(t_{\text{max}}\). Fog integration occurs exclusively between the near plane \(t_{\text{min}}\) and \(t_{\text{max}}\) to prevent fog from rendering behind opaque geometry.

Jittering and Step Discretization

Uniform sampling along ray steps causes conspicuous slice banding artifacts. To mitigate this, developers apply a pseudo-random or blue-noise offset to the starting point of the ray:

float jitter = texture(blueNoiseTex, gl_GlobalInvocationID.xy / 64.0).r;
float stepLength = (tMax - tMin) / float(NUM_STEPS);
float currentT = tMin + stepLength * jitter;

Jittering converts structured banding into high-frequency noise, which spatial or temporal denoising filters can subsequently resolve.

The Step Accumulation Loop

The raymarcher advances through the volume in fixed or exponential increments. At each step \(i\), it computes the current world position \(\vec{P} = \vec{O} + t_i \vec{D}\).

1. Density Evaluation

Density represents the concentration of particles at \(\vec{P}\). A typical density function combines an analytical height falloff with a 3D noise texture:

float getDensity(vec3 p) {
    float heightFalloff = exp(-p.y * 0.1);
    float noise = texture(densityNoiseTex, p * 0.05).r;
    return max(0.0, heightFalloff * noise);
}

2. Phase Function Evaluation

The phase function describes how strongly light scatters toward the viewer relative to the incoming light angle \(\theta\). The Henyey-Greenstein phase function models this anisotropic scattering via the asymmetry parameter \(g \in (-1, 1)\), where positive values simulate forward scattering:

float henyeyGreenstein(float cosTheta, float g) {
    float g2 = g * g;
    return (1.0 - g2) / (4.0 * 3.14159265 * pow(1.0 + g2 - 2.0 * g * cosTheta, 1.5));
}

3. Shadow and Lighting Sampling

The position \(\vec{P}\) is projected into the light's shadow space to assess visibility. If the point resides in shadow, direct in-scattering is omitted, though ambient terms may still apply.

4. Transmittance and In-Scattering Accumulation

The accumulated color and remaining transmittance are updated sequentially:

vec3 accumulatedColor = vec3(0.0);
float accumulatedTransmittance = 1.0;
float sigmaT = 0.1; // Extinction coefficient

for(int i = 0; i < NUM_STEPS; ++i) {
    vec3 currentPos = cameraPos + rayDir * currentT;
    float density = getDensity(currentPos);

    if (density > 0.001) {
        float shadow = sampleShadow(currentPos);
        float phase = henyeyGreenstein(dot(rayDir, -lightDir), 0.3);
        
        // Attenuation for this single step
        float stepOpticalDepth = density * sigmaT * stepLength;
        float stepTransmittance = exp(-stepOpticalDepth);

        // Light scattered toward the camera at this step
        vec3 inScatteredLight = lightColor * shadow * phase * density;

        // Integrated scattering over the step length
        vec3 integratedScattering = (inScatteredLight - inScatteredLight * stepTransmittance) / sigmaT;
        accumulatedColor += accumulatedTransmittance * integratedScattering;

        // Attenuate total ray transmittance
        accumulatedTransmittance *= stepTransmittance;

        if (accumulatedTransmittance < 0.01) {
            break; // Early ray termination
        }
    }

    currentT += stepLength;
}

Final Compositing

Once the loop completes or early exit triggers, the compute shader writes the results to an image texture:

imageStore(outVolumetricFog, ivec2(gl_GlobalInvocationID.xy), vec4(accumulatedColor, 1.0 - accumulatedTransmittance));

The resulting texture contains both the in-scattered light and the fog coverage factor (alpha), which a subsequent raster pass blends over the final rendered frame using standard additive and multiplicative composition:

\[\text{FinalColor} = \text{SceneColor} \times T_{\text{accumulated}} + \text{InScatteredColor}\]