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

在这里插入图片描述

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

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

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

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


1.2.1选择波形

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

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

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

  • Amplitude 振幅 A A A

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

    • 波峰每秒前进的距离。 将速度表示为相位常数 φ \varphi φ 更为方便。
    • φ = S × 2 π L \varphi = S \times\frac{2\pi}{L} φ=S×L2π
  • Direction 方向 D D D

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

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

MF_GerstnerWaves

新建一个 MF_GerstnerWaves函数,我们先把上方出现的变量在input预留:
在这里插入图片描述
以防有疑问,最后一个Node是“重路由” ,它是用来避免蓝图变意大利面的东西,无功能上的意义
在这里插入图片描述

1.2.2 法线和切线

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

Equation 6a 公式6a

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

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

格斯特纳波GerstnerWaves

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

GerstnerWaves

图1-5 GerstnerWaves

这是GerstnerWaves函数:

Equation 9 公式9

P ( x , y , t ) = ( x + ∑ ( Q i A i × D i . x × cos ⁡ ( w i D i ⋅ ( x , y ) + φ i t ) ) , y + ∑ ( Q i A i × D i . y × cos ⁡ ( w i D i ⋅ ( x , y ) + φ i t ) ) , ∑ ( A i sin ⁡ ( w i D i ⋅ ( x , y ) + φ i t ) ) ) \mathbf{P}\left(x, y, t\right) = \left(\begin{array}{l} \begin{alignedat}{3} &x&+&\sum \left(Q_{i}A_{i} \times \mathbf{D}_{i}.x \times \cos\left(w_{i}\mathbf{D}_{i} \cdot (x, y) + \varphi_{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) + \varphi_{i}t\right)\right),\\ & & &\sum \left(A_{i} \sin\left(w_{i}\mathbf{D}_{i} \cdot (x, y) + \varphi_{i}t\right)\right) \end{alignedat} \end{array}\right) P(x,y,t)= xy++(QiAi×Di.x×cos(wiDi(x,y)+φit)),(QiAi×Di.y×cos(wiDi(x,y)+φit)),(Aisin(wiDi(x,y)+φit))

∑为求和,∑( )也就是所有波加在一起,然后再在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);

能看到我们想实现这个公式还需要两个变量, x y xy xy t t t,他们分别是:
在这里插入图片描述

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

Q的计算为:
在这里插入图片描述
其中 n u m W a v e s numWaves numWaves 为波的总数,例如水面由3个GerstnerWaves组合,这里就是3
Steepness变量就是文中说的"陡峭程度"参数

可以看到公式里的( w × D ⋅ ( x , y ) + φ × t w×D·(x,y) + φ×t w×D(x,y)+φ×t )是重复的,我们先把他的SinCos算出来
(这不会影响着色器复杂度,只会让我的截图好看,避免意大利面)
注意这里的 s i n sin sin c o s cos cos 的周期为 2 p i 2pi 2pi
在这里插入图片描述
在这里插入图片描述
现在我们可以计算公式了
在这里插入图片描述
这就是顶点位移的结果,接下来我们计算法线

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

其中:
W A = w i × A i , WA=w_{i}\times A_{i}, WA=wi×Ai,
S ( ) = sin ⁡ ( w i × D i ⋅ P + φ i t ) S()=\sin\left(w_{i}\times\mathbf{D}_{i}\cdot\mathbf{P}+\varphi_{i}t\right) S()=sin(wi×DiP+φit)
C ( ) = cos ⁡ ( w i × D i ⋅ P + φ i t ) C()=\cos\bigl(w_{i}\times\mathbf{D}_{i}\cdot\mathbf{P}+\varphi_{i}t\bigr) C()=cos(wi×DiP+φit)

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

S ( ) S() S() C ( ) C() C() 我在蓝图里写成了 S i n ( ) Sin() Sin() C o s ( ) Cos() Cos()
在这里插入图片描述

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

Equation 10 公式10

B = ( 1 − ∑ ( Q i × D i . x 2 × W A × S ( ) ) , − ∑ ( Q i × D i . x × D i . y × W A × S ( ) ) , ∑ ( D i . x × W A × C ( ) ) ) \mathbf{B} = \begin{pmatrix} \begin{alignedat}{3} &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{alignedat} \end{pmatrix} B= 1(Qi×Di.x2×WA×S()),(Qi×Di.x×Di.y×WA×S()),(Di.x×WA×C())

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

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

在这里插入图片描述

Equation 11 公式11

