基于物理的实时大气散射原理与实现

基于Unity的实现 https://gitee.com/alienity/atmospheric-scattering

Paper中说明的大都是原理，在作者给出的具体实现中有较多的优化，有兴趣的朋友可以在看完原理后，直接看作者的实现及其说明。本文实现的代码也是基于作者的代码

物理模型

物理模型中大气包含空气分子和气溶胶粒子两种介质，一般是在 $R_g=6030km$ 到 $R_t=6420km$ 之间

对应这两种介质，分别会发生瑞利散射和Mie散射。瑞利散射发生的条件是粒子半径远小于入射光波长。在地球上大气的区间高度是 $R_g=6360km$ 和 $R_t=6420km$

瑞利散射

在每个点上，与入射光偏离 $\theta$ 角度后的强度等于散射系数 ${\beta }^s$ 和相函数 $P$ 的乘积。散射系数依赖于分子密度，相函数描述了角度依赖。

$$\begin{array}{rl}{\beta }_R^S(h,\lambda )& =\frac{8{\pi }^3{\left(n^2-1\right)}^2}{3N{\lambda }^4}{e}^{-\frac{h}{H_R}}\\ P_R(\mu )& =\frac{3}{16\pi }\left(1+{\mu }^2\right)\text{ where }\mu =\cos \theta \end{array}$$

其中 $h=e-R_g$ 表示的是海拔高度，$\lambda$ 是波长，$n$ 是空气折射系数，$N$ 是在海平面 $R_g$ 的分子数密度，$H_R=8km$ 是大气厚度，其中对于波长 $\lambda =(680,550,440)nm$ 的散射系数 ${\beta }_R^S=(5.8,13.5,33.1)10^{-6}{m}^{-1}$

Mie散射

气溶胶也有呈指数下降的密度，只是 height scale $H_M\simeq 1.2km$

Cornette-Shanks 相函数
$$\begin{array}{rl}{\beta }_M^S(h,\lambda )& ={\beta }_M^s(0,\lambda )e^{-\frac{h}{H_M}}\\ P_M(\mu )& =\frac{3}{8\pi }\frac{\left(1-g^2\right)\left(1+{\mu }^2\right)}{\left(2+g^2\right){\left(1+g^2-2g\mu \right)}^{3/2}}\end{array}$$

与分空气分子不一样，气溶胶吸收一定比例的光，吸收系数是 ${\beta }_M^a$，所以吸光系数是 ${\beta }_M^e={\beta }_M^s+{\beta }_M^a$

渲染等式

$L(\mathbf{x},\mathbf{v},\mathbf{s})$ 是太阳在 $\mathbf{s}$ 方向时，从方向 $\mathbf{v}$ 到达 $\mathbf{x}$ 的 radiance of light

• 点 ${\mathbf{x}}_o$ 与点 $\mathbf{x}$ 之间的 transmittance $T$；
• 在点 ${\mathbf{x}}_o$ 处反射的 radiance $\mathcal{L}$；
• 在 $\mathbf{y}$ 点朝方向 $-\mathbf{v}$ 散射的radiance $\mathcal{J}$

T ( x , x o ) = exp ⁡ ( − ∫ x x o ∑ i ∈ { R , M } β i e ( y ) d y ) I [ L ] ( x o , s ) = α ( x o ) π ∫ 2 π L ( x 0 , ω , s ) ω ⋅ n ( x o ) d ω ,  or  0 J [ L ] ( y , v , s ) = ∫ 4 π ∑ i ∈ { R , M } β i S ( y ) P i ( v . ω ) L ( y , ω , s ) d ω \begin{aligned} T\left(\mathbf{x}, \mathbf{x}_{o}\right) &=\exp \left(-\int_{\mathbf{x}}^{\mathbf{x}_{o}} \sum_{i \in\{R, M\}} \beta_{i}^{e}(\mathbf{y}) d y\right) \\ \mathcal{I}[L]\left(\mathbf{x}_{o}, \mathbf{s}\right) &=\frac{\alpha\left(\mathbf{x}_{o}\right)}{\pi} \int_{2 \pi} L\left(\mathbf{x}_{0}, \boldsymbol{\omega}, \mathbf{s}\right) \boldsymbol{\omega} \cdot \mathbf{n}\left(\mathbf{x}_{o}\right) d \omega, \text { or } 0 \\ \mathcal{J}[L](\mathbf{y}, \mathbf{v}, \mathbf{s}) &=\int_{4 \pi} \sum_{i \in\{R, M\}} \beta_{i}^{S}(\mathbf{y}) P_{i}(\mathbf{v} . \boldsymbol{\omega}) L(\mathbf{y}, \boldsymbol{\omega}, \mathbf{s}) d \omega \end{aligned} T(x,xo​)I[L](xo​,s)J[L](y,v,s)​=exp⎝⎛​−∫xxo​​i∈{R,M}∑​βie​(y)dy⎠⎞​=πα(xo​)​∫2π​L(x0​,ω,s)ω⋅n(xo​)dω, or 0=∫4π​i∈{R,M}∑​βiS​(y)Pi​(v.ω)L(y,ω,s)dω​

