Perline Noise

 

http://www.twinklingstar.cn/2015/2581/classical-perlin-noise/

 

Perlin噪声

 

本篇文章谢绝转载，也禁止用于任何商业目的

1. Perlin噪声

Ken Perlin在1983年提出了一种渐变噪声，Perlin在1985年的SIGGRAPH[1]有该算法的描述，称之为经典Perlin噪声（Classical Perlin Noise）。为了简化计算，方便使用硬件实现，并解决经典Perlin噪声中存在的错误，到2001年，Ken Perlin再次对原始的噪声算法进行了改进，称之为Simplex噪声[2]（Simplex Noise），这两种算法都可以称为Perlin噪声。但是，我们有时候也把分形噪声也称为Perlin噪声[3]，甚至在严肃的学术论文中都有这种称法[4]。为了避免歧义，本文指的Perlin噪声特指经典Perlin噪声和Simplex噪声。

Stefan Gustavson[8]指出：Simplex噪声有更小的算法复杂度，要求更少的乘法，在N维空间上，经典Perlin噪声的算法复杂度为log(2N)，但是Simplex噪声的算法复杂度为log(N2)等优点。虽然Stefan Gustavson[8]提供了对Simplex算法的注解，但是我依然不能理解Simplex噪声背后的数学原理，对Simplex噪声不作进一步阐述。

2. 经典Perlin噪声

经典Perlin噪声是Ken Perlin在1983年提出的噪声，Ken Perlin提供了一维、二维、三维算法的C实现[5]，我们无法仅从他提供的代码理解其数学原理，Matt Zucker[6]从数学原理的角度，对经典Perlin噪声进行了解读。本节将详细介绍二维、三维的Pernlin噪声的数学原理，算法的C语言实现源码如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
/* coherent noise function over 1, 2 or 3 dimensions */
/* (copyright Ken Perlin) */
 
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
 
#define B 0x100
#define BM 0xff
 
#define N 0x1000
#define NP 12   /* 2^N */
#define NM 0xfff
 
static p[B + B + 2];
static float g3[B + B + 2][3];
static float g2[B + B + 2][2];
static float g1[B + B + 2];
static start = 1;
 
static void init(void);
 
#define s_curve(t) ( t * t * (3. - 2. * t) )
 
#define lerp(t, a, b) ( a + t * (b - a) )
 
#define setup(i,b0,b1,r0,r1)\
    t = vec[i] + N;\
    b0 = ((int)t) & BM;\
    b1 = (b0+1) & BM;\
    r0 = t - (int)t;\
    r1 = r0 - 1.;
 
double noise1(double arg)
{
    int bx0, bx1;
    float rx0, rx1, sx, t, u, v, vec[1];
 
    vec[0] = arg;
    if (start) {
        start = 0;
        init();
    }
 
    setup(0, bx0,bx1, rx0,rx1);
 
    sx = s_curve(rx0);
 
    u = rx0 * g1[ p[ bx0 ] ];
    v = rx1 * g1[ p[ bx1 ] ];
 
    return lerp(sx, u, v);
}
 
float noise2(float vec[2])
{
    int bx0, bx1, by0, by1, b00, b10, b01, b11;
    float rx0, rx1, ry0, ry1, *q, sx, sy, a, b, t, u, v;
    register i, j;
 
    if (start) {
        start = 0;
        init();
    }
 
    setup(0, bx0,bx1, rx0,rx1);
    setup(1, by0,by1, ry0,ry1);
 
    i = p[ bx0 ];
    j = p[ bx1 ];
 
    b00 = p[ i + by0 ];
    b10 = p[ j + by0 ];
    b01 = p[ i + by1 ];
    b11 = p[ j + by1 ];
 
    sx = s_curve(rx0);
    sy = s_curve(ry0);
 
#define at2(rx,ry) ( rx * q[0] + ry * q[1] )
 
    q = g2[ b00 ] ; u = at2(rx0,ry0);
    q = g2[ b10 ] ; v = at2(rx1,ry0);
    a = lerp(sx, u, v);
 
    q = g2[ b01 ] ; u = at2(rx0,ry1);
    q = g2[ b11 ] ; v = at2(rx1,ry1);
    b = lerp(sx, u, v);
 
    return lerp(sy, a, b);
}
 
