虚幻引擎 Gerstner Waves -GPU Gems 从物理模型中实现有效的水体模拟

这篇文章重点在于结合GPU Gems一书中有关Gerstner Waves 的数学公式，在虚幻引擎中复现正确的Gerstner Waves和正确的法线

文中内容整理自书中，并附带我的理解，与在虚幻引擎中的实现，可以参考原文看这篇文章，原文网上很多，我这里就不转载了

内容看上去有些乱？
与虚幻无关的内容，以普通文本的形式体现，需要强调的、与虚幻有关的，我将紧挨着内容以引用块的形式附上（就像这段话）
这是因为我希望读者能随着我一起看看这本书，书讲到哪，咱们做到哪

书中还介绍了几种波，我假设你们做了课前预习，本文直接开始G波

1.2.1选择波形

我们需要一组参数来定义每个波形。这些参数包括：

• Wavelength _波长（ $L$ ）
  • 世界空间中波峰到波峰之间的距离。
  • 波长 $L$ 与角频率 $w$ 之间的关系为 $w=\frac{2\pi }{L}$。

这里将会出现第一个坑，如果你看的是英伟达的来源，他们在公式中弄丢了这个π，这篇文章中的公式将使用原书公式

• Amplitude _振幅（ $A$ ）

  • 从水平面到波峰的高度。

• Speed _速度（ $S$ ）

  • 波峰每秒前进的距离。 将速度表示为相位常数 $\phi$ 更为方便。
  • $\phi =S\times \frac{2\pi }{L}$

• Direction _方向（ $D$ ）

  • 垂直于波阵面的水平向量，波峰沿着波阵面运动。

为了在场景的动态中提供变化，我们将在一定的约束条件下随机生成这些波浪参数。随着时间的推移，我们将持续地将一个波浪渐隐，然后再以一组不同的参数将其渐显。
事实证明，这些参数是相互依赖的。必须仔细的为每个波浪生成一整套参数，这些参数需要以一种令人信服的方式组合在一起。

新建一个GerstnerWaves函数，我们先把上方出现的变量在input预留:

以防有疑问，最后一个Node是“重路由” ，它是用来避免蓝图变意大利面的东西，无功能上的意义

1.2.2 法线和切线

因为表面是显式函数，所以我们可以直接计算任何点的表面方向，而不需要依赖有限差分技术。
副法线向量 $B$ 和切线向量 $T$ 是分别 $x$ 和 $y$ 方向的偏导数。对于2D水平面中的任何点 ( $x$ , $y$ )，表面上的三维位置 $P$ 是：

Equation 6a _公式6a

$$\mathbf{N}(x,y)=\mathbf{B}(x,y)\times \mathbf{T}(x,y)$$

这个公式体现N与B、T的关系。但我们先跳过这个，最后再来处理法线，这将规避经典的法线混合问题

格斯特纳波^{GerstnerWaves}

为了有效的模拟，我们需要控制波浪的陡峭程度。如前所述，正弦波呈现出圆润的外观——这可能正是我们想要的平静、田园诗般的池塘效果。但对于粗犷的海面，我们需要形成更尖锐的波峰和更宽阔的波谷。我们可以使用公式8a和8b来实现所需的形状，但我们选择了相关的格斯特纳波。
GerstnerWaves早在有计算机图形学之前就被开发出来，用于在物理基础上模拟海水。因此GerstnerWaves提供了一些表面的微妙运动，这些变化不明显但是非常可信。（详细描述参见Tessendorf 2001）。
我们选择GerstnerWaves，因为它们有一个常被忽视的特性：将，这正是我们希望顶点集中的地方，如图1-5所示。

这是GerstnerWaves函数：

Equation 9 _公式9

$$\mathbf{P}\left(x,y,t\right)=\left(\begin{array}{l}\begin{array}{rlrl}& x& +& \sum \left(Q_i{A}_i\times {\mathbf{D}}_i.x\times \cos \left(w_i{\mathbf{D}}_i\cdot (x,y)+{\phi }_{i}t\right)\right),\\ & y& +& \sum \left(Q_i{A}_i\times {\mathbf{D}}_i.y\times \cos \left(w_i{\mathbf{D}}_i\cdot (x,y)+{\phi }_{i}t\right)\right),\\ & & & \sum \left(A_i\sin \left(w_i{\mathbf{D}}_i\cdot (x,y)+{\phi }_{i}t\right)\right)\end{array}\end{array}\right)$$

