作者：i_dovelemon
日期：2020-11-25
主题：Atmospheric Scattering, Volume Scattering, Rayleigh Scattering, Mie Scattering, Single Scattering, Multiple Scattering
引言
前文 GraphicsLab 之 Atmospheric Scattering（一）讲述了基于物理的天空渲染在 Single Scattering 情况下的基础理论知识。本篇文章将主要从代码实现的角度，详细讲解如何在 Unity 中实现一个简单粗暴版本的大气散射天空。
离散化天空渲染方程
回顾之前一篇文章中，我们最后得出的用于计算天空颜色的渲染方程，如下所示：
$$I_{final}=\int_{A}^{B}{I_{sun} * T_{cp} * S(\lambda ,\theta ,h) * T_{pa}}$$(Eq 10)
简化天空渲染方程
根据之前文章的描述，最后渲染方程积分里面，有一些数据是不会发生变化的，所以可以提取出来。
最简单直观的就是，太阳光的强度 $I_{sun}$ 不会随着积分发生变化，所以可以提取出来。
根据 Eq 2 和 Eq 3 可以得到新的 Eq 11，如下所示：
$$S(\lambda ,\theta ,h)=\beta (\lambda ,h) * \gamma (\theta)$$(Eq 2)
$$\beta (\lambda ,h)=\beta (\lambda,0) * \rho (h)$$(Eq 3)
$$S(\lambda ,\theta ,h)=\beta (\lambda,0) * \rho (h) * \gamma (\theta)$$(Eq 11)
而 Eq 11 中可以看出，$\beta (\lambda,0)$ 是海拔高度 0 处的大气散射系数，这个是常数，不会发生变化，可以提取出来。
太阳光入射方向不会发生变化，所以对于 AB 射线来说，$\theta$ 角也是固定的，所以 $\gamma (\theta)$ 也可以提取出来。
自此，我们就得到了一个新的公式：
$$I_{final}=I_{sun}*\beta (\lambda,0) * \gamma (\theta)*\int_{A}^{B}{T_{cp} * \rho (h) * T_{pa}}$$(Eq 12)
积分的离散化
上面公式的最外围是一个连续积分函数 $\int_{A}^{B}{T_{cp} * \rho (h) * T_{pa}}$ 。$T_{cp}$ 和 $T_{pa}$ 里面也存在一个 $\int_{0}^{d}\beta(\lambda,h)$ 连续积分函数。
这些函数没有解析解，所以我们需要使用数值计算的方式，计算这个积分的值。
在数值计算领域这样的计算非常的多，根据精度和待求解问题积分函数的不同，可以使用不同的方法进行。我们这里通过简单的黎曼和积分的方式来求解这个值。
黎曼和积分很简单，在积分定义域范围里面，设定一个固定的步长，计算步长中心处公式的值，然后假设整个步长范围里面所有的计算出来的值都是该大小，这样就可以得到一个离散的求和形式的公式。步长越多，计算量越大，结果越精确。关于黎曼和积分的详细描述，可以参考文献 [1] 进行了解。
所以公式 $\int_{A}^{B}{T_{cp} * \rho (h) * T_{pa}}$ 就可以被离散化为如下的形式：
$$\int_{A}^{B}{T_{cp} * \rho (h) * T_{pa}}\approx \sum_{0}^{N}{T_{cp} * \rho (h) * T_{pa}*ds}$$(Eq 13)
公式 $\int_{0}^{d}\beta(\lambda,h)$ 也可以被离散为如下形式：
$$\int_{0}^{d}\beta(\lambda,h) \approx \sum_{0}^{N}{\beta(\lambda,h)*ds}$$(Eq 14)
其中 $N$ 表示计算的次数，$ds$ 为步长大小。
将 Eq 13 带入 Eq 12 ，得到可以实际计算的离散化的公式：
$$I_{final}\approx I_{sun}*\beta (\lambda,0) * \gamma (\theta)*\sum_{0}^{N}{T_{cp} * \rho (h) * T_{pa}*ds}$$(Eq 15)
将 Eq 14 带入 Eq 6，得到如下离散化的公式：
$$T \approx e^{-{\sum_{0}^{N}{\beta(\lambda,h)*ds}}}$$(Eq 16)
后处理
大气散射天空，一般作为整个场景的背景进行渲染。所以我们这里采用后处理的方式进行天空背景的渲染。
在 Unity 中添加后处理有一套固定的方式，我这里为了简单，就直接绘制了一个面片，并且在 Vertex Shader 中不进行变换，直接投影到 ClipSpace 中，从而绘制一个覆盖全屏幕的面片出来。
Vertex Shader 代码如下所示：
1 ASOutput ASVert(ASInput input)
2 {
3 ASOutput output = (ASOutput)0;
4
5 output.positionCS = input.position * 2.0f;
6 output.positionCS.z = 1.0f;
7 output.positionCS.w = 1.0f;
8 output.uv = input.uv;
9
10 return output;
11 }
有了覆盖全屏幕的面片之后，我们就可以通过对每一个像素进行大气散射的计算，从而实现天空背景的渲染。
天空像素计算
有了覆盖全屏幕的后处理之后，就可以计算每一个像素的颜色了。如下是计算每一个像素的 PixelShader 代码：
1 float4 ASFrag(ASOutput output) : SV_TARGET
2 {
3 float3 planetPos = float3(0.0f, 0.0f, 0.0f);
4 float planetRadius = 6371e3;
5 float atmosphereRadius = 6471e3;
6 uint uViewRaySampleN = 64u;
7 uint uLightRaySampleN = 4u;
8 float3 sunIntensity = _SunIntensity;
9 float3 sunDir = normalize(_SunDirection);
10
11 float3 InSky = float3(0.0f, atmosphereRadius, atmosphereRadius);
12 float3 InGround = float3(0.0f, planetRadius + 100.0f, 0.0f);
13 float3 cameraPos = InGround;
14
15 float3 cameraView = CalculateCameraVector(output.uv, _ScreenParams.xy);
16 float3 rayleighColor = RayleighAtmosphericScatteringIntegration(
17 cameraPos, cameraView,
18 planetPos, atmosphereRadius, planetRadius,
19 uViewRaySampleN, uLightRaySampleN,
20 sunIntensity, sunDir
21 );
22 float3 mieColor = MieAtmosphericScatteringIntegration(
23 cameraPos, cameraView,
24 planetPos, atmosphereRadius, planetRadius,
25 uViewRaySampleN, uLightRaySampleN,
26 sunIntensity, sunDir
27 );
28
29 float3 color = rayleighColor + mieColor;
30 return float4(color * 1.0f, 1.0f);
31 }
前面说过，我们将大气分为了两个大类，分别使用 Rayleigh Scattering 和 Mie Scattering 进行模拟。所以上面的代码也是这样处理的，分别计算 Rayleigh Scattering 和 Mie Scattering 的颜色，然后将结果相加即可。
除了这个之外，就是一些输入参数的定义，如下所示：
- planetPos 定义了地球球体的圆心所在，这里就定义为 (0，0，0)
- planetRadius 定义了地球的半径，单位 m
- atmosphereRadius 定义了大气的半径，单位 m
- sunIntensity 定义了太阳光的光照强度，即 $I_{sun}$
- sunDir 定义了光照的入射方向
- cameraPos 定义了观察者所在的位置，这里写死了，开发者可以根据需要修改
有了这些参数之后，我们需要计算一下当前观察射线的方向，即 cameraView 。这个 cameraView 就是理论中 AB 射线的方向。每一个像素都需要计算一个独立的观察射线，如下是 在 Unity 下计算 cameraView 的代码：
1 float3 CalculateCameraVector(float2 coord, float2 screen)
2 {
3 coord.y = 1.0f - coord.y;
4 float2 uv = 2.0f * (coord.xy - float2(0.5f, 0.5f));
5 return normalize(float3(uv.x, uv.y, -1.0f));
6 }
有了这些数据之后，我们就来计算当前像素的 Rayleigh Color 和 Mie Color。
RayleighAtmosphericScatteringIntegration 和 MieAtmosphericScatteringIntegration
如下是这两个函数的代码：
1 float3 RayleighAtmosphericScatteringIntegration(
2 float3 viewPos, float3 viewDir, // View
3 float3 planetPos, float atmosphereRadius, float planetRadius, // Planet and Atmosphere
4 uint viewRaySampleN, uint lightRaySampleN, // View and Light Ray sample time
5 float3 sunIntensity, float3 sunDir // Sun
6 )
7 {
8 float la = 0.0f, lb = 0.0f;
9 bool isViewAtmosphere = IntersectAtmosphere(viewPos, viewDir, planetPos, atmosphereRadius, planetRadius, la, lb);
10 if (!isViewAtmosphere)
11 {
12 // Do not view atmoshpere, there is not scattering happen
13 return float3(0.0f, 0.0f, 0.0f);
14 }
15
16 // Split intersect ray into N segment
17 float ds = (lb - la) / viewRaySampleN;
18 float st = la;
19
20 float3 totalContribution = 0.0f;
21 for (uint i = 0u; i < viewRaySampleN; i++)
22 {
23 // Current sample position
24 float tp = (st + ds * 0.5f);
25 float3 pos = viewPos + viewDir * tp;
26
27 float height = distance(pos, planetPos) - planetRadius;
28 totalContribution = totalContribution + RayleighLightContributionIntegration(
29 height, ds, tp, st, viewPos, viewDir, sunDir, lightRaySampleN, planetPos, planetRadius, atmosphereRadius);
30
31 st = st + ds;
32 }
33
34 float3 coefficient = RayleighScatteringCoefficientAtSealevel();
35 float phase = RayleighScatteringPhase(dot(viewDir, sunDir));
36 return sunIntensity * coefficient * totalContribution * phase;
37 }
38
39 float3 MieAtmosphericScatteringIntegration(
40 float3 viewPos, float3 viewDir, // View
41 float3 planetPos, float atmosphereRadius, float planetRadius, // Planet and Atmosphere
42 uint viewRaySampleN, uint lightRaySampleN, // View and Light Ray sample time
43 float3 sunIntensity, float3 sunDir // Sun
44 )
45 {
46 float la = 0.0f, lb = 0.0f;
47 bool isViewAtmosphere = IntersectAtmosphere(viewPos, viewDir, planetPos, atmosphereRadius, planetRadius, la, lb);
48 if (!isViewAtmosphere)
49 {
50 // Do not view atmoshpere, there is not scattering happen
51 return float3(0.0f, 0.0f, 0.0f);
52 }
53
54 // Split intersect ray into N segment
55 float ds = (lb - la) / viewRaySampleN;
56 float st = la;
57
58 float3 totalContribution = 0.0f;
59 for (uint i = 0u; i < viewRaySampleN; i++)
60 {
61 // Current sample position
62 float tp = (st + ds * 0.5f);
63 float3 pos = viewPos + viewDir * tp;
64
65 float height = distance(pos, planetPos) - planetRadius;
66 totalContribution = totalContribution + MieLightContributionIntegration(
67 height, ds, tp, st, viewPos, viewDir, sunDir, lightRaySampleN, planetPos, planetRadius, atmosphereRadius);
68
69 st = st + ds;
70 }
71
72 float3 coefficient = MieScatteringCoefficientAtSealevel();
73 float phase = MieScatteringPhase(dot(viewDir, sunDir));
74 return sunIntensity * coefficient * totalContribution * phase;
75 }
可以看到这两个函数的流程是一摸一样的，它们的不同也就体现在最终计算时的一些参数不同之上。所以，下面就合在一起说明了。
第一步，我们需要计算下当前观察射线是否与大气相交。这是因为，如果我们从太空中观察地球的话，有很大可能一些观察射线是不会碰撞到大气层的。当然如果是在地球表面的话，肯定是会碰撞到的。
如果碰撞到了大气，我们就需要计算实际在大气中的线段 AB。这是因为，观察射线实际上是一个无限长的射线，而被大气所影响到的只可能是射线上的一部分，即我们需要计算的 AB。如下是 IntersectAtmosphere 函数的定义：
1 //----------------------------------------------------------------------------
2 // Intersection regin
3 //----------------------------------------------------------------------------
4 bool RayIntersectSphere(
5 float3 ro, float3 rd, // Ray
6 float3 so, float sr, // Sphere
7 out float ra, out float rb // Result
8 )
9 {
10 ra = 0.0f;
11 rb = 0.0f;
12 float a = dot(rd, rd);
13 float b = 2.0f * dot(rd, ro);
14 float c = dot(ro, ro) - (sr * sr);
15 float d = (b * b) - 4.0f * a * c;
16 if (d < 0.0f) return false;
17
18 ra = max(0.0f, (-b - sqrt(d)) / 2.0f * a); // Fuck here, ra can not be negative
19 rb = (-b + sqrt(d)) / 2.0f * a;
20 if (ra > rb) return false; // Fuck here, rb must be bigger than ra
21 return true;
22 }
23
24 bool IntersectAtmosphere(
25 float3 ro, float3 rd, // Ray
26 float3 o, float ar, float pr, // Planet and atmosphere
27 out float a, out float b // Result
28 )
29 {
30 if (!RayIntersectSphere(ro, rd, o, ar, a, b)) return false;
31
32 float pa = 0.0f, pb = 0.0f;
33 if (RayIntersectSphere(ro, rd, o, pr, pa, pb))
34 {
35 b = pa;
36 }
37
38 return true;
39 }
这个函数的实现也很简单，首先检测射线与大气相交的线段，然后在检测射线与地球相交的线段。如果射线被地球挡住的话，就使用与地球相交的线段来构成。
当我们找到了需要积分统计的 AB 线段之后，我们就可以套用上面离散化的 Eq 15，来计算最终的结果。所以，按照积分理论，我们需要确定一个步长。这里我们是指定了计算次数 $N = viewRaySampleN$，然后根据线段 AB 的长度，确定步长 $ds$，如下代码所示：
1 // Split intersect ray into N segment
2 float ds = (lb - la) / viewRaySampleN;
有了步长，有了计算次数，我们就可以循环计算 Eq 15 后面求和的部分，我这里定义为 totalContribution。
得到了 totalContribution 之后，将公式中剩下的部分一一带入，即可得到最终的结果。
下面是 Rayleigh Scattering 和 Mie Scattering 对应的 $\beta (\lambda,0)$ 函数的定义，参考前文中给出的数据：
1 float3 RayleighScatteringCoefficientAtSealevel()
2 {
3 return float3(0.00000519673f, 0.0000121427f, 0.0000296453f);
4 }
5
6 float3 MieScatteringCoefficientAtSealevel()
7 {
8 return float3(21e-6f, 21e-6f, 21e-6f);
9 }
下面是 Rayleigh Scattering 和 Mie Scattering 对应的 $\gamma (\theta)$ 函数的定义，这个函数前文也给出了：
1 float RayleighScatteringPhase(float theta)
2 {
3 return 3.0f * (1 + theta * theta) / (16.0f * 3.1415926f);
4 }
5
6 float MieScatteringPhase(float theta)
7 {
8 const float g = 0.99f;
9 const float g2 = g * g;
10 const float one_minus_g2 = 1.0f - g2;
11 const float one_add_g2 = 1.0f + g2;
12 const float two_add_g2 = 2.0f + g2;
13 float a = 3.0f * one_minus_g2 * (1.0f + theta * theta);
14 float b = 8.0f * 3.1415926f * two_add_g2 * pow(one_add_g2 - 2.0f * g * theta, 3.0f / 2.0f);
15 return a / b;
16 }
上面将公式中大部分都转化为了对应的代码，只剩下了 totalContribution 计算的部分了。
TotalContribution 的计算
totalContribution 是天空渲染公式中需要累计求和的部分 $\sum_{0}^{N}{T_{cp} * \rho (h) * T_{pa}*ds}$，单独拆分出来的代码如下所示：
1 float3 totalContribution = 0.0f;
2 for (uint i = 0u; i < viewRaySampleN; i++)
3 {
4 // Current sample position
5 float tp = (st + ds * 0.5f);
6 float3 pos = viewPos + viewDir * tp;
7
8 float height = distance(pos, planetPos) - planetRadius;
9 totalContribution = totalContribution + RayleighLightContributionIntegration(
10 height, ds, tp, st, viewPos, viewDir, sunDir, lightRaySampleN, planetPos, planetRadius, atmosphereRadius);
11
12 st = st + ds;
13 }
14
15 //------------------------------------------------------------------------------------------
16
17 float3 totalContribution = 0.0f;
18 for (uint i = 0u; i < viewRaySampleN; i++)
19 {
20 // Current sample position
21 float tp = (st + ds * 0.5f);
22 float3 pos = viewPos + viewDir * tp;
23
24 float height = distance(pos, planetPos) - planetRadius;
25 totalContribution = totalContribution + MieLightContributionIntegration(
26 height, ds, tp, st, viewPos, viewDir, sunDir, lightRaySampleN, planetPos, planetRadius, atmosphereRadius);
27
28 st = st + ds;
29 }
流程实际上也很简单，我们根据需要计算的次数 N=viewRaySampleN 进行求和。
每一次计算，我们计算当前采样点所在的位置，以及当前点距离地面的高度，如下代码所示：
1 // Current sample position
2 float tp = (st + ds * 0.5f);
3 float3 pos = viewPos + viewDir * tp;
4
5 float height = distance(pos, planetPos) - planetRadius;
6 ...
7
8 st = st + ds;
计算到了必要的参数之后，调用对应的 LightContributionIntegration 计算 contribution。
RayleighLightContributionIntegration 和 MieLightContributionIntegration
这两个函数相似，所以合在一起讨论了。
这些函数用来计算单次求和部分的公式 $T_{cp} * \rho (h) * T_{pa}*ds$，代码如下所示：
1 float3 RayleighLightContributionIntegration(float h, float ds, float tp, float ta, // Position
2 float3 viewPos, float3 viewDir, // View ray
3 float3 sunDir, // Sun direction
4 uint sampleN, // Sample
5 float3 planetPos, float planentRadius, float atmosphericRadius // Planent
6 )
7 {
8 float lightRayDepth = 0.0f;
9 float3 position = viewPos + viewDir * tp;
10 if (!RayleighOpticalDepthLightRay(position, sunDir, planetPos, planentRadius, atmosphericRadius, sampleN, lightRayDepth))
11 {
12 // Occlussion by earth
13 return float3(0.0f, 0.0f, 0.0f);
14 }
15
16 float viewRayDepth = RayleighOpticalDepthViewRay(viewPos, viewDir, ta, tp, sampleN, planetPos, planentRadius);
17 float ratio = RayleighDensityRatio(h);
18 return exp(-RayleighScatteringCoefficientAtSealevel() * (viewRayDepth + lightRayDepth)) * ratio * ds;
19 }
20
21 float3 MieLightContributionIntegration(float h, float ds, float tp, float ta, // Position
22 float3 viewPos, float3 viewDir, // View ray
23 float3 sunDir, // Sun direction
24 uint sampleN, // Sample
25 float3 planetPos, float planentRadius, float atmosphericRadius // Planent
26 )
27 {
28 float lightRayDepth = 0.0f;
29 float3 position = viewPos + viewDir * tp;
30 if (!MieOpticalDepthLightRay(position, sunDir, planetPos, planentRadius, atmosphericRadius, sampleN, lightRayDepth))
31 {
32 // Occlussion by earth
33 return float3(0.0f, 0.0f, 0.0f);
34 }
35
36 float viewRayDepth = MieOpticalDepthViewRay(viewPos, viewDir, ta, tp, sampleN, planetPos, planentRadius);
37 float ratio = MieDensityRatio(h);
38 return exp(-MieScatteringCoefficientAtSealevel() * (viewRayDepth + lightRayDepth)) * ratio * ds;
39 }
根据 Eq 3 和 Eq 16 所示，得到新 Eq 17
$$\beta(\lambda,h)=\beta(\lambda,0)*\rho(h)$$(Eq 3)
$$T \approx e^{-{\sum_{0}^{N}{\beta(\lambda,h)*ds}}}$$(Eq 16)
$$T \approx e^{-{\sum_{0}^{N}{\beta(\lambda,0)*\rho(h)*ds}}} = e^{-\beta(\lambda,0)*{\sum_{0}^{N}{\rho(h)*ds}}}$$(Eq 17)
我们将公式中的 $\sum_{0}^{N}{\rho(h)*ds}$ 定义为 Optical Depth，所以得到新的 Eq 18：
$$T \approx e^{-\beta(\lambda,0)*D}$$(Eq 18)
最终我们得到 $T_{cp} * T_{pa} = e^{-\beta(\lambda,0)*D_{cp}} * e^{-\beta(\lambda,0)*D_{pa}} = e^{-\beta(\lambda,0)*(D_{cp} + D_{pa})}$。
所以，根据上面的推导，我们的代码就很简单了。首先计算当前采样点位置 position。
接着计算 $D_{cp}$ 的值 lightRayDepth。但是需要注意的是，如果 CP 线段与地球相碰撞，就表示 P 点无法接受到太阳的光照，在地球背阴的一面，直接返回 0。
如果没有被地球遮挡，再计算 $D_{pa}$ 的值 viewRayDepth。
现在就只剩下了 $\rho(h)$，我们就根据海拔高度 $h$ 计算对应的值即可。
以上数据计算完毕之后，带入公式 $T_{cp} * \rho (h) * T_{pa}*ds$ ，就能得到单次求和的值。
如下是 Rayleigh Scattering 和 Mie Scattering 对应的 $ \rho(h)$ 的函数定义：
1 float RayleighDensityRatio(float h)
2 {
3 float H = 8e3;
4 return exp(-h / H);
5 }
6
7 float MieDensityRatio(float h)
8 {
9 float H = 1200;
10 return exp(-h / H);
11 }
如下是 Rayleigh Scattering 和 Mie Scattering 对应的 $D_{cp}$ 函数的定义：
1 bool RayleighOpticalDepthLightRay(
2 float3 p, float3 sunDir, // Light ray
3 float3 planetPos, float planentRadius, float atmosphereRadius, // Planent and Atmosphere
4 uint sampleN, // Sample
5 out float opticalDepth
6 )
7 {
8 float ta = 0.0f, tb = 0.0f;
9 RayIntersectSphere(p, sunDir, planetPos, atmosphereRadius, ta, tb);
10
11 float ds = tb / sampleN;
12 float st = 0.0f;
13
14 opticalDepth = 0.0f;
15 for (uint i = 0; i < sampleN; i++)
16 {
17 // Current sample position
18 float3 pos = p + sunDir * (st + ds * 0.5f);
19
20 // Current sample height
21 float height = distance(pos, planetPos) - planentRadius;
22 if (height < 0.0f) return false;
23
24 opticalDepth = opticalDepth + RayleighDensityRatio(height) * ds;
25
26 st = st + ds;
27 }
28
29 return true;
30 }
31
32 bool MieOpticalDepthLightRay(
33 float3 p, float3 sunDir, // Light ray
34 float3 planetPos, float planentRadius, float atmosphereRadius, // Planent and Atmosphere
35 uint sampleN, // Sample
36 out float opticalDepth
37 )
38 {
39 float ta = 0.0f, tb = 0.0f;
40 RayIntersectSphere(p, sunDir, planetPos, atmosphereRadius, ta, tb);
41
42 float ds = tb / sampleN;
43 float st = 0.0f;
44
45 opticalDepth = 0.0f;
46 for (uint i = 0; i < sampleN; i++)
47 {
48 // Current sample position
49 float3 pos = p + sunDir * (st + ds * 0.5f);
50
51 // Current sample height
52 float height = distance(pos, planetPos) - planentRadius;
53 if (height < 0.0f) return false;
54
55 opticalDepth = opticalDepth + MieDensityRatio(height) * ds;
56
57 st = st + ds;
58 }
59
60 return true;
61 }
最后是 Rayleigh Scattering 和 Mie Scattering 对应的 $D_{pa}$ 函数的定义：
1 float RayleighOpticalDepthViewRay(float3 viewPos, float3 viewDir, // View ray
2 float ta, float tp, // Position
3 uint sampleN, // Sample
4 float3 planetPos, float planentRadius // Planent
5 )
6 {
7 // Split intersect ray into N segment
8 float ds = (tp - ta) / sampleN;
9 float st = ta;
10
11 float opticalDepth = 0.0f;
12 for (uint i = 0u; i < sampleN; i++)
13 {
14 // Current sample position
15 float3 pos = viewPos + viewDir * (st + ds * 0.5f);
16
17 // Current sample height
18 float height = distance(pos, planetPos) - planentRadius;
19
20 opticalDepth = opticalDepth + RayleighDensityRatio(height) * ds;
21
22 st = st + ds;
23 }
24
25 return opticalDepth;
26 }
27
28 float MieOpticalDepthViewRay(float3 viewPos, float3 viewDir, // View ray
29 float ta, float tp, // Position
30 uint sampleN, // Sample
31 float3 planetPos, float planentRadius // Planent
32 )
33 {
34 // Split intersect ray into N segment
35 float ds = (tp - ta) / sampleN;
36 float st = ta;
37
38 float opticalDepth = 0.0f;
39 for (uint i = 0u; i < sampleN; i++)
40 {
41 // Current sample position
42 float3 pos = viewPos + viewDir * (st + ds * 0.5f);
43
44 // Current sample height
45 float height = distance(pos, planetPos) - planentRadius;
46
47 opticalDepth = opticalDepth + MieDensityRatio(height) * ds;
48
49 st = st + ds;
50 }
51
52 return opticalDepth;
53 }
结论
自此，如何计算 Single Scattering 情况下大气散射的代码就全部讲解完毕。
完整的示例工程可在 这里 查看。
需要注意的事，这里的算法十分的简单粗暴，性能并不是很好。而且只是计算了Scattering 的部分，Absortion 的部分也没有计算。文章的本意也旨在讲述大气散射的基本理论，便于后面深入了解更加复杂的，性能，效果更加出色的大气散射算法。如有不明和出错之处，请不吝指出。
