Light Pollution Ground Glow in GLSL Atmosphere Scatter

Integrating light pollution ground glow into nocturnal atmospheric scattering shaders requires expanding the standard single-scattering raymarching model to account for upward terrestrial radiance alongside celestial sources like the moon. In GLSL, this is achieved by modeling artificial city lighting as upward hemispherical emitters, sampling their irradiance along the primary view ray, and computing the accumulated in-scattered radiance attenuated by Rayleigh and Mie extinction. This article breaks down the physical modeling of ground-based light domes, the modification of the raymarching integral, and the efficient implementation of this phenomenon in real-time GLSL rendering pipelines.

The Physical Model of Urban Ground Glow

Standard daylight atmospheric models (such as Bruneton-Neyret or Nishita) assume a single directional light source located at infinity. At night, while the moon acts as a secondary directional source, urban centers act as distributed, upward-facing hemispherical emitters characterized by Garstang’s light pollution model:

\[\vec{L}_{\text{ground}}(x, z) = I_0 \cdot \rho(x, z) \cdot F(\theta)\]

Where:

  • \(I_0\) is the base luminance factor.
  • \(\rho(x, z)\) represents a 2D radiance distribution texture of urban footprints (often derived from nighttime satellite imagery).
  • \(F(\theta)\) is an angular emission function describing light emitted at zenith angle \(\theta\) (frequently modeled as a mix of Lambertian and low-angle direct upward beam emission).

In the atmosphere, this upward flux interacts predominantly with Mie particles (aerosols and pollutants located close to the planetary boundary layer, typically \(0\text{--}2\text{ km}\)) rather than pure Rayleigh molecules. This creates the localized, haze-heavy visual dome observed over cities.

Modifying the In-Scattering Integral

Along a camera ray with origin \(\mathbf{o}\) and direction \(\mathbf{d}\), the total in-scattered nocturnal radiance \(L(\mathbf{o}, \mathbf{d})\) across distance \(S\) is expressed as:

\[L(\mathbf{o}, \mathbf{d}) = \int_{0}^{S} T(\mathbf{o}, \mathbf{x}) \cdot \left[ L_{\text{moon}}(\mathbf{x}) + L_{\text{ground}}(\mathbf{x}) \right] ds\]

