This post is the first in a series to follow-up on my 2012 GPU Pro 3 article about atmospheric scattering . What I showed there was a full single-scattering solution for a planetary atmosphere running in a pixel shader, dynamic and in real time, without pre-computation or simplifying assumptions. The key to this achievement was a novel and efficient way to evaluate the Chapman function , hence the title. In the time since then I have improved on the algorithm and extended it to include aspects of multiple scattering. The latter causes horizontal diffusion (twilight situations) and vertical diffusion (deep atmospheres), and neither can be ignored for a general atmosphere renderer in a space game, for example.
I have written a Shadertoy that reflects the current state of affairs. It’s a mini flight simulator that also features clouds, and other rendering goodies. A WebGL2 capable browser is needed to run it. Under Windows, the ANGLE/Direct3D translator may take a long time to compile it (up to a minute is nothing unusual, but it runs fast afterwards). When successfully compiled it should look like this:
A little background story
The origins of my technique go back to experiments I did privately while hacking on my old Space Glider. Originally this was a demo written for DOS/DJGPP without 3-D acceleration. It featured a mini solar system with non-photorealistic graphics, in the style of Mercenary and Starglider. Sometime I dusted it off and wanted to make it more real™: larger scale, atmospheric scattering, procedural terrain, HDR, dynamic light adaptation, secondary illumination from celestial bodies, spherical harmonics, real flight dynamics, real orbital mechanics, … all that, you know, a real™ space simulator needs. Naturally, the source didn’t get prettier in the process and the project became an unreleasable mess.
Nevertheless, I came up with a working prototype that had scattering and a procedural terrain. Since it was all calculated on the CPU it had to be super efficient. I devised an analytical approximation to the airmass integral (what later became the Chapman function) so I could reduce the numerical integration from two nested loops to one loop. I made heavy use of SIMD, first by calculating the sample points in parallel, and then in a second pass over the sample points doing the color channels in parallel. This is what the result looked like (circa 2011, unreleased):
I want to emphasize that the above is rendered in software with nothing but Gouraud shading and vertex colors. This just as a reminder how nice looking visuals can be driven by low tech. It was also relatively fast (40 FPS @ 720p).
Moving it to the GPU
The article that I wrote for GPU Pro 3 describes the same algorithm as that which was running in my prototype, just with all code moved to the pixel shader, and wrapped in a ray-tracing framework (essentially “shadertoying” it). This approach is wasteful, a fact that was correctly pointed out in . One can no longer benefit from triangle interpolation and cannot parallelize threads across sample points; instead the whole thing is a big loop per pixel.
Instead, atmospheric scattering should ideally be done with a kind of map-reduce operation: First calculate all sample points in parallel (the map part) and then collect the result along rays (the reduce part), similar to how my prototype worked. This could practically be done a compute shader, but I did not yet explore this idea further. The shadertoy that was linked in the beginning of this post also does the big-loop-in-the-pixel-shader approach, for lack of other practical alternatives.
Aerial perspective and spectral sampling
The visual impact of atmospheric scattering is called aerial perspective. It is not the same as fixed-function fog, but if we want to use that as an analogy, we would have to have a fog coordinate for each color channel. To illustrate, I have taken the photo from the top of the linked Wikipedia article and separated its color channels (see the insets above). It should be evident that the distant mountains are blue not because of some blue “fog color” (that would be white, in the limit of infinite distance), but rather because each channel has a different multiplier for the fog distance, the blue channel having the highest.
Formally, aerial perspective is a transformation of some original color by a multiplicative operation, the transmittance , and an additive operation, the in-scatter (see also  for an introduction of these terms):
where , and are color vectors. In order to be universally valid and independent of the choice of color basis, the transmittance operator would have to be a tensor, as in
where , are spectral indices. If (and only if) is diagonal would it be appropriate to use the more familiar component-wise form
This opens up the question of color space. is trivially diagonal if the basis is simply the full spectrum (in infinite resolution). Otherwise the color basis must be carefully chosen to have properties of a chromatic adapation space. Following this clue I took the dominant wavelengths of the LMS cone responses as a starting point to settle for 615nm for the red wavelength, 535nm for the green wavelength and 445nm for the blue wavelength. When these wavelengths are used as sample points together with a proper transform, the resulting colors are very accurate. I used the new CIE 2012 CMFs to calculate the transform from spectral samples to sRGB colors, as follows:
Side note: The matrix from above looks red-heavy, because it contains an adaptation from white point E to D65. This was done so that the color components in the chosen 615-, 535-, 445 nm basis can be used directly as spectral point samples. For example: the color of the sunlight spectrum, normalized to unit luminance, would be in this basis, which corresponds very well to a point sampling of the smoothed sunlight spectrum at the wavelengths 615 nm, 535 nm and 445 nm. The transformation by then yields a (linear) sRGB color of which has a correlated color temperature of 5778 K.
Geometry of the aerial perspective problem
The next problem is to find the numerical values for and . Without further simplifying assumptions (such as a plane-parallel and/or homogeneous atmosphere) there is no closed form solution for this. While the transmittance is simply the antilogarithm of , the optical depth requires a path integral over the attenuation coefficient (itself the reciprocal of mean free path, so it is density dependent). The in-scatter is then a path integral over the source function attenuated by transmittance along the way. With single scattering of sunlight, and disregarding thermal emission, the source function equals incoming solar flux times a phase function times the single scattering albedo , and again attenuated by transmittance (note how the product is always dimensionless):
It is evident that finding along (semi-infinite, treating the sun as infinitely distant) rays appears as a sub-problem of evaluating . A real-time application must be able to do this in O(1). But again there is, in general, no simple solution for , since air density decreases exponentially with altitude, and the altitude above a spherical surface varies non-linearly along straight paths (and we have not even considered bending of rays at the horizon).
The Chapman function
As it is so often, problems in computer graphics have already been treated in physics literature. In this case, the optical depth along rays in a spherically symmetric, exponentially decreasing atmosphere was considered in the 1930’s by Sydney Chapman . The eponymous Chapman function
is defined as relative optical depth along a semi-infinite ray having zenith angle (read: greek letter chi), as seen by an observer at altitude (where is the radial distance to the sphere center expressed in multiples of the scale height). It is normalized to the zenith direction so for all . The following expression, distilled from a larger one in , can be taken as a virtually exact solution (although it breaks down if is near zero):
Note how this expression contains , the scaled complementary error function, which is a non-standard special function that one would need to implement on its own, so we’re back to square one. This is clearly not real-time friendly and needs to be simplified a lot.
An efficient approximation
My approximation was motivated by the fact that in the limit of a flat earth, must approach the secant function. Therefore I looked for a low-order rational function of , such that it gives unity in the zenith direction and the special value in the horizontal direction. The expression for is already much simpler since vanishes, and I aggressively simplified it even more. The extension to zenith angles was then done by some trigonometric reasoning. The final result, as originally published, is
This is the main result. Note that the expensive case for when only needs to be considered when looking downwards and the view is not already obscured by the terrain, which only applies to a small region below the horizon. A comparison between the exact and the approximate versions of is shown below:
Epilog: A literature review
At this point I want to stress how every real-time algorithm for atmospheric scattering must have an efficient way to evaluate , even when this function is not explicitly identified as such. Here is a small literature survey to illustrate that point:
|Author(s)||Type||Assumptions||appears in (disguised as)|
|Nishita et al. ’93 ||table||Parallel sunrays||Section 4.3.4|
|Preetham et al. ’99 ||analytic||Earthbound observer||Appendix A.1|
|Nielsen ’03 ||analytic||Hardcoded use case||Section 5.4.1|
|O’Neil ’04 ||table||—||Sourcecode|
|O’Neil ’05 ||analytic||Constant||Section 4.2|
|Schafhitzel et al. ’07 ||table||—||Section 4.1|
|Bruneton, Neyret ’08 ||table||—||Section 4 as|
|Schüler ’12 ||analytic||—||Listing 1|
|Yusov ’14 ||analytic||—||Section 2.4.2 as|
This was the first post in a series. I originally didn’t plan to do a series but I decided to split when the word count went over 4000, especially with the source listings. The next post will explain an overflow-free implementation of the chapman function in shader code, how transmittances can be calculated along arbitrary rays and then guide through the procedure of integrating the atmospheric in-scatter along a viewing ray.
 Naty Hoffman and Arcot J Preetham (2002): “Rendering Outdoor Light Scattering in Real Time”
 Sydney Chapman (1931): “The absorption and dissociative or ionizing effect of monochromatic radiation in an atmosphere on a rotating earth.” Proceedings of the Physical Society 43:1, 26.
 Miroslav Kocifaj (1996): “Optical Air Mass and Refraction in a Rayleigh Atmosphere.” Contributions of the Astronomical Observatory Skalnate Pleso 26:1, 23–30.
 Tomoyuki Nishita, Takao Sirai, Katsumi Tadamura and Eihachiro Nakamae (1993): “Display of the earth taking into account atmospheric scattering” Proceedings of the 20th annual conference on Computer graphics and interactive techniques, SIGGRAPH ’93, pp. 175–182.
 Arcot J Preetham, Peter Shirley and Brian Smits (1999): “A practical analytic model for daylight” Proceedings of the 26th annual conference on Computer graphics and interactive techniques, SIGGRAPH ’99, pp. 91–100
 Ralf Stockholm Nielsen (2003): “Real Time Rendering of Atmospheric Scattering Effects for Flight Simulators” Master’s thesis, Informatics and Mathematical Modelling, Technical University of Denmark
 Sean O’Neil (2004): “Real-Time Atmospheric Scattering”
 Sean O’Neil (2005): “Accurate Atmospheric Scattering” GPU Gems 2, Pearson Addison Wesley Prof, pp. 253–268
 Tobias Schafhitzel, Martin Falk and Thomas Ertl (2007): “Real-Time Rendering of Planets with Atmospheres” Journal of WSCG 15, pp. 91–98
 Éric Bruneton and Fabrice Neyret (2008): “Precomputed Atmospheric Scattering” Proceedings of the 19th Eurographics Symposium on Rendering, Comput. Graph. Forum 27:4, pp. 1079–1086
 Christian Schüler (2012): “An Approximation to the Chapman Grazing-Incidence Function for Atmospheric Scattering” GPU Pro 3, A K Peters, pp. 105–118.
 Egor Yusov (2014): “High Performance Outdoor Light Scattering using Epipolar Sampling” GPU Pro 5, A K Peters, pp. 101–126