【GDC 2015】An Introduction to Rea

作者: 离原春草 | 来源:发表于2024-05-31 00:12 被阅读0次

今天学习的是GDC 2015上分享的基于FFT的海洋渲染方案，其中通过参考知乎上的相关资料，对基于FFT的海洋模拟实现逻辑做了相对清晰的推导与总结：

通过贴图来给出某个2D的频域信号，其中贴图的uv代表两个维度的频率，数值代表该信号的振幅
基于上述频域信号，可以通过逆变换得到随时间变化的2D空域信号，这个变换过程需要在运行时完成（涉及到时间，当然也可以烘焙多张，保证首尾衔接的话即可）
频域信号用的是菲利普频谱，基于这个频谱可以在运行时计算各个频率的信号数值，也可以在离线完成计算。运行时计算的好处是：

可以根据当前环境的状态比如风向等来调制得到更拟真的效果
可以运行时调整模拟参数，实现海洋效果的动态变化

具体可以参考实现细节：

page01

page02

page03

游戏中的海洋效果：有刺客信条III以及刺客信条黑帆

page04

Crysis的海洋效果

page05

War Thunder的海洋效果

page06

海洋奇缘

page07

泰坦尼克号

page08

这些游戏跟电影有如下的一些共同点

都是基于FFT来完成模拟的
模拟的输入是一个频谱，之后对这个频谱做逆向的离散傅里叶变换就得到时域（空域）信号
时域信号是一个随时间变化的高度displacement数据

page09

为啥是FFT：DFT耗时太高了，时间复杂度是 $O(n^2)$ ，而FFT的时间复杂度则是 $O(nlog_2(n))$

在GPU中计算FFT的算法是Cooley-Tukey（蝴蝶）算法。

page10

傅里叶时域（还是频域？）贴图的分辨率应该是多少呢？正常情况下128 ~ 512就够了，泰坦尼克号跟Waterworld也就只用来2048，超过这个尺寸可能会导致浮点精度问题。

这里选用的分辨率是512。

page11

傅里叶变换的公式（可以参考[4]了解一下具体的推导过程）是：
$F(\omega) = \int_{-\infty}^{\infty} f(t) e^{-i\omega t}dt$

逆变换公式为：
$f(t) = \frac{1}{2\pi}\int_{-\infty}^{\infty} F(\omega) e^{i\omega t}d\omega$

在计算机中这个积分公式是没有解析解的，只能转换为离散的累加形式。

page12

再来看离散形式的变换与逆变换，这里的N是说将 $2\pi$ 分割成N份，N越大，我们得到的频域信号就越丰富，信号就越逼真，模拟越接近，但是计算复杂度也就越高：
$F(k) = \sum_{t=0}^{N-1} f(t) e^{-i2\pi \frac{k}{N} t}$

$f(t) = \sum_{k=0}^{N-1} F(k) e^{i2\pi \frac{k}{N} t}$

page13

参考[1]，DFT的计算其实就是一个矩阵的运算，而这个矩阵运算则可以通过排列得到如上图所示的蝴蝶累加算法，从而可以很好的实现加速计算。

page14

这里IFFT的实现是通过compute shader完成的，采用的是Radix-2算法（参考[3]），细节就不说了，网上有大把的开源实现。

page15

上面介绍的是一个一维函数的FFT实现，那么回到水体模拟上，这个公式是如何应用起来的呢？

我们知道，傅里叶变换其实是通过将信号拆解为正弦波的方式来完成对信号的表达的，每个正弦波都有一个频率，一个信号的频域表示，其实就是其拆解后得到的正弦波的振幅与频率的表示，频域信号在某个频率上的取值就是当前对应的时（空）域信号在对应频率上的振幅。

回到海洋这边来，其实海浪可以看成是由空间上无穷个2D的正弦波信号叠加而成，同样的，这个信号在频域中的表示就是海浪在不同频率的正弦波上的振幅。不同的是，这个频率是2D的，有X&Z两个方向的频率，同样的，我们的频域信号也是2维的，即在X & Z上具有不同频率的信号的振幅。

用公式来表达，先不考虑海洋波形随着时间波动的变化，我们先看下某个时间点的海洋波形，这个波形刚说了，可以看成是无穷多个2D的正弦波组成的，在空域上，可以表示为h(x, z)，即在坐标(x, z)处的水体高度。

