Integrating Optical Depth in GLSL Atmospheres

Simulating realistic planetary atmospheres in real-time requires calculating how light attenuates as it travels through varying gas and aerosol densities. In GLSL, this attenuation is governed by the Beer-Lambert law and computed by numerically integrating the optical depth along a ray passing through concentric spherical shells. This article explains the geometric foundations, density functions for multiple atmospheric layers, numerical integration via raymarching, and the GLSL implementation required to integrate optical depth accurately on the GPU.

Mathematical Foundation of Optical Depth

Optical depth (\(\tau\)) measures how much radiant power is absorbed or scattered as a ray traverses a participating medium between two points, \(\mathbf{A}\) and \(\mathbf{B}\):

\[\tau(\mathbf{A}, \mathbf{B}) = \int_{\mathbf{A}}^{\mathbf{B}} \beta_{ext}(\mathbf{x}) \, ds\]

Where \(\beta_{ext}(\mathbf{x})\) is the extinction coefficient at point \(\mathbf{x}\), which is proportional to the local atmospheric density \(\rho(h)\) at altitude \(h\):

\[\beta_{ext}(\mathbf{x}) = \beta_0 \cdot \rho(h)\]

Because real atmospheres are spherical and consist of multiple components (such as Rayleigh scattering molecules, Mie scattering aerosols, and Ozone absorption), the total optical depth is the sum of the individual layer densities integrated over the ray segment:

\[\tau_{total} = \sum_{k} \beta_{k} \int_{\mathbf{A}}^{\mathbf{B}} \rho_k(h(s)) \, ds\]

Spherical Shell Geometry and Ray-Sphere Intersection

A planetary atmosphere is geometrically bounded by the planet's radius \(R_{p}\) and the top-of-atmosphere (TOA) radius \(R_{atm}\). The altitude \(h\) of any sample point \(\mathbf{x}\) relative to the planet center \(\mathbf{C}\) is:

\[h = \|\mathbf{x} - \mathbf{C}\| - R_{p}\]

Before integrating, the ray \(\mathbf{r}(t) = \mathbf{o} + t\mathbf{d}\) must be clipped against the outer sphere (\(R_{atm}\)) and checked for occlusion against the planetary surface (\(R_{p}\)).

// Intersect ray with a sphere at origin; returns vec2(t_near, t_far)
vec2 raySphereIntersect(vec3 ro, vec3 rd, float radius) {
    float b = dot(ro, rd);
    float c = dot(ro, ro) - radius * radius;
    float d = b * b - c;
    if (d < 0.0) return vec2(-1.0);
    float sqrtD = sqrt(d);
    return vec2(-b - sqrtD, -b + sqrtD);
}

Multi-Layer Atmospheric Density Models

Different atmospheric components drop off at different rates, represented by their scale heights (\(H\)):

  1. Rayleigh Layer (Air molecules): Follows barometric exponential decay. \[\rho_R(h) = \exp(-h / H_R) \quad (H_R \approx 8.0 \text{ km})\]
  2. Mie Layer (Aerosols/Haze): Also follows exponential decay but with a lower scale height. \[\rho_M(h) = \exp(-h / H_M) \quad (H_M \approx 1.2 \text{ km})\]
  3. Ozone Absorption Layer: Typically modeled as a tent function peaking at an altitude around 25 km with a finite distribution width. \[\rho_O(h) = \max\left(0.0, 1.0 - \frac{|h - h_{peak}|}{w}\right)\]
struct DensitySample {
    float rayleigh;
    float mie;
    float ozone;
};

DensitySample sampleDensities(float altitude) {
    const float H_RAYLEIGH = 8000.0;
    const float H_MIE      = 1200.0;
    const float OZONE_PEAK = 25000.0;
    const float OZONE_WIDTH = 15000.0;

    DensitySample s;
    s.rayleigh = exp(-altitude / H_RAYLEIGH);
    s.mie      = exp(-altitude / H_MIE);
    s.ozone    = max(0.0, 1.0 - abs(altitude - OZONE_PEAK) / OZONE_WIDTH);
    return s;
}

Numerical Integration via Raymarching

Because \(h(s)\) varies non-linearly with distance \(s\) along arbitrary spherical paths, the integral is evaluated numerically using mid-point Riemann sums. The path from ray entry \(t_0\) to ray exit \(t_1\) is divided into \(N\) segments of step size \(ds = (t_1 - t_0) / N\).

For each step \(i\):

  1. Calculate the mid-point parameter: \(t = t_0 + (i + 0.5) \cdot ds\)
  2. Compute the position: \(\mathbf{x} = \mathbf{o} + t\mathbf{d}\)
  3. Evaluate altitude: \(h = \|\mathbf{x} - \mathbf{C}\| - R_{p}\)
  4. Accumulate densities multiplied by segment length \(ds\).

GLSL Implementation

The following function evaluates the integrated optical depth along a segment through the atmosphere:

// Returns integrated optical depth for Rayleigh, Mie, and Ozone layers
vec3 integrateOpticalDepth(vec3 rayOrigin, vec3 rayDir, float t0, float t1, vec3 planetCenter, float planetRadius, int sampleCount) {
    float stepSize = (t1 - t0) / float(sampleCount);
    vec3 opticalDepthAccum = vec3(0.0); // x: Rayleigh, y: Mie, z: Ozone

    for (int i = 0; i < sampleCount; ++i) {
        float t = t0 + (float(i) + 0.5) * stepSize;
        vec3 samplePos = rayOrigin + rayDir * t;
        float altitude = length(samplePos - planetCenter) - planetRadius;

        // Skip samples that dip below terrain/sea level
        if (altitude < 0.0) continue;

        DensitySample density = sampleDensities(altitude);

        opticalDepthAccum.x += density.rayleigh * stepSize;
        opticalDepthAccum.y += density.mie * stepSize;
        opticalDepthAccum.z += density.ozone * stepSize;
    }

    return opticalDepthAccum;
}

Transmittance Calculation

Once the optical depth vector \(\boldsymbol{\tau} = (\tau_R, \tau_M, \tau_O)\) is integrated, the total transmittance \(T\) across the path is computed by applying the extinction cross-sections:

vec3 computeTransmittance(vec3 opticalDepth) {
    const vec3 BETA_RAYLEIGH = vec3(5.802, 13.558, 33.100) * 1e-6; // m^-1
    const vec3 BETA_MIE      = vec3(3.996, 3.996, 3.996) * 1e-6;   // m^-1
    const vec3 BETA_OZONE    = vec3(0.650, 1.881, 0.085) * 1e-6;   // m^-1

    vec3 extinction = BETA_RAYLEIGH * opticalDepth.x 
                    + BETA_MIE      * opticalDepth.y 
                    + BETA_OZONE    * opticalDepth.z;

    return exp(-extinction);
}

This transmittance term is evaluated along both the primary view ray and secondary rays cast toward the directional light source at each step during full atmospheric scattering passes, accurately preserving the spherical curvature of the planet.