数字孪生混乱中生成世界5 柏林噪声优化 SimplexNoise

欧阳大盆栽
437 阅读6分钟

欢迎继续一起踏上 The Journey of Chaos, 关于 Shader 生成技术的一些基础性学习。噪声会分为以下几篇内容学习

  1. 随机函数与白噪声
  2. 值噪声 ValueNoise
  3. 柏林噪声 Gradient Noise
  4. 多维柏林噪声
  5. 柏林噪声优化 Simplex Noise(本篇)
  6. Worley Noise
  7. Voronoi Noise
  8. FBM深入 我在理解过程中发现并没有人给出倾斜公式的证明，可能是太简单了？ 本文会给出二维倾斜公式的证明，这应该是最具原创性的地方？

Simplex Noise好处与背景

在一篇多维柏林噪声中，发现随着噪声维度的上升，计算复杂度也来越高。 随后再 2002 年柏林有提出了一个方法 Simplex Noise。 相较于 Perlin Noise他有以下优点

  1. 计算复杂度更低 Perlin噪声：在n维空间中，计算复杂度为O(n * 2^n)，随着维度的增加，计算量呈指数级增长。 Simplex噪声：在n维空间中，计算复杂度为O(n^2)，随着维度的增加，计算量增长更为缓慢。这意味着在高维度（如4D、5D）中，Simplex噪声的计算效率更高。
  2. 更少的方向性伪影 Perlin噪声：在某些情况下，Perlin噪声会产生明显的方向性伪影，尤其是在高维度中。这些伪影会影响生成的纹理或地形的自然感。 Simplex噪声：通过使用更复杂的网格结构（如将空间划分为简单形，即n维三角形），Simplex噪声几乎消除了这些方向性伪影。这使得生成的噪声看起来更加自然和均匀。
  3. 更平滑的梯度 Perlin噪声：虽然Perlin噪声已经能够生成平滑的噪声，但在某些情况下，尤其是在高维度中，其梯度可能会出现不连续的情况。 Simplex噪声：Simplex噪声具有更平滑的梯度，并且几乎在任何地方都能计算出连续的梯度。这使得生成的噪声在视觉上更加平滑，没有明显的断点或跳跃。
  4. 硬件实现更简单 Perlin噪声：由于其计算复杂度较高，Perlin噪声在硬件实现上较为困难。 Simplex噪声：Simplex噪声的计算过程更简单，对硬件的要求更低，因此更容易在硬件上实现。
  5. 更高的频率和更复杂的纹理 Perlin噪声：在某些情况下，Perlin噪声可能会产生较平坦的区域。 Simplex噪声：Simplex噪声通常具有更高的频率和更复杂的纹理。这意味着在相同的输入条件下，Simplex噪声能够生成更丰富的细节。
  6. 更好的可扩展性 Perlin噪声：在高维度中，Perlin噪声的性能下降明显，生成速度变慢。 Simplex噪声：由于其计算复杂度较低，Simplex噪声在高维度中表现更好，能够更快地生成噪声。 总结 虽然Perlin噪声在某些应用中仍然表现良好，但Simplex噪声在高维度、需要更平滑梯度和更少方向性伪影的场景中具有明显优势。如果你的应用涉及到高维度噪声生成，或者需要更自然的纹理和地形，Simplex噪声可能是更好的选择。

但是由于 Simplex Noise 与Perlin的论文给出代码不好理解， 所以在之后的几年里，并没有大批量使用， 于是2005 Stefan Gustavson 写了 资料 2 中的论文《simplex noise demystified》，。 本篇主要内容也是对该论文进行阅读理解与翻译

简单型填充

如何在一个 N维空间里面进行用同样的一个 Tile进行填充。 例如 我们在一维空间中可以使用 长度为 1 的线段进行 1 为空间拆分， 在二维空间中可以使用 正方形进行空间拆分，在三维空间中可以使用正方体进行空间拆分。被拆分之后坐标系，我们可以叫做单型网格坐标系，单型网格的特点是所有网格单元的形状和大小相同，便于进行数值计算和优化算法。

但是除了使用正方形做网格 能不能更简单呢？ 考虑在二维空间中，我们先用三角形组成平面，然后两个三角形组成一个菱形。 菱形经过拉伸之后就变成了正方形，也就是能够组成一个平面 2501292916

