一篇文章搞懂柏林噪聲算法，附代碼講解

本文目的是以一種通俗簡單的方式介紹Ken Perlin的改進版柏林噪聲算法，講解代碼採用c#編寫，開源無償使用。若是你只是想看完整代碼，點擊這裏。

柏林噪聲是一個很是強大算法，常常用於程序生成隨機內容，在遊戲和其餘像電影等多媒體領域普遍應用。算法發明者Ken Perlin也所以算法得到奧斯卡科技成果獎（靠算法拿奧斯卡也是沒誰了666）。本文將剖析他於2002年發表的改進版柏林噪聲算法。在遊戲開發領域，柏林噪聲能夠用於生成波形，起伏不平的材質或者紋理。例如，它能用於程序生成地形（例如使用柏林噪聲來生成個人世界（Minecraft）裏的地形），火焰燃燒特效，水和雲等等。柏林噪聲絕大部分應用在2維，3維層面上，但某種意義上也能拓展到4維。柏林噪聲在1維層面上可用於卷軸地形、模擬手繪線條等。
若是將柏林噪聲拓展到4維層面，以第4維，即w軸表明時間，就能利用柏林噪聲作動畫。例如，2D柏林噪聲能夠經過插值生成地形，而3D柏林噪聲則能夠模擬海平面上起伏的波浪。下面是柏林噪聲在不一樣維度的圖像以及在遊戲中的應用場景。

正如圖所示，柏林噪聲算法能夠用來模擬許多天然中的噪聲現象。接下來讓咱們從數理上分析算法的實現原理。

---

基本原理

注意：事先聲明，本節內容大多源於this wonderful article by Matt Zucker，不過該篇文章內容也是創建在1980年所發明的柏林噪聲算法基礎上的。本文我將使用2002年發明的改進版柏林噪聲算法。所以，個人算法版本跟Zucker的版本會有些不一樣。

讓咱們從最基本的柏林噪聲函數看起：
public double perlin(double x, double y, double z);

函數接收x,y,z三個座標份量做爲輸入，並返回0.0~1.0的double值做爲輸出。那咱們應該怎麼處理輸入值？首先，咱們取3個輸入值x,y,z的小數點部分，就能夠表示爲單元空間裏的一個點了。爲了方便講解，咱們將問題降維到2維空間來討論（原理是同樣的），下圖是該點在2維空間上的表示：

圖1：小藍點表明輸入值在單元正方形裏的空間座標，其餘4個點則是單元正方形的各頂點

接着，咱們給4個頂點（在3維空間則是8個頂點）各自生成一個僞隨機的梯度向量。梯度向量表明該頂點相對單元正方形內某點的影響是正向仍是反向的（向量指向方向爲正向，相反方向爲反向）。而僞隨機是指，對於任意組相同的輸入，一定獲得相同的輸出。所以，雖然每一個頂點生成的梯度向量看似隨機，實際上並非。這也保證了在梯度向量在生成函數不變的狀況下，每一個座標的梯度向量都是肯定不變的。

舉個例子來理解僞隨機，好比咱們從圓周率π（3.14159...）的小數部分中隨機抽取某一位數字，結果看似隨機，但若是抽取小數點後1位，結果一定爲1；抽取小數點後2位，結果一定爲4。

圖2：各頂點上的梯度向量隨機選取結果

請注意，上圖所示的梯度向量並非徹底準確的。在本文所介紹的改進版柏林噪聲中，這些梯度向量並非徹底隨機的，而是由12條單位正方體（3維）的中心點到各條邊中點的向量組成：
(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)

採用這些特殊梯度向量的緣由在Ken Perlin's SIGGRAPH 2002 paper: Improving Noise這篇文章裏有具體講解。

注意：許多介紹柏林噪聲算法的文章都是根據最第一版柏林噪聲算法來說解的，預約義的梯度表不是本文所說的這12個向量。如圖2所示的梯度向量就是最第一版算法所隨機出來的梯度向量，不過這兩種算法的原理都是同樣的。

接着，咱們須要求出另外4個向量（在3維空間則是8個），它們分別從各頂點指向輸入點（藍色點）。下面有個2維空間下的例子：

圖3:各個距離向量

接着，對每一個頂點的梯度向量和距離向量作點積運算，咱們就能夠得出每一個頂點的影響值：
grad.x * dist.x + grad.y * dist.y + grad.z * dist.z