其中在大气层外 $\mathcal{L}$ 是0。最终的渲染等式如下

$$\begin{array}{rl}L(\mathbf{x},\mathbf{v},\mathbf{s})& =\left(L_0+\mathcal{R}[L]+\mathcal{S}[L]\right)(\mathbf{x},\mathbf{v},\mathbf{s})\\ L_0(\mathbf{x},\mathbf{v},\mathbf{s})& =T\left(\mathbf{x},{\mathbf{x}}_0\right)L_{sun},\text{ or }0\\ \mathcal{R}[L](\mathbf{x},\mathbf{v},\mathbf{s})& =T\left(\mathbf{x},{\mathbf{x}}_0\right)\mathcal{I}[L]\left({\mathbf{x}}_o,\mathbf{s}\right)\\ \mathcal{S}[L](\mathbf{x},\mathbf{v},\mathbf{s})& ={\int }_{\mathbf{x}}^{{\mathbf{x}}_o}T(\mathbf{x},\mathbf{y})\mathcal{J}[L](\mathbf{y},\mathbf{v},\mathbf{s})dy\end{array}$$
其中 $\mathcal{R}[L]$ 是在点 ${\mathbf{x}}_o$ 处反射并在到达$\mathbf{x}$之前被削减到的光；$\mathcal{S}[L]$ 则是 inscattered light

渲染方法

实现思路就是尽可能的预计算 $L$

因为星球是圆对称的，所以 $\mathbf{x}$ 和 $\mathbf{v}$ 可以被简化为 altitude 和 view zenith angle。则 $$

而 $L$ 可以表示为各阶散射之和
$$\begin{array}{rl}L& =L_0+(\mathcal{R}+\mathcal{S})\left[L_0\right]+(\mathcal{R}+\mathcal{S})\left[(\mathcal{R}+\mathcal{S})\left[L_0\right]\right]+\dots \\ & =L_0+L_1+L_2+\dots =L_0+L_{\ast }\end{array}$$

Zero and single scattering

通过渲染等式，我们可以准确的计算出 $L_0$ 和 $\mathcal{R}\left[{\mathcal{L}}_0\right]$

单次散射情形下，因为在阴影中是没有散射的，所有只能计算没有被阴影遮挡的空间，所以$\bar{L}$和$\overline{\mathcal{S}}\left[\bar{L}\right]$的计算被简化为4个参数(2个是$\mathbf{x}$,$\mathbf{v}$，另外2个是用于太阳方向$s$)

$$\begin{array}{rl}\mathcal{S}\left[L_0\right]& ={\int }_{\mathbf{x}}^{{\mathbf{x}}_s}T\mathcal{J}\left[{\bar{L}}_0\right]\\ & ={\int }_{\mathbf{x}}^{{\overline{\mathbf{x}}}_o}T\mathcal{J}\left[{\bar{L}}_0\right]-{\int }_{{\mathbf{x}}_s}^{{\overline{\mathbf{x}}}_o}T\mathcal{J}\left[{\bar{L}}_0\right]\end{array}$$

最终我们得到的单次散射的计算式，使用了两个预计算式，$T$是2个参数，$\overline{\mathcal{S}}\left[{\bar{L}}_0\right]$是4个参数
$$\mathcal{S}\left[L_0\right](\mathbf{x},\mathbf{v},\mathbf{s})=\overline{\mathcal{S}}\left[{\bar{L}}_0\right](\mathbf{x},\mathbf{v},\mathbf{s})-T\left(\mathbf{x},{\mathbf{x}}_s\right)\overline{\mathcal{S}}\left[{\bar{L}}_0\right]\left({\mathbf{x}}_s,\mathbf{v},\mathbf{s}\right)$$