在三维空间中， 最简单的构成空间的基本元素是 4 面体，如下图所示 2501293003

4 位空间的几何形太难想象了。如果读者中有比较好的空间想象力可以试试。 我大概能想到的一种方法就是投影。 例如我们可以通过运动的点表示线， 通过运动的线表示面，也可以通过运动的几何体表示 4 维。

在数学中，N维超立方体是一个具有N个维度的空间形状。它由N个边长相等的小正方体组成，并且每个小正方体都有一个顶点。这些顶点可以被组织成一个简单形（也称为N维金字塔），这个简单形有N+1个面。使用简单网格的优点是可以减少计算复杂度，并且更容易进行插值计算。

  1. 减少计算复杂度好理解， 因为点少了
  2. 更加容易进行插值运算怎么理解？

最后一个结论, For any given NN, N!N! simplexes can fill an N-dimensional hypercube that is flattened along the main diagonal.

问题的核心

在经典柏林噪声中，在二维平面使用 4 个点 来表示。 每个网格的结果是由 4 个点来控制。那使用简单型填充之后，两个问题需要解决

  1. 3 个点的坐标如何定义，在 正方形 tile中，其坐标就是点的坐标。
  2. 3 个点如何让噪声梯度连续呢， 我们没办法在通过双线性插值来进行每个 TIle内 P点的值了。 这两个问题的解决思路就是 Simplex Noise算法核心。本文将从二维的Simplex Noise算法分析，不会进行高维扩展，关于推导的所有基础知识点都会讲明白，因为我在查阅资料的过程中发现很多解析都默认读者有一定的基础，假定没有基础就无法理解。

顶点坐标确定

2501302938 上面这张图表示了如何进行坐标系转化，并确定顶点。 我看懂之后，觉得基本上就是这三张图。 但是他的代码是这个样子的

 double n0, n1, n2; // Noise contributions from the three corners
 // Skew the input space to determine which simplex cell we're in
 final double F2 = 0.5*(Math.sqrt(3.0)-1.0);
 double s = (xin+yin)*F2; // Hairy factor for 2D
 int i = fastfloor(xin+s);
 int j = fastfloor(yin+s);
 final double G2 = (3.0-Math.sqrt(3.0))/6.0;
 double t = (i+j)*G2;
 double X0 = i-t; // Unskew the cell origin back to (x,y) space
 double Y0 = j-t;
 double x0 = xin-X0; //The x,y distances from the cell origin
 double y0 = yin-Y0;

上来这个 F2 G2就把我干懵了。。。 这是个啥。 跟上面的图有什么关系，为什么是F2的值是0.5*(Math.sqrt(3.0)-1.0) . 这是我最大的困惑

坐标变换与逆变换

F是偏斜公式， 表示如何将上图中正方形网格倾斜成三角网格。 让我们看看上面三角网络 2501302138

其转换过程是做一个Shear让 x，y影响 15°， 然后再做一个缩放让对角线能够从 2\sqrt{2} 变成 1. 那我们知道， 对于一个给定的x轴 shear角度 θ\theta，shear变换的矩阵表达式可以写为：

(1tan(θ)01)\begin{pmatrix} 1 & tan(\theta) \\ 0 & 1 \end{pmatrix}

给定一个x轴的缩放系数 k， scale表达式写为

(k001)\begin{pmatrix} k & 0 \\ 0 & 1 \end{pmatrix}

假设变换前坐标为 (x,y)(x, y) 变换后的坐标为 (x′,y′)(x\prime, y\prime) 。那么有

(x′,y′)=(1tan(θ)tan(θ)1)⋅(k00k)⋅(x,y) (x\prime, y\prime) = \begin{pmatrix} 1 & tan(\theta) \\ tan(\theta) & 1 \end{pmatrix} \cdot \begin{pmatrix} k & 0 \\ 0 & k \end{pmatrix} \cdot (x, y)

其中 θ\theta 为 15°， 缩放系数为k。 我们单拿 x来看

