大气散射 ( Atmosphere Scattering)

大气散射介绍

当你看到蓝天白云,日落黄昏,或者从太空遥望地球,地球周围一圈蓝色光晕,有没有想过这些现象产生的原因是什么?这些物理现象用一句话解释:太阳光与大气层中的参与介质发生了大气散射现象,从而形成了天空颜色

散射

光源的辐射沿着路径前行时,碰到参与介质时会发生几种不同散射现象的行为:

  • 出射散射(Out-Scattering) 出射散射也就是光子碰到空气中的分子后,原本要射向摄像机的光线发生了偏转,光线的方向改变了

  • 入射散射(In-Scattering) 入射散射与出射散射相反,本来不指向相机的光线也可能被偏转到相机的方向,这就是入射散射。

根据粒子大小不同可分为瑞利散射(Rayleigh Scattering)米氏散射(Mie Scattering)

什么是散射?

光碰到不同大小粒子后发生散射,不同方向上的光强如下图,箭头长度表示光线强弱,箭头越长,光线越强。

大气散射模型

大气散射模型分为两种:Single Scattering和Multi Scattering。

  • Single Scattering 即只考虑太阳光经过一次散射后抵达我们的眼睛。

  • Multi Scattering,考虑光线多次散射过程的模型。显然这更符合现实世界的情况,但下一级的散射取决于与上一级散射的结果,也就是说这是一个递归方程,实现起来并不太容易。

模型简化

由于第二次散射以上的能量较少,对最终视觉的影响程度有限,所以目前Single Scattering 还是一种相对实现简单效果也不会差很多的模型,这里也主要介绍 Single Scattering。

因此单次散射大气模型求解问题如下: 我们从路径 AB 观察大气,并且求解 B 点的大气颜色,光线在大气中只发生一次散射,散射点为 P:

实际上在路径 AB 上有无数个 P 点,因此最终求解是对 AB 路径上每一个点的光照贡献进行累加(求积分)

因此整个模型可以分解成如下问题:

  • P 点光照 :光线穿过大气会衰减,到达 P 点后的光照强度。
  • 光线在 P 散射:进入相机方向的光照,即入射散射部分
  • 累加所有 P 点光照

透射函数

为了计算散射点的光照,我们先考虑如下一段光线:

  • 光线从太阳传到 C,此时是真空,所以没有损失,我们用 \(I_C\) 指代 C 点的光强
  • 在从 C 点到 P 点的过程在中,经过 CP 路径衰减,我们用 \(I_P\) 指代 P 点的光强
  • \(I_P\)\(I_C\) 的值就是 透射率(Transimittance) :\(T (CP) = \frac{I_P}{I_C}\)

P 点的光强也就是: \[I_P = I_C * \textcolor{Green}{T(CP)}\]

其中,T项衰减系数(Transmittance),它表示在某段路径上的对光照的衰减程度。该公式也可以被认为是零级散射(zero scattering),即不考虑任何散射事件、直接考虑经过衰减后光强。

散射函数:S(λ,θ,h)

P点在接受光照后,会因为散射偏转一部分,于是我们需要散射函数 \(S(\lambda,\theta,h)\) 来表示某一方向上光线的偏转情况。如下图,那些偏转了 \(\theta\) 角的光线才能被指向A点的摄像机:

于是我们可以得出从 P 点出发到 A 点的光强:

\[I_{PA}=\boxed{I_P}\textcolor{Red}{S(\lambda,\theta,h)} \textcolor{Green}{T(PA)} \]

从公式中可以看出,散射跟角度 \(\theta\)、光线波长 \(\lambda\) 以及 \(P\) 点高度 \(h\) 有关系。

散射类别

