GraphicsLab 之 Atmospheric Scattering (二)

時間 2020-11-28

標籤 html git github 算法 app ide 函數性能 spa 欄目 HTML 简体版

原文原文鏈接

做者：i_dovelemonhtml

日期：2020-11-25git

主題：Atmospheric Scattering, Volume Scattering, Rayleigh Scattering, Mie Scattering, Single Scattering, Multiple Scatteringgithub

引言

前文 GraphicsLab 之 Atmospheric Scattering（一）講述了基於物理的天空渲染在 Single Scattering 狀況下的基礎理論知識。本篇文章將主要從代碼實現的角度，詳細講解如何在 Unity 中實現一個簡單粗暴版本的大氣散射天空。算法

離散化天空渲染方程

回顧以前一篇文章中，咱們最後得出的用於計算天空顏色的渲染方程，以下所示：app

$$I_{final}=\int_{A}^{B}{I_{sun} * T_{cp} * S(\lambda ,\theta ,h) * T_{pa}}$$(Eq 10) ide

簡化天空渲染方程

根據以前文章的描述，最後渲染方程積分裏面，有一些數據是不會發生變化的，因此能夠提取出來。函數

最簡單直觀的就是，太陽光的強度 $I_{sun}$ 不會隨着積分發生變化，因此能夠提取出來。性能

根據 Eq 2 和 Eq 3 能夠獲得新的 Eq 11，以下所示：ui

$$S(\lambda ,\theta ,h)=\beta (\lambda ,h) * \gamma (\theta)$$(Eq 2)spa

$$\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 }

View Code

有了覆蓋全屏幕的面片以後，咱們就能夠經過對每個像素進行大氣散射的計算，從而實現天空背景的渲染。

天空像素計算

有了覆蓋全屏幕的後處理以後，就能夠計算每個像素的顏色了。以下是計算每個像素的 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 }

View Code

前面說過，咱們將大氣分爲了兩個大類，分別使用 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 }

View Code

有了這些數據以後，咱們就來計算當前像素的 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 }

View Code

能夠看到這兩個函數的流程是一摸同樣的，它們的不一樣也就體如今最終計算時的一些參數不一樣之上。因此，下面就合在一塊兒說明了。

第一步，咱們須要計算下當前觀察射線是否與大氣相交。這是由於，若是咱們從太空中觀察地球的話，有很大可能一些觀察射線是不會碰撞到大氣層的。固然若是是在地球表面的話，確定是會碰撞到的。

若是碰撞到了大氣，咱們就須要計算實際在大氣中的線段 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 }

View Code

這個函數的實現也很簡單，首先檢測射線與大氣相交的線段，而後在檢測射線與地球相交的線段。若是射線被地球擋住的話，就使用與地球相交的線段來構成。

當咱們找到了須要積分統計的 AB 線段以後，咱們就能夠套用上面離散化的 Eq 15，來計算最終的結果。因此，按照積分理論，咱們須要肯定一個步長。這裏咱們是指定了計算次數 $N = viewRaySampleN$，而後根據線段 AB 的長度，肯定步長 $ds$，以下代碼所示：

1 // Split intersect ray into N segment
2 float ds = (lb - la) / viewRaySampleN;

View Code

有了步長，有了計算次數，咱們就能夠循環計算 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 }

View Code

下面是 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 }

View Code

上面將公式中大部分都轉化爲了對應的代碼，只剩下了 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 }

View Code

流程實際上也很簡單，咱們根據須要計算的次數 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;

View Code

計算到了必要的參數以後，調用對應的 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 }

View Code

根據 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 }

View Code

以下是 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 }

View Code

最後是 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 }