About

Other posts

Real-time procedural sky rendering

Created in 2023, last updated in June 2024

This post demonstrates how a procedural sky can be rendered in real-time. I will describe some of the methods that I found relevant for my case, which is to render the sky being at low altitude and above the ocean. A naive implementation is firstly detailed, then the computation is gradually complexified to obtain a more convincing sky.

1. Atmospheric scattering

The sky color is due to the atmospheric scattering, that comes from the interaction of the particles in the atmosphere with the light. There are two main scattering types:

The scattering effect is locally defined by a phase function, that expresses the ratio between the incoming and the outcoming light at a certain angle. It is assumed that the Rayleigh and Mie scatterings do not interfere, their effect is therefore additive. With \( h \) being the altitude and \( \theta \) the angle between the incoming and outcoming light direction, the scattering is computed as following: $$ \sigma(h,\theta) = \beta_R(h) \; \phi_R(\theta) + \beta_M(h) \; \phi_M(\theta) $$ where

It is admitted that \( h_{R0} = 8 \) km and \( h_{M0} = 1.2 \). Both \( \beta_{R0} \) and \( \beta_{M0} \) are coefficients that control the breadth of the scattering. The scattering induces a light absortion, also called light extinction, which is the part of the light being scattered away when following a straight path: $$ \alpha(path) = \exp ( - d_{optic} ) $$ with $$ d_{optic} = \beta_{R0} \int_{path} { e^{-h/h_{R0}} dl } + \beta_{M0} \int_{path} { e^{-h/h_{M0}} dl } $$

2. Single scattering implementation

The single scattering is often detailed and implemented in online resources. It is quick to implement and it allows real-time rendering with a few tricks. The single scattering method consists in the discretization of the view ray in multiple points, and to sum the light being scattered at all of those points: $$ \overline{L} = L_{sun} \int_{P \in ray} { \alpha(V \rightarrow P) \; \alpha(P \rightarrow Sun) \; \sigma(P) \; dl } $$ where \( \alpha(path) \) is the light extinction along a straight path and \( \sigma(P) \) the scattering at the point P.

Extinction of the sun light to a point:

The extinction between the sun and a point in the atmosphere is an integral that is costly to evaluate by summation. To prevent this, many resources are using a pre-computed lookup texture. However, I decided to use a math approximation instead. I thus express the optical-distance to the sun with: $$ d_{optic,sun} = \int_{P}^{Sun} { e^{-h/h0} dl } \simeq \int_{x=0}^{+\infty} { e^{ - ( a x^2 + b x + c ) } dx } = \frac{\sqrt{\pi}}{2\sqrt{a}} e^{-c} e^{\frac{b^2}{4a}} \left( 1 - \mathrm{erf} \left( \frac{b}{2\sqrt{a}} \right) \right) $$ where \( a = \frac{1-cos^2 \theta}{2 R h0} \) , \( b = \frac{cos \theta}{h0} \) , \( c = \frac{h_{start}}{h0} \) and \( R \) the Earth radius. This approximation is fairly good while the ray does not travel towards the ground (while \(b\) is positive). I also use asymptotic approximations of the erf function to have a friendlier expression to compute. When the ray travels towards the ground, I compute a new starting location at the lowest altitude position: $$ h_{min} = (R + h_{start}) \sqrt{1 - cos \theta ^2} - R $$ If this position is below the ground, then the extinction is infinite (the light is blocked by the ground). Otherwise, the optical distance between the starting point and the lower altitude position is approximated with: $$ d_{optic,sun,pre} = \frac{-\cos \theta}{2} \left( e^{-h_{start}/h0} + e^{-h_{min}/h0} \right) (R + h_{start}) $$ Finally, the optical distance to the sun is expressed with: $$ d_{optic,sun} = d_{optic,sun,pre} + \begin{cases} e^{-c} \frac{1}{b} \left(1 - \frac{2a}{b^2} + \frac{12a^2}{b^4} \right) \; \; \; \mathrm{ when } \; b > 7 \sqrt{a} \\ e^{-c} \frac{\sqrt{\pi}}{2 \sqrt{a}} \left( 0.3480242 p - 0.0958798 p^2 + 0.7478556 p^3 \right) \; \; \; \mathrm{ with } \; \; \; p = \frac{2 \sqrt{a}}{2 \sqrt{a} + 0.47047 b} \end{cases} $$

