基於快速傅立葉轉換的水波生成(2)

• 前言

嚴格來說，這篇不是介紹如何用快速傅立葉轉換做水波生成，
而是用離散傅立葉轉換做水波生成。

為了方便理解如何用傅立葉轉換做水波生成，我將文章分為4篇，從簡而繁去做介紹。
想看主題的話可以直接看第4篇。(前提是有人看...)


本篇的文章是參考 Ocean simulation part one: using the discrete Fourier transform
(以下用作者名稱 Keith Lantz 簡稱)這個網頁的說明，製作出來的。

這個做法又稱為 spectrum-based Ocean，
Keith Lantz 根據<Simulating Ocean Water>[1] 的記載使用 Phillip's spectrum，
來模擬海浪波形。並使用二維反向傅立葉轉換將頻譜轉換到每個X-Z平面的浪高，藉此渲染出海洋形狀。

Keith Lantz 利用CPU運算實現了文獻上的做法，

本文則是挑戰修Keith Lantz的做法，改用GPU運算實現，並公開在GitHub上作為自己曾經走過遊戲業的痕跡。

不再廢話了，開始進入正題，

• 演算法分析

<Simulating Ocean Water>[1]定義了在時間t，海平面座標(x,z)的高度為


Keith Lantz 為了方便計算，將 n 和 m 正規化到 0 ~ N-1 和 0 ~ M - 1，並命名為 n '和 m'，
代入方程式得到



現在將焦點放在頻譜函式(註1)




最後，代入 n', m' 得到


就完成了整個演算法所需要的環節。

• Unity 程式碼實現

從演算法內容可以看出，初始頻率函式是個不隨時間變化的常數值，所以先計算每個(n',m')位置的初始頻率值，然後記錄在2維貼圖上，供GPU運算讀取

phillips spectrum 程式碼實現如下

float Phillips(int _n_prime, int _m_prime)
{
    Vector2 k = new Vector2(Mathf.PI * (2 * _n_prime - N) / LenN, Mathf.PI * (2 * _m_prime - M) / LenM);
    float k_length = k.magnitude;
    if (k_length < 0.000001f) return 0.0f;

    float k_length2 = k_length * k_length;
    float k_length4 = k_length2 * k_length2;

    k.Normalize();

    float k_dot_w = Vector2.Dot(k, Wind.normalized);
    float k_dot_w2 = k_dot_w * k_dot_w;

    float w_length = Wind.magnitude;
    float L = w_length * w_length / GRAVITY;
    float L2 = L * L;

    float damping = 0.001f;
    float l2 = L2 * damping * damping;

    //return Amp * Mathf.Exp(-1.0f / (k_length2 * L2)) / k_length4 * k_dot_w2;
    return Amp * Mathf.Exp(-1.0f / (k_length2 * L2)) / k_length4 * k_dot_w2 * Mathf.Exp(-k_length2 * l2);
}

本程式碼是參考 Keith Lantz 所改寫，和上個章節內容有著一些明顯的不同，

Keith Lantz 新增了 Amplitude 參數在他的演算法內，

並且回傳值多乘了  Mathf.Exp(-k_length2 * l2)

參考網路上的一些討論，似乎是為了修正在高頻區域會產生無法收束的狀況而做的一個修正

高斯隨機變數 程式碼實現如下

Vector2 GaussianRandomVariable()
{
    float x1, x2, w;
    do
    {
        x1 = 2.0f * UnityEngine.Random.value - 1.0f;
        x2 = 2.0f * UnityEngine.Random.value - 1.0f;
        w = x1 * x1 + x2 * x2;
    }
    while (w >= 1.0f);

    w = Mathf.Sqrt((-2.0f * Mathf.Log(w)) / w);
    return new Vector2(x1 * w, x2 * w);
}

得到了phillips spectrum 和 高斯隨機變數，初始頻率函式就跟著出來了

Vector2 H_Tilde0(int _n_prime, int _m_prime)
{
    Vector2 r = GaussianRandomVariable();
    return r * Mathf.Sqrt(Phillips(_n_prime, _m_prime) / 2.0f);
}

分母的根號2其實是在計算高斯隨機變數使用的，但為了節省計算，搬到外面和 phillips spectrum 一起處理

波浪的分散關係式

float Dispersion(int _n_prime, int _m_prime)
{
    // w(k) = (gk)^0.5
    float w_0 = 2.0f * Mathf.PI / 200.0f;
    float kx = Mathf.PI * (2 * _n_prime - N) / LenN;
    float kz = Mathf.PI * (2 * _m_prime - N) / LenM;

    return Mathf.Floor(Mathf.Sqrt(GRAVITY * Mathf.Sqrt(kx * kx + kz * kz)) / w_0) * w_0;
    //return Mathf.Sqrt(GRAVITY * Mathf.Sqrt(kx * kx + kz * kz));
}

在這裡也可以看出 Keith Lantz 做了部分改寫，但為了方便比較，
我以 Keith Lantz 為主

最後頻譜函數常數部分的實現方式

每個座標(n',m')都有5個常數值，將這5個常數值儲存在貼圖的5個channel(程式碼用到兩張貼圖)，剩下的就交由GPU在runtime處理。

GPU運算的實現使用shader來處理，程式碼如下

• Unity下的執行結果

專案位置

海水渲染shader是我自己寫的，可以依照專案需求更改

註1.頻譜函數是我自己取的名字，無論是在文獻[1]，或是 Keith Lantz  的網頁中，都沒有給予一個名稱，只說明是回傳一個複數振幅的函式

註2.這個也是我自己取的名字，但他的確是求出一個起始的頻譜函式

參考
[1] SIGGRAPH 2001, Tessendorf J. Simulating ocean water[J]. Simulating nature: realistic and interactive techniques.