散射系数和粒子的大小和折射率有关,这里就得解释一下大气粒子的分类。我们知道大气中有很多不同种类的大气分子/粒子,一般在大气渲染模型里把它们分成两大类。一种是小分子,例如氮气和氧气分子等,另一种是大粒子,例如各种尘埃粒子。之所以要进行分类是因为它们对光线有着不一样的行为(这里忽略对光的吸收,只考虑散射):

  • 小分子:指大小远小于光线波长的粒子。小分子对光的散射在前后方向上分布比较均匀,通常会使用 Rayleigh散射(Rayleigh Scattering) 对它们进行建模。由于这些分子大小比波长还要小很多,因此光的波长也会影响 Rayleigh 散射的程度。
  • 大粒子:指大小远大于光线波长的粒子。大粒子在发生散射的时候会把更多的光散射到前向,通常会使用 Mie散射(Mie Scattering) 对它们进行建模。由于光的波长相较于这些粒子大小来说是可以忽略的,因此我们认为 Mie 散射跟光线波长无关

Rayleigh 散射

对于足够小的粒子,光线在撞击它后会经历什么?我们通常使用瑞利散射来建模。

Rayleigh 散射曲线图如下:

散射函数如下,这个方程并不是真正的Rayleigh散射方程,你去维基百科里面查,会发现另外一个公式,这个公式来自这篇论文

\[\textcolor{Red}{S(\lambda,\theta,h)} = \frac{\pi^2(n^2-1)^2}{2}\frac{\textcolor{Gold}{\rho(h)}}{N}\frac{1}{\lambda^4}(1+cos\theta^2)\]

其中,

  • \(\lambda\) 是光的波长
  • \(n=1.00029\) 是粒子折射率
  • \(N=2.504*10^{25}\) 是海平面处的大气密度,单位是 分子数/立方米
  • \(\rho(h)\) 是高度 \(h\) 处的相对大气密度(即相对于海平面的密度,可以理解成 \(h\) 处真正的大气密度与海平面处大气密度的比值,因此它在海平面处值为 1,随着 \(h\) 增加不断减小),我们一般会使用指数函数对它的真实曲线进行数学拟合(是一种近似拟合,不要和后面提到的指数函数弄混):

\[\textcolor{Gold}{\rho(h)}=exp(-\frac{h}{H})\]

其中 H“Scale Height”,相当于是一个基准的高度。

  • 对于 Rayleigh 来说,H= 8500米
  • 对于 Mie 来说,H = 1200米

可以看出Rayleigh散射大致和波长的4次幂的成反比,波长越小(越靠近紫光)的光被散射得越厉害。所以白天的时候天空为蓝色,因为蓝光在大气里不断被散射,黄昏的时候天空会变红,因为相比于白天,阳光此时要穿越更厚得多的大气层,在到达人眼之前,大多数蓝光都被散射到其他方向,所以剩下来的就是红光了。

Rayleigh 散射系数

Rayleigh 散射的公式 \(S(\lambda,\theta,h)\) 表示在某个特定方向上的光强度,要求总散射系数就需要对散射函数做球面积分:

Tips

