真實感海洋的繪制（二）：使用快速傅裏葉變換加速波形計算
真實感海洋的繪制（二）：使用快速傅裏葉變換加速波形計算
其實上一篇博文所寫的\(H(\vec{x},t)?\)，就是二維傅裏葉變換的求和式，之前的暴力計算法屬於二維的離散傅裏葉變換（Discrete Fourier Transform, DFT），利用二維的快速傅裏葉變換（Fast Fourier Transform, FFT）可以將復雜度從\(O(n^4)?\)降低到\(O(n^2\log{n})?\)。
如果讀者不熟悉FFT，強烈建議閱讀《算法導論》相關章節，個人感覺沒有什麽資料講得比《算法導論》更清楚了。書後習題還有關於二維FFT的思考，是本文算法正確性的理論基礎。
數學推導
高度場
出於數學上推導的簡單和實現上的簡單，我們約定\(L_x=L_y=N=M=2^k\)
\[ H(\vec x,t) = \sum_{\vec{k}}h(\vec k, t)e^{i\vec{k}\cdot \vec{x}} \mbox{, where } \vec{k}=(2\pi n/L,2\pi m/L), \ \vec{x}=(x, z)\-N/2\le n, m < N/2 \]
將\(\vec{k}\)與\(\vec{x}\)帶入整理得
\[
\begin{align}
H(\vec{x},t)&=\sum_{n=-N/2}^{N/2 - 1}\sum_{m=-N/2}^{N/2-1}h(\vec{k},t)e^{2\pi nxi/N+2\pi mz/N} \&=\sum_{n=-N/2}^{N/2-1}e^{2\pi nxi/N}\sum_{m=-N/2}^{N/2-1}h(\vec{k},t)e^{2\pi mzi/N}
\end{align}
\]
為了化成便於使用FFT的形式，令\(n‘=n+N/2\), \(m‘=m+N/2\)
\[
H(\vec{x},t)=\sum_{n‘=0}^{N-1}e^{2\pi(n‘-N/2)xi/N}\sum_{m‘=0}^{N-1}h(\vec{k},t)e^{2\pi(m‘-N/2)zi/N}
\]
其中
\[
\begin{align}
e^{2\pi(n‘-N/2)xi/N}&=e^{2\pi n‘xi/N}e^{-\pi xi}\&=(-1)^x e^{2\pi n‘xi/N} (\mbox{ note that }e^{\pi i}=-1, N\ge 2)
\end{align}
\]
所以
\[
H(\vec{x},t)=(-1)^{x+z}\sum_{n‘=0}^{N-1}e^{2\pi n‘xi/N}\sum_{m‘=0}^{N-1}h(\vec{k},t)e^{2\pi m‘zi/N}
\]
這樣已經成為了可以進行FFT的形式。這兒出現了一個\((-1)^{x+z}\)，不用擔心，在全部計算完之後帶入即可。
法線
關於法線的計算有一些tricky。我們先寫出之前給出的法線公式
\[
\epsilon(\vec{x},t)=\nabla H(\vec{x},t)=\sum_{\vec{k}}i\vec{k}h(\vec{k},t)e^{i\vec{k}\cdot \vec{x}}\\begin{align}
\vec{N}(\vec{x},t)&=(0,1,0)-(\epsilon_x(\vec{x},t),0,\epsilon_z(\vec{x},t))\&=(-\epsilon_x(\vec{x},t),1,-\epsilon_z(\vec{x},t))
\end{align}
\]
問題在於計算\(\epsilon(\vec{x},t) = (\epsilon_x, \epsilon_z)\)。我們首先類似\(H\)的推導得到
\[
\epsilon(\vec{x},t)=(-1)^{x+z}\sum_{n‘=0}^{N-1}e^{2\pi n‘xi/N}\sum_{m‘=0}^{N-1}i\vec{k}h(\vec{k},t)e^{2\pi m‘zi/N}
\]
經過觀察我們發現\(\epsilon(\vec{x},t)\)的兩個方向分別只與\(\vec{k}\)的兩個方向有關，那麽可以寫成如下形式
\[
\epsilon_x(\vec{x},t)=(-1)^{x+z}\sum_{n‘=0}^{N-1}e^{2\pi n‘xi/N}\sum_{m‘=0}^{N-1}ik_xh(\vec{k},t)e^{2\pi m‘zi/N}\\epsilon_z(\vec{x},t)=(-1)^{x+z}\sum_{n‘=0}^{N-1}e^{2\pi n‘xi/N}\sum_{m‘=0}^{N-1}ik_zh(\vec{k},t)e^{2\pi m‘zi/N}\\]
然後分別計算即可。
偏置量
類似地，我們寫出偏置量的公式
\[
\vec{D}(\vec{x},t)=\sum_{\vec{k}}-i\frac{\vec{k}}{|\vec{k}|}
h(\vec{k},t)e^{i\vec{k}\cdot \vec{x}}
\]
和法線計算的方式完全一樣得到
\[
\vec{D}(\vec{x},t)=(-1)^{x+z}\sum_{n‘=0}^{N-1}e^{2\pi n‘xi/N}\sum_{m‘=0}^{N-1}-i\frac{\vec{k}}{|\vec{k}|}h(\vec{k},t)e^{2\pi m‘zi/N}
\]
分成兩個方向
\[
D_x(\vec{x},t)=(-1)^{x+z}\sum_{n‘=0}^{N-1}e^{2\pi n‘xi/N}\sum_{m‘=0}^{N-1}-i\frac{k_x}{|\vec{k}|}h(\vec{k},t)e^{2\pi m‘zi/N} \D_z(\vec{x},t)=(-1)^{x+z}\sum_{n‘=0}^{N-1}e^{2\pi n‘xi/N}\sum_{m‘=0}^{N-1}-i\frac{k_z}{|\vec{k}|}h(\vec{k},t)e^{2\pi m‘zi/N}
\]
這就完成了推導的工作。相信熟悉FFT的讀者可以對上邊的公式很容易給出一個自己的FFT計算實現。
實現結果
之後的博客將會研究如何對這樣的海面進行真實感渲染和光照計算。如果有時間的話，可能會探索如何把這個FFT模型移植到GPU上並行計算以實現更好的性能。
參考文獻
- Tessendorf, Jerry. Simulating Ocean Water. In SIGGRAPH 2002 Course Notes #9 (Simulating Nature: Realistic and Interactive Techniques), ACM Press.
- Thomas H. Cormen, Charles E. Leiserson, Ronald L. Rivest and Clifford Stein. Introduction to Algorithms, 3rd Edition, MIT Press.
真實感海洋的繪制（二）：使用快速傅裏葉變換加速波形計算