T = ( − ∑ ( Q i × D i . x × D i . y × W A × S ( ) ) , 1 − ∑ ( Q i × D i . y 2 × W A × S ( ) ) , ∑ ( D i . y × W A × C ( ) ) ) \mathbf{T}= \begin{pmatrix} \begin{alignedat}{3} & &-&\sum\bigl(Q_{i}\times\mathbf{D}_{i}.x\times\mathbf{D}_{i}.y\times WA\times S()\bigr),\\ &1&-&\sum\bigl(Q_{i}\times\mathbf{D}_{i}.y^{2}\times WA\times S()\bigr),\\ & & &\sum\bigl(\mathbf{D}_{i}.y\times WA\times C()\bigr) \end{alignedat} \end{pmatrix} T= 1(Qi×Di.x×Di.y×WA×S()),(Qi×Di.y2×WA×S()),(Di.y×WA×C())

Tangent:

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

在这里插入图片描述

Equation 12 公式12

N = ( − ∑ ( D i . x × W A × C ( ) ) , − ∑ ( D i . y × W A × C ( ) ) , 1 − ∑ ( Q i × W A × S ( ) ) ) \mathbf{N}=\begin{pmatrix} \begin{aligned} &{}- \sum\bigl(\mathbf{D}_i.x\times WA\times C()\bigr),\\ &{}- \sum\bigl(\mathbf{D}_i.y\times WA\times C()\bigr),\\ 1 &{}- \sum\bigl(Q_i\times WA\times S()\bigr) \end{aligned} \end{pmatrix} N= 1(Di.x×WA×C()),(Di.y×WA×C()),(Qi×WA×S())

这是求单个波法线的公式,但我们做多个波则不应直接使用它。
因为当你在最后试图混合波的法线时,无论你选择何种法线混合算法,其结果总是错的。
我们应该累加波的 B i n o r m a l Binormal Binormal T a n g e n t Tangent Tangent ,在最后通过叉积求解 N o r m a l Normal Normal (见公式6a) ,这将避免直接对法线进行混合,进而计算数学意义上正确的法线


我们先做两个波,以演示如何混合:

我们新建一个 ∑ ∑ 函数,MF_GerstnerWave_∑
这里先把所有波的 P o s i t i o n Position Position B i n o r m a l Binormal Binormal T a n g e n t Tangent Tangent 求和(对应公式9、10、11中的 ∑ ∑
将结果输入函数,稍后再看MF_GerstnerWave_∑ 内容
在这里插入图片描述

MF_GerstnerWave_∑


Position

公式9中,P需要为∑P的xy通道分别增加x,y:
P ( x , y , t ) = ( x + ∑ ( Q i A i × D i . x × cos ⁡ ( w i D i ⋅ ( x , y ) + φ i t ) ) , y + ∑ ( Q i A i × D i . y × cos ⁡ ( w i D i ⋅ ( x , y ) + φ i t ) ) , ∑ ( A i sin ⁡ ( w i D i ⋅ ( x , y ) + φ i t ) ) ▲ 这里 ) \mathbf{P}\left(x, y, t\right) = \left(\begin{array}{l} \begin{alignedat}{3} &x&+&\sum \left(Q_{i}A_{i} \times \mathbf{D}_{i}.x \times \cos\left(w_{i}\mathbf{D}_{i} \cdot (x, y) + \varphi_{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) + \varphi_{i}t\right)\right),\\ & & &\sum \left(A_{i} \sin\left(w_{i}\mathbf{D}_{i} \cdot (x, y) + \varphi_{i}t\right)\right)\\ &\blacktriangle& &\\ &\mathbf{这里} \end{alignedat} \end{array}\right) P(x,y,t)= xy这里++(QiAi×Di.x×cos(wiDi(x,y)+φit)),(QiAi×Di.y×cos(wiDi(x,y)+φit)),(Aisin(wiDi(x,y)+φit))

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

// Position.xy = ∑P.xy+(x,y);
Position.xyz = ∑P.xyz;

Binormal

根据公式10:

B = ( 1 − ∑ ( Q i × D i . x 2 × W A × S ( ) ) , 0 − ∑ ( Q i × D i . x × D i . y × W A × S ( ) ) , ∑ ( D i . x × W A × C ( ) ) ▲ 这里 ) \mathbf{B} = \begin{pmatrix} \begin{alignedat}{3} &1&-&\sum\left(Q_i \times \mathbf{D}_i.x^2 \times WA \times S()\right),\\ &\phantom0&-&\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)\\ &\blacktriangle& &\\ &\mathbf{这里} \end{alignedat} \end{pmatrix} B= 10这里(Qi×Di.x2×WA×S()),(Qi×Di.x×Di.y×WA×S()),(Di.x×WA×C())

Binormal.x = 1-∑B.x;
Binormal.y = 0-∑B.y;
Binormal.z =     ∑B.z;

Tangent

根据公式11:

T = ( 0 − ∑ ( Q i × D i . x × D i . y × W A × S ( ) ) , 1 − ∑ ( Q i × D i . y 2 × W A × S ( ) ) , ∑ ( D i . y × W A × C ( ) ) ▲ 这里 ) \mathbf{T}= \begin{pmatrix} \begin{alignedat}{3} &\phantom0&-&\sum\bigl(Q_{i}\times\mathbf{D}_{i}.x\times\mathbf{D}_{i}.y\times WA\times S()\bigr),\\ &1&-&\sum\bigl(Q_{i}\times\mathbf{D}_{i}.y^{2}\times WA\times S()\bigr),\\ & & &\sum\bigl(\mathbf{D}_{i}.y\times WA\times C()\bigr)\\ &\blacktriangle& &\\ &\mathbf{这里} \end{alignedat} \end{pmatrix} T= 01这里(Qi×Di.x×Di.y×WA×S()),(Qi×Di.y2×WA×S()),(Di.y×WA×C())

Tangent.x = 0-∑T.x;
Tangent.y = 1-∑T.y;
Tangent.z =     ∑T.z;

Normal

根据公式6a:
N ( x , y ) = B ( x , y ) × T ( x , y ) \mathbf{N}(x,y)=\mathbf{B}(x,y)\times\mathbf{T}(x,y) N(x,y)=B(x,y)×T(x,y)
N = B × T \mathbf{N}=\mathbf{B}\times\mathbf{T} N=B×T

Normal = Cross(Binormal,Tangent);

结合上面的计算,∑函数MF_GerstnerWave_∑ 内容为:在这里插入图片描述

对结果验证

我们先做一个使用DDXDDY实现的差分几何法线
这种方法利用了屏幕空间相邻像素的差分来重建表面法线,它受限于模型顶点精度,看上去像LowPoly,不能直接用作法线。
但是它一个很好的参考,让我们知道正确的法线看上去应该是什么样子

在这里插入图片描述
下图中,左侧为差分几何法线,右侧为我们的计算结果,可以看到其结果相一致,我们得到了准确的法线
在这里插入图片描述

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

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

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

1.2.4参数解释

波长和速度 Wavelength Speed

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

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

Equation 13 公式13

w = g × 2 π L w=\sqrt{g\times{\frac{2\pi}{L}}} w=g×L2π
其中 w w w 是角频率, g g g 是标准下的(例如980cm/s)重力常数, L L L 是波峰到波峰的长度。

海面由多个波组成,每个波的参数越少,越易于我们控制
目前为止,每个波的 S p e e d i Speed_i Speedi都是独立的参数,10个波就需要手工指定10个 S p e e d i Speed_i Speedi
但根据公式13,我们已经有了求解 S p e e d i Speed_i Speedi的条件
如果我们引入统一的全局变量 W i n d s p e e d Windspeed Windspeed 风速 G r a v i t y Gravity Gravity重力 980cm/s
就可以求解 S p e e d i Speed_i Speedi


我们先将 S p e e d Speed Speed 输入重命名为 W i n d S p e e d WindSpeed WindSpeed ,并按照公式13带入 G r a v i t y Gravity Gravity重力
Speed = WindSpeed × Gravity × w i \text{Speed} = \text{WindSpeed} \times \sqrt{\text{Gravity} \times w_i} Speed=WindSpeed×Gravity×wi

Speed = WindSpeed * Sqrt(Gravity*w);

在这里插入图片描述现在,我们的波浪可以由全局统一变量的 风向 Direction 和 风速 WindSpeed ,为每个波自动计算速度Speed

振幅 Amplitude

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

方向 Direction

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

优化数据结构

是时候整理 MF_GerstnerWaves 函数的变量了,我们先把Time也开放出来
在这里插入图片描述

在这里插入图片描述

现在函数本身已经变为4个全局变量红色,4个独立变量黄色
在这里插入图片描述
四个独立变量,也就是说,我们差一点就可以用一个Vector4就可以描述一个波了
但Direction是双通道,这4个变量实际是5数值,
在这里插入图片描述
5个值最少需要两个变量

MF_AngleToVector2

制作Shader,我们常说要用空间换时间,现在我们要反其道而行,增加计算量以减少参数输入

—— 只送大脑

我们从Direction入手,他现在是2通道值Vecter2,现在要把他反向计算为1通道值Scale
二维Direction本质就是2D方向,当我们定义一个方向为"前",任何方向就可以用角度Ange描述
在这里插入图片描述

使用-1到0的值计算分量得到Direction

RotationAngle不是标准的“Angle”,它的输入为-1到1,这样既可以方便填写,也省略了多步计算,乘π后直接求分量
默认输入为0,输出为(1,0)

回到MF_GerstnerWave修改 D i r e c t i o n Direction Direction Vector2 输入为 D i r e c t i o n A n g l e DirectionAngle DirectionAngle Scale
在这里插入图片描述
我们还省下了一个归一化,不亏
在这里插入图片描述
现在一个变量就可以定义一个波,这在数据结构上有重大意义
例如数据在材质参数集定义多个波:
在这里插入图片描述
在这里插入图片描述


施工中待续

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

  1. Navier-Stokes ↩︎

  2. Lie Transform Technique ↩︎

  3. choppiness ↩︎

  • 11
    点赞
  • 19
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值