噪声学习

关于噪声

噪声在图形学中的使用非常普遍，这里的噪声是通过程序生成一系列的随机值，通过特定的噪声算法，在给定一个输入，生成对应输出，输入和输出可以是不同维数的组合。 
噪声在图形学中常用来模拟火焰、地形、云朵等包含随机变化的对象。既然是获得随机值，为何不直接使用类似 random 这样的随机函数呢？原因是随机函数所得到的值在既定范围内分布过于均匀，这样的随机值被称为"白噪声",二维白噪声的纹理为： 
 
实际在自然界中，充满随机变化的对象，如地形，云朵等，分布不会像白噪声这样，对于地形，既有高度变化很大的山脉和峡谷，也有起伏变化很小的山丘，这些不同程度的随机变化叠加在一起形成自然的地形走向。因此才会有特定的噪声函数，来模拟不同的自然界中的噪声（随机变化）。比如，Perlin噪声用于云朵，火焰的模拟，Simplex噪声在此基础上进行优化，提高计算的效率。

噪声分类

由wiki的分类，程序噪声算法分类如下： 

以上为单噪声，与之对应的叫做 分形噪声(也可以叫做FBM噪声)，分形噪声是在将这些单噪声以不同的频率和振幅进行叠加组合，从而产生不同程度的随机细节。这里的频率指的是计算噪声时的采样距离，例如对于基于晶格的噪声们，频率越高，单位面积（特指二维）内的晶格数目越多，看起来噪声纹理“越密集”；而振幅指的就是噪声的值域。下图显示了一些基础噪声和它们fbm后的效果： 

左侧为单噪声，右侧为不同频率叠加后的效果 
【图片来自https://blog.csdn.net/candycat1992/article/details/50346469】

噪声算法的实现

主要包括常用噪声 Perlin噪声、Value噪声、Simplex噪声及简单的FBM实现。

Perlin噪声

Perlin噪声的实现有三个步骤：

1. 定义一个晶格结构，每个晶格的顶点有一个“伪随机”的梯度向量（其实就是个随机向量）。对于二维的Perlin噪声来说，晶格结构就是一个平面网格，三维的就是一个立方体网格。
2. 输入一个点（二维的话就是二维坐标，三维就是三维坐标，n维的就是n个坐标），我们找到和它相邻的那些晶格顶点（二维下有4个，三维下有8个，n维下有2^n个），计算该点到各个晶格顶点的距离向量，再分别与顶点上的梯度向量做点乘，得到2n个点乘结果。
3. 使用缓和曲线来计算它们的权重和。在原始的Perlin噪声实现中，缓和曲线是s(t)=3t^2−2t^3，在2002年，Perlin改进为s(t)=6t^5−15t^4+10t^3。这里简单解释一下，为什么不直接使用s(t)=t，即线性插值。直接使用的线性插值的话，它的一阶导在晶格顶点处（即t = 0或t = 1）不为0，会造成明显的不连续性。s(t)=3t^2−2t^3在一阶导满足连续性，s(t)=6t^5−15t^4+10t^3在二阶导上仍然满足连续性。

对应的图例说明： 

左侧为各个晶格顶点的随机向量，右侧为位于某一晶格内的输入点与晶格顶点的距离向量。

Unity下的实现效果： 

单Perlin噪声（左下角）主要代码：

float PNoise(float2 p){
float2 i=floor(p);
float2 f=frac(p);

//缓和曲线，计算权重
float2 w=f*f*(3.0-2.0*f);

//使用权重值计算插值
return lerp( lerp(dot(Hash22(i+float2(0.0,0.0)),f-float2(0.0,0.0)),
                  dot(Hash22(i+float2(1.0,0.0)),f-float2(1.0,0.0)),w.x),
             lerp(dot(Hash22(i+float2(0.0,1.0)),f-float2(0.0,1.0)),
                  dot(Hash22(i+float2(1.0,1.0)),f-float2(1.0,1.0)),w.x),w.y);
}

其中Hash22是用来生成随机向量的函数。 单Perlin噪声为左下角部分，右下角，左上角和右上角为Perlin单噪声的分形叠加效果。 
右下角部分分形叠加公式：