float noise3(float vec[3])
{
    int bx0, bx1, by0, by1, bz0, bz1, b00, b10, b01, b11;
    float rx0, rx1, ry0, ry1, rz0, rz1, *q, sy, sz, a, b, c, d, t, u, v;
    register i, j;
 
    if (start) {
        start = 0;
        init();
    }
 
    setup(0, bx0,bx1, rx0,rx1);
    setup(1, by0,by1, ry0,ry1);
    setup(2, bz0,bz1, rz0,rz1);
 
    i = p[ bx0 ];
    j = p[ bx1 ];
 
    b00 = p[ i + by0 ];
    b10 = p[ j + by0 ];
    b01 = p[ i + by1 ];
    b11 = p[ j + by1 ];
 
    t  = s_curve(rx0);
    sy = s_curve(ry0);
    sz = s_curve(rz0);
 
#define at3(rx,ry,rz) ( rx * q[0] + ry * q[1] + rz * q[2] )
 
    q = g3[ b00 + bz0 ] ; u = at3(rx0,ry0,rz0);
    q = g3[ b10 + bz0 ] ; v = at3(rx1,ry0,rz0);
    a = lerp(t, u, v);
 
    q = g3[ b01 + bz0 ] ; u = at3(rx0,ry1,rz0);
    q = g3[ b11 + bz0 ] ; v = at3(rx1,ry1,rz0);
    b = lerp(t, u, v);
 
    c = lerp(sy, a, b);
 
    q = g3[ b00 + bz1 ] ; u = at3(rx0,ry0,rz1);
    q = g3[ b10 + bz1 ] ; v = at3(rx1,ry0,rz1);
    a = lerp(t, u, v);
 
    q = g3[ b01 + bz1 ] ; u = at3(rx0,ry1,rz1);
    q = g3[ b11 + bz1 ] ; v = at3(rx1,ry1,rz1);
    b = lerp(t, u, v);
 
    d = lerp(sy, a, b);
 
    return lerp(sz, c, d);
}
 
static void normalize2(float v[2])
{
    float s;
 
    s = sqrt(v[0] * v[0] + v[1] * v[1]);
    v[0] = v[0] / s;
    v[1] = v[1] / s;
}
 
static void normalize3(float v[3])
{
    float s;
 
    s = sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
    v[0] = v[0] / s;
    v[1] = v[1] / s;
    v[2] = v[2] / s;
}
 
static void init(void)
{
    int i, j, k;
 
    for (i = 0 ; i < B ; i++) {
        p[i] = i;
 
        g1[i] = (float)((random() % (B + B)) - B) / B;
 
        for (j = 0 ; j < 2 ; j++)
            g2[i][j] = (float)((random() % (B + B)) - B) / B;
        normalize2(g2[i]);
 
        for (j = 0 ; j < 3 ; j++)
            g3[i][j] = (float)((random() % (B + B)) - B) / B;
        normalize3(g3[i]);
    }
 
    while (--i) {
        k = p[i];
        p[i] = p[j = random() % B];
        p[j] = k;
    }
 
    for (i = 0 ; i < B + 2 ; i++) {
        p[B + i] = p[i];
        g1[B + i] = g1[i];
        for (j = 0 ; j < 2 ; j++)
            g2[B + i][j] = g2[i][j];
        for (j = 0 ; j < 3 ; j++)
            g3[B + i][j] = g3[i][j];
    }
}
用简单的伪码如下所示（不理解，可以直接忽略）：
  1. 给定点输入点p；
  2. 对于每一个与点p相邻的方格端点：
    1. 挑选一个伪随机梯度向量；
    2. 计算点积；
  3. 每个维度上，均采用缓和曲线，计算出加权平均值，缓和曲线可以选用3t2−2t3。

2.1. 二维空间 

1. 找到方格点

先考虑二维的情况，噪声函数用式子Noise(x,y)表示，其中x,y是浮点数，(x,y)表示一个点，设对变量x下取整用式子floor(x)表示，即floor(4.34)=4，令：
x0=floor(x)
y0=floor(y)
x1=x0+1
y1=y0+1
用图表示点(x,y)和(x0,y0),(x1,y0),(x0,y1),(x1,y1)，如图1.所示。
pernlin_noise_1

图1. 包围点(x,y)的方格

2. 挑选伪随机梯度向量

对于方格的每个端点，分配一个伪随机梯度（Pseudorandom Gradient）向量，这是一个在二维空间上单位长度的向量。所谓的伪随机，是指相同的输入值将产生相同的输出值，并不是真正意义上的随机。可以通过下面的方法挑选伪随机梯度向量[7]，n规定表的长度：
  1. 生成一个随机排列表，P[n]；
  2. 生成一个随机梯度表，G[n]；
  3. G(i,j)=G[(i+P[j])modn]；

在实际应用中，我们不需要求模运算，设n=256，通过位运算，留下低8位的值就行，即G(i,j)=G[(i+P[j])&0xff]。由方格的四个端点，我们可以得出对应的伪随机梯度向量，如图2所示。

pernlin_noise_2

图2. 伪随机梯度向量

3. 计算点积

接着，生成一个由点(x,y)到方格四个端点的向量，如图3所示，四个向量分别为：
d⃗ 00=(x,y)−(x0,y0)
d⃗ 10=(x,y)−(x1,y0)
d⃗ 01=(x,y)−(x0,y1)
d⃗ 11=(x,y)−(x1,y1)
这四个向量分别与对应端点的伪随机梯度向量点积，得到：
s=G(x0,y0)⋅((x,y)−(x0,y0))
t=G(x1,y0)⋅((x,y)−(x1,y0))
u=G(x0,y1)⋅((x,y)−(x0,y1))
v=G(x1,y1)⋅((x,y)−(x1,y1))
它的几何意义，就是将向量d⃗ ij（i,j∈{1,2}）投影到对应的伪随机梯度向量，显然有s,t,u,v∈[0,1]。不理解其几何意义的作用，但是可以明确的一点是：s,t,u,v值由点(x,y)和对应的伪随机梯度向量唯一确定。
pernlin_noise_3