这个水波对应的频域信号，我们可以表示为 $\tilde{h}(w_x, w_z)$ ，两者可以通过傅里叶变换（空域->频域）与逆变换（频域->空域）实现转换，在海洋的实现逻辑中，我们的输入通常是频域信号，且只能用IDFT的方式来完成计算，因此按照前面的说明，这个逆变换的公式给出如下：
$h(\vec{x}) = \sum_{\vec{k}} \tilde{h}(\vec{k}) e^{i\vec{k}\vec{x}}$

其中 $\vec{k}$ 是频域上的2D向量，对应于前面的 $(w_x, w_z)$ ，或者用k做前缀，将之改写为 $\vec{k} = (k_x, k_z)$ ，根据我们FFT频域模拟的尺寸，这里有 $k_x = \frac{2\pi n}{L_x}, k_z = \frac{2\pi m}{L_z}$ ，n跟m分别是频域贴图（贴图的uv代表着两个方向的频率，对应的数值对应的是这个频率的幅度）的水平跟垂直的分辨率，也就是我们在频谱上的水平跟垂直方向上以 $\frac{2\pi}{L_x},\frac{2\pi}{L_z}$ 为单位采样的点数，其中 $L_x, L_z$ 又分别为海面X or Z方向上的patch尺寸，根据[4]中的解释，我们的采样点是分布在原点两侧的

$m \in \{-\frac{M}{2}, -\frac{M}{2} +1, ... , \frac{M}{2} - 1\} \\ n \in \{-\frac{N}{2}, -\frac{N}{2} +1, ... , \frac{N}{2} - 1\}$

在这样的采点策略下，针对海洋波形，我们就总计拿到了MxN个频域的信号（或者说正弦波），通过他们之间的叠加与逆变换得到海洋的波形，通常情况下，这两个数值越大，波形精度越高，模拟越逼真。

$\vec{x}$ 则是空域上的2D信号，对应的是当前点的xz坐标，同样的，在计算displacement的时候，我们也是针对离散点进行的，即 $\vec{x}$ 也是在海平面上以 $\frac{L_x}{N},\frac{L_z}{M}$ 为单位采样NxM个点（虽然通常N=M），同样的，在相同的覆盖范围下，这俩数值越大，顶点之间的密度就越大，水体精度就越高，当然计算消耗也同样会增加。

另外需要注意的是，指数上的两个向量的乘法是点乘，即有：
$\vec{k}\vec{x} = k_x x + k_z z$
在这个公式下我们就有当某个维度固定住的时候，我们得到的就是一个一维的信号。

在这个基础上，如果我们考虑水波随着时间变化的过程，那么我们可以将前面的频域信号改写为时间的函数，即有：

$h(\vec{x},t) = \sum_{\vec{k}} \tilde{h}(\vec{k}, t) e^{i\vec{k}\vec{x}}$

page16

接下来对这个公式做一个拆解。

page17

里面最大的问题是频域信号怎么得到，这里就没有介绍具体的背景，而是直接给出了结论，海洋波形的频域信号通常使用的是菲利普频谱，其公式为：
$\tilde{h}(\vec{k}, t) = \tilde{h_0}(\vec{k})e^{i\omega(k)t}+\tilde{h_0}^*(-\vec{k})e^{-i\omega(k)t}$

这个公式中 $\tilde{h_0}$ 跟 $\tilde{h_0}^*$ 是共轭复数（两个实部相等，虚部互为相反数的复数互为共轭复数），k则是 $\vec{k}$ 的模， $\omega(k)$ 是一个函数，其计算公式为：
$\omega(k)=\sqrt{gk}$

这里的g是重力加速度。

page18

再来看 $\tilde{h_0}(\vec{k})$ ，有如下的计算公式：
$\tilde{h_0}(\vec{k}) = \frac{1}{\sqrt 2}(\xi_r + i \xi_i)\sqrt{P_h(\vec{k})}$

$\xi_r, \xi_i$ 是相互独立的随机数，均服从均值为0，标准差为1的正态分布。

page19

$P_h(\vec{k}) = a\frac{e^{-\frac{1}{(kL)^2}}}{k^4}|\vec{k}\cdot\vec{\omega}|^2$

这个公式中的 $\vec{\omega}$ 是归一化后的风向（2D），L的计算公式为：
$L=\frac{V^2}{g}$

这里的g依然是重力加速度，V则是风速。

按照[4]的描述，这里的公式是基于统计给出的经验公式，并无过多的解析或物理含义在里面，可以不用过于纠结其来源。

