实时渲染 perlin noise

前言

本文将对图形学噪声进行介绍，其中重点介绍perlin noise，并讲解如何实现Perlin noise和优化Perlin noise

噪声

定义

噪声是一系列的随机变量，而噪声函数是给定一个输入，输出一个随机变量
噪声是将自然现象的变化加入纹理的一个工具，因为在真实世界中这些自然现象并不是完全一样的，所以噪声提供一个可控方式将随机性添加到shader中

为什么

学过计算机的都知道random()用于求得伪随机数，那么是不是噪声函数就是运用random()?但在图形学中并不是这样的，这是因为普通噪声毫无规律可言(白噪声)，而想要模拟现实中的自然噪声需要在随机中又有规律
以下是普通噪声:
噪声函数的目的：在三维空间中，提供一种高效地实现、可重复、伪随机的信号，此信号大部分能量集中在一个空间频率附近，而视觉上各向同性的(统计上不随旋转变化)

以下是柏林噪声：

分类

噪声依据自然现象进行分类，可分为两类
- 基于晶格(Lattice)
  - 梯度噪声(Gradient noise)Gradient noise - Wikipedia
    - Perlin噪声Perlin noise - Wikipedia
      
      应用范围：云朵、火焰、地形等自然环境
    - Simplex噪声Simplex noise - Wikipedia
      
      基于Perlin进行优化
  - Value noiseValue noise - Wikipedia
- 基于点
  - Worley噪声Worley noise - Wikipedia
    
    应用范围：多孔结构，如纸张木纹

Perlin noise

特点

所有视角细节都是相同的大小
可复现性。对于每个输入位置，提供一个可重复的伪随机值
实现简单快捷
有一个已知的范围(通常是[-1,1])
空间频率带宽有限
不明显的重复性
空间频率在平移下是不变的

实现

三个步骤：

初始化

在3维空间的每个整数(x,y,z)位置产生一个可重复的伪随机值，这运用哈希函数实现。哈希函数的参数是点坐标，利用点坐标得到的梯度表的索引，最后得到梯度表对应的值

排列表（Permutation Table）：乱序存放一系列索引(梯度表的索引)

梯度表（Gradient Table）：存放一系列随机梯度值
建立采样空间,对于噪声图上的像素，找到它在单位矩形/单位立方体上对应的一点，并求出它的参考点的坐标
1. 对于一维柏林噪声，采样空间为一个一维的坐标轴，轴上整数坐标位置均有一个点，参考点为两侧最近的整数点
2. 对于二维柏林噪声，采样空间为一个二维坐标系，坐标系中横纵坐标为整数的地方均有一点,参考点为组成包围该点的单位正方体的四个点
3. 三维同理,参考点为组成包围该点的单位立方体的八个点
对于不同类型的噪声，采样点在不同空间中，根据伪随机整数的梯度来计算采样空间中小数位置到参考点的距离向量和梯度向量的点积,最后进行插值

一维噪声需一次插值，二维需三次插值，三维需七次插值，呈指数增长，时间复杂度是 $O (2^{N})$

基于值的柏林噪声

初始化

决定一个点的梯度，我们需要结合哈希函数，以点坐标作为参数，以所得值作为索引，在排列表中取得对应值。也就是说，将点坐标转换为排列表的索引，最后取得索引对应的梯度表的值

//以下是二维的哈希函数
int[] perm = {...};
int[] grad = {...};

void hash( int[] gradient, int x, int y)
{
    //找到对应值在排列表中的索引
    int permIdx[] = new int[2];
    permIdx[0] = FindIndex(x);
    permIdx[1] = FindIndex(y);
    
    //通过排列表的索引值查找梯度表的索引
    int gradIdx[] = new int[2];
    gradIdx[0] = perm[permIdx[0]];
    gradIdx[0] = perm[permIdx[1]];
    
    //通过梯度表的索引找到对应的梯度值
    gradient[0] = grad[gradIdx[0]];
    gradient[1] = grad[gradIdx[1]];
}

在一维中，各个整点的梯度可以用斜率来表现,因此我们可以不使用梯度表，直接将排列表对应索引取得的值作为梯度

在经典噪声算法中，使用的排列表是由数字0到255散列组成的

int 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 };

采样、插值、计算梯度向量和其偏移量之间的点积

采样：一维柏林噪声函数的参数输入点是坐标轴上的一点，其参考点就是它左右两侧最近的整数点。因此，我们需要找到这两个点，并求出这两个点对应的梯度及相对输入点的方向向量

插值：对梯度值进行插值，公式为 $x = x 1 + t (x 2 - x 1)$ ，不如线性插值效果不太好，产生的噪声曲线很尖锐。我们使用 $t = 3 t^{2} - 2 t^{3}$ 来替换线性插值公式里的 $t$

float perlin( float x )
{
    //假设x1为x左边的索引值且x - x1的方向向量为正，x2为右边的索引值
    int x1 = floor(x);
    int x2 = x1 + 1;
    //方向向量
    float vec1 = x - x1;
    float vec2 = x - x2;
    
    //求梯度值
    float grad1 = perm[x1 % 255] * 2.0 - 255.0;
    float grad2 = perm[x2 % 255] * 2.0 - 255.0;
    
    //插值
    float t = 3 * pow( vec1, 2 ) - 2 * pow(vec1, 3);
    float product1 = grad1 * vec1;
    float product2 = grad2 * vec2;
    
    return product1 + t * (product2 - product1);
}