图3. 由点(x,y)到方格四个端点的向量

4. 计算加权平均值

接下来，考虑如何计算这几个值的加权平均，得到最终的噪声标量。
这里引入缓和曲线（Ease Curve）的概念，缓和曲线使噪声显得更加真实，采用的缓和曲线函数为：f(x)=3x2−2x3，如图4所示。
pernlin_noise_4

图4. f(x)=3x2−2x3曲线

把x−x0代入缓和函数f(x)中，得到比例系数sx=f(x−x0)，通过线性插值得到：
a=s+(t−s)⋅sx
b=u+(v−u)⋅sx
再把y−y0代入缓和函数f(x)中，得到比例系数sy=f(y−y0)，通过线性插值得到：

Noise(x,y)=a+(b−a)⋅sy

2.2. 三维空间   

在三维空间上，噪声函数用式子Noise(x,y,z)表示，与二维噪声相似，同样可以分成四个阶段来处理。三维空间上的方格是一个正方体，这里特别说明下点积的计算和加权平均值的计算。
分别把x−x0,y−y0,z−z0代入缓和函数f(x)中，得到比例系数：
sx=f(x−x0),sy=f(y−y0),sz=f(z−z0)
pernlin_noise_5

图5. 三维空间上的方格

如图5所示，计算点积
s0=G(x0,y0,z0)⋅((x,y,z)−(x0,y0,z0))
t0=G(x1,y0,z0)⋅((x,y,z)−(x1,y0,z0))
u0=G(x0,y1,z0)⋅((x,y,z)−(x0,y1,z0))
v0=G(x1,y1,z0)⋅((x,y,z)−(x1,y1,z0))
s1=G(x0,y0,z1)⋅((x,y,z)−(x0,y0,z1))
t1=G(x1,y0,z1)⋅((x,y,z)−(x1,y0,z1))
u1=G(x0,y1,z1)⋅((x,y,z)−(x0,y1,z1))
v1=G(x1,y1,z1)⋅((x,y,z)−(x1,y1,z1))
通过线性插值，计算出红色实心圆的数值：
a0=s0+(t0−s0)⋅sx
b0=u0+(v0−u0)⋅sx
a1=s1+(t1−s1)⋅sx
b1=u1+(v1−u1)⋅sx
通过线性插值，计算出紫色实心圆的数值：
z0=a0+(b0−a0)⋅sy
z1=a1+(b1−a1)⋅sy
通过线性插值，计算出绿色实心圆的数值，即是最终的结果：

Noise(x,y,z)=z0+(z1−z0)⋅sz

2.3. 实现代码

最后，提供一个经典Perlin噪声的C++实现代码，如下所示：

SrClassicalPerlinNoise.h

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
/************************************************************************      
\link   www.twinklingstar.cn
\author twinklingstar
\file   SrClassicalPerlinNoise.h
\date   2015/07/15
****************************************************************************/
#ifndef SR_CLASSIC_PERLIN_NOISE_H_
#define SR_CLASSIC_PERLIN_NOISE_H_
 
#define SR_PERLIN_N         0x100
#define SR_PERLIN_MASK      0xff
#define SR_PERLIN_PI        3.141592653589793
 
typedef float           real;
typedef unsigned char   ubyte;
 
 
class SrPerlinNoise
{
public:
    SrPerlinNoise(){}
protected:
    real    ease(real x);
 
    int     permutate(int x);
 
    real    lerp(real t, real value0, real value1);
protected:
    static ubyte mP[SR_PERLIN_N];
};
 
 
class SrClassicalPerlinNoise1D:public SrPerlinNoise
{
public:
    static SrClassicalPerlinNoise1D* create();
 
    real    noise(real x);
private:
    int     index(int ix);
 
    real    lattice(int ix, real fx);
 
    static void init();
 
    SrClassicalPerlinNoise1D(){}
    SrClassicalPerlinNoise1D(const SrClassicalPerlinNoise1D&);
    SrClassicalPerlinNoise1D& operator = (const SrClassicalPerlinNoise1D&);
private:
    static real mG[SR_PERLIN_N];
    static bool mIsInit;
};
 
class SrClassicalPerlinNoise2D:public SrPerlinNoise
{
public:
    static SrClassicalPerlinNoise2D* create();
 
    real noise(real x, real y);
private:
    int     index(int ix, int iy);
 
    real    lattice(int ix, int iy, real fx, real fy);
 
    static void init();
 
    SrClassicalPerlinNoise2D(){}
    SrClassicalPerlinNoise2D(const SrClassicalPerlinNoise2D&);
    SrClassicalPerlinNoise2D& operator = (const SrClassicalPerlinNoise2D&);
private:
    static real mG[SR_PERLIN_N][2];
    static bool mIsInit;
};
 
class SrClassicalPerlinNoise3D:public SrPerlinNoise
{
public:
    static SrClassicalPerlinNoise3D* create();
 