這正是算法所須要的值，點積運算爲兩向量長度之積，再乘以兩向量夾角餘弦：

dot(vec1,vec2) = cos(angle(vec1,vec2)) * vec1.length * vec2.length

換句話說，若是兩向量指向同一方向，點積結果爲：

1 * vec1.length * vec2.length

若是兩向量指向相反方向，則點積結果爲：

-1 * vec1.length * vec2.length

若是兩向量互相垂直，則點積結果爲0。

0 * vec1.length * vec2.length

點積也能夠理解爲向量a在向量b上的投影，當距離向量在梯度向量上的投影爲同方向，點積結果爲正數；當距離向量在梯度向量上的投影爲反方向，點積結果爲負數。所以，經過兩向量點積，咱們就知道該頂點的影響值是正仍是負的。不難看出，頂點的梯度向量直接決定了這一點。下面經過一副彩色圖，直觀地看下各頂點的影響值：

圖4：2D柏林噪聲的影響值

下一步，咱們須要對4個頂點的影響值作插值，求得加權平均值（在3維空間則是8個）。算法很是簡單（2維空間下的解法）：

// Below are 4 influence values in the arrangement:
// [g1] | [g2]
// -----------
// [g3] | [g4]
int g1, g2, g3, g4;
int u, v;   // These coordinates are the location of the input coordinate in its unit square.  
            // For example a value of (0.5,0.5) is in the exact center of its unit square.

int x1 = lerp(g1,g2,u);
int x2 = lerp(g3,h4,u);

int average = lerp(x1,x2,v);

至此，整個柏林噪聲算法還剩下最後一塊拼圖了：若是直接使用上述代碼，因爲是採用lerp線性插值計算得出的值，雖然運行效率高，但噪聲效果很差，看起來會不天然。咱們須要採用一種更爲平滑，非線性的插值函數：fade函數，一般也被稱爲ease curve(也做爲緩動函數在遊戲中普遍使用)：

圖5：ease curve

ease curve的值會用來計算前面代碼裏的u和v，這樣插值變化再也不是單調的線性變化，而是這樣一個過程：初始變化慢，中間變化快，結尾變化又慢下來（也就是在當數值趨近於整數時，變化變慢）。這個用於改善柏林噪聲算法的fade函數能夠表示爲如下數學形式：

基本上，這就是整個柏林噪聲算法的原理了！搞清了算法的各個實現關鍵步驟後，如今讓咱們着手把代碼實現出來。

---

代碼實現

在本節開始前我須要重申一遍，代碼實現是C#版本。相比Ken Perlin的Java版本實現作了小小的改動，主要是增長了代碼的整潔性和可讀性，支持噪聲重複（瓦片重複）特性。代碼徹底開源，可無償使用（考慮到這畢竟不是我原創發明的算法 - Ken Perlin纔是！）

準備工做

第一步，咱們須要先聲明一個排列表（permutation table），或者直接縮寫爲p[]數組就好了。數組長度爲256，分別隨機、無重複地存放了0-255這些數值。爲了不緩存溢出，咱們再重複填充一次數組的值，因此數組最終長度爲512：

private static readonly int[] permutation = { 151,160,137,91,90,15,                 // Hash lookup table as defined by Ken Perlin.  This is a randomly
    131,13,201,95,96,53,194,233,7,225,140,36,103,30,69,142,8,99,37,240,21,10,23,    // arranged array of all numbers from 0-255 inclusive.
    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
};

private static readonly int[] p;                                                    // Doubled permutation to avoid overflow

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

p[]數組會在算法後續的哈希計算中使用到，用於肯定一組輸入最終挑選哪一個梯度向量（從前面所列出的12個梯度向量中挑選）。後續代碼會詳細展現p[]數組的用法。

接着，咱們開始編寫柏林噪聲函數：

public 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;
    }
    
    int xi = (int)x & 255;                              // Calculate the "unit cube" that the point asked will be located in
    int yi = (int)y & 255;                              // The left bound is ( |_x_|,|_y_|,|_z_| ) and the right bound is that
    int zi = (int)z & 255;                              // plus 1.  Next we calculate the location (from 0.0 to 1.0) in that cube.
    double xf = x-(int)x;
    double yf = y-(int)y;
    double zf = z-(int)z;
    // ...
}

