渲染三维场景时经常会遇到需要渲染各种水体的情况，比如湖泊、河流、海洋等，不仅需要水体表面要有接近真实的随时间而变化的波动，还要有令人信服的颜色、反射光、透明度等细节。实时渲染水面的方法有很多，从简单的若干正弦波叠加，到《GPU Gems》中介绍的叠加Gerstner波的方法，再到如今GPU在线计算FFT得到表面高度，都是以追求效果更加逼真的同时保证计算的高效实时。我用OpenGL实现了通过合成Gestner波产生水波的方法，具体过程如下。

开发环境

Windows 7 x64，Visual Studio 2010，OpenGL版本3.0，GLSL版本1.3。

freeglut 2.8.1，GLM 0.9.5.1。GLM用于产生模型视图矩阵、透视投影矩阵和法线变换矩阵。

正弦函数波

        在一些数学书中介绍正弦函数时会提到“理想情况下的水波是正弦形状的”，但实际上，单独的水波应该是波峰尖、波谷宽的。如果用正弦波来表现这样的效果，可以选择如下变换：

        由于正弦函数的值域是[-1,1]，缩放到[0,1]区间，再做幂运算，会使函数值减小，而且距离0越小的值减小得越多。这样就能产生波峰尖、波谷宽的形状。

        下面是一组k分别等于1.0，1.5和2.0时的情况，可见k越大（k>=1），形状就越明显。

        但是只有一个参数决定这种形状过于简单，而且在CG中希望在细节多的地方（波峰）网格点较密集，在细节少的地方（波谷）网格点较稀疏。用正弦函数绘制时，如果想提高细节，只能整体提高x的细分程度，也会在波谷处增加大量的多余计算。

Gerstner波

        Gerstner波的诞生早于计算机图形学（CG），它最初在物理中用于水波的模拟。由于它的形状比较真实，而且计算量不大，所以被广泛用于CG中水波的模拟。

        Gerstner波以参数方程的形式给出：

        自变量为p，参数Q、D、A用来控制形状。Q控制波峰的尖锐度，D控制波长，A为振幅。

        Q应为较小的值，若Q很大，会在波峰处产生环，破坏波的形状。比如：

        观察x(p)的表达式可以看出，与正弦波相比，Gerstner波在波峰处的点更紧凑，在波谷处更稀疏：

波的合成

        为了产生真实的水面，需要把若干不同方向、不同参数的Gerstner波合成为一个波。即：

        在三维空间中绘制水波这样高度值频繁变化的面时，一般采用规则网格来绘制，即在x-y平面上画一张均匀的网格，对网格上的每一个点计算它的高度值（z值），这样就产生了一张高低起伏的面。随着时间的变化，每个点的高度也随之变化，就产生了动态的面。

        为了把这张网格与二维的Gerstner波结合起来，需要进行如下转换：

        假设二维Gerstner波表示为y=f(x)，三维网格表示为z=g(x,y)。则：

        (x0,y0)表示波的起点，theta角表示波传播的方向。

        初始时，网格上每点的高度设为0，每叠加一个波，就根据上面的式子计算出一个高度，加在z上。计算完所有的波后，就实现了多个波的叠加。

        由于Gerstner参数方程也在改变x（即上图的d），直接应用原式计算会增加复杂度。同时，为了尽可能地减小计算量，我采用两种固定形状的Gerstner波，每种波用11对坐标表示，计算f(x)时只需要在这11个点中计算线性内插即可。

两种波形。第一个波峰较尖，用来绘制细小的水波，第二个波峰较宽，用来绘制波长较长的水波。

线性内插函数：

用一个结构体来保存时间和每个波的波长、振幅、方向、频率和起始坐标：

计算每个网格点的z值：

法线的计算

        为了计算光照，需要知道每个顶点的法线方向。法线是根据网格上相邻4点的坐标计算的。

        由于每个点与四个网格面相邻，每个面有独立的法线，所以应该把4个面的法线的均值作为这个顶点的法线。所谓的“均值”就是把4个规范化的法线相加在除以4。CG中涉及法线的计算一般都需要在计算前规范化法线向量（使其长度为1），但因为每个网格面在x和y方向上的尺寸相同，所以法线不需要规范化就可以相加，只对加和后的法线做一次规范化即可。这样就减少了四次规范化和一次除法。

规范化函数：

        这个函数里需要注意的是，由于使用浮点数进行计算，在小数点后第8位开始容易产生误差，如果输入向量每个分量都很小，则很难保证计算结果的正确性；同时，如果输入向量的长度很小，则需要先把输入向量扩大10000倍，再进行计算，以减小误差对计算结果的影响。这个函数还支持原地计算（即输入和输出为同一个向量），适用范围是很广的。

法线计算形式如下：

网格的绘制

        最初编写代码时，我考虑过在vertex shader中计算网格坐标和投影矩阵，但在vertex shader中无法读取相邻点的坐标，不能用上面的方法计算法线。我改用根据gerstner_pt_a[]和gerstner_pt_b[]计算单个波的法线，再合成每个波的法线来近似顶点的法线（在理论上是不严谨的），但由于GLSL 1.3没有现成的对矩阵求逆的函数（inverse()在GLSL 1.5才开始支持），从而无法快捷地获得NormalMatrix，不能得到正确的法线方向。所以只好放弃，采用离线计算网格坐标、法线和几个变换矩阵。

        我使用VAO储存坐标，以GL_TRIANGLE_STRIP方式绘制的方法绘制网格。由于上面计算的坐标是按照从x到y逐行存储的，不能直接用于三角形绘制，原因如下图所示：

        之前计算得到的坐标保存在了pt_strip[]中，需要把里面的点重新排序存入vertex_data[]中，再把vertex_data[]中的数据放入VAO（储存在显存中）。变换代码如下：

        对法线的处理也是如此。产生和绑定VAO的代码如下：

绘制网格的代码如下：

使用GLM计算模型视图矩阵、透视投影矩阵和法线变换矩阵：

        为了使glm::perspective()的参数符合自己的习惯，我把glm/gtc/matrix_transform.inl的第254到264行改成如下内容：

vertex shader代码如下：

绘制的网格效果如下：

添加材质和光照

        我使用了一张512*512大小的贴图，格式为tga，图像可以在文章末尾所给的链接获取。

读取和绑定贴图：

        initTexture()函数用来读取贴图文件和生成贴图对象，这里就不详细介绍了，具体细节可以参考文章末尾给出的代码链接。为了试验法线贴图，我使用了两张贴图，在每次绘制时要交替启用两张贴图，以保证fragment shader能同时读到两张贴图：

添加材质后的效果：

        经过比较，我觉得Ward的各向异性光照效果比较好，其原理可以参考：D3DBook:(Lighting) Ward – GPWiki。

fragment shader代码如下：

添加光照后的效果：

全部代码和贴图可以在这里查看：

johnhany/OpenGLProjects/GerstnerWave

2014.7.6更新：

        为了使模型视图矩阵和透视投影矩阵更符合OpenGL的规范，调整了模型视图矩阵的计算过程►，透视投影矩阵的构成►，同时相应调整了光源的位置►。

2015.3.29更新：

          该方法在Visual Studio 2012下使用OpenGL 4.5 + GLEW + GLFW的实现可以在这里下载，所需的配置文件在这里下载。

2015.5.9更新：

          修改了网格内坐标变换计算的错误►。

2015.11.26更新：

          该方法在Visual Studio 2012下使用OpenGL 4.5 + GLEW + GLFW的实现在这里下载（项目默认目标平台是x64）。