noise(p)+1/2noise(2p)+1/4noise(4p)+...       

实现代码：

//叠加5层，初始化采样距离设置为7，可以自定义。这种噪声可以用来模拟石头、山脉这类物体。
float Noise_Sum(float2 p){
float f=0.0;
p=p*7.0;
f+=1.0000*Noise(p);p=2.0*p;
f+=0.5000*Noise(p);p=2.0*p;
f+=0.2500*Noise(p);p=2.0*p;
f+=0.1250*Noise(p);p=2.0*p;
f+=0.0625*Noise(p);p=2.0*p;
return f;
}

左上角部分分形叠加公式：

|noise(p)|+1/2|noise(2p)|+1/4|noise(4p)|+...    

实现代码：

//由于进行了绝对值操作，因此会在0值变化处出现不连续性，形成一些尖锐的效果。
//通过合适的颜色叠加，可以用这种噪声来模拟火焰、云朵这些物体。float Noise_Sum_Abs(float2 p){
float f=0.0;
p=p*7.0;
f+=1.0000*abs(Noise(p));p=2.0*p;
f+=0.5000*abs(Noise(p));p=2.0*p;
f+=0.2500*abs(Noise(p));p=2.0*p;
f+=0.1250*abs(Noise(p));p=2.0*p;
f+=0.0625*abs(Noise(p));p=2.0*p;
return f;
}

右上角部分分形叠加公式：

sin(x+|noise(p)|+1/2|noise(2p)|+1/4|noise(4p)|+...) 

实现代码：

//这个公式可以让表面沿着x方向形成一个条纹状的结构。Perlin使用这个公式模拟了一些大理石材质。
float Noise_Sum_Abs_Sin(float2 p){
float f=0.0;
p=p*7.0;
f+=1.0000*abs(Noise(p));p=2.0*p;
f+=0.5000*abs(Noise(p));p=2.0*p;
f+=0.2500*abs(Noise(p));p=2.0*p;
f+=0.1250*abs(Noise(p));p=2.0*p;
f+=0.0625*abs(Noise(p));p=2.0*p;

f=sin(f+p.x/32.0);

return f;   
}  

Value噪声

Value噪声的实现和Perlin噪声唯一的区别在于将原来的梯度向量替换成简单的伪随机值，同时不需要做点乘操作，直接在晶格点处按照缓和曲线的权重值进行叠加。 
Value噪声是一种基于晶格的操作，也需要三个步骤：

1.定义一个晶格结构，每个晶格的顶点有一个“伪随机”的值（Value）。对于二维的Value噪声来说，晶格结构就是一个平面网格，三维的就是一个立方体网格。 
2.输入一个点（二维的话就是二维坐标，三维就是三维坐标，n维的就是n个坐标），我们找到和它相邻的那些晶格顶点（二维下有4个，三维下有8个，n维下有2^n个），得到这些顶点的伪随机值。 
3.使用缓和曲线来计算它们的权重和。同样，缓和曲线可以是s(t)=3t^2−2t^3，也可以是s(t)=6t^5−15t^4+10t^3（如果二阶导不连续对效果影响较大时）。

Value噪声比 Perlin噪声实现过程更加简单，需要进行的计算也比较少。 
Unity下的实现效果： 

单Value噪声（左下角）主要代码：

float VNoise(float2 p){
float2 i=floor(p);
float2 f=frac(p);
//缓和曲线，计算插值点
float2 w=f*f*(3.0-2.0*f);

return lerp(lerp(Hash21(i+float2(0.0,0.0)),Hash21(i+float2(1.0,0.0)),w.x),
            lerp(Hash21(i+float2(0.0,1.0)),Hash21(i+float2(1.0,1.0)),w.x),w.y);
}

单Value噪声像素化明显，不过其对应的分形噪声效果可以接受。

Simplex噪声

Simplex噪声理解上相对不太友好，但其效率上要比Perlin噪声要高。 
Simplex噪声也是一种基于晶格的梯度噪声，它和Perlin噪声在实现上唯一不同的地方在于，它的晶格并不是方形（在2D下是正方形，在3D下是立方体，在更高纬度上我们称为超立方体，hypercube），而是单形。

