作者：i_dovelemon
日期：2020-04-25
主题：Perlin Noise, Curl Noise, Finite Difference Method
引言
最近在研究流体效果相关的模拟。经过一番调查，发现很多的算法都基于一定的物理原理进行模拟，计算量相对来说都比较高昂。最终寻找到一个基于噪音实现的，可在视觉上模拟流体效果的方法：Curl Noise。题图就是通过 Curl Noise 模拟的流体向量场控制的百万粒子的效果。
背景知识
在讲解什么是 Curl Noise 之前，我们需要了解一些相关背景知识。
向量场（Vector Field）
一个2D 或者 3D 的向量场，表示的是赋予空间中任意点一个 2D 或者 3D 向量的函数。公式表示如下所示：
$vec{F}left(x,y ight)=Pleft(x,y ight)vec{i}+Qleft(x,y ight)vec{j}$
$vec{F}left(x,y,z ight)=Pleft(x,y,z ight)vec{i}+Qleft(x,y,z ight)vec{j}+Rleft(x,y,z ight)vec{k}$
其中，$P$，$Q$，$R$ 各表示一个标量函数，即它们的返回值是一个标量；$vec{i}$，$vec{j}$，$vec{k}$ 各表示一个基向量。（参考文献[1]）
上面数学的解释大家可能不熟悉，但是很多人或多或少的都看过向量场的图片形式，如下所示：
散度和旋度（Curl and Divergence）
首先，我们来定义一个 $ abla$ 操作，如下所示：
$ abla=frac{partial }{partial x}vec{i}+frac{partial }{partial y}vec{j}+frac{partial }{partial z}vec{k}$
其中$partial$表示的是偏导数符号，不熟悉的读者可以去复习下微积分或者参考文献[2]。有了这个操作符之后，我们定义旋度为：
$curlvec{F}= abla imesvec{F}=(frac{partial R}{partial y}-frac{partial Q}{partial z},frac{partial P}{partial z}-frac{partial R}{partial x},frac{partial Q}{partial x}-frac{partial P}{partial y})$
其中$ imes$为叉积操作符（参考文献[3]）。
有了旋度之后，我们再来定义散度，同样的，公式如下所示：
$divvec{F}= ablacdot vec{F}=frac{partial P}{partial x}+frac{partial Q}{partial y} + frac{partial R}{partial z}$
特别的，散度和旋度之间有如下的一个关系：
$div(curlvec{F})=0$
以上内容，参考文献[4]。
根据上面的公式，我们可以知道，对于一个向量场的旋度场，它的散度为 0，即它是一个无源场（Divergence-Free）。而一个散度为 0 的向量场，表示这个场是不可压缩的流体，这对日常所见的流体来说是一个很重要的视觉性质，所以据此我们可以使用一个场的旋度场来模拟流体效果。
Curl Noise
所谓 Curl Noise，即是对一个随机向量场，进行 Curl 操作之后得到的新场。因为满足散度为 0 的特性，所以这个场看上去就具有流体的视觉特性。如果用这个场作为速度去控制粒子，即可得到开头视频中流动的效果。
2D Curl Noise
前面我们说过，需要一个随机的向量场。这里我们使用 Perlin Noise 来进行模拟，关于 Perlin Noise 网上一堆资料，这里就不再赘述。
我们假设 Perlin Noise 的函数为：
$N(x,y)$
它的返回值是一个标量值。然后据此建立一个新的向量场：
$vec{F}(x,y) = (N(x,y), N(x,y))$
然后对这个新的向量场进行 Curl 操作，即可得到旋度场。
前面只说过 3D 情况下的 Curl 操作是怎么样的，这里给出 2D 版本的 Curl 操作：
$curlvec{F}(x, y) = (frac{partial N(x,y)}{partial y}, -frac{partial N(x,y)}{partial x})$
这里就只剩下了最后一个问题，那就是形如 $frac{partial N(x,y)}{partial x}$ 这样的偏导数，该怎么计算。我们这里使用一个名为有限差分的方法(Finite Difference Method)来近似求解。
Finite Difference Method
根据文献[2]中对于偏导数的描述，我们知道 $frac{partial N(x,y)}{partial x}$ 只是一种表达方式，它的精确表示方法为：
$frac{partial N(x,y)}{partial x}= N_x(x,y) = lim_{h o0}{frac{N(x + h,y)-N(x,y)}{h}}$
而后面极限的表达方式则给了我们近似计算这个偏导数的方法，只要给定一个较小的 $h$ 值，就能够近似的得到偏导数的结果。而这种计算方法即为：有限差分方法(Finite Difference Method)。
除了上面的极限表示方法之外，还有另外一种极限表示方法，如下所示：
$frac{partial N(x,y)}{partial x}= N_x(x,y) = lim_{h o0}{frac{N(x,y)-N(x-h,y)}{h}}$
这两种差分方法分别称之为前向差分(Forward Difference)和逆向差分(Backward Difference)方法。我这里主要使用逆向差分方法。
有了计算偏导数的方法之后，我们就可以实际带到 2D Curl 操作的公式进行计算，如下是计算 2D Curl Noise 的伪代码：
vec2 computeCurl(float x, float y) { float h = 0.0001f; float n, n1, n2, a, b; n = N(x, y); n1 = N(x, y - h); n2 = N(x - h, y); a = (n - n1) / h; b = (n - n2) / h; return vec2(a, -b); }
知道怎么计算 2D Curl Noise 之后，我们用计算出来的 Curl Noise 作为速度场去控制粒子进行运动，如下是 2D Curl Noise 控制粒子运动的效果：
3D Curl Noise
有了前面 2D Curl Noise 的实现，如法炮制的实现 3D Curl Noise 的推导。
3D Perlin Noise 函数定义为：
$N(x,y,z)$
以此构造出来的 3D 向量场为：
$vec{F}(x,y,z)=(N(x,y,z),N(x,y,z)N(x,y,z))$
对这个场进行 Curl 操作，得到：
$curlvec{F}=(frac{partial N(x,y,z)}{partial y}-frac{partial N(x,y,z)}{partial z},frac{partial N(x,y,z)}{partial z}-frac{partial N(x,y,z)}{partial x},frac{partial N(x,y,z)}{partial x}-frac{partial N(x,y,z)}{partial y})$
据此，给出计算 3D Curl Noise 的伪代码：
vec3 computeCurl(float x, float y) { vec3 curl; float h = 0.0001f; float n, n1, a, b; n = N(x, y, z); n1 = N(x, y - h, z); a = (n - n1) / h; n1 = N(x, y, z - h); b = (n - n1) / h; curl.x = a - b; n1 = N(x, y, z - h); a = (n - n1) / h; n1 = N(x - h, y, z); b = (n - n1) / h; curl.y = a - b; n1 = N(x - h, y, z); a = (n - n1) / h; n1 = N(x, y - h, z); b = (n - n1) / h; curl.z = a - b; return curl; }
以下是根据得到的 3D Curl Noise，并一次控制粒子进行运动的效果：
结论
Curl Noise 在游戏中有大量的运用，Unity 的粒子系统的 Noise Module 就内置了 Curl Noise 的实现。作为游戏开发的人员，很有必要了解下这个技术的原理，便于在实际开发中灵活运用。本文的主要原理来自于参考文献[5]，感兴趣的可以深入去了解。
源代码已上传 Github：https://github.com/idovelemon/UnityProj/tree/master/CurlNoise 。
参考文献
[1] Section 5-1 : Vector Field
[2] Section 2-2：Partial Derivatives