Atmospheric Scattering Shader

Just found this rather nice Atmospheric Scattering rendering C++ code over at scratchapixel.com and thought I’d do a quick conversion to a GLSL shader as a test for the Timeline software I’m working on. Works rather nicely…

My (none optimised) fragment shader conversion is:

#version 150

#define M_PI 3.1415926535897932384626433832795

uniform float TimeOfDay; // range 0.0 -> 1.0 (0.0 = Midnight, 0.5 = Midday, etc)

const float RADIUS_EARTH = 6360e3;
const float RADIUS_ATMOSPHERE = 6420e3;
const float RAYLEIGH_SCALE_HEIGHT = 7994;
const float MIE_SCALE_HEIGHT = 1200;
const float SUN_INTENSITY = 20;

const float g = 0.76;

const vec3 betaR = vec3( 5.5e-6, 13.0e-6, 22.4e-6 );    // Rayleigh scattering coefficients at sea level
const vec3 betaM = vec3( 21e-6 );                       // Mie scattering coefficients at sea level

vec3 sunDirection = vec3( 0, 1, 0 );

const int numSamples = 16;
const int numSamplesLight = 8;

struct Ray
{
    vec3 o; //origin
    vec3 d; //direction (should always be normalized)
};

struct Sphere
{
    vec3 pos;   //center of sphere position
    float rad;  //radius
};

const Sphere SPHERE_EARTH      = Sphere( vec3( 0 ), RADIUS_EARTH );
const Sphere SPHERE_ATMOSPHERE = Sphere( vec3( 0 ), RADIUS_ATMOSPHERE );

bool intersect( in Ray ray, in Sphere sphere, out float t0, out float t1 )
{
    vec3 oc = ray.o - sphere.pos;
    float b = 2.0 * dot(ray.d, oc);
    float c = dot(oc, oc) - sphere.rad*sphere.rad;
    float disc = b * b - 4.0 * c;

    if (disc < 0.0)
        return false;

   float q;
    if (b < 0.0)         q = (-b - sqrt(disc))/2.0;     else         q = (-b + sqrt(disc))/2.0;       t0 = q;     t1 = c / q;     // make sure t0 is smaller than t1     if (t0 > t1) {
        // if t0 is bigger than t1 swap them around
        float temp = t0;
        t0 = t1;
        t1 = temp;
    }

    // if t1 is less than zero, the object is in the ray's negative direction
    // and consequently the ray misses the sphere
    if (t1 < 0.0)
        return false;

    if( t0 < 0.0 )
    {
        t0 = 0;
    }

    return( true );
}

vec3 computeIncidentLight( in Ray r )
{
    float       t0, t1;

    if( !intersect( r, SPHERE_ATMOSPHERE, t0, t1 ) )
    {
        return vec3( 1 );
    }

    float segmentLength = ( t1 - t0 ) / numSamples;
    float tCurrent = t0;

    vec3 sumR = vec3( 0 );
    vec3 sumM = vec3( 0 );

    float opticalDepthR = 0;
    float opticalDepthM = 0;

    float mu = dot( r.d, sunDirection );
    float phaseR = 3 / ( 16 * M_PI ) * ( 1 + mu * mu );
    float phaseM = 3 / (  8 * M_PI ) * ( ( 1 - g * g ) * ( 1 + mu * mu ) ) / ( ( 2 + g * g ) * pow( 1 + g * g - 2 * g * mu, 1.5 ) );

    for( int i = 0; i < numSamples ; i++ )
    {
        vec3    samplePosition = r.o + r.d * ( tCurrent + 0.5 * segmentLength );
        float   height = length( samplePosition ) - RADIUS_EARTH;

        // compute optical depth for light

        float hr = exp( -height / RAYLEIGH_SCALE_HEIGHT ) * segmentLength;
        float hm = exp( -height / MIE_SCALE_HEIGHT      ) * segmentLength;

        opticalDepthR += hr;
        opticalDepthM += hm;

        // light optical depth

        Ray lightRay = Ray( samplePosition, sunDirection );

        float lmin, lmax;

        intersect( lightRay, SPHERE_ATMOSPHERE, lmin, lmax );

        float segmentLengthLight = lmax / numSamplesLight;
        float tCurrentLight = 0;
        float opticalDepthLightR = 0;
        float opticalDepthLightM = 0;
        
        int j = 0;

        for( ; j < numSamplesLight ; j++ )
        {
            vec3 samplePositionLight = lightRay.o + lightRay.d * ( tCurrentLight + 0.5 * segmentLengthLight );

            float heightLight = length( samplePositionLight ) - RADIUS_EARTH;

            if( heightLight < 0 )
            {
                break;
            }

            opticalDepthLightR += exp( -heightLight / RAYLEIGH_SCALE_HEIGHT ) * segmentLengthLight;
            opticalDepthLightM += exp( -heightLight / MIE_SCALE_HEIGHT      ) * segmentLengthLight;

            tCurrentLight += segmentLengthLight;
        }

        if( j == numSamplesLight )
        {
            vec3 tau = betaR * ( opticalDepthR + opticalDepthLightR ) + betaM * 1.1 * ( opticalDepthM + opticalDepthLightM );
            vec3 attenuation = exp( -tau );

            sumR += hr * attenuation;
            sumM += hm * attenuation;
        }

        tCurrent += segmentLength;
    }

    return( SUN_INTENSITY * ( sumR * phaseR * betaR + sumM * phaseM * betaM ) );
}

void main()
{
    const int width = 512;
    const int height = 512;

    float a = mod( TimeOfDay - 0.5, 1 ) * 2.0 * M_PI;

    sunDirection = normalize( vec3( 0, cos( a ), sin( a ) ) );

    float x = 2 * ( gl_FragCoord.x + 0.5 ) / ( width  - 1 ) - 1;
    float y = 2 * ( gl_FragCoord.y + 0.5 ) / ( height - 1 ) - 1;

    float z2 = x * x + y * y; 

    if( z2 <= 1 )
    {
        float phi   = atan( y, x );
        float theta = acos( 1 - z2 );

        vec3 dir = vec3( sin( theta ) * cos( phi ), cos( theta ), sin( theta ) * sin( phi ) );
        vec3 pos = vec3( 0, RADIUS_EARTH + 1, 0 );

        gl_FragColor = vec4( computeIncidentLight( Ray( pos, normalize( dir ) ) ), 1 );
    }
    else
    {
        gl_FragColor = vec4( 0 );
    }
}