Followup to Atmospheric Scattering—Part 1

This post is the first in a series to fol­low-​up on my 2012 GPU Pro 3 arti­cle about atmos­pher­ic scat­ter­ing [11]. What I showed there was a full sin­gle-​scat­ter­ing solu­tion for a plan­e­tary atmos­phere run­ning in a pix­el shad­er, dynam­ic and in real time, with­out pre-​com­pu­ta­tion or sim­pli­fy­ing assump­tions. The key to this achieve­ment was a nov­el and effi­cient way to eval­u­ate the Chap­man func­tion [2], hence the title. In the time since then I have improved on the algo­rithm and extend­ed it to include aspects of mul­ti­ple scat­ter­ing. The lat­ter caus­es hor­i­zon­tal dif­fu­sion (twi­light sit­u­a­tions) and ver­ti­cal dif­fu­sion (deep atmos­pheres), and nei­ther can be ignored for a gen­er­al atmos­phere ren­der­er in a space game, for exam­ple.

I have writ­ten a Shader­toy that reflects the cur­rent state of affairs. It’s a mini flight sim­u­la­tor that also fea­tures clouds, and oth­er ren­der­ing good­ies. A WebGL2 capa­ble brows­er is need­ed to run it. Under Win­dows, the ANGLE/Direct3D trans­la­tor may take a long time to com­pile it (up to a minute is noth­ing unusu­al, but it runs fast after­wards). When suc­cess­ful­ly com­piled it should look like this:

A little background story

The ori­gins of my tech­nique go back to exper­i­ments I did pri­vate­ly while hack­ing on my old Space Glid­er. Orig­i­nal­ly this was a demo writ­ten for DOS/DJGPP with­out 3-D accel­er­a­tion. It fea­tured a mini solar sys­tem with non-pho­to­re­al­is­tic graph­ics, in the style of Mer­ce­nary and Star­glid­er. Some­time I dust­ed it off and want­ed to make it more real™: larg­er scale, atmos­pher­ic scat­ter­ing, pro­ce­dur­al ter­rain, HDR, dynam­ic light adap­ta­tion, sec­ondary illu­mi­na­tion from celes­tial bod­ies, spher­i­cal har­mon­ics, real flight dynam­ics, real orbital mechan­ics, ... all that, you know, a real™ space sim­u­la­tor needs. Nat­u­ral­ly, the source didn’t get pret­ti­er in the process and the project became an unre­leasable mess.

Nev­er­the­less, I came up with a work­ing pro­to­type that had scat­ter­ing and a pro­ce­dur­al ter­rain. Since it was all cal­cu­lat­ed on the CPU it had to be super effi­cient. I devised an ana­lyt­i­cal approx­i­ma­tion to the air­mass inte­gral (what lat­er became the Chap­man func­tion) so I could reduce the numer­i­cal inte­gra­tion from two nest­ed loops to one loop. I made heavy use of SIMD, first by cal­cu­lat­ing the sam­ple points in par­al­lel, and then in a sec­ond pass over the sam­ple points doing the col­or chan­nels in par­al­lel. This is what the result looked like (cir­ca 2011, unre­leased):


I want to empha­size that the above is ren­dered in soft­ware with noth­ing but Gouraud shad­ing and ver­tex col­ors. This just as a reminder how nice look­ing visu­als can be dri­ven by low tech. It was also rel­a­tive­ly fast (40 FPS @ 720p).

Moving it to the GPU

The arti­cle that I wrote for GPU Pro 3 describes the same algo­rithm as that which was run­ning in my pro­to­type, just with all code moved to the pix­el shad­er, and wrapped in a ray-trac­ing frame­work (essen­tial­ly “shader­toy­ing” it). This approach is waste­ful, a fact that was cor­rect­ly point­ed out in [12]. One can no longer ben­e­fit from tri­an­gle inter­po­la­tion and can­not par­al­lelize threads across sam­ple points; instead the whole thing is a big loop per pix­el.

Instead, atmos­pher­ic scat­ter­ing should ide­al­ly be done with a kind of map-reduce oper­a­tion: First cal­cu­late all sam­ple points in par­al­lel (the map part) and then col­lect the result along rays (the reduce part), sim­i­lar to how my pro­to­type worked. This could prac­ti­cal­ly be done a com­pute shad­er, but I did not yet explore this idea fur­ther. The shader­toy that was linked in the begin­ning of this post also does the big-loop-in-the-pix­el-shad­er approach, for lack of oth­er prac­ti­cal alter­na­tives.

Aer­i­al per­spec­tive — Pho­to tak­en from Wikipedia, CC-BY-SA Joaquim Alves Gas­par

Aerial perspective and spectral sampling

