GraphicsLab 之 Atmospheric Scattering (二)
作者：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 的部分也没有计算。文章的本意也旨在讲述大气散射的基本理论，便于后面深入了解更加复杂的，性能，效果更加出色的大气散射算法。如有不明和出错之处，请不吝指出。
