Perlin噪声

1. Perlin噪声
2. 经典Perlin噪声
3. 地形的应用
参考

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

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$ ；
对于每一个与点相邻的方格端点：
1. 挑选一个伪随机梯度向量；
2. 计算点积；
每个维度上，均采用缓和曲线，计算出加权平均值，缓和曲线可以选用 $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”或者按住鼠标左键并移动鼠标，可以控制视角。

2015-7-15 14-16-49

图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/>.

3 comments on “Perlin噪声”

Pingback: 不只是噪声，更是数学美 ---浅谈Perlin Noise - 算法网
caiying Commented on: April 19, 2018

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

Reply to this Comment
1. Anonymous Commented on: December 25, 2020
  
  个人以为换成s-u那么公式要改成a = u + (t -u)*sx，b也要做同样处理。这里的操作是先在平行的两条边上映射出两个点，然后在将这两个点映射出（x,y）点的向量，那么平行的两条边是垂直方向或者水平方向都是可以的
  
  Reply to this Comment

TwinklingStar

1. Perlin噪声

2. 经典Perlin噪声

2.1. 二维空间

2.2. 三维空间

2.3. 实现代码

3. 地形的应用

参考

3 comments on “Perlin噪声”

Leave a reply Cancel reply

TwinklingStar

Categories

Recent Post

Archive