Sky Simulation

For my game, I would like to have a realistic sky in the background. I have a preference for a procedurally-generated sky over a skybox. But then, I need to learn about different ways of modeling the sky.

I first considered reusing the computations in the sky example from Three.js’s documentation, as it is liberally licensed. However, I decided I’d like to understand more of the theory behind it and use my own implementation instead of copying the Three.js’s example’s code.

The Three.js example’s implementation references the Preetham model. Introduced in 1999, the Preetham model seems to have been the baseline standard for a procedural sky for a long time. The original paper does its best to explain the scientific terminology to graphics practitioners, but there is still enough left unsaid that I found the paper hard to interpret in isolation. The most useful reference I found providing necessary background knowledge is Timothy Kol’s technical report. Kol also implements and compares the newer 2012 Hošek–Wilkie model. The remainder of this post mostly distills information from these three papers.

Both the Preetham model and the Hošek–Wilkie model are based on the Perez model. The Perez model is a single formula that produces the relative luminosity of some element of the sky based on the direction and some constants.

Y_{\mathit{rel}} = \mathcal F(\theta, \gamma) = (1 + Ae^{B / \cos\theta})(1 + Ce^{D\gamma} + E\cos^2 \gamma)

Here 𝜃 is the “latitude” (angle down from the zenith, aka the top of the sky) and 𝛾 is the angle between the viewing direction and the sun (note: great circle distance-style, not angle at equator aka “longitude”). The constants A, B, C, D, and E have interpretations (you can read about it in Kol’s paper), but you don’t have to understand them and we’ll come back to this later. This alone doesn’t get you a color though – a couple more steps are necessary:

  • First, you need to run this three different times with different sets of constants – so from AY, BY, CY, DY, and EY, you compute YY; from Ax, Bx, Cx, Dx, and Ex, you compute Yx; and from Ay (careful! Y and y are different!), By, Cy, Dy, and Ey, you compute Yy.
  • These Y values are relative luminances that need to be turned into absolute luminance. It’s not described this way in any of the papers, but let’s say that’s multiplying it by another constant F. So we have Y = FYY.
  • Now we have three values Yx, Yy, and Y. These are colors in the CIE xyY color space. These need to be converted to sRGB. I’ll omit details, but basically you first need to convert to CIE XYZ, then to linear RGB, then to sRGB.

In the Hošek–Wilkie model, the Perez formula is modified to this:

\textstyle\mathcal F(\theta, \gamma) = (1 + Ae^{B/(0.01 + \cos\theta)})(C + De^{E\gamma} + F\cos^2 \gamma + G\,(\frac{1 + \cos^2 \gamma}{(1 + H^2 - 2H\cos\gamma)^{3/2}}) + I \cos^{1/2} \theta)

You can see an intuitive diagram explaining the reason behind the changes on slide 12 of Hošek and Wilkie’s SIGGRAPH presentation. Hošek–Wilkie’s model has many more constants – 9 instead of 5 (or 10 instead of 6, if you count what we earlier called F in Y = FYY).

That’s really all you need inside of your shader! Find the angles, evaluate the formulas (three times with different sets of constants) with those angles, convert relative to absolute luminance, and convert from CIE xyY to sRGB. Bingo, right?

There’s one major piece missing – how do you get those constants? Aye, there’s the rub. This is the other major difference between the Preetham model and the Hošek–Wilkie model. But not to get ahead of ourselves – there’s a reason I’m burying the lede here. If you do not care about a dynamic time of day, atmospheric conditions, or anything like that – that is, if your sky is static – then you don’t need to write code for this. You’ll need to do some work here to compute these constants ahead of time, but then you can plop those constants right down into your shader, hard-coded. I won’t frown. If you don’t need a dynamic sky, don’t pay for the complexity! YAGNI.

There’s also a happy medium point too – if your sky isn’t completely static in that the time of day changes but nothing else changes, there’s one “constant” you’ll have to write code for to recompute, but all the others stay the same. (Specifically, the “F” in the Y = FYY mentioned earlier will change, but constants A–E in Preetham and A–I in Hošek–Wilkie remain constant.)