单形可以认为是在N维空间里，选出一个最简单最紧凑的多边形，让它可以平铺整个N维空间。比如一维空间下的单形是等长的线段（1-单形），把这些线段收尾相连即可铺满整个一维空间。在二维空间下，单形是三角形（2-单形），可以把等腰三角形连接起来铺满整个平面。三维空间下的单形，即3-单形就是四面体。更高维空间的单形也是存在的。

使用单行的好处在于——它的顶点数很少，要远小于超立方体（hypercube）的顶点个数。总结起来，在N维空间下，超立方体的顶点数目是2n，而单形的顶点数目是n+1，这样在计算梯度噪声时可以大大减少需要计算的顶点权重数目。

在理解了单形后，Simplex噪声的计算过程其实和Perlin噪声基本一样。以二维空间下的为例。二维空间下的单形即是等边三角形，如下图所示。这些单形组成了一个单形网格结构，和Perlin噪声类似，这些网格顶点处也存储了伪随机梯度向量。 

【图片来自https://blog.csdn.net/candycat1992/article/details/50346469】

当输入一点后，找到该点所在的三角形（图中红色三角形），再找到该三角形三个顶点的梯度向量和每个顶点到输入点的差值向量，把每个顶点的梯度向量和插值向量做点乘，得到三个点乘结果。最后，把它们按权重进行叠加混合，这个权重与输入点到每个顶点的有关，即每个顶点的噪声贡献度为：

(r^2−|dist|^2)^4×dot(dist,grad)

其中，dist是输入点到该顶点的距离向量，grad是该点存储的伪随机梯度向量，r^2的取值是0.5或0.6。取0.5时可以保证没有不连续的间断点，在连续性并不那么明显时可以取0.6得到更好的视觉效果。在Perlin原始的论文中，r的取值是0.6（后面会讲到0.5的值是如何得到的）。当得到单形每个顶点的噪声贡献度后，就可以把它们相加起来。为了把结果归一到-1到1的范围，还需要在最后乘以一个系数。

目前需要解决的是如何找到输入点所在的单形？在计算Perlin噪声时，判断输入点所在的正方形是非常容易的，只需要对输入点下取整即可找到，这里可以把单形进行坐标偏斜（skewing），把平铺空间的单形变成一个新的网格结构，这个网格结构是由超立方体组成的，而每个超立方体又由一定数量的单形构成。对应的图解： 

【图片来自https://blog.csdn.net/candycat1992/article/details/50346469】

之前讲到的单形网格如上图中的红色网格所示，它们有一些等边三角形组成（注意到这些等边三角形是沿空间对角线排列的）。经过坐标倾斜后，它们变成了后面的黑色网格，这些网格由正方形组成，每个正方形是由之前两个等边三角形变形而来的三角形组成。这个把N维空间下的单形网格变形成新网格的公式如下：

在二维空间下，取n为2即可。这样变换之后，可以按照之前方法判断该点所在的超立方体，在二维下即为正方形。这样可以得到Simplex噪声的第一步：

1.坐标偏斜：把输入点坐标进行坐标偏斜，对坐标下取整得到输入点所在的超立方体xi=floor(x′),yi=floor(y′),...，同时可以得到小数部分xf=x′−xi,yf=y′−yi,...，这些小数部分可以进一步判断输入点所在的单形以及计算权重。

最终是要得到输入点所在的单形，而不是超立方体。因此需要继续做判断。还是如之前的图所示，经过坐标偏斜后，一个正方形由两个三角形组成，可以判断xf和yf之间的关系来判断输入点位于哪个三角形内，并得到该三角形的三个顶点。由此，Simplex噪声的第二步：