球面积分示意图,这里假定 r 等于 1,对球面积分时,因为 \(d\theta\)\(d\phi\) 非常小,因此球面可以看成是矩形,因此蓝色框长度为 \(w=rsin\theta d\phi\) (弦长公式:蓝色面片上边缘半径 \(rsin\theta\),角度 \(d\phi\),高为 \(d\theta\),则面积 \(S=rsin\theta d\phi d\theta\)

因此:

\[\beta(\lambda, h) = \int_{0}^{2\pi}{\int_{0}^{\pi}{S(\lambda,\theta,h)}sin\theta d\theta d\phi}\]

最后得出:

\[\beta(\lambda, h) = \frac{8\pi^3(n^2-1)^2}{3}\frac{\textcolor{Gold}{\rho(h)}}{N}\frac{1}{\lambda^4}\]

特别地,当 \(h=0\) 时,即表示海平面的散射系数:

\[ \beta(\lambda, 0) = \frac{8\pi^3(n^2-1)^2}{3} \frac{ \textcolor{Gold}{1} }{N} \frac{1}{\lambda^4} \]

举例来说就,对于红绿蓝的波长,在海平面时,可以得出这些结果:

\[\begin{aligned} \beta(\textcolor{Red}{680nm}) &= 5.2*10^{-6} \\ \beta(\textcolor{Green}{550nm}) &= 12.1*10^{-6} \\ \beta(\textcolor{DodgerBlue}{440nm}) &= 29.6*10^{-6} \end{aligned}\]

有些博客跟论文里计算结果如下: \(\beta_rrgb=(5.8, 13.5, 33.1) \cdot 10^{-6}\)\(n=1.0003,N=2.545\cdot10^{25}\)

>>> 8*math.pi**3*(1.00029**2-1)**2/3.0/(2.504e25)/((680e-9)**4)
5.196731735928312e-06
>>> 8*math.pi**3*(1.00029**2-1)**2/3.0/(2.504e25)/((550e-9)**4)
1.2142697926864656e-05
>>> 8*math.pi**3*(1.00029**2-1)**2/3.0/(2.504e25)/((440e-9)**4)
2.964525861050941e-05

总散射系数公式也可以拆解为:

\[\beta(\lambda, h) = \frac{8\pi^3(n^2-1)^2}{3}\frac{\textcolor{Gold}{\rho(h)}}{N}\frac{1}{\lambda^4}=\beta(\lambda) \textcolor{Gold}{\rho(h)}\]

其中:

\[\beta(\lambda) = \frac{8\pi^3(n^2-1)^2}{3}\frac{1}{N}\frac{1}{\lambda^4}\]

Rayleigh 相位函数

Rayleigh 散射的原方程可以分解成两个分量,一个是上面的有关能量强度的散射系数,另一个分量则和它几何形状有关

\[S(\lambda,\theta,h)=\beta(\lambda, h)\underbrace{\gamma(\theta)}_{geometry}\]

\(\gamma(\theta)\) 可以用前两者之比得到。

\[\begin{aligned} \gamma(\theta) &= \frac{S(\lambda,\theta,h)}{\beta(\lambda, h)}\\ &=\frac{3}{16\pi}(1+cos^2\theta) \end{aligned}\]

该表达式不依赖波长,\(\gamma(\theta)\) 的函数图像就是之前两瓣的偶极子形状

PS: 整个散射函数可以拆分成 \(\textcolor{Red}{S(\lambda,\theta,h)}= \textcolor{DodgerBlue}{ \beta(\lambda)} \textcolor{Gold}{\rho(h)} \textcolor{Purple}{\gamma(\theta)}\)

Mie 散射

光线的散射是很复杂的,瑞利散射对原子和分子可以有效模拟,但是对光线和大物体(气溶胶)的交互没什么作用。地球大气布满了气溶胶,所以仅是瑞利散射是不足以复现的,为此常用的就是 米氏散射(Mie scattering),这种散射类型倾向于扩散光线,让光源看起来比实际大,但不会改变光的颜色

Mie 散射系数

米散射系数很难去推导计算,基于统计数据的结果,米散射系数可以常量进行估计。论文10中给出海平面的米散射系数:

\[\beta(\textcolor{Red}{680nm},\textcolor{Green}{550nm} , \textcolor{DodgerBlue}{440nm} = (2.2*10^{-5}, 2.2*10^{-5}, 2.2*10^{-5})\]

Mie 散射相位函数

真实的 Mie 散射相位函数图像如下:

因此 Mie 相位函数也是近似方式得出,我们可以利用Henyey-Greenstein 函数,该函数在1941年提出:

\[\gamma_{HG}(\theta)=\frac{1}{4\pi}\frac{1-g^2}{(1+g^2-2gcos\theta)^\frac{3}{2}}\]

在 1992 年 Cornette 对其进行了改进:

\[\gamma_{M}(\theta)=\frac{3}{8\pi}\frac{1-g^2}{(2+g^2)}\frac{1+cos^2\theta}{(1+g^2-2gcos\theta)^\frac{3}{2}}\]

该函数引入了一个参数 \(g(-1<g<1)\),决定了前后散射的相对强度,代表散射方向的平均余弦。正的 \(g\) 值会增加正向散射的光量,负的 \(g\) 值表示向后散射的光量。 \(g=0\) 的值导致各向同性散射。

对于地球大气层,一般可以取 \(g=0.758\)。然后由于上述两个函数计算成本比较高,快速应用可以使用 Schlick 相位函数:

\[\lambda_S(\theta)=\frac{1}{4\pi}\frac{1-g^2}{(1+gcos\theta)^2}\]

衰减系数 \(\textcolor{Green}{T}\)

\(T(a,b)\) 表示 \(a,b\) 之间的透射率(transmittance),其定义为9

\[T(a,b)=exp(-\int_a^b\sigma_t(p)dl_p\]

其中 \(\sigma_t\) 为介质衰减系数,

  • Rayleigh 在大气中,Rayleigh 散射可以用 Rayleigh theory 近似描述:

    \(\sigma_RT(h,\lambda)=\beta(\lambda, h)\)

即:

\[\textcolor{Green}{T(PA)}=exp\{-\int_P^A{\textcolor{DodgerBlue}{\beta(\lambda, h)}}ds\}\]

带入之前求解的总散射系数公式: \(\beta(\lambda,h)\)

\[\textcolor{Green}{T(PA)}=exp\{-\textcolor{DodgerBlue}{\beta(\lambda)}\int_P^A{\textcolor{Gold}{\rho(h)}}ds\}\]

其中的积分项对应的是路径长度 \(x\),可以理解为在光传播路径上,对大气密度函数 \(\textcolor{Gold}{\rho(h)}\) 求积分,这部分被称为 光学距离——Optical Depth

\[\textcolor{Darkorange}{D(PA)}=\int_P^A{\textcolor{Gold}{\rho(h)}}ds\}\]

总结

至此,我们简单推导了 \(S\)\(T\),这样一开始的公式可以写成:

\[\begin{aligned} I_{PA} &= \boxed{I_P}\textcolor{Red}{S(\lambda,\theta,h)} \underline {\textcolor{Green}{T(PA)} } \\ &=\boxed{I_S\textcolor{Green}{T(CP)}} \textcolor{Red}{S(\lambda,\theta,h)} \underline { \textcolor{Green}{T(PA)} } \\ &= \boxed { I_S \{exp\{-\textcolor{DodgerBlue}{\beta(\lambda)}\textcolor{Darkorange}{D(CP)}\}\} } \textcolor{Red}{ \{ \textcolor{DodgerBlue}{ \beta(\lambda)} \textcolor{Gold}{\rho(h)} \textcolor{Purple}{\gamma(\theta)} \} } \underline { \{exp\{-\textcolor{DodgerBlue}{\beta(\lambda)}\textcolor{Darkorange}{D(PA)}\}\} }\\ &= I_S \textcolor{Red}{ \{ \textcolor{DodgerBlue}{\beta(\lambda)} \textcolor{Gold}{\rho(h)} \textcolor{Purple}{\gamma(\theta)} \} } \{exp\{ -\textcolor{DodgerBlue}{\beta(\lambda)} \{ \textcolor{Darkorange}{D(CP)} + \textcolor{Darkorange}{D(PA)} \} \} \} \end{aligned} \]

其中,一般情况下,太阳光照角度固定,可以看成是平行光,\(\textcolor{Purple}{\gamma(\theta)}\) 可以看成常量,\(\textcolor{DodgerBlue}{\beta(\lambda)}\) 也是常量,也就是海平面附近的散射系数,在 Shader 中作为输入即可,\(\textcolor{Darkorange}{D(PA)}\) 等在 Shader 中计算。

总结起来,公式表示的含义就是,计算观察路径AB上某一点P的光照贡献:

  • 光线从太阳出发,到达大气边缘的C点
  • 经过路径CP上的衰减到达P点(0级散射)
  • 在P点发生一次散射,将一部分光散射到了观察方向上(1级散射)
  • 这部分光又经过AP路径上的衰减最终到达了我们的眼睛

数值积分

P只是AB上的某一点,我们需要对路径AB上的每一点进行积分:

\[\begin{aligned} I_A &= \sum_{i=0}^{n}I_{P_n} \\ &= \int_{A}^{B} I_{PA} ds \\ &= \int_{A}^{B} I_S \textcolor{Red}{ \{ \textcolor{DodgerBlue}{\beta(\lambda)} \textcolor{Gold}{\rho(h)} \textcolor{Purple}{\gamma(\theta)} \} } \{ exp\{ -\textcolor{DodgerBlue}{\beta(\lambda)} \{ \textcolor{Darkorange}{D(CP)} + \textcolor{Darkorange}{D(PA)} \} \} \} ds \\ &= I_S \textcolor{DodgerBlue}{\beta(\lambda)} \textcolor{Purple}{\gamma(\theta)} \int_{A}^{B} {exp\{ -\textcolor{DodgerBlue}{\beta(\lambda)} \{ \textcolor{Darkorange}{D(CP)} + \textcolor{Darkorange}{D(PA)} } \}\} \textcolor{Gold}{\rho(h)} ds \end{aligned} \]

Tips

在大气shader中,我们将 \(AB\) 分割成等长的若干段,然后取这些段的中点记为 \(P_i\),然后分别计算每个点的 \(I_{P_i}\),最后将这些点的光照累加起来。

大气球体shader框架思路

光线与大气相交

从之前得出的公式中可以得知,我们需要在 \(AB\) 路径上求光学距离,我们第一步需要求出 \(AB\) 的长度:

如下图所示,灰色部分是大气,我们需要求出观察路径 \(AB\),其中已知:

  • \(O\):相机位置
  • \(C\):大气球心位置
  • \(L\):相机到大气球心长度
  • \(R\):大气球半径
  • \(\vec{OD}\) :观察方向,单位向量

如图可以依次求出,

  • \(OT = \vec{OC} \cdot \vec{OD}\)
  • \(TC = \sqrt{L^2 - OC^2}\)
  • \(AT = \sqrt{R^2 - TC^2}\)
  • \(OA = OT-AT\)
  • \(OB = OT+BT=OT+AT\)
bool RayIntersect(
    // Ray
    float3 O, // Origin
    float3 D, // Direction

    // Sphere
    float3 C, // Centre
    float R,    // Radius
    out float AO, // First intersection time
    out float BO  // Second intersection time
)
{
    float3 L = C - O;
    float DT = dot (L, D);
    float R2 = R * R;

    float CT2 = dot(L,L) - DT*DT;
    
    // CT 长度超过了求半径 R,此时视线跟大气不相交
    // Intersection point outside the circle
    if (CT2 > R2)
        return false;

    float AT = sqrt(R2 - CT2);
    float BT = AT;

    AO = DT - AT;
    BO = DT + BT;
    return true;
}

光线与行星相交

视线有可能会被行星阻挡,因此需要做两次相交函数,一次判断大气,一次判断行星。

// Intersections with the atmospheric sphere
//     D    : 观察方向向量
// worldPos : 相机坐标 
float tA;    // Atmosphere entry point (worldPos + D * tA)
float tB;    // Atmosphere exit point  (worldPos + D * tB)
if (!RayIntersect(O, D, _PlanetCentre, _AtmosphereRadius, tA, tB))
    return fixed4(0,0,0,0); // The view rays is looking into deep space

// Is the ray passing through the planet core?
float pA, pB;
if (RayIntersect(O, D, _PlanetCentre, _PlanetRadius, pA, pB))
    tB = pA;

采样视线

之前我们得到了 \(AB\) 路径上每个点 \(P\) 的光照公式:

\[\begin{aligned} I_A &= I_S \textcolor{DodgerBlue}{\beta(\lambda)} \textcolor{Purple}{\gamma(\theta)} \int_{A}^{B} {exp\{ -\textcolor{DodgerBlue}{\beta(\lambda)} \{ \textcolor{Darkorange}{D(CP)} + \textcolor{Darkorange}{D(PA)} } \}\} \textcolor{Gold}{\rho(h)} ds \\ &= I_S \textcolor{DodgerBlue}{\beta(\lambda)} \textcolor{Purple}{\gamma(\theta)} \sum_{P \in AB} \textcolor{Green}{T(CP)} \textcolor{Green}{T(PA)} \textcolor{Gold}{\rho(h)} ds \end{aligned}\]

线段上有无数个 \(P\) 点,为了近似求解 \(I\),我们需要把线段划分为数个长度为 \(ds\) 的小段

\(AB\) 被划分的段数,就是视线采样数(View Samples),在 Shader 可以这么处理:

// Numerical integration to calculate
// the light contribution of each point P in AB
float3 totalViewSamples = 0;
float time = tA;
float ds = (tB - tA) / (float)(_ViewSamples);
for (int i = 0; i < _ViewSamples; i ++)
{
    // Point position
    // (sampling in the middle of the view sample segment)
    // O : 相机位置
    // D : 相机方向向量
    float3 P = O + D * (time + ds * 0.5);

    // T(CP) * T(PA) * ρ(h) * ds
    totalViewSamples += ViewSampling(P, ds);

    time += ds;
}

// I = I_S * β(λ) * γ(θ) * totalViewSamples
float3 I = _SunIntensity *  _ScatteringCoefficient * phase * totalViewSamples;

接下来就是求解每个 \(P\) 点上的 \(\textcolor{Green}{T(CP)} \textcolor{Green}{T(PA)} \textcolor{Gold}{\rho(h)}ds\)

光学深度 \(PA\)

回顾一下之前的公式:

\[ \textcolor{Green}{T(CP)} \textcolor{Green}{T(PA)} = exp\{ -\textcolor{DodgerBlue}{\beta(\lambda)} \{ \textcolor{Darkorange}{D(CP)} + \textcolor{Darkorange}{D(PA)} \} \}\]

其中:

\[\textcolor{Darkorange}{D(PA)}=\sum_{Q \in PA} exp \{ -\frac{h_Q}{H} \}ds\]

\(h_Q\) 表示当前点的高度,有一点注意 \(\textcolor{Gold}{\rho(h)}=exp(-\frac{h_Q}{H})\),这里顺便可以求解 \(AB\) 上的积分 \(\textcolor{Gold}{\rho(h)}ds\)

Shader 实现如下:

// + Accumulator for the optical depth
float opticalDepthPA = 0;

float3 totalViewSamples = 0;
float time = tA;
float ds = (tB - tA) / (float)(_ViewSamples);
for (int i = 0; i < _ViewSamples; i ++)
{
    float3 P = O + D * (time + ds * 0.5);

    ///////// Begin ViewSampling()
    // T(CP) * T(PA) * ρ(h) * ds

    // 1: ρ(h) * ds

    //            C : 行星球形坐标
    // _ScaleHeight : Rayleigh/Mie 散射里的基准高度

    float height = distance(C, P) - _PlanetRadius;
    float opticalDepthSegment = exp(-height / _ScaleHeight) * ds;


    // 2: D(PA)
    opticalDepthPA += opticalDepthSegment;

    ////////    End ViewSampling()

    time += ds;
}

// I = I_S * β(λ) * γ(θ) * totalViewSamples
float3 I = _SunIntensity *  _ScatteringCoefficient * phase * totalViewSamples;

光学深度 \(CP\)

累加符号里 \(P\) 点的光量贡献还剩下 \(CP\) 的光学深度,我们在求 \(CP\) 的光学深度时,建立一个 LightSampling 方法,也就从 \(P\) 点出发指向太阳进行采样,把太阳的出射点叫做 \(C\) (注意跟之前的球心坐标 \(C\) 区分),然后注意下面 \(C_0\) 的情况,光线被星球遮挡,此时要忽略 \(P_0\) 点的贡献。

Shader 实现如下:

bool LightSampling(
    float3 P,    // Current point within the atmospheric sphere
    float3 S,    // S 是光照方向,从 P 点到 C 点的平行光
    out float opticalDepthCA
)
{
    float _; // don't care about this one

    // 这里 C 是局部变量,表示光线在大气表面的入射点
    float C;
    RayInstersect(P, S, _PlanetCentre, _AtmosphereRadius, _, C);

    // Samples on the segment PC
    // _LightSamples : CP 上的分段次数

    float time = 0;
    float ds = distance(P, P + S * C) / (float)(_LightSamples);
    for (int i = 0; i < _LightSamples; i ++)
    {
        // Q 点如上图所示
        // S : 光照方向
        float3 Q = P + S * (time + lightSampleSize * 0.5);
        float height = distance(_PlanetCentre, Q) - _PlanetRadius;
        // Inside the planet
        if (height < 0)
            return false;

        // Optical depth for the light ray
        opticalDepthCA += exp(-height / _RayScaleHeight) * ds;

        time += ds;
    }

    return true;
}

然后之前的代码加入 \(CP\) 的计算如下:

float opticalDepthPA = 0;

float3 totalViewSamples = 0;
float time = tA;
float ds = (tB - tA) / (float)(_ViewSamples);
for (int i = 0; i < _ViewSamples; i ++)
{
    float3 P = O + D * (time + ds * 0.5);

    ///////// Begin ViewSampling()
    // T(CP) * T(PA) * ρ(h) * ds

    // ρ(h) * ds
    float height = distance(C, P) - _PlanetRadius;
    float opticalDepthSegment = exp(-height / _ScaleHeight) * ds;

    // D(PA)
    opticalDepthPA += opticalDepthSegment;

    // D(CP)
    float opticalDepthCP = 0;
    bool overground = LightSampling(P, S, opticalDepthCP);

    if (overground)
    {
        // Combined transmittance
        // T(CP) * T(PA) = T(CPA) = exp{ -β(λ) [D(CP) + D(PA)]}
        float3 transmittance = exp(
            -_ScatteringCoefficient *
            (opticalDepthCP + opticalDepthPA)
        );

        // Light contribution
        // T(CPA) * ρ(h) * ds
        totalViewSamples += transmittance * opticalDepthSegment;
    }

    ////////    End ViewSampling()

    time += ds;
}

// I = I_S * β(λ) * γ(θ) * totalViewSamples
float3 I = _SunIntensity *  _ScatteringCoefficient * phase * totalViewSamples;

散射系数

之前我们重点以 Rayleigh 散射来介绍如何计算单次大气散射,接下来加上 Mie 散射,之前的方程:

\[\begin{aligned} I &= I_S \textcolor{DodgerBlue}{\beta(\lambda)} \textcolor{Purple}{\gamma(\theta)} \sum_{P \in AB} L(P) \\ \end{aligned}\]

加上 Mie 反射后:

\[\begin{aligned} I &= \overbrace{ I_S \textcolor{DodgerBlue}{\beta_R(\lambda)} \textcolor{Purple}{\gamma_R(\theta)} \sum_{P \in AB} L_R(P) }^{Rayleigh \ Scattering} + \overbrace{ I_S \textcolor{DodgerBlue}{\beta_M(\lambda)} \textcolor{Purple}{\gamma_M(\theta)} \sum_{P \in AB} L_M(P) }^{Mie \ Scattering} \\ &=I_S \left( \overbrace{ \textcolor{DodgerBlue}{\beta_R(\lambda)} \textcolor{Purple}{\gamma_R(\theta)} \sum_{P \in AB} L_R(P) }^{Rayleigh \ Scattering} + \overbrace{ \textcolor{DodgerBlue}{\beta_M(\lambda)} \textcolor{Purple}{\gamma_M(\theta)} \sum_{P \in AB} L_M(P) }^{Mie \ Scattering} \right) \end{aligned}\]

因此代码实现也需要做修改:

bool LightSampling(
    float3 P,
    float3 S,
    // Modify : float ->  float2  
    // 分别计算光学深度 .xy  ->  Rayleight Mie
    out float2 opticalDepthCA
)
{
    float _;
    float C;
    RayInstersect(P, S, _PlanetCentre, _AtmosphereRadius, _, C);

    float time = 0;
    float ds = distance(P, P + S * C) / (float)(_LightSamples);
    for (int i = 0; i < _LightSamples; i ++)
    {
        float3 Q = P + S * (time + lightSampleSize * 0.5);
        float height = distance(_PlanetCentre, Q) - _PlanetRadius;
        if (height < 0)
            return false;

        // 分别计算 Rayleigh、Mie 的光学深度,基准的高度值不一样
        opticalDepthCA.x += exp(-height / _RayScaleHeight.x) * ds;
        opticalDepthCA.y += exp(-height / _RayScaleHeight.y) * ds;

        time += ds;
    }

    return true;
}
// Modify : float -> float2
float2 opticalDepthPA = 0;

float3 totalRayViewSamples = 0;
float3 totalMieViewSamples = 0;

float time = tA;
float ds = (tB - tA) / (float)(_ViewSamples);
for (int i = 0; i < _ViewSamples; i ++)
{
    float3 P = O + D * (time + ds * 0.5);

    ///////// Begin ViewSampling()
    // T(CP) * T(PA) * ρ(h) * ds

    // ρ(h) * ds
    float height = distance(C, P) - _PlanetRadius;
    // float -> float2
    float2 opticalDepthSegment = 0;
    opticalDepthSegment.x = exp(-height / _ScaleHeight.x) * ds;
    opticalDepthSegment.y = exp(-height / _ScaleHeight.y) * ds;

    // D(PA)
    opticalDepthPA += opticalDepthSegment;

    // D(CP)
    // Modify : float -> float2
    float2 opticalDepthCP = 0;
    bool overground = LightSampling(P, S, opticalDepthCP);

    if (overground)
    {
        // Combined transmittance
        // T(CP) * T(PA) = T(CPA) = exp{ -β(λ) [D(CP) + D(PA)]}
        float3 transmittanceRay = exp(
            -_RayScatteringCoefficient  *
            (opticalDepthCP.x + opticalDepthPA.x)
        );

        float3 transmittranceMie = exp(
            -_MieScatteringCoefficient  *
            (opticalDepthCP.x + opticalDepthPA.y)
        );

        // Light contribution
        // T(CPA) * ρ(h) * ds
        totalRayViewSamples  += transmittanceRay * opticalDepthSegment;
        totalMieViewSamples  += transmittanceMie * opticalDepthSegment;
    }

    ////////    End ViewSampling()

    time += ds;
}

// I = I_S [βr(λ) * γr(θ) * totalRayViewSamples
//     + βm(λ) * γm(θ) * totalMieViewSamples)
float3 I = _SunIntensity * (
 _RayScatteringCoefficient * rayPhase * totalRayViewSamples +
 _MieScatteringCoefficient * miePhase * totalMieViewSamples
);

其中常量部分:

\(\_RayScatteringCoefficient(\textcolor{Red}{680nm},\textcolor{Green}{550nm},\textcolor{DodgerBlue}{440nm}) = (5.8, 13.5, 33.1)*10^{-6}\) \(\_MieScatteringCoefficient(\textcolor{Red}{680nm},\textcolor{Green}{550nm},\textcolor{DodgerBlue}{440nm})=(2.2, 2.2, 2.2)*10^{-5}\)

\(rayPhase\)\(miePhase\) 提前计算好的相位函数值,如果光源会改变方向,可以在 Shader 中增加函数来计算数值。

\(\_ScaleHeight=(8500, 1200)\)

参考资料

1.Volumetric Atmospheric Scattering

2.Unity基于体绘制的大气散射shader

3.【实战】从零实现一套完整单次大气散射

4.基于物理的大气渲染 -- 冯乐乐

5.Extended Disney "principled" Shader

6.Atmospheric Scattering for Unity 5

7.Sébastien Hillaire Website

8.基于物理的大气散射

9.预计算大气散射模型:原理与实现

10.Precomputed Atmospheric Scattering