This is the second of three technical articles detailing my recent work on improving the look and customizability of the skies in Cloud Party. Part 1 covered background and research on sky dome rendering; part 2 (this one) will go into detail on the implementation and optimization of the sky dome rendering technique; and part 3 will cover the details of the cloud system.

Recall from last time that our goal is to create a sky renderer that looks good (not necessarily correct, but believable), is customizable by users, reacts to changes in sun angle, and runs very fast. Speed is important because we want to run and look good on a wide variety of GPUs, and ideally the sky rendering should not take up very much of our frame time.

If you missed it last time, go check out the demo!

The code in this article is licensed under a Creative Commons Attribution 3.0 Unported License.

Initial Implementation

The first simplification to the scattering problem we are going to make is to only compute single-scattering (assuming light bounces of particles in the atmosphere at most one time) instead of multiple-scattering (where light bounces around many times). We are also going to assume that the density of the atmosphere is constant.

With these simplifications, using the ATI paper I referenced last time I wrote the following (completely unoptimized) GLSL functions:

vec3 beta_rayleigh, beta_mie;
float g_mie;
float PI = 3.14159265359;

vec3 calcExtinction(float dist) {
  return exp(-dist * (beta_rayleigh + beta_mie));
}

vec3 calcScattering(float cos_theta) {
  // phase for rayleigh scattering
  float phase_rayleigh = (3.0 / (16.0 * PI)) * (1.0 + cos_theta * cos_theta);

  // phase for mie scattering, using Henyey/Greenstein approximation
  float phase_mie = (1.0 / (4.0 * PI)) * (1.0 - g_mie) * (1.0 - g_mie);
  phase_mie /= pow(1.0 + g_mie * g_mie - 2.0 * g_mie * cos_theta, 3.0 / 2.0);

  return (phase_rayleigh * beta_rayleigh + phase_mie * beta_mie) / (beta_rayleigh +
beta_mie);
}

Additionally we need to calculate the distance through the atmosphere for a particular position and vector. To do this I use a ray-sphere intersection test, which I simplified by assuming the position is located inside the atmosphere:

float PLANET_RADIUS = 6371000;
float ATMO_RADIUS = PLANET_RADIUS + 100000;

float opticalDepth(in vec3 position, in vec3 ray) {
  // the incoming position is in a space with the origin on the surface of the planet
  // convert to a space with the origin at the center of the planet
  position.z += PLANET_RADIUS;

  // ray-sphere intersection, assuming position is inside the sphere
  float a0 = ATMO_RADIUS * ATMO_RADIUS - dot(position, position);
  float a1 = dot(position, ray);
  return sqrt(a1 * a1 + a0) - a1;
}

Now all that is left is to use these functions to calculate single-scattering for the sun. The standard method for doing this is to step along the view ray and then cast a ray towards the sun, calculating the amount of light that makes it through the atmosphere to that point. Then we compute the amount of in-scattered light from each of these rays and accumulate:

int NUM_STEPS = 20;

vec3 current_position = vec3(0.0);
float ray_dist = opticalDepth(current_position, view_vector);
float step_size = ray_dist / float(NUM_STEPS);

// calculate extinction for a single step
vec3 step_extinction = calcExtinction(step_size);

// calculate the in-scattering for a single step
vec3 step_scattering = calcScattering(dot(view_vector, sun_vector));
step_scattering *= 1.0 - step_extinction;

vec3 extinction = vec3(1.0);
vec3 in_scatter = vec3(0.0);

for (int ii = 0; ii <= NUM_STEPS; ii++) {
  // cast a ray towards the sun and calculate the incoming extincted light
  float light_ray_dist = opticalDepth(current_position, sun_vector);
  vec3 incoming_light = sun_light * calcExtinction(light_ray_dist);

  // multiply by scattering to determine how much light
  // bounces into the view direction
  in_scatter += incoming_light * step_scattering * extinction;

  // accumulate extinction for this step
  extinction *= step_extinction;
  current_position += step_size * view_vector;
}

gl_FragColor.xyz = in_scatter;
 

Improved Implementation

The iterative steps in this method can get very expensive to calculate, and since I did not want to precompute a solution I needed to simplify it even further. Instead of calculating all those incoming light values for a particular view ray, what if we could calculate the average light value hitting the view ray? Then we could use that light value in place of the constant light in the initial solution.