上面的代碼很直觀。首先，對輸入座標使用求餘運算符%，求出[0,repeat)範圍內的餘數。緊接着聲明xi, yi, zi三個變量。它們表明了輸入座標落在了哪一個單元正方形裏。咱們還要限制座標在[0,255]這個範圍內，避免訪問數組p[]時，出現數組越界錯誤。這也產生了一個反作用：柏林噪聲每隔256個整數就會再次重複。但這不是太大的問題，由於算法不只能處理整數，還能處理小數。最後，咱們經過xf, yf, zf三個變量（也就是x,y,z的小數部分值），肯定了輸入座標在單元正方形裏的空間位置（就是前面所示的小藍點）。

Fade函數

如今咱們須要用代碼表示前面所提到的fade函數（圖5）。正如上文所提，函數的數學表示：

代碼實現以下：

public static double fade(double t) {
                                                        // Fade function as defined by Ken Perlin.  This eases coordinate values
                                                        // so that they will ease towards integral values.  This ends up smoothing
                                                        // the final output.
    return t * t * t * (t * (t * 6 - 15) + 10);         // 6t^5 - 15t^4 + 10t^3
}

public double perlin(double x, double y, double z) {
    // ...

    double u = fade(xf);
    double v = fade(yf);
    double w = fade(zf);

    // ...
}

代碼所計算得出的u / v / w變量將在後面的插值計算中使用到。

哈希函數

柏林噪聲哈希函數用於給每組輸入計算返回一個惟1、肯定值。哈希函數在維基百科的定義以下：

哈希函數是一種從任何一種數據中建立小的數字「指紋」的方法，輸入數據有任何細微的不一樣，都會令輸出結果徹底不同

下面代碼就是柏林噪聲算法所使用的哈希函數。它使用了早前咱們聲明的p[]數組：

public double perlin(double x, double y, double z) {
    // ...

    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)];

    // ...
}

public int inc(int num) {
    num++;
    if (repeat > 0) num %= repeat;
    
    return num;
}

代碼的哈希函數，對包圍着輸入座標（小藍點）的周圍8個單元正方形的索引座標進行了哈希計算。inc()函數用於將輸入值增長1，同時保證範圍在[0,repeat)內。若是不須要噪聲重複，inc()函數能夠簡化成單純將輸入值增長1。因爲哈希結果值是從p[]數組中獲得的，因此哈希函數的返回值範圍限定在[0,255]內。

梯度函數

我時常認爲Ken Perlin的最第一版算法裏的grad()函數寫法過於複雜，使人費解。咱們只要明白grad()函數的做用在於計算隨機選取的梯度向量以及頂點位置向量的點積。Ken Perlin巧妙地使用了位翻轉(bit-flipping)技巧來實現：

public static double grad(int hash, double x, double y, double z) {
    int h = hash & 15;                                    // Take the hashed value and take the first 4 bits of it (15 == 0b1111)
    double u = h < 8 /* 0b1000 */ ? x : y;                // If the most significant bit (MSB) of the hash is 0 then set u = x.  Otherwise y.
    
    double v;                                             // In Ken Perlin's original implementation this was another conditional operator (?:).  I
                                                          // expanded it for readability.
    
    if(h < 4 /* 0b0100 */)                                // If the first and second significant bits are 0 set v = y
        v = y;
    else if(h == 12 /* 0b1100 */ || h == 14 /* 0b1110*/)  // If the first and second significant bits are 1 set v = x
        v = x;
    else                                                  // If the first and second significant bits are not equal (0/1, 1/0) set v = z
        v = z;
    
    return ((h&1) == 0 ? u : -u)+((h&2) == 0 ? v : -v); // Use the last 2 bits to decide if u and v are positive or negative.  Then return their addition.
}

下面代碼則是以另外一種使人容易理解的方式完成了這個任務（並且在不少語言版本的運行效率都優於前面一種）：

// Source: http://riven8192.blogspot.com/2010/08/calculate-perlinnoise-twice-as-fast.html
public static double grad(int hash, double x, double y, double z)
{
    switch(hash & 0xF)
    {
        case 0x0: return  x + y;
        case 0x1: return -x + y;
        case 0x2: return  x - y;
        case 0x3: return -x - y;
        case 0x4: return  x + z;
        case 0x5: return -x + z;
        case 0x6: return  x - z;
        case 0x7: return -x - z;
        case 0x8: return  y + z;
        case 0x9: return -y + z;
        case 0xA: return  y - z;
        case 0xB: return -y - z;
        case 0xC: return  y + x;
        case 0xD: return -y + z;
        case 0xE: return  y - x;
        case 0xF: return -y - z;
        default: return 0; // never happens
    }
}