    real    noise(real x, real y, real z);
;
private:
    int     index(int ix, int iy, int iz);
 
    real    lattice(int ix, int iy, int iz, float fx, float fy, float fz);
 
    static void init();
 
    SrClassicalPerlinNoise3D(){}
    SrClassicalPerlinNoise3D(const SrClassicalPerlinNoise3D&);
    SrClassicalPerlinNoise3D& operator = (const SrClassicalPerlinNoise3D&);
private:
    static real mG[SR_PERLIN_N][3];
    static bool mIsInit;
};
 
 
#endif

SrClassicalPerlinNoise.cpp

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
/************************************************************************      
\link   www.twinklingstar.cn
\author twinklingstar
\file   SrClassicalPerlinNoise.cpp
\date   2015/07/15
****************************************************************************/
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "SrClassicalPerlinNoise.h"
 
 
/************************************************************************/
/*                     class SrPerlinNoise                              */
/************************************************************************/
 
real SrPerlinNoise::ease(real x)
{
    return (real)(x * x * (3 - 2 * x));
}
 
int SrPerlinNoise::permutate(int x)
{
    return mP[x & SR_PERLIN_MASK];
}
 
real SrPerlinNoise::lerp(real t, real value0, real value1)
{
    return (real)(value0 + t * (value1 - value0));
}
 
ubyte SrPerlinNoise::mP[SR_PERLIN_N] = {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
};
 
/************************************************************************/
/*               class SrClassicalPerlinNoise1D                         */
/************************************************************************/
SrClassicalPerlinNoise1D* SrClassicalPerlinNoise1D::create()
{
    if(!mIsInit)
    {
        mIsInit = true;
        init();
    }
    return (new SrClassicalPerlinNoise1D());
}
 
real SrClassicalPerlinNoise1D::noise(real x)
{
    int x0 = (int)floor(x);
    real fx0 = x - x0;
    real fx1 = fx0 - 1;
 
    real sx = ease(fx0);
 
    real u = lattice(x0    , fx0);
    real v = lattice(x0 + 1, fx1);
 
    return lerp(sx, u, v);
}
int SrClassicalPerlinNoise1D::index(int ix)
{
    return permutate(ix);
}
 
real SrClassicalPerlinNoise1D::lattice(int ix, real fx)
{
    int indx = index(ix);
    return (real)(mG[indx] * fx );
}
 
void SrClassicalPerlinNoise1D::init()
{
    int i;
    for (i = 0 ; i < SR_PERLIN_N ; i++)
    {
        mG[i] = (float)((rand() % (SR_PERLIN_N + SR_PERLIN_N)) - SR_PERLIN_N) / SR_PERLIN_N;
    }
}
 
bool SrClassicalPerlinNoise1D::mIsInit = false;
 
real SrClassicalPerlinNoise1D::mG[SR_PERLIN_N] = {};
 
/************************************************************************/
/*              class SrClassicalPerlinNoise2D                          */
/************************************************************************/
 
SrClassicalPerlinNoise2D* SrClassicalPerlinNoise2D::create()
{
    if(!mIsInit)
    {
        mIsInit = true;
        init();
    }
    return (new SrClassicalPerlinNoise2D());
}
 
real SrClassicalPerlinNoise2D::noise(real x, real y)
{
    int x0 = (int)floor(x);
    real fx0 = x - x0;
    real fx1 = fx0 - 1;
    real sx = ease(fx0);
 
    int y0 = (int)floor(y);
    real fy0 = y - y0;
    real fy1 = fy0 - 1;
    real sy = ease(fy0);
 
    real s = lattice(x0    ,y0    , fx0, fy0);
    real t = lattice(x0 + 1,y0    , fx1, fy0);
    real u = lattice(x0    ,y0 + 1, fx0, fy1);
    real v = lattice(x0 + 1,y0 + 1, fx1, fy1);
 
    real a = lerp(sx, s, t);
    real b = lerp(sx, u, v);
    return lerp(sy, a, b);
}
 
int SrClassicalPerlinNoise2D::index(int ix, int iy)
{
    return permutate(ix + permutate(iy));
}
 
real SrClassicalPerlinNoise2D::lattice(int ix, int iy, real fx, real fy)
{
    int indx = index(ix, iy);
    return (real)(mG[indx][0] * fx + mG[indx][1] * fy);
}
 
 
 
void SrClassicalPerlinNoise2D::init()
{
    int i, j;
    real s;
    for (i = 0 ; i < SR_PERLIN_N ; i++)
    {
        for (j = 0 ; j < 2 ; j++)
            mG[i][j] = (float)((rand() % (SR_PERLIN_N + SR_PERLIN_N)) - SR_PERLIN_N) / SR_PERLIN_N;
 
        s = sqrt((real)(mG[i][0] * mG[i][0] + mG[i][1] * mG[i][1]));
        mG[i][0] = mG[i][0] / s;
        mG[i][1] = mG[i][1] / s;
    }
}
 
bool SrClassicalPerlinNoise2D::mIsInit = false;
 