The method I devised works by picking a single point along the view ray to use to calculate this “average” incoming light. This is not correct, but with some experimentation and analysis of a few example scenarios I was able to come up with a result that looks reasonable. When the sun is high in the sky I select a point farther from the viewer; this means there is less atmosphere for the light to go through, so less light is extincted (this makes the sky more blue). When the sun is close to the horizon I select a point closer to the viewer. In this case there is more atmosphere for the light to travel through, so more light is extincted (this makes the sky more red).

This actually works extremely well, producing accurate-looking skies for different times of day. It is also very efficient to render, leaving more GPU cycles free to draw the rest of the scene! Here is the unoptimized GLSL code:

// optical depth along view ray
float ray_dist = opticalDepth(vec3(0.0), view_vector);
vec3 extinction = calcExtinction(ray_dist);

// optical depth for incoming light hitting the view ray
float ray_scale = 0.15 + 0.75 * sun_vector.z;
vec3 light_ray_pos = view_vector * (ray_dist * ray_scale);
float light_ray_dist = opticalDepth(light_ray_pos, sun_vector);

// optical depth for edge of atmosphere:
// this handles the case where the sun is low in the sky and
// the view is facing away from the sun; in this case the distance
// the light needs to travel is much greater
float light_ray_dist_full = opticalDepth(view_vector * ray_dist, sun_vector);
light_ray_dist = max(light_ray_dist, light_ray_dist_full);

// cast a ray towards the sun and calculate the incoming extincted light
vec3 incoming_light = sun_light * calcExtinction(light_ray_dist);

// calculate the in-scattering
vec3 scattering = calcScattering(dot(view_vector, sun_vector));
scattering *= 1.0 - extinction;

// combine
vec3 in_scatter = incoming_light * scattering;

gl_FragColor.xyz = in_scatter;

Finishing Touches

In addition to the sky dome itself (which is modeling the light bounced around from the sun), I also render the sun. This is merely a disk (calculated by angle between view vector and light vector) drawn with the extincted sunlight color.

I also apply extinction to the sunlight color used when lighting the rest of the scene. This way the lighting responds to the change of sun angle and matches the sky. Also I calculate in-scattering for the light direction and add that to the lighting (if you don’t do this you will notice the light is over-saturated at noon).

Note that these calculations are all in linear space, so make sure you are using linear space lighting with HDR tonemapping.

For customization I expose various parameters to the user. The planet size and atmosphere thickness are exposed as multipliers on Earth values (so 1.0 for each is Earth normal). Atmosphere density and pollution sliders affect Rayleigh and Mie scattering strength and phase. Finally, the sun disk size and intensity are tunable as well.

Optimizations

Optimization is key to achieving good performance for any graphics technique. The first pass of optimizations I did on the sky rendering are standard shader optimizations. To start, I moved all operations on values that are constant for the frame to the CPU. This includes calculating the scattering values from the input parameters, precomputing the square of the atmosphere radius for the raycast code, and various others. After that I reformulated many of the calculations in the shader to match the typical low-level shader instruction set. See Emil Persson’s Low-Level Thinking In High-Level Languages for more details.

For my second optimization pass, I realized that there is no reason to draw the sky using a dome or sphere mesh. Instead, a single fullscreen quad drawn with depth forced to the far end of the depth range works just as well with far less vertex and pixel-quad overhead. In fact, you could take this optimization even further by drawing only a single triangle that goes beyond the edges of the screen but covers it entirely; that way you don’t have the repeated pixel-quads along the screen diagonal.

Finally, I draw the sky after drawing opaque geometry, with depth testing turned on. That way we only pay the pixel shader cost for the skydome where it is actually visible.

Final Code

Here is the final optimized shader code:

vec3 calcExtinction(float dist) {
  return exp(dist * sky_params6.xyz);
}

vec3 calcScattering(float cos_theta) {
  float r_phase = (cos_theta * cos_theta) * sky_params6.w + sky_params6.w;
  float m_phase = sky_params1.w * pow(sky_params2.w * cos_theta + sky_params3.w, -1.5);
  return sky_params2.xyz * r_phase + (sky_params3.xyz * m_phase);
}

float baseOpticalDepth(in vec3 ray) {
  float a1 = sky_params4.x * ray.z;
  return sqrt(a1 * a1 + sky_params4.w) - a1;
}

float opticalDepth(in vec3 pos, in vec3 ray) {
  pos.z += sky_params4.x;
  float a0 = sky_params4.y - dot(pos, pos);
  float a1 = dot(pos, ray);
  return sqrt(a1 * a1 + a0) - a1;
}