Where:

  • \(\mathbf{x} = \mathbf{o} + s\mathbf{d}\) is the sample point along the ray.
  • \(T(\mathbf{o}, \mathbf{x}) = \exp\left(-\int_{0}^{s} (\beta_R(s') + \beta_M(s')) ds'\right)\) is the optical transmittance.
  • \(L_{\text{ground}}(\mathbf{x})\) is the in-scattering contribution from all surface light sources to point \(\mathbf{x}\).

Because evaluating an area integral over the entire ground footprint at every raymarch step is prohibitive for real-time applications, \(L_{\text{ground}}(\mathbf{x})\) is approximated by sampling an array of dominant point-like light centers or by querying a downsampled, mipmapped 2D ground irradiance map.

GLSL Implementation

Below is a GLSL fragment implementation demonstrating the evaluation of nocturnal scatter with an integrated city light dome.

#version 450 core

out vec4 FragColor;
in vec2 uv;

uniform vec3 u_CameraPos;
uniform vec3 u_RayDir;
uniform vec3 u_MoonDir;
uniform vec3 u_MoonColor;

// Light pollution parameters
uniform vec3 u_CityPosition;      // World-space ground position (e.g., vec3(x, 0.0, z))
uniform vec3 u_CityRadiance;      // Color and intensity of light emission (e.g., sodium vapor orange)
uniform float u_CityRadius;       // Effective radius of the urban footprint

const int STEPS = 32;
const float EARTH_RADIUS = 6371000.0;
const float ATMO_RADIUS  = 6471000.0;
const float HR = 8000.0; // Rayleigh scale height
const float HM = 1200.0; // Mie scale height (aerosol trap)

const vec3 BETA_R = vec3(5.8e-6, 13.5e-6, 33.1e-6); // Rayleigh scattering coefficients
const vec3 BETA_M = vec3(21.0e-6);                   // Mie scattering coefficients
const vec3 BETA_EXT_M = BETA_M * 1.11;

// Cornette-Shanks or Henyey-Greenstein phase function for Mie
float PhaseMie(float cosTheta, float g) {
    float g2 = g * g;
    return (3.0 * (1.0 - g2) * (1.0 + cosTheta * cosTheta)) / 
           (2.0 * (2.0 + g2) * pow(1.0 + g2 - 2.0 * g * cosTheta, 1.5)) / (4.0 * 3.14159265);
}

float PhaseRayleigh(float cosTheta) {
    return 3.0 / (16.0 * 3.14159265) * (1.0 + cosTheta * cosTheta);
}

// Intersect ray with atmosphere shell
bool RaySphereIntersect(vec3 ro, vec3 rd, float rad, out float t0, out float t1) {
    float b = dot(ro, rd);
    float c = dot(ro, ro) - rad * rad;
    float d = b * b - c;
    if (d < 0.0) return false;
    t0 = -b - sqrt(d);
    t1 = -b + sqrt(d);
    return true;
}

void main() {
    float tNear, tFar;
    if (!RaySphereIntersect(u_CameraPos, u_RayDir, ATMO_RADIUS, tNear, tFar)) {
        FragColor = vec4(0.0);
        return;
    }
    tNear = max(tNear, 0.0);

    float stepSize = (tFar - tNear) / float(STEPS);
    vec3 rayleighAccum = vec3(0.0);
    vec3 mieAccum = vec3(0.0);
    vec3 opticalDepth = vec3(0.0);

    float cosMoon = dot(u_RayDir, u_MoonDir);
    float phaseR_Moon = PhaseRayleigh(cosMoon);
    float phaseM_Moon = PhaseMie(cosMoon, 0.76);

    for (int i = 0; i < STEPS; ++i) {
        float s = tNear + (float(i) + 0.5) * stepSize;
        vec3 p = u_CameraPos + u_RayDir * s;
        float h = length(p) - EARTH_RADIUS;

        if (h < 0.0) break; // Occluded by terrain/planet

        float densityR = exp(-h / HR) * stepSize;
        float densityM = exp(-h / HM) * stepSize;

        // Cumulative extinction from camera to current sample point
        vec3 stepExtinction = BETA_R * densityR + BETA_EXT_M * densityM;
        opticalDepth += stepExtinction;
        vec3 tau = exp(-opticalDepth);

        // 1. Lunar directional scattering
        rayleighAccum += densityR * tau;
        mieAccum      += densityM * tau;

        // 2. Terrestrial Light Pollution Scattering
        vec3 toCity = u_CityPosition - p;
        float distToCitySq = dot(toCity, toCity);
        float distToCity = sqrt(distToCitySq);
        vec3 cityLightDir = toCity / distToCity;

        // Phase angle between view ray and incoming ground light vector
        float cosGround = dot(u_RayDir, -cityLightDir);
        float phaseM_Ground = PhaseMie(cosGround, 0.65);

        // Ground-to-sample point extinction approximation (Beer-Lambert over distToCity)
        float cityTransmittance = exp(-BETA_EXT_M.x * densityM * (distToCity / stepSize));
        
        // Geometric attenuation and cosine distribution for upward emission
        float geometricDropoff = 1.0 / max(distToCitySq, 1.0);
        float upwardFactor = clamp(dot(normalize(u_CityPosition), -cityLightDir), 0.0, 1.0);

        vec3 groundIrradiance = u_CityRadiance * upwardFactor * geometricDropoff * cityTransmittance;

        // Light pollution heavily affects the Mie phase (aerosol haze)
        mieAccum += densityM * groundIrradiance * phaseM_Ground * tau;
    }

    vec3 moonRadiance = u_MoonColor * (rayleighAccum * BETA_R * phaseR_Moon + mieAccum * BETA_M * phaseM_Moon);
    vec3 finalScatter = moonRadiance + (mieAccum * BETA_M);

    FragColor = vec4(finalScatter, 1.0);
}

Performance Optimization and LUT Integration

Executing multi-source ground sampling within a high-step raymarch loop introduces substantial overhead. Production rendering pipelines typically decouple and accelerate the pass using three primary strategies:

  1. Dual-Resolution Passes: The celestial nocturnal pass (moonlight) is computed at full resolution or read from precomputed Look-Up Tables (LUTs), while the ground glow Mie scatter is raymarched at quarter resolution (\(0.25\times\)) with \(8\text{--}16\) steps and bilaterally upsampled using depth-aware weights.
  2. Ground Luminance Shadow Maps: For complex terrain, occlusions of the city glow behind mountain ranges are resolved by projecting a hemispherical shadow map centered on the city centroid, preventing unnatural glow on opposite mountain faces.
  3. Aerosol Stratification Separation: Because light pollution primarily scatters below \(2\text{ km}\), the integration distance along the ray can be clamped to the planetary boundary layer rather than marching through the full atmosphere thickness.