真实感海洋的绘制（二）：使用快速傅里叶变换加速波形计算
其实上一篇博文所写的\(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.