float cos_theta = dot(view_vec, sun_vector);

// optical depth along view ray
float ray_dist = baseOpticalDepth(view_vec);

// extinction of light along view ray
vec3 extinction = calcExtinction(ray_dist);

// optical depth for incoming light hitting the view ray
vec3 light_ray_pos = view_vec * (ray_dist * sky_params4.z);
float light_ray_dist = opticalDepth(light_ray_pos, sun_vector);

// optical depth for edge of atmosphere:
// this handles the case where the sun is low in the sky and
// the view is facing away from the sun; in this case the distance
// the light needs to travel is much greater
float light_ray_dist_full = opticalDepth(view_vec * ray_dist, sun_vector);

light_ray_dist = max(light_ray_dist, light_ray_dist_full);

// cast a ray towards the sun and calculate the incoming extincted light
vec3 incoming_light = calcExtinction((light_ray_dist);

// calculate the in-scattering
vec3 scattering = calcScattering(cos_theta);
scattering *= 1.0 - extinction;

// combine
vec3 in_scatter = incoming_light * scattering;

// sun disk
float sun_strength = clamp(cos_theta * sky_params1.x + sky_params1.y, 0.0, 1.0);
sun_strength *= sun_strength;
vec3 sun_disk = extinction * sun_strength;

gl_FragColor.xyz = sky_params5.xyz * (sky_params5.w * sun_disk + in_scatter);
gl_FragColor.w = 1.0;
 

And the CPU setup code:

function sqr(val) {
  return val * val;
}

function pow4(val) {
  val = val * val;
  return val * val;
}

var sky_lambda = Vec3.create(680e-9, 550e-9, 450e-9);
var sky_k = Vec3.create(0.686, 0.678, 0.666);

var earth_radius = 6.371e6;
var earth_atmo_thickness = 0.1e6;

var clarity = 1 + sky_props.clarity;
var two_pi = 2 * Math.PI;

// compute betaR
var factor = 1.86e-31 / (clarity * Math.max(sky_props.density, 0.001));
this.sky_params2.x = factor / pow4(sky_lambda.x);
this.sky_params2.y = factor / pow4(sky_lambda.y);
this.sky_params2.z = factor / pow4(sky_lambda.z);

// combute betaM
factor = 1.36e-19 * Math.max(sky_props.pollution, 0.001);
this.sky_params3.x = factor * sky_k.x * sqr(two_pi / sky_lambda.x);
this.sky_params3.y = factor * sky_k.y * sqr(two_pi / sky_lambda.y);
this.sky_params3.z = factor * sky_k.z * sqr(two_pi / sky_lambda.z);

// betaR + betaM, -(betaR + betaM), betaR / (betaR + betaM), betaM / (betaR + betaM)
this.sky_params1.fromAdd(this.sky_params2, this.sky_params3);
this.sky_params6.fromAdd(this.sky_params2, this.sky_params3).mulScalar(-1);
this.sky_params2.div(this.sky_params1);
this.sky_params3.div(this.sky_params1);

// mie scattering phase constants
var g = (1 - sky_props.pollution) * 0.2 + 0.75;
this.sky_params1.w = sqr(1 - g) / (4 * Math.PI);
this.sky_params2.w = -2 * g;
this.sky_params3.w = 1 + sqr(g);

var planet_radius = earth_radius * sky_props.planet_scale;
var atmo_radius = planet_radius + earth_atmo_thickness * sky_props.atmosphere_scale;
this.sky_params4.x = planet_radius;
this.sky_params4.y = atmo_radius * atmo_radius;
this.sky_params4.z = 0.15 + 0.75 * light_props.sun_vector.z;
this.sky_params4.w = atmo_radius * atmo_radius - planet_radius * planet_radius;

// sun disk cutoff
this.sky_params1.y = -(1 - 0.015 * sky_props.sun_disk_radius);
this.sky_params1.x = 1 / (1 + this.sky_params1.y);
this.sky_params1.y *= this.sky_params1.x;

this.sky_params5.fromVec3(this.sun_color).mulScalar(sky_props.brightness);
this.sky_params5.w = sky_props.sun_disk_intensity;

this.sky_params6.w = clarity * 3 / (16 * Math.PI);
 

In the final article of this series I will explain how we render clouds in Cloud Party, including details of shape generation and lighting.