extinction visualization,
at surface (0 km) 		extinction visualization,
at 40 km altitude

Extinction between the viewer and a point:

Since the light is already accumulated along the view ray by summation, the extinction will be compute in the same loop. While the scattered light is evaluated on a single point per segment, the extinction can be evaluated with a more accurate formula:

Results:

This implementation of the sky scattering gives a plausible result with low effort, however the resulting sky isn't realistic enough. In short terms, at midday, the sky looks too dark compared with the horizon light's intensity. To increase the scattering strength may appear as a solution, but the sky will appear yellowed even at midday. The overall light intensity behaves weirdly, with a maximal intensity not at midday. The light intensity is lower on the sides, producing spurious pattern. This comes from the Rayleigh phase function. With Mie scattering, the horizon appears dark (not shown here). This is because the poor sampling of the Mie scattering effect, where the light is highly absorbed at low altitude. On the picture below, the sky is rendered with Rayleigh scattering only, from midday to sunset.

Shader code:

(click to expand). 
uniform vec3  sunDirectionIncoming;
uniform float factorRayleigh;
uniform float factorMie;
uniform float gMie;
const float rE   = 6.3e3f; // radius of Earth (km)
const vec2  h0   = vec2(8., 1.2); // ch. altitude (Rayleigh, Mie) (km)
const float hMax = 50.f; // max altitude for computations
const vec3  scatteringCoefs_R = vec3(0.11f, 0.32f, 0.68f);
const vec3  scatteringCoefs_M = vec3(0.50f, 0.50f, 0.50f);

vec3 localScatteringLight(float h, float c)
{
  vec3 fR = exp(-h / h0.x) * scatteringCoefs_R * factorRayleigh;
  vec3 fM = exp(-h / h0.y) * scatteringCoefs_M * factorMie;
  float cc = c * c;
  float gg = gMie * gMie;
  float den = 1.f + gg - 2.f * c * gMie;
  float phaseR = 0.05968310365946075f * (1.f + cc);
  float phaseM = 0.1193662073189215f * (1.f - gg) * (1.f + cc) / (sqrt(den * den * den) * (2.f + gg));
  return fR * phaseR + fM * phaseM;
}

vec2 opticalDistanceToSun(float csz, float hStart, float h0)
{
  float optD = 0.f;
  if (csz < 0.f)
  {
    float hMin = (rE + hStart) * sqrt(1.f - csz * csz) - rE;
    if (hMin < 0.f) return 1.e30f; // + infinite
    optD += 0.5f * (exp(-hMin/h0) + exp(-hStart/h0)) * (rE + hStart) * (-csz);
    hStart = hMin;
    csz = 0.f;
  }
  float a = 0.5f * (1.f - csz * csz) / rE / h0;
  float b = csz / h0;
  float c = hStart / h0;
  float sqrta = sqrt(a);
  float expc = exp(-c);
  if (b > 7.f*sqrta)
  {
    float ainvbb = a / (b * b);
    optD += expc * (1.f - 2.f*ainvbb + 12.f*ainvbb*ainvbb) / b;
  }
  else
  {
    float p = (2.f * sqrta) / (2.f * sqrta + 0.47047f * b);
    float ret = 0.5f * expc * sqrt(3.14159f / a) * p * (0.3480242f - 0.0958798f * p + 0.7478556f * p * p);
    optD += min(ret, 1.e30f);
  }
  return optD;
}

