今天学习的是GDC 2015上分享的基于FFT的海洋渲染方案，其中通过参考知乎上的相关资料，对基于FFT的海洋模拟实现逻辑做了相对清晰的推导与总结：
- 通过贴图来给出某个2D的频域信号，其中贴图的uv代表两个维度的频率，数值代表该信号的振幅
- 基于上述频域信号，可以通过逆变换得到随时间变化的2D空域信号，这个变换过程需要在运行时完成（涉及到时间，当然也可以烘焙多张，保证首尾衔接的话即可）
- 频域信号用的是菲利普频谱，基于这个频谱可以在运行时计算各个频率的信号数值，也可以在离线完成计算。运行时计算的好处是：
- 可以根据当前环境的状态比如风向等来调制得到更拟真的效果
- 可以运行时调整模拟参数，实现海洋效果的动态变化
具体可以参考实现细节：
游戏中的海洋效果：有刺客信条III以及刺客信条黑帆
Crysis的海洋效果
War Thunder的海洋效果
海洋奇缘
泰坦尼克号
这些游戏跟电影有如下的一些共同点
- 都是基于FFT来完成模拟的
- 模拟的输入是一个频谱，之后对这个频谱做逆向的离散傅里叶变换就得到时域（空域）信号
- 时域信号是一个随时间变化的高度displacement数据
为啥是FFT：DFT耗时太高了，时间复杂度是，而FFT的时间复杂度则是
在GPU中计算FFT的算法是Cooley-Tukey（蝴蝶）算法。
傅里叶时域（还是频域？）贴图的分辨率应该是多少呢？正常情况下128 ~ 512就够了，泰坦尼克号跟Waterworld也就只用来2048，超过这个尺寸可能会导致浮点精度问题。
这里选用的分辨率是512。
傅里叶变换的公式（可以参考[4]了解一下具体的推导过程）是：
逆变换公式为：
在计算机中这个积分公式是没有解析解的，只能转换为离散的累加形式。
再来看离散形式的变换与逆变换，这里的N是说将分割成N份，N越大，我们得到的频域信号就越丰富，信号就越逼真，模拟越接近，但是计算复杂度也就越高：
参考[1]，DFT的计算其实就是一个矩阵的运算，而这个矩阵运算则可以通过排列得到如上图所示的蝴蝶累加算法，从而可以很好的实现加速计算。
这里IFFT的实现是通过compute shader完成的，采用的是Radix-2算法（参考[3]），细节就不说了，网上有大把的开源实现。
上面介绍的是一个一维函数的FFT实现，那么回到水体模拟上，这个公式是如何应用起来的呢？
我们知道，傅里叶变换其实是通过将信号拆解为正弦波的方式来完成对信号的表达的，每个正弦波都有一个频率，一个信号的频域表示，其实就是其拆解后得到的正弦波的振幅与频率的表示，频域信号在某个频率上的取值就是当前对应的时（空）域信号在对应频率上的振幅。
回到海洋这边来，其实海浪可以看成是由空间上无穷个2D的正弦波信号叠加而成，同样的，这个信号在频域中的表示就是海浪在不同频率的正弦波上的振幅。不同的是，这个频率是2D的，有X&Z两个方向的频率，同样的，我们的频域信号也是2维的，即在X & Z上具有不同频率的信号的振幅。
用公式来表达，先不考虑海洋波形随着时间波动的变化，我们先看下某个时间点的海洋波形，这个波形刚说了，可以看成是无穷多个2D的正弦波组成的，在空域上，可以表示为h(x, z)，即在坐标(x, z)处的水体高度。
这个水波对应的频域信号，我们可以表示为，两者可以通过傅里叶变换（空域->频域）与逆变换（频域->空域）实现转换，在海洋的实现逻辑中，我们的输入通常是频域信号，且只能用IDFT的方式来完成计算，因此按照前面的说明，这个逆变换的公式给出如下：
其中是频域上的2D向量，对应于前面的
，或者用k做前缀，将之改写为
，根据我们FFT频域模拟的尺寸，这里有
，n跟m分别是频域贴图（贴图的uv代表着两个方向的频率，对应的数值对应的是这个频率的幅度）的水平跟垂直的分辨率，也就是我们在频谱上的水平跟垂直方向上以
为单位采样的点数，其中
又分别为海面X or Z方向上的patch尺寸，根据[4]中的解释，我们的采样点是分布在原点两侧的
在这样的采点策略下，针对海洋波形，我们就总计拿到了MxN个频域的信号（或者说正弦波），通过他们之间的叠加与逆变换得到海洋的波形，通常情况下，这两个数值越大，波形精度越高，模拟越逼真。
则是空域上的2D信号，对应的是当前点的xz坐标，同样的，在计算displacement的时候，我们也是针对离散点进行的，即
也是在海平面上以
为单位采样NxM个点（虽然通常N=M），同样的，在相同的覆盖范围下，这俩数值越大，顶点之间的密度就越大，水体精度就越高，当然计算消耗也同样会增加。
另外需要注意的是，指数上的两个向量的乘法是点乘，即有：
在这个公式下我们就有当某个维度固定住的时候，我们得到的就是一个一维的信号。
在这个基础上，如果我们考虑水波随着时间变化的过程，那么我们可以将前面的频域信号改写为时间的函数，即有：
接下来对这个公式做一个拆解。
里面最大的问题是频域信号怎么得到，这里就没有介绍具体的背景，而是直接给出了结论，海洋波形的频域信号通常使用的是菲利普频谱，其公式为：
这个公式中跟
是共轭复数（两个实部相等，虚部互为相反数的复数互为共轭复数），k则是
的模，
是一个函数，其计算公式为：
这里的g是重力加速度。
再来看，有如下的计算公式：
是相互独立的随机数，均服从均值为0，标准差为1的正态分布。
这个公式中的是归一化后的风向（2D），L的计算公式为：
这里的g依然是重力加速度，V则是风速。
按照[4]的描述，这里的公式是基于统计给出的经验公式，并无过多的解析或物理含义在里面，可以不用过于纠结其来源。
这里对公式中各个变量的含义做了基本的介绍，这里补上前面遗漏的项，前面介绍过的这里就不再介绍：
- a是全局的海浪的振幅
基于高度偏移的公式，我们可以通过解析的方式得到其法线。
先来计算一下梯度向量（某个函数即这里的高度变化率最大的方向，因为这里的变量还都是水平方向的两个轴，因此这里的向量还是一个水平方向的向量），因为不包含梯度计算相关的
信息，因此可以直接移出来：
继续推导，我们有：
代入前面的公式，我们就可以得到梯度向量的计算公式：
针对垂直向上的向量、梯度向量与法线向量的关系[4]:
注意梯度向量本身是没有归一化的，我们可以有：，将前面的公式代入，则我们有：
前面得到的是高度方向的偏移，但是我们实际上需要的是三个方向的偏移，因此还需要再额外处理一下另外两个方向的。
在Gerstner Wave的实现中，为了得到尖锐的波峰效果（choppy wave），需要基于如下公式[4]在XZ方向上做挤压：
而同样的，FFT实现的海浪效果也需要做挤压，挤压公式为：
如下图[4]所示，当挤压力度过大的时候，就会出现水波的相互穿刺，即可以理解为水波碰撞在一起，此时会产生浮沫：
这个过程其实就是某个形状（比如上图中的近似四边形，通常称之为一个图元）发生了翻转，正面翻到了背面，用数学概念来描述，就是该图元有向面积变为负数，而这个在数学公式上来表达就是对应的雅可比矩阵计算结果为负数，具体可以参考[4]的推导。
挤压强度假设为，我们就有：
而这个公式跟前面Gerstner Wave的挤压公式整体形式上其实是一样的，只是换了一种写法。
经过上面计算，我们就得到了空域中各个点的偏移向量，这个可以用一张3通道的贴图来存储，我们就有了上图所示的Displacement Map
参考NVIDIA的解说图，我们回顾一下上面的计算逻辑：
- 我们基于频谱函数或者说频谱贴图
就能得到高度偏移
- 同样基于频谱函数跟挤压公式，我们可以得到在水平方向上的偏移
- 将上述三者合到一起，我们就有了displacement map
- 基于Displacement map，我们就有计算出法线贴图，同时基于雅可比矩阵计算结果，就能够计算出foam intensity map（也就是上图中的folding）
经过上述计算之后，我们就能得到上述相对真实的海面模拟效果
这里采用的是Crysis Engine最开始实现海洋方案时的mesh方案，有点类似于Projected Grid，不同的是，这里选择的不是基于场景相机的视角做的投影，而是基于场景相机为中心，俯视角做的投影（那不就是挂在相机上的一个固定的mesh吗，只是不需要存储mesh的数据而已）。
这个方案实现简单，可以自动实现近密远疏效果，而且据说不会有瑕疵
从原理上来分析，Projected Grid那种随着相机移动或旋转而导致顶点跳变的问题应该还是会有，但这里又单独提了Projected grid的瑕疵，意思就是这个问题确实是被消灭了，不知道这里具体干了啥。。。
这是线框模式下的网格展示
模拟逻辑用C++实现，着色则通过材质实现。
水体着色这里主要关注反射跟折射效果。
而菲涅尔项则用于给出反射跟折射光强的比例，但是完整的菲涅尔计算逻辑太过复杂，这里采用的是Schlick Fresnel近似计算公式。
也就是我们高光BRDF计算中常用的那个，如果避开复杂的指数计算（如何避开？），实测这个近似计算公式要比原始公式快30%。
得到菲涅尔项之后，就可以通过lerp实现折射跟反射的混合。
除了这个计算逻辑之外，还借用了UE的SSR来添加带有遮挡信息的环境反射效果。
[1]. 傅立叶变换简记
[2]. 离散傅里叶变换DFT详解及应用
[3]. FFT原理及实现(Radix-2)
[4]. fft海面模拟(一)