x′=k(x+ytan(15∘))∵tan(15∘)=2−3∴x′=k⋅(x+y(2−3))=k⋅((3+1)⋅(3−1)2x+(3−1)22y)=k⋅(3−1)⋅(3+12x+3−12y)=k⋅(3−1)⋅(x+3−12⋅(x+y))\begin{aligned} x\prime &= k(x+ytan(15^\circ)) \\ \because \quad tan(15^\circ) &= 2 - \sqrt{3} \\ \therefore x\prime &= k \cdot (x+y( 2 - \sqrt{3})) \\ &= k \cdot ( \frac{(\sqrt3 + 1) \cdot (\sqrt3 - 1)}{2} x + \frac{(\sqrt3-1)^2} {2}y ) \\ &= k \cdot (\sqrt3 - 1) \cdot (\frac{\sqrt3 + 1}{2} x + \frac{\sqrt3 - 1}{2}y) \\ & = k \cdot (\sqrt3 - 1) \cdot ( x + \frac{\sqrt3 - 1}{2} \cdot (x+y)) \end{aligned}

看到没有。。 *为什么是F2的值是0.5(Math.sqrt(3.0)-1.0)**，这个问题的答案就在这里...

x′=x+(x+y+...)∗Fy′=y+(x+y+...)∗FF=(n+1)−1n\begin{aligned} x\prime &= x + (x + y + ...) * F \\ y\prime &= y + (x + y + ...) * F \\ F &= \frac{\sqrt(n+1) - 1}{n} \end{aligned}

那偏斜方程回到正常坐标系其公式，可以通过解公式得

x=x′−(x′+y′+...)∗Gy=y′−(x′+y′+...)∗GG=1−1/(n+1)n\begin{aligned} x &= x\prime - (x\prime + y\prime + ...) * G \\ y &= y\prime - (x\prime + y\prime + ...) * G \\ G &= \frac{1 - 1 / \sqrt(n+1) }{n} \end{aligned}

至于为什么他需要写成这种形式， 我个人感觉是因为这种计算形式的 GPU非常方便... . 只有加法和乘法，没有三角函数一类复杂的运算。 纯属个人瞎感觉，不负责任。。。。反正我最后的实现一定不用这种方法，，我要用最好理解的方式写代码。。。

三个顶点坐标计算

有了变换与逆变换， 接下来就是计算三个顶点。 先设定变换前坐标系为coordcoord,变换后坐标系为coord′coord^\prime. 在经典柏林噪声 4 个顶点的计算方式为

p0=floor(uv)+vec2(0,0)p1=floor(uv)+vec2(0,1)p2=floor(uv)+vec2(1,0)p3=floor(uv)+vec2(1,1)\begin{aligned} p_0 &= floor(uv) + vec2(0, 0) \\ p_1 &= floor(uv) + vec2(0, 1) \\ p_2 &= floor(uv) + vec2(1, 0) \\ p_3 &= floor(uv) + vec2(1, 1) \end{aligned}

一样我们首先找出 vec2(0, 0) 的这个 P0点。 先求出 P点在变换之后的点P′P^\prime，通过 frac获取网格内坐标(i′,j′)(i^\prime, j^\prime)， 这个点其实也是一个相对于 P0点的移动向量i′j′⃗\vec{i^\prime j^\prime}。再将这个 Vector进行坐标逆变换获得新的ij⃗\vec{i j}。通过将原有的 P点与该向量进行做差得到在正常坐标下的 P0点坐标。

F2=3−12G2=3−36P′=P+(P.x+P.y)∗F2i′j′⃗=floor(P′)ij⃗=i′j′⃗−(i′+j′)∗G2P0=P−ij⃗\begin{aligned} F2 &= \frac{\sqrt3 - 1}{2} \\ G2 &= \frac{3 - \sqrt3}{6} \\ P^\prime &= P + (P.x + P.y) * F2 \\ \vec{i^\prime j^\prime} &= floor(P^\prime) \\ \vec{i j} &= \vec{i^\prime j^\prime} - (i^\prime + j^\prime) * G2 \\ P_0 &= P - \vec{i j} \end{aligned}

找到第一个点之后，我们继续找与 P0相连的点P1. 由于我们讲一个正方形分为两个三角形，所以三角形的与 P0相连的点，需要通过与（1,1）向量做叉积判断左右， 如下图所示 2501305641

由于仿射变换保留了直线平行线的线性变换。所以我们可以直接通过坐标空间 Coord的 ij来判定