The visu­al impact of atmos­pher­ic scat­ter­ing is called aer­i­al per­spec­tive. It is not the same as fixed-func­tion fog, but if we want to use that as an anal­o­gy, we would have to have a fog coor­di­nate for each col­or chan­nel. To illus­trate, I have tak­en the pho­to from the top of the linked Wikipedia arti­cle and sep­a­rat­ed its col­or chan­nels (see the insets above). It should be evi­dent that the dis­tant moun­tains are blue not because of some blue “fog col­or” (that would be white, in the lim­it of infi­nite dis­tance), but rather because each chan­nel has a dif­fer­ent mul­ti­pli­er for the fog dis­tance, the blue chan­nel hav­ing the high­est.

For­mal­ly, aer­i­al per­spec­tive is a trans­for­ma­tion of some orig­i­nal col­or L_0 by a mul­ti­plica­tive oper­a­tion, the trans­mit­tance T, and an addi­tive oper­a­tion, the in-scat­ter I (see also [1] for an intro­duc­tion of these terms):

    \[L = T L_0 + I,\]

where L, L_0 and I are col­or vec­tors. In order to be uni­ver­sal­ly valid and inde­pen­dent of the choice of col­or basis, the trans­mit­tance oper­a­tor would have to be a ten­sor, as in

    \[L_\nu = T^\mu_\nu L_0_\mu + I_\nu,\]

where \mu, \nu are spec­tral indices. If (and only if) T is diag­o­nal would it be appro­pri­ate to use the more famil­iar com­po­nent-wise form

    \[L_\nu = T_\nu L_0_\nu + I_\nu.\]

This opens up the ques­tion of col­or space. T is triv­ial­ly diag­o­nal if the basis is sim­ply the full spec­trum (in infi­nite res­o­lu­tion). Oth­er­wise the col­or basis must be care­ful­ly cho­sen to have prop­er­ties of a chro­mat­ic ada­p­a­tion space. Fol­low­ing this clue I took the dom­i­nant wave­lengths of the LMS cone respons­es as a start­ing point to set­tle for 615nm for the red wave­length, 535nm for the green wave­length and 445nm for the blue wave­length. When these wave­lengths are used as sam­ple points togeth­er with a prop­er trans­form, the result­ing col­ors are very accu­rate. I used the new CIE 2012 CMFs to cal­cu­late the trans­form from spec­tral sam­ples to sRGB col­ors, as fol­lows:

    \[M_{\mathrm{615,535,445}\rightarrow\mathrm{sRGB}} = \begin{pmatrix} 1.6218 & -0.4493 & 0.0325 \\ -0.0374 & 1.0598 & -0.0742 \\ -0.0283 & -0.1119 & 1.0491 \end{pmatrix}.\]

Side note: The matrix above con­tains an adap­ta­tion from white point E to D65. This is done so that col­or com­po­nents in the mono­chro­mat­ic basis can direct­ly rep­re­sent spec­tral point sam­ples. For exam­ple, the col­or of the sun­light spec­trum, nor­mal­ized to unit lumi­nance, would be \scriptstyle \left( \begin{smallmatrix}0.9420 & 1.0269 & 1.0241\end{smallmatrix} \right)^T} in this basis, which cor­re­sponds very well to a point sam­pling of the smoothed sun­light spec­trum at the wave­lengths 615 nm, 535 nm and 445 nm. The trans­for­ma­tion then yields a (lin­ear) sRGB col­or of \scriptstyle \left( \begin{smallmatrix}1.0997 & 0.9771 & 0.9329\end{smallmatrix} \right)^T which has a cor­re­lat­ed col­or tem­per­a­ture of 5778 K.

Geometry of the aerial perspective problem

The next prob­lem is to find the numer­i­cal val­ues for T and I. With­out fur­ther sim­pli­fy­ing assump­tions (such as a plane-par­al­lel and/or homo­ge­neous atmos­phere) there is no closed form solu­tion for this. While the trans­mit­tance is sim­ply the antilog­a­rithm of \tau, the opti­cal depth \tau requires a path inte­gral over the atten­u­a­tion coef­fi­cient \alpha (itself the rec­i­p­ro­cal of mean free path, so it is den­si­ty depen­dent). The in-scat­ter is then a path inte­gral over the source func­tion J atten­u­at­ed by trans­mit­tance along the way. With sin­gle scat­ter­ing of sun­light, and dis­re­gard­ing ther­mal emis­sion, the source func­tion equals incom­ing solar flux F times a phase func­tion f times the sin­gle scat­ter­ing albe­do \omega_0, and again atten­u­at­ed by trans­mit­tance (note how the prod­uct \alpha \, \mathrm{d}u is always dimen­sion­less):

    \begin{align*} T &=e^{-\tau}, & \tau &= \int \alpha(u) \mathrm{d}u^{\,}, & I &= \int J(u) e^{-\tau(u)} \alpha(u) \mathrm{d}u, \\ & & & & J &= \omega_0 f(\theta) F_\mathrm{sun} e^{-\tau_\mathrm{sun}}. \end{align*}