real SrClassicalPerlinNoise2D::mG[SR_PERLIN_N][2] = {};
 
/************************************************************************/
/*                class SrClassicalPerlinNoise3D                        */
/************************************************************************/
 
SrClassicalPerlinNoise3D* SrClassicalPerlinNoise3D::create()
{
    if(!mIsInit)
    {
        mIsInit = true;
        init();
    }
    return (new SrClassicalPerlinNoise3D());
}
 
real SrClassicalPerlinNoise3D::noise(real x, real y, real z)
{
    int x0 = (int)floor(x);
    real fx0 = x - x0;
    real fx1 = fx0 - 1;
    real wx = ease(fx0);
 
    int y0 = (int)floor(y);
    real fy0 = y - y0;
    real fy1 = fy0 - 1;
    real wy = ease(fy0);
 
    int z0 = (int)floor(z);
    real fz0 = z - z0;
    real fz1 = fz0 - 1;
    real wz = ease(fz0);
 
    real vx0 = lattice(x0    , y0, z0, fx0, fy0, fz0);
    real vx1 = lattice(x0 + 1, y0, z0, fx1, fy0, fz0);
    real vy0 = lerp(wx, vx0, vx1);
 
    vx0 = lattice(x0    , y0 + 1, z0, fx0, fy1, fz0);
    vx1 = lattice(x0 + 1, y0 + 1, z0, fx1, fy1, fz0);
    real vy1 = lerp(wx, vx0, vx1);
 
    real vz0 = lerp(wy, vy0, vy1);
 
    vx0 = lattice(x0    , y0, z0 + 1, fx0, fy0, fz1);
    vx1 = lattice(x0 + 1, y0, z0 + 1, fx1, fy0, fz1);
    vy0 = lerp(wx, vx0, vx1);
 
    vx0 = lattice(x0    , y0 + 1, z0 + 1, fx0, fy1, fz1);
    vx1 = lattice(x0 + 1, y0 + 1, z0 + 1, fx1, fy1, fz1);
    vy1 = lerp(wx, vx0, vx1);
 
    real vz1 = lerp(wy, vy0, vy1);
 
    return lerp(wz, vz0, vz1);
}
 
int SrClassicalPerlinNoise3D::index(int ix, int iy, int iz)
{
    return permutate(ix + permutate(iy + permutate(iz)));
}
 
real SrClassicalPerlinNoise3D::lattice(int ix, int iy, int iz, float fx, float fy, float fz)
{
    int indx = index(ix, iy, iz);
    return (real)(mG[indx][0] * fx + mG[indx][1] * fy + mG[indx][2] * fz);
}
 
void SrClassicalPerlinNoise3D::init()
{
    int i;
    for (i = 0 ; i < SR_PERLIN_N ; i++)
    {
        real z = 1.0f - 2.0f * ((real)::rand()/((real)RAND_MAX+1));
        real r = (real)sqrt(1.0f - z * z);
        real theta = 2 * (real)SR_PERLIN_PI * ((real)::rand()/((real)RAND_MAX+1));
 
        mG[i][0]    = r * (real)cos(theta);
        mG[i][1]    = r * (real)sin(theta);
        mG[i][2]    = z;
    }
}
 
 
bool SrClassicalPerlinNoise3D::mIsInit = false;
 
real SrClassicalPerlinNoise3D::mG[SR_PERLIN_N][3] = {};

3. 地形的应用

float3x3[9]描述了一种基于经典Perlin噪声的地形高度生成算法，具体参见链接对该算法的描述，这里提供DirectX的地形实现类的代码，整个工程代码参见github链接（https://github.com/twinklingstar20/twinklingstar_cn_PerlinNoise）。

实现的效果图6所示，滚动鼠标或者按键“A”, “W”, “S”, “D” , “U” , “J”或者按住鼠标左键并移动鼠标，可以控制视角。

2015-7-15 14-16-49

图6. 基于经典Perlin噪声的地形

SrTerrain.h

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
/************************************************************************      
\link   www.twinklingstar.cn
\author twinklingstar
\file   SrTerrain.h
\date   2015/05/20
****************************************************************************/
#ifndef SR_TERRAIN_H_
#define SR_TERRAIN_H_
 
#include <math.h>
#include <stdlib.h>
 
#include "SrClassicalPerlinNoise.h"
#include "SrModel.h"
 
 
//Terrain parameter.
struct SrTerrainParam
{
    float       mXSize;
    float       mYSize;
    float       mZSize;
    D3DXVECTOR3 mCenter;
    int         mNumRows;
    int         mNumCols;
    float       mF;
    float       mD;
    float       mErode;
 
};
 
/**
\brief Generate terrain based on Perlin noise.
 */
class SrTerrain:public SrModel
{
private:
    struct TERRAINVERTEX
    {
        D3DXVECTOR3   mPos;
        D3DXVECTOR3   mNormal;
        D3DXVECTOR2   mTexture;
    };
 