[ij]×[11]>0→i−j>0−>i>j\begin{bmatrix} i \\ j \\ \end{bmatrix} \times \begin{bmatrix} 1 \\ 1 \\ \end{bmatrix} > 0 \to i - j > 0 -> i > j

说明 P在右边，那么 P1点就是 增加（1，0）。同时这里的(1, 0)是在 coord′coord^\prime的 (1, 0)。同样把它当做一个平移向量，可以代入逆变换公式得到

Vector(1,0)=Vector(1−c,−c)Vector(0,1)=Vector(−c,1−c)P1=P0+Vector(1,0)P2=P1+Vector(0,1)\begin{aligned} Vector(1, 0) &= Vector(1 - c, -c) \\ Vector(0, 1) &= Vector(-c ,1- c) \\ P_1 &= P0 + Vector(1, 0) \\ P_2 &= P1 + Vector(0, 1) \end{aligned}

终于我们求出了三个点的坐标。。。。到这里看了一天了，， 好笨好笨啊。。

网格内噪声插值

还记得我们在经典柏林噪声中需要求P点到端点距离向量，并通过这个向量与随机梯度做 Dot。随后进行双线性插值。 就像贝塞尔曲线有几个控制点，经典柏林噪声就是有角上的 4 个控制点。 2501275409 同时在经典柏林噪声，通过 ease函数控制在边缘处导数为 0，保证在边界处值导数相等不发生突变。 在 Simplex也要解决这两个问题

  1. 如何通过角上的三个点，控制 P点的值
  2. 网格边界左右两边导数值相等

而回到SImpleNoise中. 如果在边界上的点只受到边界两个端点的影响是不是就能保证边界处导数相等啦 2501313855 考虑上D,B两点圆的范围内，会收到D,B两点影响，而 E与 E'两个点只受到 AC两个点影响。 2501314557 当上图 EG距离趋近于 0 的时候 u,v两个向量相等， 两边导数值也相对相等。 这解决了边界问题。

所以我们可以通过P点与三个点的距离作为控制变量， 距离越远受到的影响越小。 P点到角的线，叫做径向。 希望衰减的程度更加剧烈，于是做了一个平方。那么三角形的高是多少呢？这有助于我们确定上图圆的大小。 回到第一节的坐标变换与逆变换，我们的 K值与 Shear之后的边长，可以得到。三角形的高为0.52

shearedEdgeLength=63k=3−1heightOfTriangle=63⋅(3−1)⋅cos(30∘)=0.52\begin{aligned} shearedEdgeLength &= \frac{\sqrt6}{3} \\ k &= \sqrt3 - 1 \\ heightOfTriangle &= \frac{\sqrt6}{3} \cdot (\sqrt3 - 1) \cdot cos(30^\circ) = 0.52 \end{aligned}

当然也可以去比这个 height稍微小一点的值对于边为 1 的等边三角形，其圆最小半径为 33\frac{\sqrt3}{3} 可以满足三角形内任何一个点都在其中一个圆中。也就是 r的取值范围 [0.3, 0.52] 于是有公式推理出最终点 P的值

dn=P−Pngradientn=hash(id+offsetn)r=0.52value=∑n=13max(0,r2−dn2)⋅dot(gradientn,dn)normalized=70.∗value\begin{aligned} d_n &= P - P_n \\ gradient_n &= hash(id + offset_n) \\ r &= 0.52 \\ value &= \sum_{n=1}^{3} max(0, r^2 - d_n^2) \cdot dot(gradient_n, d_n) \\ normalized &= 70. * value \end{aligned}

最后那个 70 保障随机值的最大值能够到 1，这并非一个准确值。 最大值多了一些少了一些其实没关系，我看到资料中有知乎文章给出了 70 的结果是

134⋅16⋅2⋅2≈170\frac{1}{3}^4 \cdot \frac{1}{\sqrt{6}} \cdot \sqrt2 \cdot 2 \approx \frac{1}{70}

没太能 get到

资料

  1. the book of shaders noise thebookofshaders.com/11/
  2. simplex noise解惑： www.researchgate.net/publication…
  3. zhuanlan.zhihu.com/p/240763739
  4. lygia snoise源码： lygia.xyz/generative/…
avatar
文章
阅读
粉丝