以上的源碼能夠點擊這裏查看。不管如何，上面的兩種實現並無實質差異。他們都是從如下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)

隨機挑選結果其實取決於前一步所計算得出的哈希值（grad()函數的第一個參數）。後面3個參數則表明由輸入點指向頂點的距離向量（最終拿來與梯度向量進行點積）。

插值整合

通過前面的幾步計算，咱們得出了8個頂點的影響值，並將它們進行平滑插值，得出了最終結果：

public double perlin(double x, double y, double z) {
    // ...

    double x1, x2, y1, y2;
    x1 = lerp(    grad (aaa, xf  , yf  , zf),           // The gradient function calculates the dot product between a pseudorandom
                grad (baa, xf-1, yf  , zf),             // gradient vector and the vector from the input coordinate to the 8
                u);                                     // surrounding points in its unit cube.
    x2 = lerp(    grad (aba, xf  , yf-1, zf),           // This is all then lerped together as a sort of weighted average based on the faded (u,v,w)
                grad (bba, xf-1, yf-1, zf),             // values we made earlier.
                  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);
    
    return (lerp (y1, y2, w)+1)/2;                      // For convenience we bind the result to 0 - 1 (theoretical min/max before is [-1, 1])
}

// Linear Interpolate
public static double lerp(double a, double b, double x) {
    return a + x * (b - a);
}

---

利用倍頻實現更天然的噪聲

最後讓咱們再思考下，除了前面所講的計算，還有其餘辦法能夠令噪聲結果更加天然嗎？雖然柏林噪聲算法必定程度上模擬了天然噪聲，但仍沒有徹底表現出天然噪聲的不規律性。舉個現實例子，現實地形會有大段連綿、高聳的山地，也會有丘陵和蝕坑，更小點的有大塊岩石，甚至更小的鵝卵石塊，這都屬於地形的一部分。那如何讓柏林噪聲算法模擬出這樣的天然噪聲特性，解決方法也很簡單：咱們可使用不一樣的頻率（frequencies）和振幅（amplitudes）參數進行多幾回柏林噪聲計算，而後將結果疊加在一塊兒。頻率是指採樣數據的間隔，振幅是指返回值的幅度範圍。

圖6：不一樣頻率和振幅參數下的柏林噪聲結果

將全部結果疊加在一塊兒，咱們就能獲得如下結果：

圖7：圖6全部噪聲的疊加結果

很明顯，這樣的噪聲結果更加使人信服。上面的6組噪聲被稱之爲噪聲的不一樣倍頻(Octave)。隨着倍頻增大，噪聲對於最終疊加噪聲的影響程度變小。固然，倍頻組數的增長，會線性地增長代碼執行時間，在遊戲運行時使用噪聲算法，再好不要使用超過幾組倍頻（好比，當你想在60fps下模擬火焰特效時，最好不要這麼幹）。然而，作數據預處理時，就很適合使用多組倍頻疊加來模擬更天然的噪聲（好比用於提早生成遊戲地形等）。

那咱們應該分別挑選多大的頻率和振幅來進行噪聲計算呢？這個能夠經過persistence參數肯定。Hugo Elias對persistence的定義使用以下：

以上公式i的值取決於倍頻數量，代碼實現也很簡單：

public double OctavePerlin(double x, double y, double z, int octaves, double persistence) {
    double total = 0;
    double frequency = 1;
    double amplitude = 1;
    double maxValue = 0;  // Used for normalizing result to 0.0 - 1.0
    for(int i=0;i<octaves;i++) {
        total += perlin(x * frequency, y * frequency, z * frequency) * amplitude;
        
        maxValue += amplitude;
        
        amplitude *= persistence;
        frequency *= 2;
    }
    
    return total/maxValue;
}

---

小結

以上就是算法的所有內容，咱們如今可使用算法制造噪聲了。再次聲明，你能夠點擊這裏查看完整源碼。

---

本文原文出自Understanding Perlin Noise，翻譯不易，若是你已經看到這裏了，順手點個讚唄！