Getting the constants

Both models have an input parameter of turbidity. Think of it as a scale from completely clear to completely overcast and foggy. Excerpted from Figure 3 in Preetham’s paper (itself derived from a paper by E.J. McCartney, see paper for details; some numbers eyeballed from plot):

Turbidity Meteorological range Characterization
1 256 km Pure air
1.2 128 km Exceptionally clear
1.8 48 km Very clear
3 12 km Clear
7 5 km Light haze
20 2.5 km Haze
40 1 km Thin fog

This parameter is denoted T. In the Preetham model (ignoring the F I added earlier), that’s the only parameter. Each constant (by which I mean e.g. AY) is a linear function of T, shown in the paper as a set of matrix expressions. Note that the Preetham model is only formulated for T ∈ [2,6].

Now about that F. There’s a separate model for determining the luminance of the zenith, briefly mentioned in section A.2 of the Preetham paper:

Y_z = (4.0453T - 4.9710)\tan((\tfrac49 - \tfrac T{120})\cdot(\pi-2\theta_s)) - 0.2155T + 2.4192

Where 𝜃s is the sun angle. Then F_Y = Y_z / \mathcal F(0, \theta_s). Fx and Fy can be computed similarly given xz and yz. Those are computed as a polynomial over T and 𝜃s; the factors are also in section A.2.

In the Hošek–Wilkie model, the range of T has been extended to T ∈ [1,10], and there is an additional parameter of ground albedo. This reflects (haha, get it) how the appearance of the sky may differ depending on the color of the ground – for example, in a snowy area, the snow reflects much of the light back into the sky, changing its color slightly (see Figure 10 in Hošek–Wilkie’s paper for an example).

The computation of constants is more complicated in Hošek–Wilkie. For each integer turbidity value and for albedo 0 and 1, there is a Bézier function for each constant. This Bézier function is evaluated. For non-integer turbidity and albedo, the values for adjacent integers are computed and the results are linearly interpolated.

Wait, where’s the sun?

Still with me? You’ve tried all this, and your sky looks great – but there seems to be something missing. It’s like we’ve gone to great lengths to draw an iconic scene of Paris from the Champ de Mars, with all the park’s trees and straight ahead, nothing but the Jardins du Trocadéro and Palais de Chaillot. Wait, shouldn’t the Eiffel Tower be there?

Up ’till now, we’ve been talking about modeling the sky and how the sky looks when illuminated by the sun. But we haven’t actually modeled what the sun itself should look like!

Section 3.1 of the Preetham paper briefly describes modeling sunlight. I think it should be possible to use this to include a sun disc with an appropriate color, but none of the figures in the paper actually show the sun, so I’m not sure how realistic it is.

Hošek and Wilkie have a follow-up 2013 paper on this, Adding a Solar-Radiance Function to the Hošek–Wilkie Skylight Model. I haven’t read this paper yet, but it seems complicated enough that I won’t be covering it in this post. However, do note that code implementing the solar radiance function is available in the code tarball provided at the authors’ web site.

Limitations and newer work

Neither of these models deal with night time. And the Hošek–Wilkie model explicitly excludes dusk – that is, 𝜃s < 10°. One commercial implementation explains that they fall back to Preetham for dusk, blending at the transition, and with a proprietary model for night illumination.

The same lab bringing us the Hošek–Wilkie model has published two further papers, and now calls it the Prague Sky Model:

These newer models better support dusk conditions. However, they start to get more complicated, and the data files start to get quite large. As an aside, it seems like the data files for the 2021 model have disappeared. But it seems like at least for the ground-level subset, the file format has stayed the same in the 2021 and 2022 papers, so the 2022 data can be used with the 2021 code if desired.

The newer Prague Sky Models seem to be mostly targeted at non-real-time applications. For real-time applications, they mention Hillaire 2020, which is used for sky simulation in Unreal Engine. Bruneton 2017 also seems commonly referenced, including by Hillaire 2020.

I have not looked into these papers in detail. I may later. For now, consider this post to be living in 2012.

Tags:

Comments are closed.