1. Introduction

As a ray of light travels through a medium, the radiance, commonly called the intensity, may change depending on the contents of the medium. A medium which affects the radiance of light is called a participating medium and such a medium can influence light in several ways. Scattering of light happens in a vast variety of environments and is used in creating visual effects such as: atmospheric haze, sky-light computation, light passing through clouds or smoke, light passing through water and sub-surface scattering. These effects are very complex and usually require large amounts of computation due to the mathematical functions that need to be evaluated.

For the purpose of real-time computer graphics, scattering algorithms are often simplified and adapted to run on graphics hardware. Another option is to allow artists to setup functions and parameters that mimic real light scattering, thus saving computation time.

2. Harris' Model for Cloud Rendering and Lighting

The solution proposed by Harris is to approximate the scattering integral over the volume of the cloud using graphics hardware to speed up the process. The basic ideas are

1. The radiance absorbed in each cloud volume unit, modeled as a metaball, is stored in a texture which is used for splatting.
2. The product which approximates the integral is calculated by reading from the draw buffer the previous splat result and using it for the current splat which is blended back into the buffer.
3. Two scattering directions are used: from the sun to the cloud center and from the eye point to the cloud center. This accounts for most of the light scattered in the cloud.

To improve the system, the clouds can be rendered in impostor textures at a resolution based on distance from the camera. The impostors are updated when the change in angle between the camera and the cloud center increases above a certain threshold.

3. Cloud Lighting

First of all, the clouds are stored as an array of particles which represent the metaballs. The shapes can be modeled in various ways – by an artist, using a fluid motion equation solver, using procedural noise techniques – but that is beyond the scope of this article. The particles have position, color, size, albedo (the ratio of reflected light to the incoming light), extinction (reduction of the intensity of light) and an alpha component can be added for simulating cloud formation and extinction.

When lighting the clouds we are all approximating light scattered in the light's direction (from the cloud center to the sun). To do this we will sort our particles away from the light position thus the closest particle will be rendered first. This is fairly obvious since the particle which is not occluded in any way should receive the most amount of light. We can use the square of the distance to the point for sorting and this is stored in the DistanceToCam member. We will use a function that can sort an array of particles both away and towards a point:

switch (mode)
{
  case SORT_AWAY:
    sort(Cloud->Puffs.begin(), Cloud->Puffs.end(), SortAway);
    break;
  case SORT_TOWARD:
    sort(Cloud->Puffs.begin(), Cloud->Puffs.end(), SortToward);
    break;
}

class SortAwayComparison
{
public:
  bool operator () (CloudPuff puff1, CloudPuff puff2)
  {
    return puff1.DistanceToCam < puff2.DistanceToCam;
  }
} SortAway;

class SortTowardComparison
{
public:
  bool operator () (CloudPuff puff1, CloudPuff puff2)
  {
    return puff1.DistanceToCam > puff2.DistanceToCam;
  }
} SortToward;

First of all after sorting our particles we need to setup the camera to be placed in the sun's position, viewing the cloud center and the projection should map the cloud onto the whole viewport. The size of the viewport can be chosen arbitrarily but a value of 32, proposed by Harris is fine. We will use an orthographic projection because it will not deform far away particles.

glMatrixMode(GL_PROJECTION);
glPushMatrix();
glLoadIdentity();
glOrtho(-Cloud->Radius-pr, Cloud->Radius+pr,
        -Cloud->Radius-pr, Cloud->Radius+pr, d - r, d + r);

glMatrixMode(GL_MODELVIEW);
glPushMatrix();
glLoadIdentity();
gluLookAt(Sun.x, Sun.y, Sun.z,
          Cloud->Center.x, Cloud->Center.y, Cloud->Center.z,
          0, 1, 0);

glPushAttrib(GL_VIEWPORT_BIT);
glViewport(0, 0, SplatBufferSize, SplatBufferSize);

The lighting equation given by Harris is:

This equation relates the intensity of the current particle I_k with the intensity of the previous particle I_k-1 and the transparency of the previous particle T_k-1. The other term is g_k-1 where g_k = a_k x τ_k x p(l, -l) x I_k x γ/4π. As we can see g_k-1 is related to the intensity of the k-1 fragment so we will have to read this value from the frame buffer. All the other values in the equation (albedo, extinction, solid angle) are constants and the p(l, -l) element is a phase function which will be discussed later. The original intensity I₀ is fully bright thus the buffer is first cleared to white.

This equation can be encoded for graphics hardware through blending. Blending calculates the sum of the incoming fragment multiplied with the "source factor" and the existing fragment in the buffer multiplied with a "destination factor". In our case the source factor is 1 (as g_k-1 doesn't have any coefficient) and the destination factor isT_k-1, the transmittance of the fragment in the splat texture. The result of the blend operation is a color which is then read back for the next particle. Since opacity is stored in the splat texture, transmittance T_k-1 will equal one minus opacity. This gives the blend function with parameters GL_ONE, GL_ONE_MINUS_SRC_ALPHA.

A rough description of the lighting process can be formulated: we start from a fully bright buffer and use the splat textures to decrease the luminance thus "darkening" the color of the particles as they get farther from the light.

We start by looping over the particles and calculating the screen position of it's center using projection:

double CenterX, CenterY, CenterZ;
gluProject(Cloud->Puffs[i].Position.x, 
           Cloud->Puffs[i].Position.y, 
           Cloud->Puffs[i].Position.z, 
           mm, mp, vp, &CenterX, &CenterY, &CenterZ);

(note: here the puffs were in world space)

The, using the solid angle over which we will read back the pixels from the buffer, the size of the splat buffer and the cloud radius we compute the area which will be read:

Area = Cloud->Puffs[i].DistanceToCam * SolidAngle; //squared distance
Pixels = (int)(sqrt(Area) * SplatBufferSize / (2 * Cloud->Radius));
if (Pixels < 1) Pixels = 1;

ReadX = (int)(CenterX-Pixels/2);
if (ReadX < 0) ReadX = 0;
ReadY = (int)(CenterY-Pixels/2);
if (ReadY < 0) ReadY = 0;

buf = new float[Pixels * Pixels];
//we only need the red component since this is greyscale
glReadBuffer(GL_BACK);
glReadPixels(ReadX, ReadY, Pixels, Pixels, GL_RED, GL_FLOAT, buf);

Finally we compute the average intensity in the area and calculate the color of the current particle which will be splatted, following the equation above.

avg = 0.0f;
for (j = 0; j < Pixels * Pixels; j++) avg += buf[j];
avg /= (Pixels * Pixels);

delete [] buf;

//Light color * 
// average color from solid angle (sum * solidangle / (pixels^2 * 4pi))
// * albedo * extinction
// * rayleigh scattering in the direction of the sun (1.5f)
// (only for rendering, don't store)

factor = SolidAngle / (4 * PI);

ParticleColor.R = LightColor.R * Albedo * Extinction * avg * factor;
ParticleColor.G = LightColor.G * Albedo * Extinction * avg * factor;
ParticleColor.B = LightColor.B * Albedo * Extinction * avg * factor;
ParticleColor.A = 1.0f - exp(-Extinction);

This color, stored in ParticleColor will be stored as the color which will be used for rendering later on, but when we are splatting the particles for lighting we need to include the phase function. This is always equal to 1.5 in the direction of the light so we scale up our color by this value before rendering this particle as a billboard.

The phase function