    #define D3DFVF_TERRAIN_VERTEX (D3DFVF_XYZ | D3DFVF_NORMAL| D3DFVF_TEX1)
 
public:
    SrTerrain(const TCHAR* file, const SrTerrainParam& para);
 
    virtual bool init(IDirect3DDevice9* pd3dDevice);
 
    virtual void render(IDirect3DDevice9* pd3dDevice) const;
 
    ~SrTerrain();
 
protected:
    void calHeightMap(float *h, D3DXVECTOR3* m);
    /**
    \brief Calculate the normal of each vertex.
     */
    void calNormal(D3DXVECTOR3* m, D3DXVECTOR3* norm);
    /**
    \brief Calculate the texture coordinates of each vertex.
     */
    void calTextureCoord(D3DXVECTOR2* coords);
 
    void calVertes(TERRAINVERTEX* v, D3DXVECTOR3 * m, D3DXVECTOR3* n, D3DXVECTOR2* coords);
 
    void buildMesh(IDirect3DDevice9* pd3dDevice);
    /**
    \brief Calculate the height.
     */
    void calHeight(float* height);
 
    void addPerlinNoise(float* height, float f, float base = 10.0f);
 
    void perturb(float* height, float f, float d);
 
    void erode(float * height, float smoothness);
 
    void smoothen(float * height);
private:
    SrClassicalPerlinNoise3D*   mPerlin;
    LPDIRECT3DVERTEXBUFFER9     mVB;
    LPDIRECT3DTEXTURE9          mTexture;
    SrTerrainParam              mParam;
    TCHAR                       mTextureFile[MAX_PATH];
 
};
 
#endif

SrTerrain.cpp

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
/************************************************************************      
\link   www.twinklingstar.cn
\author twinklingstar
\file   SrTerrain.cpp
\date   2015/05/20
****************************************************************************/
#include "SrTerrain.h"
 
 
 
SrTerrain::SrTerrain(const TCHAR* file,const SrTerrainParam& para):SrModel()
{
    mPerlin = SrClassicalPerlinNoise3D::create();
    mTexture         = NULL;
    mParam  = para;
    float x = mParam.mCenter.x - mParam.mXSize * 0.5f;
    float z = mParam.mCenter.z - mParam.mZSize * 0.5f;
    float y = mParam.mCenter.y - mParam.mYSize * 0.5f;
    D3DXMatrixTranslation(&mMatrix, x, y, z);
 
    lstrcpy(mTextureFile, file);
 
}
 
bool SrTerrain::init(IDirect3DDevice9* pd3dDevice)
{
    buildMesh(pd3dDevice);
    D3DXCreateTextureFromFile( pd3dDevice,mTextureFile, &mTexture);
 
    return true;
}
 
void SrTerrain::calHeightMap(float *h, D3DXVECTOR3* m)
{
    int r, c, indx = 0;
 
    float cstep = (float)mParam.mXSize / (float)mParam.mNumCols;
    float rstep = (float)mParam.mZSize / (float)mParam.mNumRows;
 
    for(r = 0; r < mParam.mNumRows ; r += 1)
    {
        for(c = 0 ; c < mParam.mNumCols ; c += 1)
        {
            m[indx].x = c * cstep;
            m[indx].y = h[indx];
            m[indx].z = r * rstep;
 
            indx += 1;
        }
    }
}
 
void SrTerrain::calNormal(D3DXVECTOR3* m, D3DXVECTOR3* norm)
{
    int r, c, count = 0;
    float len;
    D3DXVECTOR3 v1, v2, v3, s1, s2, n, sum;
    D3DXVECTOR3* normals = new D3DXVECTOR3[(mParam.mNumRows - 1) * (mParam.mNumCols - 1)];
 
    for(r = 0; r < mParam.mNumRows - 1 ; r += 1)
    {
        for(c = 0 ; c < mParam.mNumCols - 1 ; c += 1)
        {
            v1 = m[(r * mParam.mNumCols) + c];
            v2 = m[(r * mParam.mNumCols) + c + 1];
            v3 = m[((r + 1) * mParam.mNumCols) + c];
 
            s1 = v3 - v1;
            s2 = v2 - v1;
 
            n.x = (s1.y * s2.z) - (s1.z * s2.y);
            n.y = (s1.z * s2.x) - (s1.x * s2.z);
            n.z = (s1.x * s2.y) - (s1.y * s2.x);
 
            normals[r * (mParam.mNumCols - 1) + c] = n;
        }
    }
 
    for(r = 0; r < mParam.mNumRows ; r += 1)
    {
        for(c = 0 ; c < mParam.mNumCols ; c += 1)
        {
            sum = D3DXVECTOR3(0,0,0);
            count = 0;
            if(((r - 1) >= 0) && (c - 1) >= 0)
            {
                sum += normals[(r - 1) * (mParam.mNumCols - 1) + c - 1];
                count += 1;
            }
            if((r < (mParam.mNumRows - 1)) && ((c - 1) >= 0))
            {
                sum += normals[r * (mParam.mNumCols - 1) + c - 1];
                count += 1;
            }
            if(((r - 1) >= 0) && (c < (mParam.mNumCols - 1)))
            {
                sum += normals[(r - 1) * (mParam.mNumCols - 1) + c];
                count += 1;
            }
            if((r < (mParam.mNumRows - 1)) && (c < (mParam.mNumCols - 1)))
            {
                sum += normals[r * (mParam.mNumCols - 1) + c];
                count += 1;
            }
            sum /= (float) count;
            len = sqrt(sum.x * sum.x + sum.y * sum.y + sum.z * sum.z);
            norm[r * mParam.mNumCols + c]= (sum / len);
        }
    }
    delete []normals;
}
 
 
void SrTerrain::calTextureCoord(D3DXVECTOR2* coords)
{
    int r, c, indx = 0;
 
    for(r = 0; r < mParam.mNumRows ; r += 1)
    {
        for(c = 0 ; c < mParam.mNumCols ; c += 1)
        {
            coords[indx] = D3DXVECTOR2((float)c/(float)mParam.mNumCols * 10.0f, (1 - (float)r/(float)mParam.mNumRows) * 10.0f);
            indx += 1;
        }
    }
}
 