∑为求和，∑( )也就是所有波加在一起，然后再在xy通道分别加x和y。
这里我们想要先实现一个波，也就是∑( )的内容，那么公式可以写为：

Px = Q×A×Dx×cos(w×D·(x,y) + φ×t);
Py = Q×A×Dy×cos(w×D·(x,y) + φ×t);
Pz =      A       ×sin(w×D·(x,y) + φ×t);

能看到我们想实现这个公式还需要两个变量， $xy$ 和 $t$，他们分别是：

公式中的 $Q_i$ 是一个控制波浪陡峭程度的参数。
对于单个波浪 $i$，$Q_i=0$ 产生常见的滚动正弦波，而 $Q_i=\frac{1}{w_i{A}_i}$ 产生尖锐的波峰。应避免使用较大的 $Q_i$ 值，因为它们会在波峰上方形成"环\卷"。
实际上，我们可以将 $Q$ 作为“陡峭程度”参数留给制作艺术家来指定，允许范围是 $0$ 到 $1$，并使用 $Q_i=\frac{Q}{w_i{A}_i\times \text{numWaves}}$ 来变化，从完全平滑的波浪到我们能产生的最尖锐的波浪。

Q的计算为：

其中 $numWaves$ 为波的总数，例如水面由3个GerstnerWaves组合，这里就是3
Steepness变量就是文中说的"陡峭程度"参数

可以看到公式里的( $w\times D\cdot (x,y)+\phi \times t$ )是重复的,我们先把他的SinCos算出来
(这不会影响着色器复杂度，只会让我的截图好看，避免意大利面)
注意这里的 $sin$ 和 $cos$ 的周期为 $2pi$


现在我们可以计算公式了

这就是顶点位移的结果，接下来我们计算法线

值得注意的是，公式3和公式9之间唯一的区别是顶点的横向移动。他们的高度是相同的。这意味着我们不再有一个严格的高度函数。即，$\mathbf{P}(x,y,t).x\ne x$ 。然而，该函数仍然容易求导，并且有一些项可以方便的消去。

其中:
$WA=w_i\times A_i,$
$S()=\sin \left(w_i\times {\mathbf{D}}_i\cdot \mathbf{P}+{\phi }_{i}t\right)$
$C()=\cos (w_i\times {\mathbf{D}}_i\cdot \mathbf{P}+{\phi }_{i}t)$

又一个万人坑，这里的 $P$ 实际上指的是 ( $x$，$y$ )

$S()$ 和 $C()$ 我在蓝图里写成了 $Sin()$ 和 $Cos()$

求导后得到切线空间的基础向量是：

Equation 10 _公式10

$$\mathbf{B}=\left(\begin{array}{c}\begin{array}{rlrl}& 1& -& \sum \left(Q_i\times {\mathbf{D}}_i.x^2\times WA\times S()\right),\\ & & -& \sum \left(Q_i\times {\mathbf{D}}_i.x\times {\mathbf{D}}_i.y\times WA\times S()\right),\\ & & & \sum \left({\mathbf{D}}_i.x\times WA\times C()\right)\end{array}\end{array}\right)$$

Binormal计算可以写公式为：(注意这里仍然先实现一个波，因此计算公式的∑的内容)

Bx = Q×D.x^2        ×WA×sin();
By = Q×D.x    ×D.y×WA×sin();
Bz =       Dx             ×WA×cos();

Equation 11 _公式11

$$\mathbf{T}=\left(\begin{array}{c}\begin{array}{rlrl}& & -& \sum (Q_i\times {\mathbf{D}}_i.x\times {\mathbf{D}}_i.y\times WA\times S()),\\ & 1& -& \sum (Q_i\times {\mathbf{D}}_i.y^2\times WA\times S()),\\ & & & \sum ({\mathbf{D}}_i.y\times WA\times C())\end{array}\end{array}\right)$$

Tangent:

Tx = Q×D.x×D.y×WA×sin();
Ty = Q×D.y^2     ×WA×sin();
Tz =      Dy          ×WA×cos();

Equation 12 _公式12

$$\mathbf{N}=\left(\begin{array}{c}\begin{array}{rl}& -\sum ({\mathbf{D}}_i.x\times WA\times C()),\\ & -\sum ({\mathbf{D}}_i.y\times WA\times C()),\\ 1& -\sum (Q_i\times WA\times S())\end{array}\end{array}\right)$$

