基于Unity的实现 https://gitee.com/alienity/atmospheric-scattering
Paper中说明的大都是原理,在作者给出的具体实现中有较多的优化,有兴趣的朋友可以在看完原理后,直接看作者的实现及其说明。本文实现的代码也是基于作者的代码
物理模型
物理模型中大气包含空气分子和气溶胶粒子两种介质,一般是在 R g = 6030 k m R_{g}=6030 km Rg=6030km 到 R t = 6420 k m R_{t}=6420km Rt=6420km 之间
对应这两种介质,分别会发生瑞利散射和Mie散射。瑞利散射发生的条件是粒子半径远小于入射光波长。在地球上大气的区间高度是
R
g
=
6360
k
m
R_{g}=6360km
Rg=6360km 和
R
t
=
6420
k
m
R_{t}=6420km
Rt=6420km
瑞利散射
在每个点上,与入射光偏离 θ \theta θ 角度后的强度等于散射系数 β s \beta^{s} βs 和相函数 P P P 的乘积。散射系数依赖于分子密度,相函数描述了角度依赖。
β R S ( h , λ ) = 8 π 3 ( n 2 − 1 ) 2 3 N λ 4 e − h H R P R ( μ ) = 3 16 π ( 1 + μ 2 ) where μ = cos θ \begin{aligned} \beta_{R}^{S}(h, \lambda) &=\frac{8 \pi^{3}\left(n^{2}-1\right)^{2}}{3 N \lambda^{4}} e^{-\frac{h}{H_{R}}} \\ P_{R}(\mu) &=\frac{3}{16 \pi}\left(1+\mu^{2}\right) \quad \text { where } \mu=\cos \theta \end{aligned} βRS(h,λ)PR(μ)=3Nλ48π3(n2−1)2e−HRh=16π3(1+μ2) where μ=cosθ
其中 h = e − R g h = e - R_{g} h=e−Rg 表示的是海拔高度, λ \lambda λ 是波长, n n n 是空气折射系数, N N N 是在海平面 R g R_{g} Rg 的分子数密度, H R = 8 k m H_{R}=8km HR=8km 是大气厚度,其中对于波长 λ = ( 680 , 550 , 440 ) n m \lambda=(680,550,440)nm λ=(680,550,440)nm 的散射系数 β R S = ( 5.8 , 13.5 , 33.1 ) 1 0 − 6 m − 1 \beta_{R}^{S}=(5.8,13.5,33.1)10^{-6}m^{-1} βRS=(5.8,13.5,33.1)10−6m−1
Mie散射
气溶胶也有呈指数下降的密度,只是 height scale H M ≃ 1.2 k m H_{M} \simeq 1.2 km HM≃1.2km
Cornette-Shanks 相函数
β
M
S
(
h
,
λ
)
=
β
M
s
(
0
,
λ
)
e
−
h
H
M
P
M
(
μ
)
=
3
8
π
(
1
−
g
2
)
(
1
+
μ
2
)
(
2
+
g
2
)
(
1
+
g
2
−
2
g
μ
)
3
/
2
\begin{aligned} \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}-2 g \mu\right)^{3 / 2}} \end{aligned}
βMS(h,λ)PM(μ)=βMs(0,λ)e−HMh=8π3(2+g2)(1+g2−2gμ)3/2(1−g2)(1+μ2)
与分空气分子不一样,气溶胶吸收一定比例的光,吸收系数是 β M a \beta_{M}^{a} βMa,所以吸光系数是 β M e = β M s + β M a \beta_{M}^{e}=\beta_{M}^{s}+\beta_{M}^{a} βMe=βMs+βMa
渲染等式
L ( x , v , s ) L(\mathbf{x}, \mathbf{v}, \mathbf{s}) L(x,v,s) 是太阳在 s \mathbf{s} s 方向时,从方向 v \mathbf{v} v 到达 x \mathbf{x} x 的 radiance of light
- 点 x o \mathbf{x}_{o} xo 与点 x \mathbf{x} x 之间的 transmittance T T T;
- 在点 x o \mathbf{x}_{o} xo 处反射的 radiance L \mathcal{L} L;
- 在 y \mathbf{y} y 点朝方向 − v -\mathbf{v} −v 散射的radiance J \mathcal{J} 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⎝⎛−∫xxoi∈{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ω
其中在大气层外 L \mathcal{L} L 是0。最终的渲染等式如下
L
(
x
,
v
,
s
)
=
(
L
0
+
R
[
L
]
+
S
[
L
]
)
(
x
,
v
,
s
)
L
0
(
x
,
v
,
s
)
=
T
(
x
,
x
0
)
L
s
u
n
,
or
0
R
[
L
]
(
x
,
v
,
s
)
=
T
(
x
,
x
0
)
I
[
L
]
(
x
o
,
s
)
S
[
L
]
(
x
,
v
,
s
)
=
∫
x
x
o
T
(
x
,
y
)
J
[
L
]
(
y
,
v
,
s
)
d
y
\begin{aligned} 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_{s u n}, \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}) d y \end{aligned}
L(x,v,s)L0(x,v,s)R[L](x,v,s)S[L](x,v,s)=(L0+R[L]+S[L])(x,v,s)=T(x,x0)Lsun, or 0=T(x,x0)I[L](xo,s)=∫xxoT(x,y)J[L](y,v,s)dy
其中
R
[
L
]
\mathcal{R}[L]
R[L] 是在点
x
o
\mathbf{x}_{o}
xo 处反射并在到达
x
\mathbf{x}
x之前被削减到的光;
S
[
L
]
\mathcal{S}[L]
S[L] 则是 inscattered light
渲染方法
实现思路就是尽可能的预计算 L L L
因为星球是圆对称的,所以 x \mathbf{x} x 和 v \mathbf{v} v 可以被简化为 altitude 和 view zenith angle。则 $$
而
L
L
L 可以表示为各阶散射之和
L
=
L
0
+
(
R
+
S
)
[
L
0
]
+
(
R
+
S
)
[
(
R
+
S
)
[
L
0
]
]
+
…
=
L
0
+
L
1
+
L
2
+
…
=
L
0
+
L
∗
\begin{aligned} 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]+\ldots \\ &=L_{0}+L_{1}+L_{2}+\ldots=L_{0}+L_{*} \end{aligned}
L=L0+(R+S)[L0]+(R+S)[(R+S)[L0]]+…=L0+L1+L2+…=L0+L∗
Zero and single scattering
通过渲染等式,我们可以准确的计算出 L 0 L_{0} L0 和 R [ L 0 ] \mathcal{R\left[L_{0}\right]} R[L0]
单次散射情形下,因为在阴影中是没有散射的,所有只能计算没有被阴影遮挡的空间,所以 L ˉ \bar{L} Lˉ和 S ‾ [ L ˉ ] \overline{\mathcal{S}}\left[\bar{L}\right] S[Lˉ]的计算被简化为4个参数(2个是 x \mathbf{x} x, v \mathbf{v} v,另外2个是用于太阳方向 s s s)
S [ L 0 ] = ∫ x x s T J [ L ˉ 0 ] = ∫ x x ‾ o T J [ L ˉ 0 ] − ∫ x s x ‾ o T J [ L ˉ 0 ] \begin{aligned} \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{aligned} S[L0]=∫xxsTJ[Lˉ0]=∫xxoTJ[Lˉ0]−∫xsxoTJ[Lˉ0]
最终我们得到的单次散射的计算式,使用了两个预计算式,
T
T
T是2个参数,
S
‾
[
L
ˉ
0
]
\overline{\mathcal{S}}\left[\bar{L}_{0}\right]
S[Lˉ0]是4个参数
S
[
L
0
]
(
x
,
v
,
s
)
=
S
‾
[
L
ˉ
0
]
(
x
,
v
,
s
)
−
T
(
x
,
x
s
)
S
‾
[
L
ˉ
0
]
(
x
s
,
v
,
s
)
\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)
S[L0](x,v,s)=S[Lˉ0](x,v,s)−T(x,xs)S[Lˉ0](xs,v,s)
Multiple scattering
L 2 + … = R [ L ∗ ] + S [ L ∗ ] L_{2}+\ldots=\mathcal{R}\left[L_{*}\right]+\mathcal{S}\left[L_{*}\right] L2+…=R[L∗]+S[L∗]
因为相对于单次散射,多次散射的值更小,所以地表反射的贡献就直接忽略掉,所以在计算积分的时候就没有考虑遮蔽了
S [ L ∗ ] ≃ ∫ x x s T J [ L ˉ ∗ ] \mathcal{S}\left[L_{*}\right] \simeq \int_{\mathbf{x}}^{\mathbf{x}_{s}} T \mathcal{J}\left[\bar{L}_{*}\right] S[L∗]≃∫xxsTJ[Lˉ∗]
我们在
R
[
L
∗
]
\mathcal{R}\left[L_{*}\right]
R[L∗] 使用水平半球的地表切面
1
+
n
⋅
n
2
\frac{1+\mathbf{n}\cdot\mathbf{n}}{2}
21+n⋅n 导致的遮蔽效果
R
^
[
L
ˉ
∗
]
=
T
(
x
,
x
0
)
α
(
x
0
)
π
1
+
n
(
x
o
)
⋅
n
‾
(
x
o
)
2
E
‾
[
L
ˉ
∗
]
(
x
o
,
s
)
E
‾
[
L
ˉ
∗
]
(
x
o
,
s
)
=
∫
2
π
L
ˉ
∗
(
x
o
,
ω
,
s
)
ω
⋅
n
‾
(
x
0
)
d
ω
,
or
0
\begin{array}{l} \hat{\mathcal{R}}\left[\bar{L}_{*}\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}_{*}\right]\left(\mathbf{x}_{o}, \mathbf{s}\right) \\\\ \overline{\mathcal{E}}\left[\bar{L}_{*}\right]\left(\mathbf{x}_{o}, \mathbf{s}\right)=\int_{2 \pi} \bar{L}_{*}\left(\mathbf{x}_{o}, \boldsymbol{\omega}, \mathbf{s}\right) \boldsymbol{\omega} \cdot \overline{\mathbf{n}}\left(\mathbf{x}_{0}\right) d \omega, \text { or } 0 \end{array}
R^[Lˉ∗]=T(x,x0)πα(x0)21+n(xo)⋅n(xo)E[Lˉ∗](xo,s)E[Lˉ∗](xo,s)=∫2πLˉ∗(xo,ω,s)ω⋅n(x0)dω, or 0
按照单次散射的形式重写得到
L ≃ L 0 + R [ L 0 ] + R ^ [ L ˉ ∗ ] + S ‾ [ L ˉ ] ∣ x − T ( x , x s ) S ‾ [ L ˉ ] ∣ x s L \simeq L_{0}+\mathcal{R}\left[L_{0}\right]+\hat{\mathcal{R}}\left[\bar{L}_{*}\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}} L≃L0+R[L0]+R^[Lˉ∗]+S[Lˉ]∣x−T(x,xs)S[Lˉ]∣xs
预计算
- 预计算 Transmittance 的 2D 纹理,对于所有的 x \mathbf{x} x和 v \mathbf{v} v计算出 T ( x , x ‾ o ( x , v ) ) T\left(\mathbf{x}, \overline{\mathbf{x}}_{o}(\mathbf{x}, \mathbf{v})\right) T(x,xo(x,v)),即 T ( x , v ) \mathbb{T}(\mathbf{x}, \mathbf{v}) T(x,v) 。
- 预计算 E ‾ [ L ˉ ∗ ] \overline{\mathcal{E}}\left[\bar{L}_{*}\right] E[Lˉ∗] 和 S ‾ [ L ˉ ] \overline{\mathcal{S}}[\bar{L}] S[Lˉ] 到两个纹理 E \mathbb{E} E 和 S \mathbb{S} S 中,并在每次循环中计算散射 L ˉ i \bar{L}_{i} Lˉi,把 Δ E \Delta\mathbb{E} ΔE 和 Δ S \Delta\mathbb{S} ΔS 加到 E \mathbb{E} E 和 S \mathbb{S} S 中,最后地表反射 R ‾ \overline{\mathcal{R}} R 可以通过 R ‾ [ L ] ( x , v , s ) = T ( x , x ‾ o ) α ˉ π ˉ E ‾ [ L ] ( x ‾ o , s ) \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) R[L](x,v,s)=T(x,xo)πˉαˉE[L](xo,s) 计算得到。
- 这个算法使用到了三个临时纹理 Δ E \Delta \mathbb{E} ΔE, Δ S \Delta \mathbb{S} ΔS 和 Δ J \Delta \mathbb{J} ΔJ,分别包含的是第 i i i 次迭代的 E ‾ [ L ˉ i ] \overline{\mathcal{E}}\left[\bar{L}_{i}\right] E[Lˉi], S ‾ [ L ˉ i ] \overline{\mathcal{S}}\left[\bar{L}_{i}\right] S[Lˉi] 以及 J [ L ˉ i ] \mathcal{J}\left[\bar{L}_{i}\right] J[Lˉi]
关于参数化和把4D的Scattering Texture转换为3D的纹理,建议直接看作者给出的代码实现和说明
[1] https://ebruneton.github.io/precomputed_atmospheric_scattering/