void SrTerrain::calVertes(TERRAINVERTEX* v, D3DXVECTOR3 * m, D3DXVECTOR3* n, D3DXVECTOR2* coords)
{
    int r, c, indx1, indx2, indx3, indx4;
    int indx = 0;
    for( r = 0; r < mParam.mNumRows - 1; r += 1)
    {
        for(c = 0 ; c < mParam.mNumCols - 1 ; c += 1)
        {
            indx1 = r * mParam.mNumCols + c;
            indx2 = (r + 1) * mParam.mNumCols + c;
            indx3 = r * mParam.mNumCols + (c + 1);
            indx4 = (r + 1) * mParam.mNumCols + (c + 1);
 
            //Triangle 1.
            v[indx].mPos     = m[indx1];
            v[indx].mNormal  = n[indx1];
            v[indx].mTexture = coords[indx1];
            indx += 1;
 
            v[indx].mPos     = m[indx2];
            v[indx].mNormal  = n[indx2];
            v[indx].mTexture = coords[indx2];
            indx += 1;
 
            v[indx].mPos     = m[indx4];
            v[indx].mNormal  = n[indx4];
            v[indx].mTexture = coords[indx4];
            indx += 1;
            //Triangle 2.
            v[indx].mPos     = m[indx1];
            v[indx].mNormal  = n[indx1];
            v[indx].mTexture = coords[indx1];
            indx += 1;
 
            v[indx].mPos     = m[indx4];
            v[indx].mNormal  = n[indx4];
            v[indx].mTexture = coords[indx4];
            indx += 1;
 
            v[indx].mPos     = m[indx3];
            v[indx].mNormal  = n[indx3];
            v[indx].mTexture = coords[indx3];
            indx += 1;
        }
    }
}
 
void SrTerrain::buildMesh(IDirect3DDevice9* pd3dDevice) 
{
    int count = mParam.mNumCols * mParam.mNumRows;
    float *h            = new float[count];
    D3DXVECTOR3* n      = new D3DXVECTOR3[count];
    D3DXVECTOR3* m      = new D3DXVECTOR3[count];
    D3DXVECTOR2* coords = new D3DXVECTOR2[count];
    memset(h, 0, sizeof(float) * count);
 
    calHeight(h);
 
    calHeightMap(h, m);
 
    calNormal(m, n);
 
    calTextureCoord(coords);
 
    TERRAINVERTEX* v;
    int numVerts = (mParam.mNumRows - 1) * (mParam.mNumCols - 1);
    pd3dDevice->CreateVertexBuffer(6 * numVerts * sizeof(TERRAINVERTEX), 0, D3DFVF_TERRAIN_VERTEX, D3DPOOL_MANAGED, &mVB, NULL );
    mVB->Lock( 0, 6 * numVerts * sizeof(TERRAINVERTEX), (void**)&v, 0 );
    calVertes(v, m, n, coords);
    mVB->Unlock();
 
    delete []h;
    delete []n;
    delete []m;
    delete []coords;
}
void SrTerrain::render(IDirect3DDevice9* pd3dDevice)const
{
    pd3dDevice->SetTransform( D3DTS_WORLD, &mMatrix );
    pd3dDevice->SetTexture(0, mTexture);
    pd3dDevice->SetStreamSource( 0, mVB, 0, sizeof(TERRAINVERTEX) );
    pd3dDevice->SetFVF( D3DFVF_TERRAIN_VERTEX );
    pd3dDevice->DrawPrimitive( D3DPT_TRIANGLELIST, 0, 2 * (mParam.mNumCols - 1) * (mParam.mNumRows - 1));
 
}
 
SrTerrain::~SrTerrain()
{
    SAFE_RELEASE(mTexture);
    SAFE_RELEASE(mVB);
    delete mPerlin;
}
 
void SrTerrain::calHeight(float* height)
{
    memset(height, 0, sizeof(float) * mParam.mNumCols * mParam.mNumRows);
 
    addPerlinNoise(height, mParam.mYSize);
    perturb(height, mParam.mF, mParam.mD);
    for (int i = 0; i < 10; i++ )
        erode(height, mParam.mErode);
    smoothen(height);
}
 