Multiple scattering

$L_2+\dots =\mathcal{R}\left[L_{\ast }\right]+\mathcal{S}\left[L_{\ast }\right]$

因为相对于单次散射，多次散射的值更小，所以地表反射的贡献就直接忽略掉，所以在计算积分的时候就没有考虑遮蔽了

$$\mathcal{S}\left[L_{\ast }\right]\simeq {\int }_{\mathbf{x}}^{{\mathbf{x}}_s}T\mathcal{J}\left[{\bar{L}}_{\ast }\right]$$

我们在 $\mathcal{R}\left[L_{\ast }\right]$ 使用水平半球的地表切面 $\frac{1+\mathbf{n}\cdot \mathbf{n}}{2}$ 导致的遮蔽效果
$$\begin{array}{l}\hat{\mathcal{R}}\left[{\bar{L}}_{\ast }\right]=T\left(\mathbf{x},{\mathbf{x}}_0\right)\frac{\alpha \left({\mathbf{x}}_0\right)}{\pi }\frac{1+\mathbf{n}\left({\mathbf{x}}_o\right)\cdot \overline{\mathbf{n}}\left({\mathbf{x}}_o\right)}{2}\overline{\mathcal{E}}\left[{\bar{L}}_{\ast }\right]\left({\mathbf{x}}_o,\mathbf{s}\right)\\ \\ \overline{\mathcal{E}}\left[{\bar{L}}_{\ast }\right]\left({\mathbf{x}}_o,\mathbf{s}\right)={\int }_{2\pi }{\bar{L}}_{\ast }\left({\mathbf{x}}_o,\mathbf{\omega },\mathbf{s}\right)\mathbf{\omega }\cdot \overline{\mathbf{n}}\left({\mathbf{x}}_0\right)d\omega ,\text{ or }0\end{array}$$

按照单次散射的形式重写得到

$$L\simeq L_0+\mathcal{R}\left[L_0\right]+\hat{\mathcal{R}}\left[{\bar{L}}_{\ast }\right]+\overline{\mathcal{S}}[\bar{L}{]}_{\mid \mathbf{x}}-T\left(\mathbf{x},{\mathbf{x}}_s\right)\overline{\mathcal{S}}[\bar{L}{]}_{\mid {\mathbf{x}}_s}$$

预计算

1. 预计算 Transmittance 的 2D 纹理，对于所有的$\mathbf{x}$和$\mathbf{v}$计算出 $T\left(\mathbf{x},{\overline{\mathbf{x}}}_o(\mathbf{x},\mathbf{v})\right)$，即 $\mathbb{T}(\mathbf{x},\mathbf{v})$ 。
2. 预计算 $\overline{\mathcal{E}}\left[{\bar{L}}_{\ast }\right]$ 和 $\overline{\mathcal{S}}[\bar{L}]$ 到两个纹理 $\mathbb{E}$ 和 $\mathbb{S}$ 中，并在每次循环中计算散射 ${\bar{L}}_i$，把 $\Delta \mathbb{E}$ 和 $\Delta \mathbb{S}$ 加到 $\mathbb{E}$ 和 $\mathbb{S}$ 中，最后地表反射 $\overline{\mathcal{R}}$ 可以通过 $\overline{\mathcal{R}}[L](\mathbf{x},\mathbf{v},\mathbf{s})=T\left(\mathbf{x},{\overline{\mathbf{x}}}_o\right)\frac{\bar{\alpha }}{\bar{\pi }}\overline{\mathcal{E}}[L]\left({\overline{\mathbf{x}}}_o,\mathbf{s}\right)$ 计算得到。
3. 这个算法使用到了三个临时纹理 $\Delta \mathbb{E}$，$\Delta \mathbb{S}$ 和 $\Delta \mathbb{J}$，分别包含的是第 $i$ 次迭代的 $\overline{\mathcal{E}}\left[{\bar{L}}_i\right]$，$\overline{\mathcal{S}}\left[{\bar{L}}_i\right]$ 以及 $\mathcal{J}\left[{\bar{L}}_i\right]$

关于参数化和把4D的Scattering Texture转换为3D的纹理，建议直接看作者给出的代码实现和说明

[1] https://ebruneton.github.io/precomputed_atmospheric_scattering/

[2] https://hal.inria.fr/inria-00288758/document