从军北征
[唐代][李益]
天山雪后海风寒, 横笛遍吹行路难。
碛里征人三十万, 一时回首月中看。
GPU Gems 2
GPU Gems 2 is now available, right here, online. You can purchase a beautifully printed version of this book, and others in the series, at a 30% discount courtesy of InformIT and Addison-Wesley.
Please visit our Recent Documents page to see all the latest whitepapers and conference presentations that can help you with your projects.
Chapter 16. Accurate Atmospheric Scattering
Sean O'Neil
16.1 Introduction
Generating realistic atmospheric scattering for computer graphics has always been a difficult problem, but it is very important for rendering realistic outdoor environments. The equations that describe atmospheric scattering are so complex that entire books have been dedicated to the subject. Computer graphics models generally use simplified equations, and very few of them run at interactive frame rates. This chapter explains how to implement a real-time atmospheric scattering algorithm entirely on the GPU using the methods described in Nishita et al. 1993. Figure 16-1 shows screenshots from the scattering demo included on this book's CD.
本文将使用Nishita在1993年提出的大气散射算法实现一种实时交互的大气散射表现，效果如下图。
Figure 16-1 Screenshots from the Scattering Demo
Many atmospheric scattering models assume that the camera is always on or very close to the ground. This makes it easier to assume that the atmosphere has a constant density at all altitudes, which simplifies the scattering equations in Nishita et al. 1993 tremendously. Refer to Hoffman and Preetham 2002 for an explanation of how to implement these simplified equations in a GPU shader. This implementation produces an attractive scattering effect that is very fast on DirectX 8.0 shaders. Unfortunately, it doesn't always produce very accurate results, and it doesn't work well for a flight or space simulator, in which the camera can be located in space or very high above the ground.
很多游戏中实现的大气散射都是基于相机离地不高的假设来实现的，在这个假设下，大气的密度将不再随着海拔高度而变化，从而可以对散射公式做出极大的简化，不过这种方式对于拥有飞行视角的游戏以及太空视角的游戏而言，就会存在问题了。
This chapter explains how to implement the full equations from Nishita et al. 1993 in a GPU shader that runs at interactive frame rates. These equations model the atmosphere more accurately, with density falling off exponentially as altitude increases. O'Neil 2004 describes a similar algorithm that ran on the CPU, but that algorithm was too CPU-intensive. It was based on a precalculated 2D lookup table with four floating-point channels. Calculating the color for each vertex required several lookups into the table, with extra calculations around each lookup. At the time that article was written, no GPU could support such operations in a shader in one pass.
O'Neil曾在2004年给过一个用CPU实现的大气散射算法，这个算法通过LUT来实现各个顶点的color的计算，不过这种算法在一次color计算中需要用到多次LUT操作，另外还伴随有其他计算消耗，对于CPU来说是非常不友好的，不过在04年，还不能通过GPU来完成相关计算操作（这也是为什么当时使用CPU来计算的原因），因此将这个算法放到现在，通过GPU来实现的话，性能将会有进一步的提升空间。本文将给出一种用GPU来实现的大气散射算法，这个算法在保证显示质量的前提下去掉了LUT操作，所有的shader操作只需要Shader Model 2.0即可完成。
In this chapter, we eliminate the lookup table without sacrificing image quality, allowing us to implement the entire algorithm in a GPU shader. These shaders are small and fast enough to run in real time on most GPUs that support DirectX Shader Model 2.0.
16.2 Solving the Scattering Equations
The scattering equations have nested integrals that are impossible to solve analytically; fortunately, it's easy to numerically compute the value of an integral with techniques such as the trapezoid rule. Approximating an integral in this manner boils down to a weighted sum calculated in a loop. Imagine a line segment on a graph: break up the segment into n sample segments and evaluate the integrand at the center point of each sample segment. Multiply each result by the length of the sample segment and add them all up. Taking more samples makes the result more accurate, but it also makes the integral more expensive to calculate.
散射公式由多个嵌套的积分运算组成，想要实时完成这个公式的计算会需要很高的消耗，不过积分运算可以通过诸如不规则四边形方法将之简化成求和运算，比如对于一条线段上的积分，可以通过将这条线段分割成多个segment，之后将对应segment端点的数值乘以segment长度并求和来拟合积分运算结果，通过增加segment数目可以提升精度，但是同时也会导致消耗增加。
In our case, the line segment is a ray from the camera through the atmosphere to a vertex. This vertex can be part of the terrain, part of the sky dome, part of a cloud, or even part of an object in space such as the moon. If the ray passes through the atmosphere to get to the vertex, scattering needs to be calculated. Every ray should have two points defined that mark where the ray starts passing through the atmosphere and where it stops passing through the atmosphere. We'll call these points A and B, and they are shown in Figure 16-2. When the camera is inside the atmosphere, A is the camera's position. When the vertex is inside the atmosphere, B is the vertex's position. When either point is in space, we perform a sphere-intersection check to find out where the ray intersects the outer atmosphere, and then we make the intersection point A or B.
从相机到屏幕像素对应的采样点的这一条路径中，只要有穿过大气的部分，就需要考虑其中的散射，假设这条线段与大气层的交点分别为A,B，那么示意图可以如下所示：
Figure 16-2 The Geometry of Atmospheric Scattering
Now we have a line segment defined from point A to point B, and we want to approximate the integral that describes the atmospheric scattering across its length. For now let's take five sample positions and name their points P 1 through P 5. Each point P 1 through P 5 represents a point in the atmosphere at which light scatters; light comes into the atmosphere from the sun, scatters at that point, and is reflected toward the camera. Consider the point P 5, for example. Sunlight goes directly from the sun to P 5 in a straight line. Along that line the atmosphere scatters some of the light away from P 5. At P 5, some of this light is scattered directly toward the camera. As the light from P 5 travels to the camera, some of it gets scattered away again.
在这个示意图中，给出的散射范围线段上指定了五个采样点，分别用P1~P5标识，相机所捕捉到的散射光强就由五个点的散射光强累加而来，注意，这里给的散射是从阳光直射到采样点，之后经由采样点反射的inscattering。
16.2.1 Rayleigh Scattering vs. Mie Scattering
Another important detail is related to how the light scattering at the point P is modeled. Different particles in the atmosphere scatter light in different ways. The two most common forms of scattering in the atmosphere are Rayleigh scattering and Mie scattering. Rayleigh scattering is caused by small molecules in the air, and it scatters light more heavily at the shorter wavelengths (blue first, then green, and then red). The sky is blue because the blue light bounces all over the place, and ultimately reaches your eyes from every direction. The sun's light turns yellow/orange/red at sunset because as light travels far through the atmosphere, almost all of the blue and much of the green light is scattered away before it reaches you, leaving just the reddish colors.
在散射计算之前，先要选定大气散射的模型。总的来说，目前主流的散射有两种：瑞利散射与米氏散射。其中瑞利散射对应的是大气中的小尺寸粒子molecule，而米氏散射对应的则是尺寸稍大的粒子aerosols比如灰尘以及雾霾粒子等。瑞利散射的散射系数与光波波长的四次方成反比，因此其对于小波长的光波（比如蓝光）的散射强度更大，这就是为什么太阳光经由大气散射后在天空上呈现蓝色的主要原因，也是傍晚天空呈现橘红色的原因（傍晚阳光经行路径长，小波长光波都已经散射掉了）。米氏散射的散射系数与波长关系较小，在阴天的时候，天空呈现灰色，且太阳周边出现一圈光晕就是米氏散射的作用，此外由于冰水粒子的散射而导致的彩虹也是米氏散射的作用。
Mie scattering is caused by larger particles in the air called aerosols (such as dust and pollution), and it tends to scatter all wavelengths of light equally. On a hazy day, Mie scattering causes the sky to look a bit gray and causes the sun to have a large white halo around it. Mie scattering can also be used to simulate light scattered from small particles of water and ice in the air, to produce effects like rainbows, but that is beyond the scope of this chapter. (Refer to Brewer 2004 for more information.)
16.2.2 The Phase Function
The phase function describes how much light is scattered toward the direction of the camera based on the angle (the angle between the two green rays in Figure 16-2) and a constant g that affects the symmetry of the scattering. There are many different versions of the phase function. This one is an adaptation of the Henyey-Greenstein function used in Nishita et al. 1993.
相函数描述的是光线方向与散射方向夹角对于散射强度的影响，有很多不同的实现版本，这里给出的是Nishita推导出来的Henyey-Greenstein的改编版本，其中g是一个影响散射对称性的常量。
Rayleigh scattering can be approximated by setting g to 0, which greatly simplifies the equation and makes it symmetrical for positive and negative angles. Negative values of g scatter more light in a forward direction, and positive values of g scatter more light back toward the light source. For Mie aerosol scattering, g is usually set between -0.75 and -0.999. Never set g to 1 or -1, as it makes the equation reduce to 0.
瑞利散射可以通过将g设置为0来实现简化模拟，当g设置为0的时候，就表明散射在正方向跟反方向的强度是对称的，正方向指的是沿着光照前进的方向，而反方向指的是背离光线前进的方向。米氏散射通常会将g设置为-0.75或者-0.999，千万不要将g设置为1或者-1，这会导致整个相函数的值为0.
16.2.3 The Out-Scattering Equation
The out-scattering equation is the inner integral mentioned earlier. The integral part determines the "optical depth," or the average atmospheric density across the ray from point Pa to point Pb multiplied by the length of the ray. This is also called the "optical length" or "optical thickness." Think of it as a weighting factor based on how many air particles are in the path of the light along the ray. The rest of the equation is made up of constants, and they determine how much of the light those particles scatter away from the ray.
outscattering计算公式给出如上，公式中的积分项指的是光线传播路径上的大气浓度（其中的exp指的是对应位置的大气浓度），这一项有一个专有名词用于指代：光学厚度或者光学长度。用于表明光线传播过程中的损耗或者衰减，而其余部分则是一个常量，用于修正大气散射的强度。
To compute the value of this integral, the ray from Pa to Pb will be broken up into segments and the exponential term will be evaluated at each sample point. The variable h is the height of the sample point. In my implementation, the height is scaled so that 0 represents sea level and 1 is at the top of the atmosphere. In theory, the atmosphere has no fixed top, but for practical purposes, we have to choose some height at which to render the sky dome. H 0 is the scale height, which is the height at which the atmosphere's average density is found. My implementation uses 0.25, so the average density is found 25 percent of the way up from the ground to the sky dome.
如前所述，为了简化计算，积分运算会被拆解成多个指数的求和运算。积分公式中的h指的是对应采样点的海拔（指数高度雾），在本文的视线中，h的取值范围为[0,1]，0指的是海平面，而1指的是大气的顶层的高度。理论上来说，大气范围是上不封顶的，不过出实际计算需要，这里需要指定一个用于进行天空渲染的高度，H0指的是平均大气浓度所对应的高度，本文中这一项取值为0.25，表明从海平面到天穹垂直路径中25%的位置的浓度为平均的大气浓度。
The constant is the wavelength (or color) of light and K() is the scattering constant, which is dependent on the density of the atmosphere at sea level. Rayleigh and Mie scattering each have their own scattering constants, including the scale height (H 0). Rayleigh and Mie scattering also differ in how they depend on wavelength. The Rayleigh scattering constant is usually divided by 4. In most computer graphics models, the Mie scattering constant is not dependent on wavelength, but at least one implementation divides it by 0.84. Wherever this equation depends on wavelength, it must be solved separately for each of the three color channels and separately for each type of scattering.
积分之外的常量K指的是散射常量，这个常量跟海平面处的大气密度有关。瑞利散射跟米氏散射的散射常量是不一样的，且锐利散射的常量跟波长有关，而米氏散射则与之无关。瑞利散射的常量通常要除以4，米氏散射则除以0.84.对于跟波长相关的散射而言，这个公式就需要写成RGB版本的。
16.2.4 The In-Scattering Equation
The in-scattering equation describes how much light is added to a ray through the atmosphere due to light scattering from the sun. For each point P along the ray from Pa to Pb , PPc is the ray from the point to the sun and PPa is the ray from the sample point to the camera. The out-scattering function determines how much light is scattered away along the two green rays in Figure 16-2. The remaining light is scaled by the phase function, the scattering constant, and the intensity of the sunlight, Is (). The sunlight intensity does not have to be dependent on wavelength, but this is where you would apply the color if you wanted to create an alien world revolving around a purple star.
inscattering公式给出如上，对于每个从Pa到Pb的点P而言，PPc指的是P到太阳的射线，而PPa指的是从P到相机的射线。积分项计算完成后还需要乘上相函数F，散射常量K以及太阳光强I（这里的太阳光可以不需要指定颜色，表明太阳光是白色的，不过如果需要设置阳光颜色，那么只需要修正I即可）
16.2.5 The Surface-Scattering Equation
To scatter light reflected from a surface, such as the surface of a planet, you must take into account the fact that some of the reflected light will be scattered away on its way to the camera. In addition, extra light is scattered in from the atmosphere. Ie () is the amount of light emitted or reflected from a surface, and it is attenuated by an out-scattering factor. The sky is not a surface that can reflect or emit light, so only Iv () is needed to render the sky. Determining how much light is reflected or emitted by a surface is application-specific, but for reflected sunlight, you need to account for the out-scattering that takes place before the sunlight strikes the surface (that is, Is () x exp(-t(Pc Pb ,))), and use that as the color of the light when determining how much light the surface reflects.
上面给出考虑了inscattering跟outscattering两种散射后的计算公式，其中Iv表示太阳的inscattering光强，后面的部分给出的是物体反射光强考虑outscattering导致的传播衰减后的反射光强。对于天空而言，第二部分是0，因此可以直接根据Iv计算，此外，对于物体表面反射的太阳光光强，在这部分光强抵达表面之前，也需要考虑传播过程中的衰减。
16.3 Making It Real-Time
Let's find out how poorly these equations will perform if they're implemented as explained in the preceding section, with five sample points for the in-scattering equation and five sample points for each of the integrals to compute the out-scattering equations. This gives 5 x (5 + 5) samples at which to evaluate the functions for each vertex. We also have two types of scattering, Rayleigh and Mie, and we have to calculate each one for the different wavelengths of each of the three color channels (RGB). So now we're up to approximately 2 x 3 x 5 x (5 + 5), or 300 computations per vertex. To make matters worse, the visual quality suffers noticeably when using only five samples for each integral. O'Neil 2004 used 50 samples for the inner integrals and five for the outer integral, which pushes the number of calculations up to 3,000 per vertex!
按照前面所说的，每个积分只使用五个采样点计算，把各种简化操作带上，每个顶点也要面临接近300个计算单位的消耗，而且在这些简化条件的作用下，其输出的效果也难以令人满意。
We don't need to go any further to know that this will not run very fast. Nishita et al. 1993 used a precalculated 2D lookup table to cut the number of calculations in half, but it still won't run in real time. Their lookup table took advantage of the fact that the sun is so far away that its rays can be considered parallel. (This idea is reflected in the in-scattering equation, in Section 16.2.4.) This makes it possible to calculate a lookup table that contains the amount of out-scattering for rays going from the sun to any point in the atmosphere. This table replaces one of the out-scattering integrals with a lookup table whose variables are altitude and angle to the sun. Because the rays to the camera are not parallel, the out-scattering integral for camera rays still had to be solved at runtime.
因此不需要再进一步考虑就知道这种实现方式不可行，此外，即使用上Nishita的LUT方法，也不过是将计算量减掉一半，还是无法满足需求。
In O'Neil 2004, I proposed a better 2D lookup table that allows us to avoid both out-scattering integrals. The first dimension, x, takes a sample point at a specific altitude above the planet, with 0.0 being on the ground and 1.0 being at the top of the atmosphere. The second dimension, y, represents a vertical angle, with 0.0 being straight up and 1.0 being straight down. At each (x, y) pair in the table, a ray is fired from a point at altitude x to the top of the atmosphere along angle y. The lookup table had four channels, two reserved for Rayleigh scattering and two reserved for Mie scattering. One channel for each simply contained the atmospheric density at that altitude, or exp(-h/H 0). The other channel contained the optical depth of the ray just described. Because the lookup table was precomputed, I chose to use 50 samples to approximate the optical depth integral, resulting in very good accuracy.
前面提到O'Neil给出的CPU LUT方法计算消耗比较小，本文准备按照这种思路来设计算法。这里的LUT二维贴图中横坐标x表示的是采样点相对于海平面的高度，0表示海平面，1表示天穹高度；纵坐标y表示采样点所在的射线相对于垂直上方向的夹角，0表示垂直向上，1表示垂直向下。根据这两个维度就可以从LUT中读取一个float4，其中前面两个通道对应的是瑞利散射的数据，后面两个通道对应的是米氏散射的数据。而在某种散射的两个通道的数据中，第一个通道对应的是当前采样点的密度：exp(-h/H0)；而第二个通道对应的是光学厚度。采用这种方式对复杂的公式项进行提前预计算可以降低很多计算消耗，可以使得即使每个积分计算采用50个采样点，依然能够得到实时的帧率。
As with the lookup table proposed in Nishita et al. 1993, we can get the optical depth for the ray to the sun from any sample point in the atmosphere. All we need is the height of the sample point (x) and the angle from vertical to the sun (y), and we look up (x, y) in the table. This eliminates the need to calculate one of the out-scattering integrals. In addition, the optical depth for the ray to the camera can be figured out in the same way, right? Well, almost. It works the same way when the camera is in space, but not when the camera is in the atmosphere. That's because the sample rays used in the lookup table go from some point at height x all the way to the top of the atmosphere. They don't stop at some point in the middle of the atmosphere, as they would need to when the camera is inside the atmosphere.
借助这个LUT，只需要提供正确的横纵坐标，即可从表中读取对应的光学厚度以及大气密度，而不需要进行任何的积分运算，事实上，这个说法是有条件的，那就是相机需要处于太空中，而不能位于大气层中，因为LUT中存储的数据是从大气中某一点到大气层顶部某点的直连线段的数据，而没有从大气层某点到大气层另一点的直连线段的数据，对于这种情况就需要进行两次LUT查询，并将之结果相减即可。
Fortunately, the solution to this is very simple. First we do a lookup from sample point P to the camera to get the optical depth of the ray passing through the camera to the top of the atmosphere. Then we do a second lookup for the same ray, but starting at the camera instead of starting at P. This will give us the optical depth for the part of the ray that we don't want, and we can subtract it from the result of the first lookup. Examine the rays starting from the ground vertex (B 1) in Figure 16-3 for a graphical representation of this.
Figure 16-3 Problems with the Improved Lookup Table
There's only one problem left. When a vertex is above the camera in the atmosphere, the ray from sample point P through the camera can pass through the planet itself. The height variable is not expected to go so far negative, and it can cause the numbers in the lookup table to go extremely high, losing precision and sometimes even encountering a floating-point overflow. The way to avoid this problem is to reverse the direction of the rays when the vertex is above the camera. Examine the rays passing through the sky vertex (B 2) in Figure 16-3 for a graphical representation of this.
这里还有一个问题需要解决，那就是如果采样点位置高于相机位置，那么从采样点向着相机引线就会跟地平面相交，按照这种方式来进行LUT采样得到的结果可能会存在溢出以及精度受损等异常，为了避免这个问题，可以考虑将连线方向进行翻转，如上图所示，将P2到相机A点的连线，改成从P2到B2的方向。
So now we've reduced 3,000 calculations per vertex to something closer to 2 x 3 x 5 x (1 + 1), or 60 (assuming five samples for the out-scattering integral for the eye ray). Implemented in software on an Athlon 2500+ with an inefficient brute-force rendering method, I was able to get this method to run between 50 and 100 frames per second.
经过上述处理后，原来可能需要3000个单元计算的消耗就被降低到60个了，在Athlon 2500+机器上可以跑到50~100帧每秒。
16.4 Squeezing It into a Shader
At this point, I felt that this algorithm had been squeezed about as much as it could, but I knew that I had to squeeze it even smaller to fit it into a shader. I didn't want this algorithm to require Shader Model 3.0, so having lookup tables in textures used by the vertex shaders wasn't possible. I decided to take a different approach and started to mathematically analyze the results of the lookup table. Even though I knew I couldn't come up with a way to simplify the integral equations, I had hoped that I might be able to simulate the results closely enough with a completely different equation.
到此为止，这个算法看起来像是无可优化了，不过由于这里不打算使用SM 3.0，因此无法在VS中进行LUT采样，因此考虑能不能给出一种方式对这个算法的实现进行优化。
16.4.1 Eliminating One Dimension
I started by plotting the results of the lookup table on a graph. I plotted height from 0.0 to 1.0 along the x axis and the lookup table result (optical depth) along the y axis. For various angles sampled from 0 to 1, I plotted a separate line on the graph. I noticed right away that each line dropped off exponentially as x went from 0 to 1, but the scale of each line on the graph varied dramatically. This makes sense, because as the angle of the ray goes from straight up to straight down, the length of the ray increases dramatically.
这里先尝试保持一维不变，只考虑另一维的变化所对应的数据的变化。即沿着x方向考虑数据的变化，以及沿着y方向考虑数据的变化。得到的结论是，对于y保持不变的情况下，沿着x方向，数据的变化规律是按照指数形式进行衰减的（x方向代表的是高度，大气密度跟光学厚度随着高度指数衰减，很好理解），不过对于不同的y，衰减的缩放系数是不同的。
To make it easier to compare the shapes of the curves of each line, I decided to normalize them. I took the optical depth value at x = 0 (or height = 0) for each line and divided all of the values on that line by that value. This scaled all lines on the graph to start at (x = 0, y = 1) and work their way down toward (x = 1, y = 0). To my surprise, almost all of the normalized lines fell right on top of each other on the graph! The curve was almost exactly the same for every one, and that curve turned out to be exp(-4x).
为了能够直观的对不同的y指的x方向变化曲线进行比对，这里考虑对每条x方向的变化曲线进行归一化，也就是将每条曲线的数据的值除以曲线的起始端点（x = 0，也就是高度为0）的数值，令人震颤的事情发生了，这里观察到所有的曲线的形状都接近一致（大概规律服从exp(-4x)），且每一条曲线都在上一条曲线的上方（是指在y方向上有个偏移吗？）
This makes some sense, because the optical depth equation is the integral of exp(-h/H 0). I chose H 0 to be 0.25, so exp(-4h) is a common factor. It is still a bit puzzling, however, as the h inside the integral is not the same as the height x outside the integral, which is only the height at the start of the ray. The h value inside the integral does not vary linearly, and it has more to do with how it passes through a spherical space than with the starting height. There is some variation in the lines, and the variation increases as the angle increases. The variation gets worse exponentially as the angle increases over 90 degrees. Because we don't care about angles that are much larger than 90 degrees (because the ray passes through the planet), exp(-4x) works very well for eliminating the x axis of the lookup table.
之所以会有exp(-4x)的规律，是因为我们前面选取了H0位0.25，那么exp(-h/H0)就等同于exp(-4h)，不过这里还是有点让人疑惑，因为从定义上来说，h并不等同于x（没太看懂这一段的论述）。实际上对于不同的y值，对应的数据曲线之间还是存在区别的，且y值越大，区别也就越明显，最大值在y=0.5（90度角的时候，大于90度的情况这里是不需要考虑的，因为其方向就是朝着地平面了，这种情况在前面说过，会用反方向的线段来替代），最终实验表明，使用exp(-4x)的公式能够很好的消除LUT中的x的维度，减少查表的消耗。
16.4.2 Eliminating the Other Dimension
Now that the x dimension (height) of the lookup table is being handled by exp(-4x), we need to eliminate the y dimension (angle). The only part of the optical depth that is not handled by exp(-4x) is the scale used to normalize the lines on the graph explained previously, which is the value of the optical depth at x = 0. So now we create a new graph by plotting the angle from 0 to 1 on the x axis and the scale of each angle on the y axis. For lack of a better name, I call this the scale function.
下面尝试对y轴数据进行模拟简化，前面说到x轴经过归一化处理后得到exp(-4x)的简化模拟公式，这里对于不同的y值，其归一化除去的数值是不同的。这里尝试在保持x不变的情况下，绘制y从0到1的变化规律，这个变化规律，我们起了个名字叫做scale function
The first time I looked at the scale function, I noticed that it started at 0.25 (the scale height) and increased on some sort of accelerating curve. Thinking that it might be exponential, I divided the scales by the scale height (to make the graph start at 1) and took the natural logarithm of the result. The result was another accelerating curve that I didn't recognize. I tried a number of curves, but nothing fit well on all parts of the curve. I ended up using graphical analysis software to find a "best fit" equation for the curve, and it came back with a polynomial equation that was not pretty but fit the values well.
最开始观察scale function的走势的时候，我们发现其起点为0.25（为什么不是0？），这个值我们称之为scale height，之后一路往上走，我们起初怀疑是指数变化曲线，不过在将之除以scale height（归一化）并求取对数后发现，并没有得到想要的线性规律，而是另一个无法识别的曲线，最终我们尝试借助曲线模拟软件来对这个曲线的走势进行模拟，得到了一个不太直观但是模拟表现良好的多项式公式。
One significant drawback to this implementation is that the scale function is dependent on the scale height and the ratio between the atmosphere's thickness and the planet's radius. If either value changes, you need to calculate a new scale function. In the demo included on this book's CD, the atmosphere's thickness (the distance from the ground to the top of the atmosphere) is 2.5 percent of the planet's radius, and the scale height is 25 percent of the atmosphere's thickness. The radius of the planet doesn't matter as long as those two values stay the same.、
scale function的一个问题在于，其规律是受scale height以及大气厚度（从地平面到天穹顶部的距离）还有星球半径所影响的，其中任何一个数值发生变化，都需要重新对scale function进行计算，不过对于地球而言，其半径是不变的，且这里我们可以取地球半径的2.5%作为大气层厚度，因此，这两个值就可以认为是不变的（那还有scale height呢？）。
16.5 Implementing the Scattering Shaders
Now that the problem has been solved mathematically, let's look at how the demo was implemented. The C++ code for the demo is fairly simple. The gluSphere() function is called to render both the ground and the sky dome. The front faces are reversed for the sky dome so that the inside of its sphere is rendered. It uses a simple rectangular Earth texture to make it possible to see how the ground scattering affects colors on the ground, and it uses a simple glow texture billboard to render the moon. No distinct sun is rendered, but the Mie scattering creates a glow in the sky dome that looks like the sun (only when seen through the atmosphere).
LUT已经可以通过数学公式完整表达了，下面来看一下具体实现，这里使用gluSphere来进行地球与天穹的渲染，其中天穹的渲染需要反过来（避免被裁减掉），另外使用了一个贴图公告板来绘制月亮，这里并没有直接渲染太阳，不过由于米氏散射的原因，通过大气散射可以观察到太阳形状的glow效果。
I have provided shader implementations in both Cg and GLSL on the book's CD. The ground, the sky, and objects in space each have two scattering shaders, one for when the camera is in space and one for when the camera is in the atmosphere (this avoids conditional branching in the shaders). The ground shaders can be used for the terrain, as well as for objects that are beneath the camera. The sky shaders can be used for the sky dome, as well as for objects that are above the camera. The space shaders can be used for any objects outside the atmosphere, such as the moon.
这里提供了两种散射shader，一种对应于太空中的效果，另一种对应于大气层中的效果（避免shader 分支）。
The naming convention for the shaders is "render_object From camera_position". So the SkyFromSpace shader is used to render the sky dome when the camera is in space. There is also a common shader that contains some common constants and helper functions used throughout the shaders. Let's use SkyFromSpace as an example.
16.5.1 The Vertex Shader
As you can see in Listing 16-1, SkyFromSpace.vert is a fairly complex vertex shader, but hopefully it's easy enough to follow with the comments in the code and the explanations provided here. Kr is the Rayleigh scattering constant, Km is the Mie scattering constant, and ESun is the brightness of the sun. Rayleigh scatters different wavelengths of light at different rates, and the ratio is 1/pow(wavelength, 4). Referring back to Figure 16-2, v3Start is point A from the previous examples and v3Start + fFar * v3Ray is point B. The variable v3SamplePoint goes from P 1 to Pn with each iteration of the loop.
The variable fStartOffset is actually the value of the lookup table from point A going toward the camera. Why would we need to calculate this when it's at the outer edge of the atmosphere? Because the density is not truly zero at the outer edge. The density falls off exponentially and it is close to zero, but if we do not calculate this value and use it as an offset, there may be a visible "jump" in color when the camera enters the atmosphere.
大气层边缘的大气密度并不等于0，而是按照指数形式无穷趋近于0，如果在计算的时候不考虑这一点，在相机从外太空进入大气的时候，就会发生一个颜色的跳变。
You may have noticed that the phase function is missing from this shader. The phase function depends on the angle toward the light source, and it suffers from tessellation artifacts if it is calculated per vertex. To avoid these artifacts, the phase function is implemented in the fragment shader.
在VS中没有考虑相函数的影响，是因为如果在VS中计算这一项的话，会导致结果存在较大的瑕疵（类似在VS中采样贴图导致color计算的色块一样），因此将之放到了PS中。
Example 16-1. SkyFromSpace.vert, Which Renders the Sky Dome When the Camera Is in Space
#include "Common.cg"
vertout main(float4 gl_Vertex : POSITION,
uniform float4x4 gl_ModelViewProjectionMatrix,
uniform float3 v3CameraPos, // The camera's current position
uniform float3 v3LightDir, // Direction vector to the light source
uniform float3 v3InvWavelength, // 1 / pow(wavelength, 4) for RGB
uniform float fCameraHeight, // The camera's current height
uniform float fCameraHeight2, // fCameraHeight^2
uniform float fOuterRadius, // The outer (atmosphere) radius
uniform float fOuterRadius2, // fOuterRadius^2
uniform float fInnerRadius, // The inner (planetary) radius
uniform float fInnerRadius2, // fInnerRadius^2
uniform float fKrESun, // Kr * ESun
uniform float fKmESun, // Km * ESun
uniform float fKr4PI, // Kr * 4 * PI
uniform float fKm4PI, // Km * 4 * PI
uniform float fScale, // 1 / (fOuterRadius - fInnerRadius)
uniform float fScaleOverScaleDepth) // fScale / fScaleDepth
{
// Get the ray from the camera to the vertex and its length (which
// is the far point of the ray passing through the atmosphere)
float3 v3Pos = gl_Vertex.xyz;
float3 v3Ray = v3Pos - v3CameraPos;
float fFar = length(v3Ray);
v3Ray /= fFar;
// Calculate the closest intersection of the ray with
// the outer atmosphere (point A in Figure 16-3)
float fNear = getNearIntersection(v3CameraPos, v3Ray, fCameraHeight2,
fOuterRadius2);
// Calculate the ray's start and end positions in the atmosphere,
// then calculate its scattering offset
float3 v3Start = v3CameraPos + v3Ray * fNear;
fFar -= fNear;
float fStartAngle = dot(v3Ray, v3Start) / fOuterRadius;
float fStartDepth = exp(-fInvScaleDepth);
float fStartOffset = fStartDepth * scale(fStartAngle);
// Initialize the scattering loop variables
float fSampleLength = fFar / fSamples;
float fScaledLength = fSampleLength * fScale;
float3 v3SampleRay = v3Ray * fSampleLength;
float3 v3SamplePoint = v3Start + v3SampleRay * 0.5;
// Now loop through the sample points
float3 v3FrontColor = float3(0.0, 0.0, 0.0);
for(int i=0; i<nSamples; i++) {
float fHeight = length(v3SamplePoint);
float fDepth = exp(fScaleOverScaleDepth * (fInnerRadius - fHeight));
float fLightAngle = dot(v3LightDir, v3SamplePoint) / fHeight;
float fCameraAngle = dot(v3Ray, v3SamplePoint) / fHeight;
float fScatter = (fStartOffset + fDepth * (scale(fLightAngle) Ð
scale(fCameraAngle)));
float3 v3Attenuate = exp(-fScatter *
(v3InvWavelength * fKr4PI + fKm4PI));
v3FrontColor += v3Attenuate * (fDepth * fScaledLength);
v3SamplePoint += v3SampleRay;
}
// Finally, scale the Mie and Rayleigh colors
vertout OUT;
OUT.pos = mul(gl_ModelViewProjectionMatrix, gl_Vertex);
OUT.c0.rgb = v3FrontColor * (v3InvWavelength * fKrESun);
OUT.c1.rgb = v3FrontColor * fKmESun;
OUT.t0 = v3CameraPos - v3Pos;
return OUT;
}
16.5.2 The Fragment Shader
This shader should be fairly self-explanatory. One thing that is not shown in Listing 16-2 is that I commented out the math in the getRayleighPhase() function. The Rayleigh phase function causes the blue sky to be the brightest at 0 degrees and 180 degrees and the darkest at 90 degrees. However, I feel it makes the sky look too dark around 90 degrees. In theory, this problem arises because we are not implementing multiple scattering. Single scattering only calculates how much light is scattered coming directly from the sun. Light that is scattered out of the path of a ray just vanishes. Multiple scattering attempts to figure out where that light went, and often it is redirected back toward the camera from another angle. It is a lot like radiosity lighting, and it is not considered feasible to calculate in real time. Nishita et al. 1993 uses an ambient factor to brighten up the darker areas. I feel that it looks better to just leave out the Rayleigh phase function entirely. Feel free to play with it and see what you like best.
下面给出PS的实现代码，这里需要注意的是getRayleighPhase()函数，瑞利散射相函数是导致天空呈现蓝色的原因，且其颜色亮度在0度以及180度达到最大值，在90度是最小的，这样就会导致天空显得比较阴暗，这是因为我们的计算只考虑了一次散射的原因，即只考虑从太阳直接散射过来的光照，将散射到其他方向的光照就直接认定为消失不见了，而在实际情况下，大气是存在多次散射的，即散射向其他方向的光照会以另一个方向另一个角度再次散射进入人眼，不过这就有点类似于辐射度照明（radiosity lighting）了，其计算过程过于复杂，不太适合实时渲染。Nishita通过增加一个环境光因子的方式对来加强天空阴暗部分的亮度实现对多次散射效果的伪造。
Example 16-2. SkyFromSpace.frag, Which Renders the Sky Dome When the Camera Is in Space
#include "Common.cg"
float4 main(float4 c0 : COLOR0,
float4 c1 : COLOR1,
float3 v3Direction : TEXCOORD0,
uniform float3 v3LightDirection,
uniform float g,
uniform float g2) : COLOR
{
float fCos = dot(v3LightDirection, v3Direction) / length(v3Direction);
float fCos2 = fCos * fCos;
float4 color = getRayleighPhase(fCos2) * c0 +
getMiePhase(fCos, fCos2, g, g2) * c1;
color.a = color.b;
return color;
}
16.6 Adding High-Dynamic-Range Rendering
Atmospheric scattering doesn't look very good without high-dynamic-range (HDR) rendering. This is because these equations can very easily generate images that are too bright or too dark. See Figure 16-4. It looks much better when you render the image to a floating-point buffer and scale the colors to fall within the range 0..1 using an exponential curve with an adjustable exposure constant.
在大气渲染的计算中，如果不开启HDR的话，其渲染效果会比较奇怪，因为实际的运算中经常会出现结果由于表达范围的限制而变得过曝或者过暗，如下图所示：左边开启了HDR，右边么有。
Figure 16-4 The Importance of High Dynamic Range
The exposure equation used in this demo is very simple: 1.0 – exp(-fExposure x color). The exposure constant works like the aperture of a camera or the pupils of your eyes. When there is too much light, the exposure constant is reduced to let in less light. When there is not enough light, the exposure constant is increased to let in more light. Without reasonable limits set on the exposure constant, daylight can look like moonlight and vice versa. Regardless, the relative brightness of each color is always preserved.
这里给出的曝光度公式比较简单： 1.0 – exp(-fExposure x color)。曝光度用于控制场景的整体亮度，如果想要在渲染中模拟瞳孔或者相机的成像小孔的自动曝光功能，就需要调整这个值。
The HDR rendering implemented in this demo is also very simple. An OpenGL pbuffer is created with a floating-point pixel format and render-to-texture enabled. The same rendering code as before is used, but everything is rendered to the pbuffer's rendering context. Afterward, the main rendering context is selected and a single quad is rendered to fill the screen. The floating-point buffer is selected as a texture, and a simple exponential scaling shader is used to scale the colors to the normal 0–255 range. The exposure constant is set manually and can be changed at runtime in the demo, along with several of the scattering constants used.
这里HDR的实现使用的也是非常简单的方法，从HDR到LDR的转换用的是一个指数衰减公式，衰减常量是手动设置的。
Fortunately, modern GPUs support floating-point render targets. Some critical features such as alpha blending aren't available until you move up to GeForce 6 Series GPUs, so blending the sky dome with objects rendered in space won't work for older GPUs.
16.7 Conclusion
Although it took a fair number of simplifications to be able to interactively render scenes with the model described, the demo looks pretty sharp, especially with HDR rendering. Tweaking the parameters to the system can be a source of much entertainment; I had never imagined that the sky at the horizon would be orange at high noon if our atmosphere were a little thicker, or that sunsets would look that much more incredible. Red and orange skies produce sunsets of different shades of blue, and purple skies produce green sunsets. Once I even managed to find settings that produced a rainbow sunset.
In any case, the work is not done here. There are a number of improvements that could be made to this algorithm:
关于这个算法，还有很多可以改进的地方：
Create an atmospheric-density scale function that is not hard-coded to a specific scale height and the ratio between the atmosphere thickness and the planetary radius.--给出一个不是使用硬编码的方式计算的scale function
Split up the Rayleigh and Mie scattering so that they can have different scale depths.--为瑞利散射跟米氏散射指定不同的scale depth
Provide a better way to simulate multiple scattering and start using the Rayleigh phase function again.
Allow light from other sources, such as the moon, to be scattered, as well as sunlight. This would be necessary to create a moonlight effect with a halo surrounding the moon.--考虑多次散射的影响
Atmospheric scattering also makes distant objects look larger and blurrier than they are. The moon would look much smaller to us if it were not viewed through an atmosphere. It may be possible to achieve this in the HDR pass with some sort of blurring filter that is based on optical depth.--大气散射通常会导致远处的物体显得大而模糊，这种效果可以通过在HDR中使用某种基于光学厚度的模糊算法来实现。
Change the sample points taken along the ray from linear to exponential based on altitude. Nishita et al. 1993 explains why this is beneficial.--将采样点的分布从线性修正成指数
Add other atmospheric effects, such as clouds, precipitation, lightning, and rainbows. See Harris 2001, Brewer 2004, and Dobashi et al. 2001.--添加其他大气效果，云层，闪电，彩虹以及雨雪
16.8 References
Brewer, Clint. 2004. "Rainbows and Fogbows: Adding Natural Phenomena." NVIDIA Corporation. SDK white paper. http://download.nvidia.com/developer/SDK/Individual_Samples/DEMOS/Direct3D9/src/HLSL_RainbowFogbow/docs/RainbowFogbow.pdf
Dobashi, Y., T. Yamamoto, and T. Nishita. 2001. "Efficient Rendering of Lightning Taking into Account Scattering Effects due to Clouds and Atmospheric Particles." http://nis-lab.is.s.u-tokyo.ac.jp/~nis/cdrom/pg/pg01_lightning.pdf
Harris, Mark J., and Anselmo Lastra. 2001. "Real-Time Cloud Rendering." Eurographics 2001 20(3), pp. 76–84. http://www.markmark.net/PDFs/RTClouds_HarrisEG2001.pdf
Hoffman, N., and A. J. Preetham. 2002. "Rendering Outdoor Scattering in Real Time." ATI Corporation. http://www.ati.com/developer/dx9/ATI-LightScattering.pdf
Nishita, T., T. Sirai, K. Tadamura, and E. Nakamae. 1993. "Display of the Earth Taking into Account Atmospheric Scattering." In Proceedings of SIGGRAPH 93, pp. 175–182. http://nis-lab.is.s.u-tokyo.ac.jp/~nis/cdrom/sig93_nis.pdf
O'Neil, Sean. 2004. "Real-Time Atmospheric Scattering." GameDev.net. http://www.gamedev.net/columns/hardcore/atmscattering/
