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