2 单形分割：将之前得到的(xf,yf,...)中的数值按降序排序，来决定输入点位于变形后的哪个单形内。这个单形的顶点是由按序排列的(0, 0, …, 0)到(1, 1, …, 1)中的n+1个顶点组成，共有n!种可能性。可以按下面的过程来得到这n+1个顶点：从零坐标(0, 0, …, 0)开始，找到当前最大的分量，在该分量位置加1，直至添加了所有分量。例如，对于二维空间来说，如果xf,yf满足xf>yf，那么对应的3个单形坐标为：首先找到(0, 0)，由于x分量比较大，因此下一个坐标是(1, 0)，接下来是y分量，坐标为(1, 1)；对于三维空间来说，如果xf,yf,zf满足xf>zf>yf，那么对应的4个单形坐标位：首先从(0, 0, 0)开始，接下来在x分量上加1得(1, 0, 0)，再在z分量上加1得(1, 0, 1)，最后在y分量上加1得(1, 1, 1)。这一步的算法复杂度即为排序复杂度O(n^2)。

找到了对应的单形后，后面的工作就比较简单了。首先找到该单形各个顶点上的伪随机梯度向量，这就是第三步：

3 梯度选取：在偏斜后的超立方体网格上获取该单形的各个顶点的伪随机梯度向量。

现在需要的东西基本都准备好了，最后一步就是计算所有单形顶点对输出的噪声贡献度。

4 贡献度取和：首先需要把单形顶点变回到之前由单形组成的单形网格。这一步需要使用第一步公式的逆函数来求得： 

我们由此可以得到输入点到这些单形顶点的位移向量。这些向量有两个用途，一个是为了和顶点梯度向量点乘，另一个是为了得到之前提到的距离值dist，来据此求得每个顶点对结果的贡献度：

(r^2−|dist|^2)^4×dot(dist,grad)

这里r^2为什么取0.5，由于要求经过第一步坐标偏斜后得到的网格宽度为1，因此可以倒推出在变形前单形网格中每个单形边的边长为sqrt(2/3)，这样一来单形每个顶点到对面边的距离（即高）的长度为sqrt(2)/2，它的平方即为0.5。不仅是二维，在其他维度下，每个单形顶点到对面边/面的距离都是0.5。

在把它们相加并返回最终结果前，需要保证返回值的范围在-1到1，这就需要计算n每个分量相加后的最大值。这个最大值可以在输入点在某一边中点时取得，此时|a|=sqrt(2)/2,|b|=|c|=1/sqrt(6)，得h=(0,1/3,1/3)。取最大值时，梯度和距离向量的点乘结果可以认为是两者模的乘积，而梯度模的最大值为sqrt(2)，因此最后和的最大值为：

(1/3)^4(1/sqrt(6))sqrt(2)*2≈1/70

因此，最后把结果乘以70。那么，如果r2取0.6，大约乘24.51。利用这个想法，可以在任意维度下计算最后的伸缩值，例如在三维下，单形，即正四面体的边长为sqrt(3)/2，当r^2取0.5时，最后大概需要乘以31.32。

虽然理解上Simplex噪声相比于Perlin噪声更难理解，但由于它的效果更好、速度更优，因此很多情况下会替代Perlin噪声。

Unity下的实现效果: 

单Simplex噪声（左下角）主要代码：

float SNoise(float2 p){
const float K1=0.366025404;  //(sqrt(3)-1)/2
const float K2=0.211324685; //(3-sqrt(3))/6

//将输入点进行坐标偏移，向下取整得到原点，转换到超立方体空间
float2 i=floor(p+(p.x+p.y)*K1);
//得到转换前输入点到原点距离（单形空间下）
float2 a=p-(i-(i.x+i.y)*K2);
//确定顶点在哪个三角形内
float2 o=(a.x<a.y)?float2(0.0,1.0):float2(1.0,0.0);
//得到转换前输入点到第二个顶点的距离
float2 b=a-o+K2;
//得到转换前输入点到第三个顶点的距离
float2 c=a-1+2*K2;  

//根据权重计算每个顶点的贡献度
float3 h=max((0.5-float3(dot(a,a),dot(b,b),dot(c,c))),0.0);
float3 n=h*h*h*h*float3(dot(a,Hash22(i)),dot(b,Hash22(i+o)),dot(c,Hash22(i+1.0)));

//乘以系数，做归一化处理
return dot(float3(70.0,70.0,70.0),n);
} 

参考： 
https://blog.csdn.net/candycat1992/article/details/50346469