前文学习了1D、2D Perlin噪声的生成算法,本文将更进一步,学习3D、4D Perlin噪声的生成。我们不仅需要2D的噪声面纹理,也需要3D的体纹理(后面将进一步学习),同时,我们可以利用3D体纹理来生成2D动态面纹理,也可以利用4D的超纹理来生成3D的动态体纹理,这就是我们为什么要学习4D Perlin噪声的原因。4D Perlin噪声的生成需要用到超晶格的概念,比较难以理解,不过原理是一样的。当然,也可以生成5D、6D Perlin噪声,但一般很少用到4D以上的噪声,这也不是我们的重点。3D,4D Perlin噪声计算复杂度是O(2^n),这对3D,4D来说是很大的计算量了,同时,为便于理解,我们的代码不做任何优化,对大尺寸的4D噪声来说,纯CPU计算开销将是非常可观的。通过对比,我们也可以看到GPU对这种计算密集型的体系优势还是很明显的。
(一)、3D Perlin Noise
跟1D,2D Perlin噪声一样,我们也采用三步来逐步实现之。
第一步,:生成伪随机梯度值PRNG 。在2D中,我们采用在单位圆上分布均匀的单位向量作为PRNG值。在3D中,我们采用立方体中心点到晶格12条边的中点向量作为PRNG值,这将产生12个向量。
在计算中,为了方便与操作,我们另外加了4个向量,使总的向量个数达到16个,但需要明白的是,加的这4个向量与原12个向量是重复的,如下代码,4个重复的向量我用空行隔开了。
new Vector3(1,1,0),
new Vector3(-1,1,0),
new Vector3(1,-1,0),
new Vector3(-1,-1,0),
new Vector3(1,0,1),
new Vector3(-1,0,1),
new Vector3(1,0,-1),
new Vector3(-1,0,-1),
new Vector3(0,1,1),
new Vector3(0,-1,1),
new Vector3(0,1,-1),
new Vector3(0,-1,-1),
new Vector3(1,1,0),
new Vector3(0,-1,1),
new Vector3(-1,1,0),
new Vector3(0,-1,-1)
我们在这里也没有使用DirectxMath库,为了使用Vector3,我们自定义了Vector3类。排列表permutation与2D是一样的,唯一不同的时,在使用时,我们又建了一个512长度的P数组,这样是为了防止在计算选择向量时索引越界溢出,稍后的代码中可以看到p[p[p[Xi] + Yi] + Zi]的索引值是可能大于255的。与2D算法一样,现在有了排列表,我们就可以根据输入值所在晶格选择相应的梯度值了。选择梯度值的代码如下:
aaa = p[p[p[Xi] + Yi] + Zi];
aba = p[p[p[Xi] + Yi + 1] + Zi];
aab = p[p[p[Xi] + Yi] + Zi + 1];
abb = p[p[p[Xi] + Yi + 1] + Zi + 1];
baa = p[p[p[Xi + 1] + Yi] + Zi];
bba = p[p[p[Xi + 1] + Yi + 1] + Zi];
bab = p[p[p[Xi + 1] + Yi] + Zi + 1];
bbb = p[p[p[Xi + 1] + Yi + 1] + Zi + 1];
这里的梯度值索引选择也是有规律的,我们用a代表为原值,而b代表为加1的值,这里的选择与之后的点乘要能对应上,不然,在晶格边界处将会出现不连续。这样我们就得到了晶格的8个顶点梯度值索引,我们还要让他与0xF相与,保证梯度值索引小于等15(因为我们只有16个梯度值)。有了梯度值索引后,我们通过gradients就可以得到梯度值了。
第二步:生成距离向量。与2D一样,距离向量的生成方法完全一样,这里就不在赘述,代码出如:
Xi = GetFloor((float)x * Frequency);
Yi = GetFloor((float)y * Frequency);
Zi = GetFloor((float)z * Frequency);
fracX = (float)x * Frequency - (float)Xi;
fracY = (float)y * Frequency - (float)Yi;
fracZ = (float)z * Frequency - (float)Zi;
得到了距离向量和梯度值,我们将这两个值作点积,这与2D算法也是一样的。
第三步:插值。因为是3D,总共有8个顶点,所以点积后得到8个值,得到8个值后我们对这8个值进行三次线性插值,这就是第三步了。
x1 = Interpolate(DotProduct(gradients[aaa & 0xF], fracX, fracY, fracZ), DotProduct(gradients[baa & 0xF], fracX - 1, fracY, fracZ), u);
x2 = Interpolate(DotProduct(gradients[aba & 0xF], fracX, fracY - 1, fracZ), DotProduct(gradients[bba & 0xF], fracX - 1, fracY - 1, fracZ), u);
y1 = Interpolate(x1, x2, v);
x1 = Interpolate(DotProduct(gradients[aab & 0xF], fracX, fracY, fracZ - 1), DotProduct(gradients[bab & 0xF], fracX - 1, fracY, fracZ - 1), u);
x2 = Interpolate(DotProduct(gradients[abb & 0xF], fracX, fracY - 1, fracZ - 1), DotProduct(gradients[bbb & 0xF], fracX - 1, fracY - 1, fracZ - 1), u);
y2 = Interpolate(x1, x2, v);
z = Interpolate(y1, y2, w);
至此,3D Perlin噪声生成完毕,我们在平面上我们可以将第三维作为时间梯度代入,这样就可以生成动态的二维图像了。CPU与GPU生成的图像分别为:
利用CPU生成的3D Perlin噪声
利用GPU生成的3D Perlin噪声
(二)、4D Perlin Noise
3D Perlin 噪声学习完了,4D Perlin噪声就简单了,算法都一样,但是还需要注意的是对原理的理解。4D Perlin噪声涉及到超立方体,超正方体有8个胞(立方体)、24个面(正方形)、32条棱和16个顶点,所以在扩展到4D的时候,梯度值采用超立方体到32条棱中点向量将会产生32个向量。只要理解这个了,其他的都比较好理解了。还是遵照我们前面的三步法来处理。
第一步,:生成伪随机梯度值PRNG 。在3D中,我们采用立方体中心点到晶格12条边的中点向量作为PRNG值,在4DK ,我们也采用超立方体中心点到晶格32条棱的中点向量作为PRNG,这将产生32个梯度向量值。
new Vector4(0,-1,-1,-1),
new Vector4(0,-1,-1,1),
new Vector4(0,-1,1,-1),
new Vector4(0,-1,1,1),
new Vector4(0,1,-1,-1),
new Vector4(0,1,-1,1),
new Vector4(0,1,1,-1),
new Vector4(0,1,1,1),
new Vector4(-1,-1,0,-1),
new Vector4(-1,1,0,-1),
new Vector4(1,-1,0,-1),
new Vector4(1,1,0,-1),
new Vector4(-1,-1,0,1),
new Vector4(-1,1,0,1),
new Vector4(1,-1,0,1),
new Vector4(1,1,0,1),
new Vector4(-1,0,-1,-1),
new Vector4(1,0,-1,-1),
new Vector4(-1,0,-1,1),
new Vector4(1,0,-1,1),
new Vector4(-1,0,1,-1),
new Vector4(1,0,1,-1),
new Vector4(-1,0,1,1),
new Vector4(1,0,1,1),
new Vector4(0,-1,-1,0),
new Vector4(0,-1,-1,0),
new Vector4(0,-1,1,0),
new Vector4(0,-1,1,0),
new Vector4(0,1,-1,0),
new Vector4(0,1,-1,0),
new Vector4(0,1,1,0),
new Vector4(0,1,1,0)
我们自定义了Vector4类。排列表permutation和算法处理与3D一样,现在有了排列表,我们就可以根据输入值所在晶格选择相应的梯度值了。选择梯度值的代码如下:
aaaa = p[p[p[p[Xi] + Yi] + Zi] + Wi];
aaba = p[p[p[p[Xi] + Yi] + Zi + 1] + Wi];
aaab = p[p[p[p[Xi] + Yi] + Zi] + Wi + 1];
aabb = p[p[p[p[Xi] + Yi] + Zi + 1] + Wi + 1];
abaa = p[p[p[p[Xi] + Yi + 1] + Zi] + Wi];
abba = p[p[p[p[Xi] + Yi + 1] + Zi + 1] + Wi];
abab = p[p[p[p[Xi] + Yi + 1] + Zi] + Wi + 1];
abbb = p[p[p[p[Xi] + Yi + 1] + Zi + 1] + Wi + 1];
baaa = p[p[p[p[Xi + 1] + Yi] + Zi] + Wi];
baba = p[p[p[p[Xi + 1] + Yi] + Zi + 1] + Wi];
baab = p[p[p[p[Xi + 1] + Yi] + Zi] + Wi + 1];
babb = p[p[p[p[Xi + 1] + Yi] + Zi + 1] + Wi + 1];
bbaa = p[p[p[p[Xi + 1] + Yi + 1] + Zi] + Wi];
bbba = p[p[p[p[Xi + 1] + Yi + 1] + Zi + 1] + Wi];
bbab = p[p[p[p[Xi + 1] + Yi + 1] + Zi] + Wi + 1];
bbbb = p[p[p[p[Xi + 1] + Yi + 1] + Zi + 1] + Wi + 1];
这样我们就得到了晶格的16个顶点梯度值索引,我们还要让他与0x1F(即31)相与,保证梯度值索引小于等31(因为我们只有32个梯度值)。有了梯度值索引后,我们通过gradients就可以得到梯度值了。
第二步:生成距离向量。与3D一样,距离向量的生成方法完全一样,这里就不在赘述,代码出如:
Xi = GetFloor((float)x * Frequency);
Yi = GetFloor((float)y * Frequency);
Zi = GetFloor((float)z * Frequency);
Wi = GetFloor((float)w * Frequency);
fracX = (float)x * Frequency - (float)Xi;
fracY = (float)y * Frequency - (float)Yi;
fracZ = (float)z * Frequency - (float)Zi;
fracW = (float)w * Frequency - (float)Wi;
得到了距离向量和梯度值,我们将这两个值作点积,这与3D算法也是一样的。
第三步:插值。因为是3D,总共有16个顶点,所以点积后得到16个值,得到16个值后我们对这16个值进行四次线性插值,这就是第三步了,代码很简单。
x1 = Interpolate(DotProduct(gradients[aaaa & 0x1F], fracX, fracY, fracZ, fracW), DotProduct(gradients[baaa & 0x1F], fracX - 1, fracY, fracZ, fracW), r);
x2 = Interpolate(DotProduct(gradients[aaba & 0x1F], fracX, fracY, fracZ-1,fracW), DotProduct(gradients[baba & 0x1F], fracX - 1, fracY, fracZ -1 ,fracW), r);
y1 = Interpolate(x1, x2, g);
x1 = Interpolate(DotProduct(gradients[aaab & 0x1F], fracX, fracY, fracZ,fracW -1), DotProduct(gradients[baab & 0x1F], fracX - 1, fracY, fracZ,fracW -1), r);
x2 = Interpolate(DotProduct(gradients[aabb & 0x1F], fracX, fracY, fracZ - 1,fracW -1), DotProduct(gradients[babb & 0x1F], fracX - 1, fracY, fracZ - 1, fracW -1), r);
y2 = Interpolate(x1, x2, g);
x1 = Interpolate(DotProduct(gradients[abaa & 0x1F], fracX, fracY-1, fracZ, fracW), DotProduct(gradients[bbaa & 0x1F], fracX - 1, fracY-1, fracZ, fracW), r);
x2 = Interpolate(DotProduct(gradients[abba & 0x1F], fracX, fracY-1, fracZ - 1, fracW), DotProduct(gradients[bbba & 0x1F], fracX - 1, fracY-1, fracZ - 1, fracW), r);
y3 = Interpolate(x1, x2, g);
x1 = Interpolate(DotProduct(gradients[abab & 0x1F], fracX, fracY-1, fracZ, fracW - 1), DotProduct(gradients[bbab & 0x1F], fracX - 1, fracY-1, fracZ, fracW - 1), r);
x2 = Interpolate(DotProduct(gradients[abbb & 0x1F], fracX, fracY-1, fracZ - 1, fracW - 1), DotProduct(gradients[bbbb & 0x1F], fracX - 1, fracY-1, fracZ - 1, fracW - 1), r);
y4 = Interpolate(x1, x2, g);
z1 = Interpolate(y1, y2, b);
z2 = Interpolate(y3, y4, b);
noise = Interpolate(z1, z2, a);
至此,4D Perlin噪声生成完毕,在平面上很难表现出4D,我们只能将第三维用作时间轴,第四维给取定值。但在3D空间中,我们可以将物体的三维空间座标作为参数,外加将第四维用作时间轴,这样就可以生成三维的图像了。CPU与GPU生成的图像分别为:
这是将第三维作时间,第四维了定值的CPU图像
这是将第四维作时间的GPU图像
(三)、小结
Perlin噪声3D,4D实现方法与2D差别不大,只要理解了其原理,算法就出来了,从上面我们实现的方法也可以看到,3D,4D Perlin噪声沿用了1D,2D是的步骤。Perlin噪声的算法复杂是O(2^n),到了4D,明显可以感觉到CPU计算的压力。
(四)、Perlin写的Java参考代码
在3D,4D Perlin噪声的实现中,我们采用了最直观的方式,没有采用优化技巧,对比Perlin本人写的Java参考代码,Ken Perlin在处理点乘的时候采用了一点技巧,使代码难于理解了,这里我们稍作解释,我们对Perlin 3D 代码中的grad函数分析一下。我们只要明白grad()函数的作用在于计算随机选取的梯度向量与顶点位置向量的点积。Ken Perlin巧妙地使用了位翻转(bit-flipping)技巧来实现,他的原代码如下:
public static double grad(int hash, double x, double y, double z) {
int h = hash & 15;
double u = h < 8 /* 0b1000 */ ? x : y;
double v;
if(h < 4 /* 0b0100 */)
v = y;
else if(h == 12 /* 0b1100 */ || h == 14 /* 0b1110*/)
v = x;
else
v = z;
return ((h&1) == 0 ? u : -u)+((h&2) == 0 ? v : -v);
}
这段代码可以改写成如下:
public static double grad(int hash, double x, double y, double z)
{
switch(hash & 0xF)
{
case 0x0: return x + y;
case 0x1: return -x + y;
case 0x2: return x - y;
case 0x3: return -x - y;
case 0x4: return x + z;
case 0x5: return -x + z;
case 0x6: return x - z;
case 0x7: return -x - z;
case 0x8: return y + z;
case 0x9: return -y + z;
case 0xA: return y - z;
case 0xB: return -y - z;
case 0xC: return y + x;
case 0xD: return -y + z;
case 0xE: return y - x;
case 0xF: return -y - z;
default: return 0; // never happens
}
}
在代码中,hash & 0xF 即是将得到的梯度向量索引值转到0-15范围内,返回值其实是个点乘值, x + y =dot( (x,y,z),(1,1,0)) ,(x,y,z)可以看成是(x,y,z)-(0,0,0),即是输入点到最低点(0,0,0)的距离向量,(1,1,0)就是hash & 0xF=0x0时选取的梯度向量,所以最后的结果为x+y。再分析一下调用函数的方式:
grad (aaa, xf , yf , zf),
grad (baa, xf-1, yf , zf),
grad (aba, xf , yf-1, zf),
grad (bba, xf-1, yf-1, zf),
grad (aab, xf , yf , zf-1),
grad (bab, xf-1, yf , zf-1),
grad (abb, xf , yf-1, zf-1),
grad (bbb, xf-1, yf-1, zf-1),
这切就都很清楚了,剩下的就可以分析理解了。4D也一样,Ken Perlin采用了同样的办法来处理以达到减少计算量的目的。
3D Perlin噪声,Ken Perlin Java 参考代码
public final class ImprovedNoise {
static public double noise(double x, double y, double z) {
int X = (int)Math.floor(x) & 255, // FIND UNIT CUBE THAT
Y = (int)Math.floor(y) & 255, // CONTAINS POINT.
Z = (int)Math.floor(z) & 255;
x -= Math.floor(x); // FIND RELATIVE X,Y,Z
y -= Math.floor(y); // OF POINT IN CUBE.
z -= Math.floor(z);
double u = fade(x), // COMPUTE FADE CURVES
v = fade(y), // FOR EACH OF X,Y,Z.
w = fade(z);
int A = p[X ]+Y, AA = p[A]+Z, AB = p[A+1]+Z, // HASH COORDINATES OF
B = p[X+1]+Y, BA = p[B]+Z, BB = p[B+1]+Z; // THE 8 CUBE CORNERS,
return lerp(w, lerp(v, lerp(u, grad(p[AA ], x , y , z ), // AND ADD
grad(p[BA ], x-1, y , z )), // BLENDED
lerp(u, grad(p[AB ], x , y-1, z ), // RESULTS
grad(p[BB ], x-1, y-1, z ))),// FROM 8
lerp(v, lerp(u, grad(p[AA+1], x , y , z-1 ), // CORNERS
grad(p[BA+1], x-1, y , z-1 )), // OF CUBE
lerp(u, grad(p[AB+1], x , y-1, z-1 ),
grad(p[BB+1], x-1, y-1, z-1 ))));
}
static double fade(double t) { return t * t * t * (t * (t * 6 - 15) + 10); }
static double lerp(double t, double a, double b) { return a + t * (b - a); }
static double grad(int hash, double x, double y, double z) {
int h = hash & 15; // CONVERT LO 4 BITS OF HASH CODE
double u = h<8 ? x : y, // INTO 12 GRADIENT DIRECTIONS.
v = h<4 ? y : h==12||h==14 ? x : z;
return ((h&1) == 0 ? u : -u) + ((h&2) == 0 ? v : -v);
}
static final int p[] = new int[512], permutation[] = { 151,160,137,91,90,15,
131,13,201,95,96,53,194,233,7,225,140,36,103,30,69,142,8,99,37,240,21,10,23,
190, 6,148,247,120,234,75,0,26,197,62,94,252,219,203,117,35,11,32,57,177,33,
88,237,149,56,87,174,20,125,136,171,168, 68,175,74,165,71,134,139,48,27,166,
77,146,158,231,83,111,229,122,60,211,133,230,220,105,92,41,55,46,245,40,244,
102,143,54, 65,25,63,161, 1,216,80,73,209,76,132,187,208, 89,18,169,200,196,
135,130,116,188,159,86,164,100,109,198,173,186, 3,64,52,217,226,250,124,123,
5,202,38,147,118,126,255,82,85,212,207,206,59,227,47,16,58,17,182,189,28,42,
223,183,170,213,119,248,152, 2,44,154,163, 70,221,153,101,155,167, 43,172,9,
129,22,39,253, 19,98,108,110,79,113,224,232,178,185, 112,104,218,246,97,228,
251,34,242,193,238,210,144,12,191,179,162,241, 81,51,145,235,249,14,239,107,
49,192,214, 31,181,199,106,157,184, 84,204,176,115,121,50,45,127, 4,150,254,
138,236,205,93,222,114,67,29,24,72,243,141,128,195,78,66,215,61,156,180
};
static { for (int i=0; i < 256 ; i++) p[256+i] = p[i] = permutation[i]; }
}
4D Perlin噪声,Ken Perlin Java 参考代码
// JAVA REFERENCE IMPLEMENTATION OF IMPROVED NOISE IN 4D - COPYRIGHT 2002 KEN PERLIN.
public final class ImprovedNoise4D {
static public double noise(double x, double y, double z, double w) {
int X = (int)Math.floor(x) & 255, // FIND UNIT HYPERCUBE
Y = (int)Math.floor(y) & 255, // THAT CONTAINS POINT.
Z = (int)Math.floor(z) & 255,
W = (int)Math.floor(w) & 255;
x -= Math.floor(x); // FIND RELATIVE X,Y,Z,W
y -= Math.floor(y); // OF POINT IN CUBE.
z -= Math.floor(z);
w -= Math.floor(w);
double a = fade(x), // COMPUTE FADE CURVES
b = fade(y), // FOR EACH OF X,Y,Z,W.
c = fade(z),
d = fade(w);
int A = p[X ]+Y, AA = p[A]+Z, AB = p[A+1]+Z, // HASH COORDINATES OF
B = p[X+1]+Y, BA = p[B]+Z, BB = p[B+1]+Z, // THE 16 CORNERS OF
AAA = p[AA]+W, AAB = p[AA+1]+W, // THE HYPERCUBE.
ABA = p[AB]+W, ABB = p[AB+1]+W,
BAA = p[BA]+W, BAB = p[BA+1]+W,
BBA = p[BB]+W, BBB = p[BB+1]+W;
return lerp(d, // INTERPOLATE DOWN.
lerp(c,lerp(b,lerp(a,grad(p[AAA ], x , y , z , w),
grad(p[BAA ], x-1, y , z , w)),
lerp(a,grad(p[ABA ], x , y-1, z , w),
grad(p[BBA ], x-1, y-1, z , w))),
lerp(b,lerp(a,grad(p[AAB ], x , y , z-1, w),
grad(p[BAB ], x-1, y , z-1, w)),
lerp(a,grad(p[ABB ], x , y-1, z-1, w),
grad(p[BBB ], x-1, y-1, z-1, w)))),
lerp(c,lerp(b,lerp(a,grad(p[AAA+1], x , y , z , w-1),
grad(p[BAA+1], x-1, y , z , w-1)),
lerp(a,grad(p[ABA+1], x , y-1, z , w-1),
grad(p[BBA+1], x-1, y-1, z , w-1))),
lerp(b,lerp(a,grad(p[AAB+1], x , y , z-1, w-1),
grad(p[BAB+1], x-1, y , z-1, w-1)),
lerp(a,grad(p[ABB+1], x , y-1, z-1, w-1),
grad(p[BBB+1], x-1, y-1, z-1, w-1)))));
}
static double fade(double t) { return t * t * t * (t * (t * 6 - 15) + 10); }
static double lerp(double t, double a, double b) { return a + t * (b - a); }
static double grad(int hash, double x, double y, double z, double w) {
int h = hash & 31; // CONVERT LO 5 BITS OF HASH TO 32 GRAD DIRECTIONS.
double a=y,b=z,c=w; // X,Y,Z
switch (h >> 3) { // OR, DEPENDING ON HIGH ORDER 2 BITS:
case 1: a=w;b=x;c=y;break; // W,X,Y
case 2: a=z;b=w;c=x;break; // Z,W,X
case 3: a=y;b=z;c=w;break; // Y,Z,W
}
return ((h&4)==0 ? -a:a) + ((h&2)==0 ? -b:b) + ((h&1)==0 ? -c:c);
}
static final int p[] = new int[512], permutation[] = { 151,160,137,91,90,15,
131,13,201,95,96,53,194,233,7,225,140,36,103,30,69,142,8,99,37,240,21,10,23,
190, 6,148,247,120,234,75,0,26,197,62,94,252,219,203,117,35,11,32,57,177,33,
88,237,149,56,87,174,20,125,136,171,168, 68,175,74,165,71,134,139,48,27,166,
77,146,158,231,83,111,229,122,60,211,133,230,220,105,92,41,55,46,245,40,244,
102,143,54, 65,25,63,161, 1,216,80,73,209,76,132,187,208, 89,18,169,200,196,
135,130,116,188,159,86,164,100,109,198,173,186, 3,64,52,217,226,250,124,123,
5,202,38,147,118,126,255,82,85,212,207,206,59,227,47,16,58,17,182,189,28,42,
223,183,170,213,119,248,152, 2,44,154,163, 70,221,153,101,155,167, 43,172,9,
129,22,39,253, 19,98,108,110,79,113,224,232,178,185, 112,104,218,246,97,228,
251,34,242,193,238,210,144,12,191,179,162,241, 81,51,145,235,249,14,239,107,
49,192,214, 31,181,199,106,157,184, 84,204,176,115,121,50,45,127, 4,150,254,
138,236,205,93,222,114,67,29,24,72,243,141,128,195,78,66,215,61,156,180
};
static { for (int i=0; i < 256 ; i++) p[256+i] = p[i] = permutation[i]; }
}
(四)、代码下载
VS2015+C#版
链接: http://download.csdn.net/download/yolon3000/9977386
Unity2017.1.1f1+cg版
链接: http://download.csdn.net/download/yolon3000/9977389
参考文献:
https://mrl.nyu.edu/~perlin/noise/ImprovedNoise4D.java
http://mrl.nyu.edu/~perlin/noise/