page20

page21

page22

这里对公式中各个变量的含义做了基本的介绍，这里补上前面遗漏的项，前面介绍过的这里就不再介绍：

a是全局的海浪的振幅

page23

基于高度偏移的公式 $h(\vec{x},t) = \sum_{\vec{k}} \tilde{h}(\vec{k}, t) e^{i\vec{k}\vec{x}}$ ，我们可以通过解析的方式得到其法线。

先来计算一下梯度向量（某个函数即这里的高度变化率最大的方向，因为这里的变量还都是水平方向的两个轴，因此这里的向量还是一个水平方向的向量），因为 $\tilde{h}$ 不包含梯度计算相关的 $\vec{x}$ 信息，因此可以直接移出来：
$\nabla h(\vec{x},t) = \sum_{\vec{k}} \tilde{h}(\vec{k}, t) \nabla e^{i\vec{k}\vec{x}}$

继续推导，我们有：
$\nabla e^{i\vec{k}\vec{x}} = (\frac{\partial {e^{i(k_x x+k_z z)}}}{\partial x}, \frac{\partial {e^{i(k_x x+k_z z)}}}{\partial z}) \\ =e^{i(k_x x+k_z z)} * i(k_x, k_z) \\ = i\vec k e^{i \vec k \vec x}$

代入前面的公式，我们就可以得到梯度向量的计算公式：
$\nabla h(\vec{x},t) = \sum_{\vec{k}} i\vec k \tilde{h}(\vec{k}, t) e^{i \vec k \vec x}$

针对垂直向上的向量、梯度向量与法线向量的关系[4]:

注意梯度向量本身是没有归一化的，我们可以有： $\vec N = \vec up - \vec {grad}$ ，将前面的公式代入，则我们有：
$\vec N = (0, 1, 0) - (\nabla h_x(\vec{x},t), 0, \nabla h_z(\vec{x},t)) \\ = ( - \nabla h_x(\vec{x},t), 1, - \nabla h_z(\vec{x},t))$

page24

前面得到的是高度方向的偏移，但是我们实际上需要的是三个方向的偏移，因此还需要再额外处理一下另外两个方向的。

在Gerstner Wave的实现中，为了得到尖锐的波峰效果（choppy wave），需要基于如下公式[4]在XZ方向上做挤压：
$P(x,y,t)= \begin{cases} x+\sum Q_i A_i \times D_i.x \times cos(w_iD_i(x, y)+\phi_i t) \\ y+\sum Q_i A_i \times D_i.y \times cos(w_iD_i(x, y)+\phi_i t) \\ \sum A_i sin(w_iD_i \cdot (x, y)+\phi_i t) \end{cases}$

而同样的，FFT实现的海浪效果也需要做挤压，挤压公式为：
$\vec D(\vec x, t) = \sum_{\vec k} - i \frac{\vec k}{k}\tilde h(\vec k, t)e^{i \vec k \vec x}$

page31

page32

如下图[4]所示，当挤压力度过大的时候，就会出现水波的相互穿刺，即可以理解为水波碰撞在一起，此时会产生浮沫：

这个过程其实就是某个形状（比如上图中的近似四边形，通常称之为一个图元）发生了翻转，正面翻到了背面，用数学概念来描述，就是该图元有向面积变为负数，而这个在数学公式上来表达就是对应的雅可比矩阵计算结果为负数，具体可以参考[4]的推导。

page25

挤压强度假设为 $\lambda$ ，我们就有：
$\vec x^{'} = \vec x + \lambda *\vec D(\vec x, t)$

而这个公式跟前面Gerstner Wave的挤压公式整体形式上其实是一样的，只是换了一种写法。

经过上面计算，我们就得到了空域中各个点的偏移向量 $\vec D(x, y, z)$ ，这个可以用一张3通道的贴图来存储，我们就有了上图所示的Displacement Map

page26

参考NVIDIA的解说图，我们回顾一下上面的计算逻辑：

我们基于频谱函数或者说频谱贴图 $P_h(\vec k)$ 就能得到高度偏移 $D_y$
同样基于频谱函数跟挤压公式，我们可以得到在水平方向上的偏移 $D_x, D_z$
将上述三者合到一起，我们就有了displacement map
基于Displacement map，我们就有计算出法线贴图，同时基于雅可比矩阵计算结果，就能够计算出foam intensity map（也就是上图中的folding）