本篇文章谢绝转载，也禁止用于任何商业目的
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 \left( {{2}^{N}} \right)\)，但是Simplex噪声的算法复杂度为\(\log \left( {{N}^{2}} \right)\)等优点。虽然Stefan Gustavson[8]提供了对Simplex算法的注解，但是我依然不能理解Simplex噪声背后的数学原理，对Simplex噪声不作进一步阐述。
2. 经典Perlin噪声
经典Perlin噪声是Ken Perlin在1983年提出的噪声，Ken Perlin提供了一维、二维、三维算法的C实现[5]，我们无法仅从他提供的代码理解其数学原理，Matt Zucker[6]从数学原理的角度，对经典Perlin噪声进行了解读。本节将详细介绍二维、三维的Pernlin噪声的数学原理，算法的C语言实现源码如下所示：
/* 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];
}
}
用简单的伪码如下所示（不理解，可以直接忽略）：
- 给定点输入点\(p\)；
- 对于每一个与点\(p\)相邻的方格端点：
- 挑选一个伪随机梯度向量；
- 计算点积；
- 每个维度上，均采用缓和曲线，计算出加权平均值，缓和曲线可以选用\(3{{t}^{2}}-2{{t}^{3}}\)。
2.1. 二维空间
1. 找到方格点
先考虑二维的情况，噪声函数用式子\(Noise\left( x,y \right)\)表示，其中\(x,y\)是浮点数，\(\left( x,y \right)\)表示一个点，设对变量\(x\)下取整用式子\(floor(x)\)表示，即\(floor\left( 4.34 \right)=4\)，令：
\({{x}_{0}}=floor\left( x \right)\)
\({{y}_{0}}=floor\left( y \right)\)
\({{x}_{1}}={{x}_{0}}+1\)
\({{y}_{1}}={{y}_{0}}+1\)
用图表示点\(\left( x,y \right)\)和\(\left( {{x}_{0}},{{y}_{0}} \right),\left( {{x}_{1}},{{y}_{0}} \right),\left( {{x}_{0}},{{y}_{1}} \right),\left( {{x}_{1}},{{y}_{1}} \right)\)，如图1.所示。
图1. 包围点\(\left( x,y \right)\)的方格
2. 挑选伪随机梯度向量
对于方格的每个端点，分配一个伪随机梯度（Pseudorandom Gradient）向量，这是一个在二维空间上单位长度的向量。所谓的伪随机，是指相同的输入值将产生相同的输出值，并不是真正意义上的随机。可以通过下面的方法挑选伪随机梯度向量[7]，\(n\)规定表的长度：
- 生成一个随机排列表，\(P[n]\)；
- 生成一个随机梯度表，\(G[n]\)；
- \(G\left( i,j \right)=G\left[ \left( i+P\left[ j \right] \right)\bmod n \right]\)；
在实际应用中，我们不需要求模运算，设\(n=256\)，通过位运算，留下低8位的值就行，即\(G\left( i,j \right)=G\left[ \left( i+P\left[ j \right] \right)\And 0xff \right]\)。由方格的四个端点，我们可以得出对应的伪随机梯度向量，如图2所示。
图2. 伪随机梯度向量
3. 计算点积
接着，生成一个由点\(\left( x,y \right)\)到方格四个端点的向量，如图3所示，四个向量分别为：
\({{\vec{d}}_{00}}=\left( x,y \right)-\left( {{x}_{0}},{{y}_{0}} \right)\)
\({{\vec{d}}_{10}}=\left( x,y \right)-\left( {{x}_{1}},{{y}_{0}} \right)\)
\({{\vec{d}}_{01}}=\left( x,y \right)-\left( {{x}_{0}},{{y}_{1}} \right)\)
\({{\vec{d}}_{11}}=\left( x,y \right)-\left( {{x}_{1}},{{y}_{1}} \right)\)
这四个向量分别与对应端点的伪随机梯度向量点积，得到：
\(s=G\left( {{x}_{0}},{{y}_{0}} \right)\cdot \left( \left( x,y \right)-\left( {{x}_{0}},{{y}_{0}} \right) \right)\)
\(t=G\left( {{x}_{1}},{{y}_{0}} \right)\cdot \left( \left( x,y \right)-\left( {{x}_{1}},{{y}_{0}} \right) \right)\)
\(u=G\left( {{x}_{0}},{{y}_{1}} \right)\cdot \left( \left( x,y \right)-\left( {{x}_{0}},{{y}_{1}} \right) \right)\)
\(v=G\left( {{x}_{1}},{{y}_{1}} \right)\cdot \left( \left( x,y \right)-\left( {{x}_{1}},{{y}_{1}} \right) \right)\)
它的几何意义，就是将向量\({{\vec{d}}_{ij}}\)（\(i,j\in \left\{ 1,2 \right\}\)）投影到对应的伪随机梯度向量，显然有\(s,t,u,v\in \left[ 0,1 \right]\)。不理解其几何意义的作用，但是可以明确的一点是：\(s,t,u,v\)值由点\(\left( x,y \right)\)和对应的伪随机梯度向量唯一确定。
图3. 由点\(\left( x,y \right)\)到方格四个端点的向量
4. 计算加权平均值
接下来，考虑如何计算这几个值的加权平均，得到最终的噪声标量。
这里引入缓和曲线（Ease Curve）的概念，缓和曲线使噪声显得更加真实，采用的缓和曲线函数为：\(f\left( x \right)=3{{x}^{2}}-2{{x}^{3}}\)，如图4所示。
图4. \(f\left( x \right)=3{{x}^{2}}-2{{x}^{3}}\)曲线
把\(x-{{x}_{0}}\)代入缓和函数\(f\left( x \right)\)中，得到比例系数\({{s}_{x}}=f\left( x-{{x}_{0}} \right)\)，通过线性插值得到：
\(a=s+(t-s)\cdot {{s}_{x}}\)
\(b=u+(v-u)\cdot {{s}_{x}}\)
再把\(y-{{y}_{0}}\)代入缓和函数\(f\left( x \right)\)中，得到比例系数\({{s}_{y}}=f\left( y-{{y}_{0}} \right)\)，通过线性插值得到：
\(Noise\left( x,y \right)=a+(b-a)\cdot {{s}_{y}}\)
2.2. 三维空间
在三维空间上，噪声函数用式子\(Noise\left( x,y,z \right)\)表示，与二维噪声相似，同样可以分成四个阶段来处理。三维空间上的方格是一个正方体，这里特别说明下点积的计算和加权平均值的计算。
分别把\(x-{{x}_{0}},y-{{y}_{0}},z-{{z}_{0}}\)代入缓和函数\(f\left( x \right)\)中，得到比例系数：
\({{s}_{x}}=f\left( x-{{x}_{0}} \right),{{s}_{y}}=f\left( y-{{y}_{0}} \right),{{s}_{z}}=f\left( z-{{z}_{0}} \right)\)
图5. 三维空间上的方格
如图5所示，计算点积
\({{s}_{0}}=G\left( {{x}_{0}},{{y}_{0}},{{z}_{0}} \right)\cdot \left( \left( x,y,z \right)-\left( {{x}_{0}},{{y}_{0}},{{z}_{0}} \right) \right)\)
\({{t}_{0}}=G\left( {{x}_{1}},{{y}_{0}},{{z}_{0}} \right)\cdot \left( \left( x,y,z \right)-\left( {{x}_{1}},{{y}_{0}},{{z}_{0}} \right) \right)\)
\({{u}_{0}}=G\left( {{x}_{0}},{{y}_{1}},{{z}_{0}} \right)\cdot \left( \left( x,y,z \right)-\left( {{x}_{0}},{{y}_{1}},{{z}_{0}} \right) \right)\)
\({{v}_{0}}=G\left( {{x}_{1}},{{y}_{1}},{{z}_{0}} \right)\cdot \left( \left( x,y,z \right)-\left( {{x}_{1}},{{y}_{1}},{{z}_{0}} \right) \right)\)
\({{s}_{1}}=G\left( {{x}_{0}},{{y}_{0}},{{z}_{1}} \right)\cdot \left( \left( x,y,z \right)-\left( {{x}_{0}},{{y}_{0}},{{z}_{1}} \right) \right)\)
\({{t}_{1}}=G\left( {{x}_{1}},{{y}_{0}},{{z}_{1}} \right)\cdot \left( \left( x,y,z \right)-\left( {{x}_{1}},{{y}_{0}},{{z}_{1}} \right) \right)\)
\({{u}_{1}}=G\left( {{x}_{0}},{{y}_{1}},{{z}_{1}} \right)\cdot \left( \left( x,y,z \right)-\left( {{x}_{0}},{{y}_{1}},{{z}_{1}} \right) \right)\)
\({{v}_{1}}=G\left( {{x}_{1}},{{y}_{1}},{{z}_{1}} \right)\cdot \left( \left( x,y,z \right)-\left( {{x}_{1}},{{y}_{1}},{{z}_{1}} \right) \right)\)
通过线性插值，计算出红色实心圆的数值：
\({{a}_{0}}={{s}_{0}}+\left( {{t}_{0}}-{{s}_{0}} \right)\cdot {{s}_{x}}\)
\({{b}_{0}}={{u}_{0}}+\left( {{v}_{0}}-{{u}_{0}} \right)\cdot {{s}_{x}}\)
\({{a}_{1}}={{s}_{1}}+\left( {{t}_{1}}-{{s}_{1}} \right)\cdot {{s}_{x}}\)
\({{b}_{1}}={{u}_{1}}+\left( {{v}_{1}}-{{u}_{1}} \right)\cdot {{s}_{x}}\)
通过线性插值，计算出紫色实心圆的数值：
\({{z}_{0}}={{a}_{0}}+\left( {{b}_{0}}-{{a}_{0}} \right)\cdot {{s}_{y}}\)
\({{z}_{1}}={{a}_{1}}+\left( {{b}_{1}}-{{a}_{1}} \right)\cdot {{s}_{y}}\)
通过线性插值，计算出绿色实心圆的数值，即是最终的结果：
\(Noise(x,y,z)={{z}_{0}}+\left( {{z}_{1}}-{{z}_{0}} \right)\cdot {{s}_{z}}\)
2.3. 实现代码
最后，提供一个经典Perlin噪声的C++实现代码，如下所示：
SrClassicalPerlinNoise.h
/************************************************************************
\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
/************************************************************************
\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”或者按住鼠标左键并移动鼠标，可以控制视角。
图6. 基于经典Perlin噪声的地形
SrTerrain.h
/************************************************************************
\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
/************************************************************************
\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;
}
}
}
参考
- Ken Perlin. "An image synthesizer." ACM Siggraph Computer Graphics, vol.19, no.3, pp.287-296, 1985.
- Ken Perlin. "Improving noise." ACM Transactions on Graphics (TOG), vol.21, no.3, 2002.
- Hugo Elias. “Perlin Noise”. website <http://freespace.virgin.net/hugo.elias/models/m_perlin.htm>.
- “Perlin 噪声.” website <https://zh.wikipedia.org/wiki/Perlin%E5%99%AA%E5%A3%B0>.
- Ken Perlin. “Noise and Turbulence.” website <http://www.mrl.nyu.edu/~perlin/doc/oscar.html>.
- Matt Zucker. “The Perlin noise math FAQ.” website <http://webstaff.itn.liu.se/~stegu/TNM022-2005/perlinnoiselinks/perlin-noise-math-faq.html>.
- Ken Perlin. “Making Noise.” website <http://www.noisemachine.com/talk1/>.
- Stefan Gustavson. "Simplex noise demystified." Linköping University, Linköping, Sweden, Research Report, 2005.
- “Generating realistic and playable terrain height”. website <http://www.float4x4.net/index.php/2010/06/generating-realistic-and-playable-terrain-height-maps/>.
感谢博主的演算过程！想问一个问题：为什么在二维加权的时候，插值计算的结果是 a=s+(t−s)⋅sx; b=u+(v−u)⋅sx呢？为什么是t-s, v-u，而不是s-u, t-v之类呢？是有方向或者其他的限制么？