void SrTerrain::addPerlinNoise(float* height, float f, float base)
{
    int r, c;
    for (r = 0; r < mParam.mNumRows; r += 1)
    {
        for (c = 0; c < mParam.mNumCols; c += 1)
        {
            height[r * mParam.mNumCols + c] += base * mPerlin->noise(f * r / (float)mParam.mNumRows, f * c / (float)mParam.mNumCols, 0);
        }
    }
}
void SrTerrain::perturb(float* height, float f, float d)
{
    int u, v, i, j;
    float* tmp = new float[mParam.mNumCols * mParam.mNumRows];
    for (i = 0; i < mParam.mNumRows; i += 1)
    {
        for (j = 0; j < mParam.mNumCols; j += 1)
        {
            u = i + (int)(mPerlin->noise(f * i / (float)mParam.mNumRows, f * j / (float)mParam.mNumCols, 0) * d);
            v = j + (int)(mPerlin->noise(f * i / (float)mParam.mNumRows, f * j / (float)mParam.mNumCols, 1) * d);
            if (u < 0) u = 0; if (u >= mParam.mNumRows) u = mParam.mNumRows - 1;
            if (v < 0) v = 0; if (v >= mParam.mNumCols) v = mParam.mNumCols - 1;
            tmp[i * mParam.mNumCols + j] = height[u * mParam.mNumCols + v];
        }
    }
    memcpy(height, tmp, sizeof(float) * mParam.mNumCols * mParam.mNumRows);
 
    delete []tmp;
}
void SrTerrain::erode(float * height, float smoothness)
{
    int i, j;
    for (i = 1; i < mParam.mNumRows - 1; i++)
    {
        for (j = 1; j < mParam.mNumCols - 1; j++)
        {
            float d_max = 0.0f;
            int match[] = {0, 0};
 
            for (int u = -1; u <= 1; u++)
            {
                for (int v = -1; v <= 1; v++)
                {
                    if(abs(u) +abs(v) > 0)
                    {
                        float d_i = height[i * mParam.mNumCols + j] - height[(i + u) * mParam.mNumCols + (j + v)];
                        if (d_i > d_max)
                        {
                            d_max = d_i;
                            match[0] = u; match[1] = v;
                        }
                    }
                }
            }
            if(0 < d_max && d_max <= smoothness)
            {
                float d_h = 0.5f * d_max;
                height[i * mParam.mNumCols + j] -= d_h;
                height[(i + match[0]) * mParam.mNumCols + (j + match[1])] += d_h;
            }
        }
    }
}
 
void SrTerrain::smoothen(float * height)
{
    for (int i = 1; i < mParam.mNumRows - 1; ++i)
    {
        for (int j = 1; j < mParam.mNumCols - 1; ++j)
        {
            float total = 0.0f;
            for (int u = -1; u <= 1; u++)
            {
                for (int v = -1; v <= 1; v++)
                {
                    total += height[(i + u) * mParam.mNumCols + (j + v)];
                }
            }
            height[i * mParam.mNumCols + j] = total / 9.0f;
        }
    }
}

参考

  1. Ken Perlin. "An image synthesizer." ACM Siggraph Computer Graphics, vol.19, no.3, pp.287-296, 1985.
  2. Ken Perlin. "Improving noise." ACM Transactions on Graphics (TOG), vol.21, no.3, 2002.
  3. Hugo Elias. “Perlin Noise”. website <http://freespace.virgin.net/hugo.elias/models/m_perlin.htm>.
  4. “Perlin 噪声.” website <https://zh.wikipedia.org/wiki/Perlin%E5%99%AA%E5%A3%B0>.
  5. Ken Perlin. “Noise and Turbulence.” website <http://www.mrl.nyu.edu/~perlin/doc/oscar.html>.
  6. Matt Zucker. “The Perlin noise math FAQ.” website <http://webstaff.itn.liu.se/~stegu/TNM022-2005/perlinnoiselinks/perlin-noise-math-faq.html>.
  7. Ken Perlin. “Making Noise.” website <http://www.noisemachine.com/talk1/>.
  8. Stefan Gustavson. "Simplex noise demystified." Linköping University, Linköping, Sweden, Research Report, 2005.
  9. “Generating realistic and playable terrain height”. website <http://www.float4x4.net/index.php/2010/06/generating-realistic-and-playable-terrain-height-maps/>.
Posted in: Computer Graphics 
spacer
More Articles
← 实时阴影技术
凹凸纹理 →

One comment on “Perlin噪声”

    1.  
      caiyingCommented on: April 19, 2018

      感谢博主的演算过程！想问一个问题：为什么在二维加权的时候，插值计算的结果是 a=s+(t−s)⋅sx; b=u+(v−u)⋅sx呢？为什么是t-s, v-u，而不是s-u, t-v之类呢？是有方向或者其他的限制么？

posted on 2018-04-23 12:52  fanbird2008  阅读(445)  评论(0)    收藏  举报

导航