这是求法线的公式，但我们做多个波时不应直接使用它。当试图直接混合多个法线时，无论你选择何种法线混合算法来近似,其结果总是错误的。
我们应该使用混合多个波后，使用累加的Binormal，Tangent计算数学意义上正确的法线
这里我们先做两个波，以演示如何混合
我们新建一个∑ 函数，来处理波的混合

这里先把所有波的P,B,T相加,（也就是公式中 $\sum$ 的意义）
然后在函数中计算 $\sum$ 之外的，我们刚才没有处理的内容，
我们回过头再看P、B、T三个公式是如何处理$\sum$ 的

Position
原公式中，P需要为∑P的xy通道分别增加x，y：

但由于虚幻中顶点空间为0，所以并不需要做这个偏移，直接使用相加后的Position

Binormal

newB.x = 1-∑B.x;
newB.y =   -∑B.y;
newB.z =     ∑Bz;

Tangent

newT.x =   -∑T.x;
newT.y = 1-∑T.y;
newT.z =     ∑T.z;

法线$N=Cross(B,T)$,结合上面的计算，函数内容为：

我们先做一个使用DDXDDY实现的顶点法线：
这种做法可以计算出顶点准确的法线，但它受限于模型顶点精度，不能直接用作法线。但是它一个很好的参考，让我们知道正确的法线看上去应该是什么样子

左侧为ddxddy计算的法线，右侧为我们的计算结果，可以看到其结果相一致，我们得到了顶点准确的法线

当我们完成了这些，能注意到，Gerstner Waves函数有一些参数是可以通过物理计算求出.
接下来演示使用真实的物理常数，对算法做一些优化，去掉一些手工的变量，
待续

这些公式不像4b、5b和6b方程那样简洁明了，但它们计算起来非常高效。

在形成波峰环的背景下，仔细观察法线的 $z$ 分量证明了这一点非常有趣。虽然 Tessendorf (2001) 从流体动力学的纳维-斯托克斯¹描述和“李变换技术²”中推导出他的“切碎效应³”，最终结果是在频率域中表达的格斯特纳波的一个变体。
在频率域中，可以避免并检测到波顶的环形，但在空间域中，我们可以清楚地看到正在发生的事情。
当 $Q_i\times w_i\times A_i$ 的和大于 $1$ 时，我们法线的 $z$ 分量在峰值处可能变为负值，因为我们的波浪会在自身上方形成环。
只要我们选择的 $Q_i$ 使得这个和始终小于或等于 $1$，我们将形成尖锐的峰值但永远不会形成环。

1.2.4参数解释

波长和速度 ^{WavelengthSpeed}

我们首先选择合适的波长。与其追求现实世界的分布，不如最大化我们能承担的少数波浪的效果。
相似长度的波浪的叠加突出了水面的动态性。
因此，我们选择一个中值波长，并在该长度的一半到两倍之间生成随机波长。中值波长在创作过程中被编写，它可以随时间变化。
例如，在暴风雨期间，波浪可能比晴朗平静时显著更大。
注意，我们不能改变正在活动的波的波长。即便是逐渐改变，波浪的波峰也会向原点扩展或收缩，这看起来非常不自然。
因此，我们改变当前的平均波长，随着时间的推移，当波浪淡出以后，它们将基于新的长度重新生成。方向也是如此。

根据波长，我们可以轻松计算它在表面上的传播速度。Tessendorf 2001中给出忽略高次项的传播关系：

Equation 13 _公式13

$$w=\sqrt{g\times \frac{2\pi }{L}}$$
其中 $w$ 是角频率，$g$ 是标准下的（例如980cm/s）重力常数，$L$ 是波峰到波峰的长度。

振幅 ^Amplitude

如何处理振幅是一个见仁见智的问题。
虽然振幅可以看做波长和当前天气条件的函数来求得振幅导数，但我们还是要使用在编写时指定的常数（或脚本化的）比率。
更准确地说，与中值波长一起，艺术家指定了一个中值振幅。
对于任何大小的波浪，其振幅与波长的比率将匹配中值振幅与中值波长的比率。

方向 ^Direction

波浪行进的方向与其他参数完全独立，因此我们可以根据自己选择的任何标准为每个波浪选择一个方向。
如前所述，我们从大致是风向的恒定向量开始。然后我们从风向的恒定角度内随机选择方向。
这个恒定角度在内容创建时被指定，或者可能被脚本化。

施工中待续

• 参考来自Gerstner Waves圣经—— Mark Finch - Cyan Worlds GPU Gems