vec3 skyColor(vec3 ray)
{
  float dl = 1.f;
  vec3  pt = vec3(0.f, 0.f, rE) + 0.5f * dl * ray; // mid-point integration
  vec3  lightAcc = vec3(0.f);
  vec2  optD_VP = vec2(0.f, 0.f); // optical-distance: viewer - point
  while (true)
  {
    float rP      = length(pt);
    float hP      = rP - rE; // elevation
    if (hP < 0.f)  return vec3(0.f); // ground reflexion ?
    if (hP > hMax) break; // we're done
    float csz     = clamp(dot(pt, -sunDirectionIncoming) / rP, -1.f, 1.f);
    vec2  optD_PS = vec2(opticalDistanceToSun(csz, hP, h0.x) * factorRayleigh,
                         opticalDistanceToSun(csz, hP, h0.y) * factorMie);
    vec2  rhoP    = exp(-hP/h0) * vec2(factorRayleigh, factorMie);
    // direct scattering
    {
      vec3 globalOptD = (optD_VP.x + optD_PS.x) * scatteringCoefs_R +
                        (optD_VP.y + optD_PS.y) * scatteringCoefs_M;
      vec3 tr0 = glm::exp(-globalOptD);
      vec3 _localExtinct = rhoP.x * scatteringCoefs_R + rhoP.y * scatteringCoefs_M;
      vec3 trSelf = (1.f - exp(-dl * _localExtinct)) / _localExtinct;
      vec3 localScatteringPS = localScatteringLight(hP, cosRS);
      lightAcc += tr0 * trSelf * localScatteringPS;
    }
    pt += dl * ray;
    optD_VP += step * rhoP;
  }
  {
    float cosRS = dot(-sunDirectionIncoming, ray); // cos of angle between the ray and the sunlight
    float sunDisk = min(exp(6.e4f * (cosRS-0.999962f)), 1.f); // the sun-disk is 0.25 degree radius
    lightAcc += sunDisk * vec3(1.f, 1.f, 1.f);
  }
  return lightAcc;
}

3. Sky illumination with multiple atmospheric scattering

To improve the results, the indirect scattering needs to be taken into account. Badfully, methods based on ray-tracing would be too computationally intense to stay in real time. One of the possible methods is to model the entire atmosphere, and iteratively refine the atmosphere global illumination. This way, the indirect scattering will natively be taken into account.

The atmosphere modeling consists in partitioning the atmosphere into pieces. Each cell is defined by an altitude range and an angular range. Because that the Earth - Sun system has a radial symmetry, a cell is characterized only by its altitude and its angular position relatively on the sun direction. In fact, a single cell represents a whole ring defined by its altitude and the local sun elevation. Therefore, atmosphere model can be described in a 2D space.

Each cell would in turn receive and emit light. The light is emitted in all directions. I have modeled the emitted light with two contributions. The first is the scattering from the direct sun light. It is kept as an analytic expression. It is computed once because it is constant alongside the computations. The second is the scattering from the light emitted by the atmosphere. It is interpolated with some spherical harmonics. Finally, the emitted light by a cell is defined by: $$ L_{emit}(\vec{r}) = c_R \phi_R(\theta(\vec{r})) + c_M \phi_M(\theta(\vec{r})), g) + \sum_i c_i \psi_i(\vec{r})) $$ where \( \phi_R \) and \( \phi_M \) are the Rayleigh and Mie phase functions, \( c_R \), \( c_M \) the coefficients of the direct scattering, \( \psi_i \) the spherical harmonics, and \( c_i \) their coefficients. The computations looks like: 

atmModel.reset();
for (cell in atmModel)
  cell.emitDirect += scatterLight(receivedLightSun(), sunDirection());
for (int step = 0; step < 6; step += 1)
{
  for (cell in atmModel)
  {
    emitIndirectNew = 0;
    for (direction in sphere)
      emitIndirectNew += scatterLight(atmModel.receivedLightFromCells(direction), direction);
    cell.emitIndirect = emitIndirectNew;
  }
}

These computations can be issued asynchronously on the CPU or on the GPU. Furthermore, the atmosphere global illumination does not depend on the sun elevation angle, so its computation won't be issued on every frame. In fact, it only depends on the atmosphere properties, which are:

Results:


 Received light by the cell
(without the direct sun light) 		Emitted light by the cell Cell map
Cell at the 45° sun-inclination,
in the air,
with Rayleigh scattering only
Cell at the 0° sun-inclination,
near the ground,
with Rayleight and Mie scattering

Sky-dome generation:

The sky dome illumination by ray-marching on the atmosphere cells. This process can be issued asynchronously too, and the real-time shader would just sample a sky map. The sky domes can be generated at any step for the atmosphere global illumination computation. Here 2 examples, generated from the step 0, taking account of the direct scattering only, and from the step 3, taking account of 3 levels of indirect scattering.

Sky dome at midday,
with Rayleigh scattering 		Sky dome at sunset,
with Rayleigh and Mie scattering

I can observe that the indirect scattering does not contribute significantly to the sky illumination. It still improve the realism. Some work are planned in the future to improve the global feeling, without excluding to artificially exaggerate the amount of indirect scattering.

Sampling:

In the final aplication, the sky shader consists in sampling a sky dome map. Let define the sky dome coordinate:

The shader code folds into few lines of code: 

uniform sampler2D skyLookUp;
uniform vec3      skyLookUpFactor;
const int         skyLoopUpResolution = 24;
vec3 skyColor(vec3 ray)
{
  float sP = abs(ray.z); // sin(p)
  float cY = dot(-sunDirectionIncoming.xy, ray.xy) / (sqrt(1.f - sP*sP) + 1.e-6f); // cos(y)
  // non-uniform sampling, based on resolution of 24x24 pixels:
  vec2  uv = 0.021f + 0.958f * vec2(sP, 1.f - sqrt(0.5f * (1.f - cY)));
  return skyLookUpFactor * texture(skyLookUp, uv).xyz;
}

Note that the sky-domes can be tabulated, based on multiple sky configurations (sun evelation, scattering parameters). For example, I decided to tabulate with 6 * 16 table-points when targeting mobile devices or web platform.

Sky-dome regression:

We can go a step further by fitting the sky-dome illumination in an arbitrary model, with only few coefficients. Finally, only few floats will have to be sent to the GPU shader. The GPU shader now needs to reconstruct the sky-dome illumination, with help of few arithmetic functions.

Work in progress ...

The shader code folds into few lines of code: 

uniform vec3 sunDirectionIncoming;
uniform vec3 skyCoefs[5];
vec3 skyColor(vec3 ray)
{
  float sP = abs(ray.z);
  float cYcP = dot(-sunDirectionIncoming.xy, ray.xy);
  float cS = dot(-sunDirectionIncoming, ray);
  vec3  color = ...;
  return max(vec3(0.f, 0.f, 0.f), color);
}

4. Rendering of the clouds

The clouds are the visible effect of the scattering induced by the water droplets. Then we only need to consider the Mie scattering. I will use an even simpler scattering phase function: a constant value. Indeed, the multiple scattering tends to scatter the light as if a smoother phase function was initially used. The scattered light is: $$ \sigma_{cloud}(\rho) = \frac{1}{4 \pi} \rho $$ where \( \rho \) is the local scattering strength, based on the cloud's particle local density. The light extinction is computed the same way: $$ \alpha(path) = \exp \left( - \int_{path} { \rho dl } \right) $$ A simple direct ray-cast method, enhanced with some pre-computations, is enough to have convincing results.

Work in progress ...

2D-based clouds:

I use 2D-textures to model the shape of the clouds, in addition of 2D-noise. On the view ray, the local scattering strengh of the cloud is considered constant. Let's use \( D \) as a fake depth and \( \rho \) as the scattering strengh. Then the received scattered light can be expressed as: $$ L = \int_{0}^{D} { \left( \frac{1}{4 \pi} L_{sun} + L_{ambiant} \right) \rho e^{-\rho l} dl } = \left( \frac{1}{4 \pi} L_{sun} + L_{ambiant} \right) \left( 1 - e^{-\rho D} \right) $$ We can observe that the scattered sky color is constant (leaving out the density-dependent multiplier). To obtain better results, we can fake the cloud volume by applying a faked normal. As such the intensity of the scattered sun light would be modulated.

Shader code (click to expand). 
TODO

3D-based clouds:

I use 3D-textures to model the shape of the clouds, in addition of 3D-noise. On the texture, additional pre-computed geometry parameters are embeded.

Shader code (click to expand). 
TODO