柏林梯度噪声

思路

为方便讲解此处假定为2d平面。

取输入值(x, y)(下图中的蓝点)的小数部分(如此将问题抽象到一个单位立方体内)
为输入点周围四个点生成一个伪随机的梯度向量(在后面有讲)
计算输入点到周围四个点的距离向量
对每个顶点进行梯度向量和距离向量的点积运算，求得该点的影响值
对四个顶点的影响值进行lerp

初始化

int repeat;

Perlin(int repeat = -1) 
{
	this->repeat = repeat;
}

double perlin(double x, double y, double z) 
{
    if(repeat > 0) 
    {                   // If we have any repeat on, change the coordinates to their "local" repetitions
        x = x % repeat;
        y = y % repeat;
        z = z % repeat;
    }
   	
    // clamp[0,255]
    int xi = (int)x & 255;            
    int yi = (int)y & 255;             
    int zi = (int)z & 255; 
    // 小数部分
    double xf = x-(int)x;
    double yf = y-(int)y;
    double zf = z-(int)z;
    // ...
}

排列表

作用

用于哈希函数，确定梯度向量

实现

为避免缓存溢出，需要重复填充一次排列表

int[] permutation = { 151,160,137,91,90,15,                 // Hash lookup table as defined by Ken Perlin
    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
};

int []p;
Perlin()
{
    p = new int[512];
    for(int x=0;x<512;x++) 
    {
        p[x] = permutation[x%256];
    }
}

哈希函数

目的

为每个输入计算一个确定值

实现

// 将输入值增加1，同时保证范围[0,repeat)
int inc(int num) 
{
    num++;
    if (repeat > 0) num %= repeat;
    
    return num;
}

// 对围着输入坐标的8个索引坐标(单元正方形)进行哈希运算
{
    int aaa, aba, aab, abb, baa, bba, bab, bbb;
    aaa = p[p[p[    xi ]+    yi ]+    zi ];
    aba = p[p[p[    xi ]+inc(yi)]+    zi ];
    aab = p[p[p[    xi ]+    yi ]+inc(zi)];
    abb = p[p[p[    xi ]+inc(yi)]+inc(zi)];
    baa = p[p[p[inc(xi)]+    yi ]+    zi ];
    bba = p[p[p[inc(xi)]+inc(yi)]+    zi ];
    bab = p[p[p[inc(xi)]+    yi ]+inc(zi)];
    bbb = p[p[p[inc(xi)]+inc(yi)]+inc(zi)];
}

梯度函数(Grad)

目的

基于值的perlin noise的噪声图会有明显的块状，而基于梯度的perlin noise正是用于解决该问题。这里不再是在每个点随机生成值，而是使用向量。
梯度向量
下图中，正方形四个点上的向量即为梯度向量(大致是这个并不完全是，这个是优化前采用的随机梯度向量)。优化后的柏林噪声不再用随机的梯度向量，而是从由单位正方体的中心点指向各条边中点的12个向量中随机选择：(1,1,0),(-1,1,0),(1,-1,0),(-1,-1,0), (1,0,1),(-1,0,1),(1,0,-1),(-1,0,-1), (0,1,1),(0,-1,1),(0,1,-1),(0,-1,-1)

实现
该实现的含义：从由单位正方体的中心点指向各条边中点的12个向量中随机挑选一个梯度向量

// hash是由哈希函数计算得到的值
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
    }
}

Fade函数(lerp因子)

优化后的版本lerp因子t不再是 $3 t^{2} - 2 t^{3}$ ，而是 $6 t^{5} - 15 t^{4} + 10 t^{3}$
3t^2 - 2t^3的缺点

二次导数含有非0值，在使用噪声函数的导数时会出现失真效果

实现

double fade(double t) 
{                                                      
    return t * t * t * (t * (t * 6 - 15) + 10);         // 6t^5 - 15t^4 + 10t^3
}

三线性插值

实现

对求得的8个顶点值进行三线性插值

double lerp(double a, double b, double x) 
{
    return a + x * (b - a);
};

{
    double x1, x2, y1, y2;
    x1 = lerp(grad(aaa, xf, yf, zf), grad(baa, xf-1, yf, zf), u);                     
    x2 = lerp(grad (aba, xf, yf-1, zf), grad (bba, xf-1, yf-1, zf), u);
    y1 = lerp(x1, x2, v);

    x1 = lerp(grad(aab, xf, yf, zf-1), grad(bab, xf-1, yf, zf-1), u);
    x2 = lerp(grad (abb, xf, yf-1, zf-1), grad (bbb, xf-1, yf-1, zf-1), u);
    y2 = lerp (x1, x2, v);
    
    //// For convenience we bind the result to 0 - 1 (theoretical min/max before is [-1, 1])
    return (lerp (y1, y2, w)+1)/2;  
}