柏林噪聲是一個非常強大演算法，經常用於程式生成隨機內容，在遊戲和其他像電影等多媒體領域廣泛應用。演算法發明者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

圖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)

注意：許多介紹柏林噪聲演算法的文章都是根據最初版柏林噪聲演算法來講解的，預定義的梯度表不是本文所說的這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變數將在後面的插值計算中使用到。

雜湊函式

柏林噪聲雜湊函式用於給每組輸入計算返回一個唯一、確定值。雜湊函式在維基百科的定義如下：

雜湊函式是一種從任何一種資料中建立小的數字“指紋”的方法，輸入資料有任何細微的不同，都會令輸出結果完全不一樣

下面程式碼就是柏林噪聲演算法所使用的雜湊函式。它使用了早前我們宣告的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;
}

小結

以上就是演算法的全部內容，我們現在可以使用算法制造噪聲了。再次宣告，你可以點選這裡檢視完整原始碼。