It is evi­dent that find­ing \tau along (semi-infi­nite, treat­ing the sun as infi­nite­ly dis­tant) rays appears as a sub-prob­lem of eval­u­at­ing J. A real-time appli­ca­tion must be able to do this in O(1). But again there is, in gen­er­al, no sim­ple solu­tion for \tau, since air den­si­ty decreas­es expo­nen­tial­ly with alti­tude, and the alti­tude above a spher­i­cal sur­face varies non-lin­ear­ly along straight paths (and we have not even con­sid­ered bend­ing of rays at the hori­zon).

The Chapman function

As it is so often, prob­lems in com­put­er graph­ics have already been treat­ed in physics lit­er­a­ture. In this case, the opti­cal depth along rays in a spher­i­cal­ly sym­met­ric, expo­nen­tial­ly decreas­ing atmos­phere was con­sid­ered in the 1930’s by Syd­ney Chap­man [2]. The epony­mous Chap­man func­tion


is defined as rel­a­tive opti­cal depth along a semi-infi­nite ray hav­ing zenith angle \chi (read: greek let­ter chi), as seen by an observ­er at alti­tude x (where x is the radi­al dis­tance to the sphere cen­ter expressed in mul­ti­ples of the scale height). It is nor­mal­ized to the zenith direc­tion so \operatorname{Ch}(x,0)=1 for all x. The fol­low­ing expres­sion, dis­tilled from a larg­er one in [3], can be tak­en as a vir­tu­al­ly exact solu­tion (although it breaks down if x is near zero):

    \begin{align*} \operatorname{Ch}&(x,\chi) = \frac{1}{2} \Bigg[ \cos(\chi) + \operatorname{erfcx} \left( \sqrt{ \frac{x \cos^2\chi}{2} } \right) \left( \frac{1}{x} + 2 - \cos^2\chi \right) \sqrt{\frac{\pi x}{2}} \Bigg] \\ \end{align*}

Note how this expres­sion con­tains \operatorname{erfcx}, the scaled com­ple­men­tary error func­tion, which is a non-stan­dard spe­cial func­tion that one would need to imple­ment on its own, so we’re back to square one. This is clear­ly not real-time friend­ly and needs to be sim­pli­fied a lot.

An efficient approximation

My approx­i­ma­tion was moti­vat­ed by the fact that in the lim­it of a flat earth, \lim_{x\rightarrow \infty}\operatorname{Ch}(x,\chi) must approach the secant func­tion. There­fore I looked for a low-order ratio­nal func­tion of \cos \chi, such that it gives uni­ty in the zenith direc­tion and the spe­cial val­ue c=\operatorname{Ch}(x,90^\circ) in the hor­i­zon­tal direc­tion. The expres­sion for c is already much sim­pler since \cos\chi van­ish­es, and I aggres­sive­ly sim­pli­fied it even more. The exten­sion to zenith angles \chi>90^\circ was then done by some trigono­met­ric rea­son­ing. The final result, as orig­i­nal­ly pub­lished, is

    \begin{align*} \operatorname{Ch}(x,\chi) &\approx \begin{dcases} { c \over ( c - 1 ) \cos \chi + 1 } , & \text{if } \chi \leq 90^\circ \\ { c \over ( c - 1 ) \cos \chi - 1 } + 2 \, c \, e^{ x - x \sin \chi } \sqrt{ \sin \chi } , & \text{if } \chi > 90^\circ \\ \end{dcases} \\ \mathrm{with} \qquad & \\ c &= \sqrt{ \pi x \over 2 } \qquad \text{for } x > 10. \end{align*}

This is the main result. Note that the expen­sive case for when \chi > 90^\circ only needs to be con­sid­ered when look­ing down­wards and the view is not already obscured by the ter­rain, which only applies to a small region below the hori­zon. A com­par­i­son between the exact and the approx­i­mate ver­sions of \operatorname{Ch}(x,\chi) is shown below:

Epilog: A literature review

At this point I want to stress how every real-time algo­rithm for atmos­pher­ic scat­ter­ing must have an effi­cient way to eval­u­ate \operatorname{Ch}(x,\chi), even when this func­tion is not explic­it­ly iden­ti­fied as such. Here is a small lit­er­a­ture sur­vey to illus­trate that point:

Author(s) Type Assump­tions appears in (dis­guised as)
Nishi­ta et al. ’93 [4] table Par­al­lel sun­rays Sec­tion 4.3.4
Preetham et al. ’99 [5] ana­lyt­ic Earth­bound observ­er Appen­dix A.1
Nielsen ’03 [6] ana­lyt­ic Hard­cod­ed use case Sec­tion 5.4.1
O’Neil ’04 [7] table Source­code
O’Neil ’05 [8] ana­lyt­ic Con­stant \operatorname{Ch}_{\parallel} Sec­tion 4.2
Schafhitzel et al. ’07 [9] table Sec­tion 4.1
Brune­ton, Neyret ’08 [10] table Sec­tion 4 as \mathbb{T}(x,v)
Schüler ’12 [11] ana­lyt­ic List­ing 1
Yusov ’14 [12] ana­lyt­ic Sec­tion 2.4.2 as \operatorname{T}(A\rightarrow B)


This was the first post in a series. I orig­i­nal­ly didn’t plan to do a series but I decid­ed to split when the word count went over 4000, espe­cial­ly with the source list­ings. The next post will explain an over­flow-free imple­men­ta­tion of the chap­man func­tion in shad­er code, how trans­mit­tances can be cal­cu­lat­ed along arbi­trary rays and then guide through the pro­ce­dure of inte­grat­ing the atmos­pher­ic in-scat­ter along a view­ing ray.

[1] Naty Hoff­man and Arcot J Preetham (2002): “Ren­der­ing Out­door Light Scat­ter­ing in Real Time”…

[2] Syd­ney Chap­man (1931): “The absorp­tion and dis­so­cia­tive or ion­iz­ing effect of mono­chro­mat­ic radi­a­tion in an atmos­phere on a rotat­ing earth.” Pro­ceed­ings of the Phys­i­cal Soci­ety 43:1, 26.…

[3] Miroslav Koci­f­aj (1996): “Opti­cal Air Mass and Refrac­tion in a Rayleigh Atmos­phere.” Con­tri­bu­tions of the Astro­nom­i­cal Obser­va­to­ry Skalnate Ple­so 26:1, 23–30.…

[4] Tomoyu­ki Nishi­ta, Takao Sir­ai, Kat­su­mi Tadamu­ra and Eihachi­ro Naka­mae (1993): “Dis­play of the earth tak­ing into account atmos­pher­ic scat­ter­ing” Pro­ceed­ings of the 20th annu­al con­fer­ence on Com­put­er graph­ics and inter­ac­tive tech­niques, SIGGRAPH ’93, pp. 175–182.…

[5] Arcot J Preetham, Peter Shirley and Bri­an Smits (1999): “A prac­ti­cal ana­lyt­ic mod­el for day­light” Pro­ceed­ings of the 26th annu­al con­fer­ence on Com­put­er graph­ics and inter­ac­tive tech­niques, SIGGRAPH ’99, pp. 91–100…

[6] Ralf Stock­holm Nielsen (2003): “Real Time Ren­der­ing of Atmos­pher­ic Scat­ter­ing Effects for Flight Sim­u­la­tors” Master’s the­sis, Infor­mat­ics and Math­e­mat­i­cal Mod­el­ling, Tech­ni­cal Uni­ver­si­ty of Den­mark…

[7] Sean O’Neil (2004): “Real-Time Atmos­pher­ic Scat­ter­ing”…

[8] Sean O’Neil (2005): “Accu­rate Atmos­pher­ic Scat­ter­ing” GPU Gems 2, Pear­son Addi­son Wes­ley Prof, pp. 253–268…

[9] Tobias Schafhitzel, Mar­tin Falk and Thomas Ertl (2007): “Real-Time Ren­der­ing of Plan­ets with Atmos­pheres” Jour­nal of WSCG 15, pp. 91–98…

[10] Éric Brune­ton and Fab­rice Neyret (2008): “Pre­com­put­ed Atmos­pher­ic Scat­ter­ing” Pro­ceed­ings of the 19th Euro­graph­ics Sym­po­sium on Ren­der­ing, Com­put. Graph. Forum 27:4, pp. 1079–1086…

[11] Chris­t­ian Schüler (2012): “An Approx­i­ma­tion to the Chap­man Graz­ing-Inci­dence Func­tion for Atmos­pher­ic Scat­ter­ing” GPU Pro 3, A K Peters, pp. 105–118.…

[12] Egor Yusov (2014): “High Per­for­mance Out­door Light Scat­ter­ing using Epipo­lar Sam­pling” GPU Pro 5, A K Peters, pp. 101–126…

2 thoughts on “Followup to Atmospheric Scattering—Part 1

  1. Pingback: Followup to GPU Pro3 Atmospheric Scattering—Part 1 via @KostasAAA | hanecci's blog : はねっちブログ

Leave a Reply

Your email address will not be published.