使用OpenGL模拟水面

本文会对一个基于波动方程来模拟水面实现的OpenGL程序进行分析。会给出程序基本组成，绘制流程示意图。并且对着色器中使用的原理进行推导并给出比较详细的注释。

前言

关于程序源代码，其中使用的着色器其实都是来自于一个别人的js项目，我把它们应用到了Qt中，只做了少量修改。你可以在我之前的文章Qt使用OpenGL实现水波特效中获得源代码和一些基本的解释。
即使你不使用Qt，只要你对这个项目实现的效果感兴趣，想做一些修改或者是移植，也可以阅读本文章，因为尽管不同的图形界面开发框架和语言使用OpenGL的方法不尽相同，但是着色器代码和绘制流程是一致或者至少是相近的，此文应该能帮你节省一些自己分析研究源码的时间。

程序可能会有一点文不对题的嫌疑，尽管模拟水面高度确确实实是程序的一部分，但此程序最终所做的并不是将水面高度可视化。

如果你想更直观的感受到此算法下的水面高度是如何变化的，你可以使用webgl-water项目中的demo，这个项目模拟了很多的物理现象，相应的也就复杂得多。我主要所参考的js项目看起来就是把这个项目中的一部分简化并提取出来了。你可以下载webgl-water的源码并修改参数来获得更为明显的效果。
这里给一张webgl-water的截图吧。

举个具体的例子，如果你下载了此项目的源码，那么在webgl-water项目源码的main.js文件约140行处可以修改生成水波函数的参数。

water.addDrop(pointOnPlane.x, pointOnPlane.z, 0.03, 0.01);

参数0.03影响点击鼠标时生成水波的半径，0.01影响生成水波的强度（高度）。

OpenGL相关资源示意

一共三个着色器程序，各包含一个顶点着色器和一个片段着色器，作用已在图片中说明，代码后文均会给出并解释。

关于两个帧缓冲，并非固定的1用来写，2用来读，在循环绘制的过程中需要不断的交换帧缓冲，使update_program总能读取到最新的水面高度数据，接下来的绘制流程中也会进一步说明。由于我们是使用纹理附件的方式存储水面高度，但水面高度并非是一个RGBA颜色值，而是一个浮点数，所以在设置帧缓冲的纹理附件格式时，与使用普通图片纹理的设置会有一些不同。c++中可以使用下方这行代码中的参数进行设置。

glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA32F_ARB, 宽度,高度, 0, GL_RGBA, GL_FLOAT, NULL);

此外，我们还需要一组顶点，并将他们存入一个VBO和VAO中。
由于此项目中所有的绘制都相当于是一个纹理图片的绘制，所以需要的顶点数组比较简单，就是一个正方形的四个顶点，使用GL_TRIANGLES方式进行绘制的话，则需绘制两个三角形共六个顶点。

static GLfloat vertArray[]={-1.0,1.0,-1.0,-1.0,1.0,-1.0,1.0,1.0,-1.0,1.0,1.0,-1.0
};

绘制流程

循环绘制示意图

通过这个循环绘制过程，便可以在屏幕展示出水面的动态效果，此过程在程序运行期间不断的执行。

如果我没理解错，只使用一个帧缓冲也能完成绘制，但本项目包括本文所提到的参考项目中的绘制流程就是如上图所示，而绘制完一帧后交换帧缓冲的操作其实只交换了帧缓冲的序号而并非是内容。

添加水波过程示意图

此过程只在添加水波时被执行。

着色器代码详解

方便起见我使用常量字符串的形式将着色器代码放到了程序主体中

drop_program

顶点着色器：

static const char* globVert="attribute vec2 vertex;\n""varying vec2 coord;\n""void main() {\n""  coord = vertex * 0.5 + 0.5;\n"" gl_Position = vec4(vertex, 0.0, 1.0);\n""}\n";

由于绘制空间坐标范围是[-1,1]，而纹理坐标范围是[0,1]。coord = vertex * 0.5 + 0.5;这一句用来实现由绘制坐标到纹理坐标的转换。

片段着色器：

static const char* dropFrag="precision highp float;\n""const float PI = 3.141592653589793;\n""uniform sampler2D texture;\n"//此纹理来自于当前读帧缓冲"uniform vec2 center;\n"//波源中心"uniform float radius;\n"//波源半径"uniform float strength;\n"//波源强度"uniform float ratio;\n"//绘制窗口长宽比，用来转换坐标，以在长方形的窗口内绘制圆形的水波"varying vec2 coord;\n"//在顶点着色器中计算的纹理坐标"void main() {\n"" vec4 info = texture2D(texture, coord);\n"//取得对应坐标处的水面高度值"    float x=center.x * 0.5 + 0.5 - coord.x;\n"" float y=(center.y * 0.5 + 0.5 - coord.y)*ratio;\n"" float drop = max(0.0, 1.0 - length(vec2(x,y)) / radius);\n"//"   drop = 0.5 - cos(drop * PI) * 0.5;\n"//这几句会以center为中心，添加一个拱形的波源"  info.r += drop * strength;\n"//纹理的r分量中储存的即是水面高度数据"  gl_FragColor = info;\n""}\n";

update_program

update_program以波动方程为数学依据，负责以帧为单位演算水面高度的变化。想搞懂这个算法是如何工作的，还是需要一些篇幅来进行展开说明的。

在其他人的一篇博客，动态水面的模拟中，我发现这篇文章使用的公式与我参考的项目有非常相近的形式。
他还对该原理的原文进行了翻译，实时流体模拟中的数学方法。
你也可以选择在线阅读Chapter 15 Fluid and Cloth Simulation原文，文档页码462页。

但是此项目使用的算法与该书中提供的形式并不完全相同，我又自己进行了一些推导，最终我认为update_programe使用的算法的确是基于波动方程的，但是在上述文章的基础上进一步做了近似和简化，下面会给出一部分推导过程，请结合上方的文章学习。

关键的不同之处是对于水波衰减的处理。上方文章中波动方程的无阻尼形式是下方推导的基础。

那么首先来看看二维波动方程的无阻尼形式：

∂2z∂t2=c2⟮∂2z∂x2+∂2z∂y2⟯(1)\tag 1\frac {\partial^2 z} {\partial t^2}=c^2 \big\lgroup\frac {\partial^2 z} {\partial x^2}+\frac {\partial^2 z} {\partial y^2} \big\rgroup ∂t2∂2z=c2⎩⎧∂x2∂2z+∂y2∂2z⎭⎫(1)

二维波动方程中，变量z即可代表波平面的高度，接下来，我们使用z(i,j,k)来表示二维平面波在坐标(i,j)处的点在时间k所处的高度。
变量t表示时间
变量x,y表示二维平面坐标。
变量c是波的传播速度。

根据导数的定义，最终式(1)可以变形为：

z(i,j,k+1)−2z(i,j,k)+z(i,j,k−1)t2=c2⟮z(i+1,j,k)−2z(i,j,k)+z(i−1,j,k)d2+z(i,j+1,k)−2z(i,j,k)+z(i,j−1,k)d2⟯(2)\tag 2\frac {z(i,j,k+1)-2z(i,j,k)+z(i,j,k-1)} {t^2}=\\c^2 \big\lgroup\frac {z(i+1,j,k)-2z(i,j,k)+z(i-1,j,k)} {d^2}+\\\frac {z(i,j+1,k)-2z(i,j,k)+z(i,j-1,k)} {d^2} \big\rgroup t2z(i,j,k+1)−2z(i,j,k)+z(i,j,k−1)=c2⎩⎧d2z(i+1,j,k)−2z(i,j,k)+z(i−1,j,k)+d2z(i,j+1,k)−2z(i,j,k)+z(i,j−1,k)⎭⎫(2)
其中d表示我们相邻采样点之间的距离。
这个过程的详细推导可以在上方文章中找到，去掉阻尼项即可。
接下来我们来分析一下(2)式中等式左边部分的物理意义。
结论是，在微小的时间间隔内，我们把波在某一点的z轴方向的运动看做匀加速直线运动，那么等式左边的物理意义即是该点在这段时间内的加速度。
横轴为时间，纵轴为高度。上图表示的都是坐标点(i,j)处的信息，xy平面坐标是不变的。
v1=z(i,j,k)−z(i,j,k−1)tv2=z(i,j,k+1)−z(i,j,k)ta(i,j,k)=v2−v1tv1=\frac {z(i,j,k)-z(i,j,k-1)} t \\ v2=\frac {z(i,j,k+1)-z(i,j,k)} t \\ a(i,j,k)=\frac {v2-v1} t v1=tz(i,j,k)−z(i,j,k−1)v2=tz(i,j,k+1)−z(i,j,k)a(i,j,k)=tv2−v1
v1和v2表示速度,a(i,j,k)则表示加速度。

方便起见，我们令帧之间的间隔t为单位时间T=1，直接在推导中忽略这个参数，如果你带着t推导结果应该也一样，大概会麻烦一些。
文章指出，传播速度c的取值有一个范围，超出此范围的取值会使迭代方程发散。这个范围是：
0<c<d2tut+20<c<\frac d {2t}\sqrt{ut+2} 0<c<2tdut+2
经过我们的简化，则有
0<c2<d220<c^2<\frac {d^2} 2 0<c2<2d2
文章中做出了合理假设并严格推导出了这个范围。抛开推导的话，这个取值范围的限制你也可以类比理解为，当我们计算微积分（比如计算图形面积）时，微元dx需要取的足够小才能得到正确的结果。

虽然上式中是小于号，但后续处理中使用了一个衰减系数，所以此处我们可以取右侧的边界值，最终(2)式被我们化简为：
a(i,j,k)=12⟮z(i+1,j,k)+z(i−1,j,k)+z(i,j+1,k)+z(i,j−1,k)−4z(i,j,k)⟯(3)\tag 3a(i,j,k)=\frac 1 2\lgroup{z(i+1,j,k)+z(i-1,j,k)+z(i,j+1,k)+z(i,j-1,k)}\\-{4z(i,j,k)}\rgroup a(i,j,k)=21⟮z(i+1,j,k)+z(i−1,j,k)+z(i,j+1,k)+z(i,j−1,k)−4z(i,j,k)⟯(3)

观察等式右边，我们可以发现z的第三个坐标均是k，即在我们的假设条件下，仅由当前帧的水面高度数据即可算出此时的加速度,由此进一步我们可以算出下一帧的水面高度
v2=v1+a(i,j,k)∗T=v1+a(i,j,k)z(i,j,k+1)=z(i,j,k)+v2∗T=z(i,j,k)+v2v2=v1+a(i,j,k)*T =v1+a(i,j,k)\\ z(i,j,k+1)=z(i,j,k)+v2*T=z(i,j,k)+v2 v2=v1+a(i,j,k)∗T=v1+a(i,j,k)z(i,j,k+1)=z(i,j,k)+v2∗T=z(i,j,k)+v2
由上式可以看出，引入物理意义之后，我们只凭一帧的数据z(i,j,k)即可计算出下一帧z(i,j,k+1)。不过此时的水波是没有衰减的。如果在每次计算出v2后，乘以一个衰减值，就可实现衰减效果了。
v2′=v2∗dampingv2'=v2*damping v2′=v2∗damping
这个程序中最复杂的部分到此就解释完毕了，现在我们可以阅读着色器的代码并尝试将其与上方的推导对应上了。

顶点着色器：

static const char* globVert="attribute vec2 vertex;\n""varying vec2 coord;\n""void main() {\n""  coord = vertex * 0.5 + 0.5;\n"" gl_Position = vec4(vertex, 0.0, 1.0);\n""}\n";

同drop_program

片段着色器：

static const char* updateFrag="precision highp float;\n""uniform sampler2D texture;\n""uniform vec2 delta;\n"//采样间距，对应上文中的d，程序中表现为水波效果的分辨率"uniform float damping;\n"//衰减系数"varying vec2 coord;\n""void main() {\n""  vec4 info = texture2D(texture, coord);\n""   vec2 dx = vec2(delta.x, 0.0);\n""    vec2 dy = vec2(0.0, delta.y);\n""    float average = (\n""        texture2D(texture, coord - dx).r +\n""       texture2D(texture, coord - dy).r +\n""       texture2D(texture, coord + dx).r +\n""      texture2D(texture, coord + dy).r\n"" ) * 0.25;\n"//average-info.r=0.5*a(i,j,k)"   info.g += (average - info.r) * 2.0;\n"//v2=v1+a(i,j,k)"   info.g *= damping;\n"//v2'=v2*damping" info.r += info.g;\n"//z(i,j,k+1)=z(i,j,k)+v2'"  gl_FragColor = info;\n"//写入帧缓冲"}\n";

根据上方的分析和推导，帧缓冲的r分量存储的是水面高度数据，而g分量存储的则是速度数据。

render_program

顶点着色器：

static const char* renderVert="precision highp float;\n""attribute vec2 vertex;\n""varying vec2 ripplesCoord;\n"//这里我与原项目坐标转换的方式不同，此参数实际上属于弃用"varying vec2 backgroundCoord;\n""void main() {\n""  backgroundCoord=vertex*0.5+0.5;\n"" ripplesCoord=backgroundCoord;\n"//不过为了修改简便起见，保留了该参数并直接使用上面的值" gl_Position = vec4(vertex.x, vertex.y, 0.0, 1.0);\n""}\n";

片段着色器：

static const char* renderFrag="precision highp float;\n""uniform sampler2D samplerBackground;\n"//背景图片纹理"uniform sampler2D samplerRipples;\n"//帧缓冲纹理"uniform vec2 delta;\n"//采样间距，程序中表现为水波效果的分辨率"uniform float perturbance;\n""varying vec2 ripplesCoord;\n""varying vec2 backgroundCoord;\n""void main() {\n"" float height = texture2D(samplerRipples, ripplesCoord).r;\n""    float heightX = texture2D(samplerRipples, vec2(ripplesCoord.x + delta.x, ripplesCoord.y)).r;\n""    float heightY = texture2D(samplerRipples, vec2(ripplesCoord.x, ripplesCoord.y + delta.y)).r;\n""    vec3 dx = vec3(delta.x, heightX - height, 0.0);\n"//x轴方向切线向量"    vec3 dy = vec3(0.0, heightY - height, delta.y);\n"//y轴方向切线向量//这里原项目有点令人疑惑,不知为何把z分量作为第二个参数//下方这一行要取得xy方向上的偏移分量也因此需要改为xz"  vec2 offset = -normalize(cross(dy, dx)).xz;\n"//cross为叉乘，可以得到水面该点处的法线向量//下面这行假设固定方向上有一个光源并计算了反射光，水波的白色高亮即是由specular的值决定的"  float specular = pow(max(0.0, dot(offset, normalize(vec2(-0.6, 1.0)))), 4.0);\n""    gl_FragColor = texture2D(samplerBackground, backgroundCoord + offset * perturbance) + specular;\n"//这里使用近似值简化了计算"}\n";

这里解释一下这个着色器的最后一行

视线经过水面后会发生折射，本来应该偏移至offset’的位置，但此时我们认为，无论入射角如何，视线一旦经过水面，便会折射到水面的法线方向，在物理上则是认为此时视线经过的液体折射率为无穷大，这显然是不可能的。但是当入射角很小或者是水面到水底的高度很低时（即参数perturbance很小），offset’与offset之间的差值也很小，便可以直接使用offset作为近似值。

出于验证的想法，我也动用了自己贫瘠的物理和数学知识写了一个更精确的计算光线折射的着色器：

static const char* renderFrag="precision highp float;\n""uniform sampler2D samplerBackground;\n""uniform sampler2D samplerRipples;\n""uniform vec2 delta;\n""uniform float perturbance;\n""varying vec2 ripplesCoord;\n""varying vec2 backgroundCoord;\n""void main() {\n""    float height = texture2D(samplerRipples, ripplesCoord).r;\n""    float heightX = texture2D(samplerRipples, vec2(ripplesCoord.x + delta.x, ripplesCoord.y)).r;\n""    float heightY = texture2D(samplerRipples, vec2(ripplesCoord.x, ripplesCoord.y + delta.y)).r;\n""    vec3 dx = vec3(delta.x,  0.0,heightX - height);\n""  vec3 dy = vec3(0.0, delta.y,heightY - height);\n""   vec3 normal = normalize(cross(dy,dx));\n"//单位化的法线向量"   vec3 ref = refract(vec3(0.0,0.0,-1.0),-normal,1.0/1.333);\n"//opengl有内建函数refract可以计算经过折射的光线方向//1.333为水的折射率"   vec3 ofs = ref*(0.05+height)/dot(vec3(0.0,0.0,-1.0),ref);\n"//ofs的xy分量即为上图中更加精确的偏移offset'"  vec2 offset = normalize(cross(dy, dx)).xy;\n""   float specular = pow(max(0.0, dot(offset, normalize(vec2(-0.6, 1.0)))), 4.0);\n""    gl_FragColor = texture2D(samplerBackground, backgroundCoord + ofs.xy) + specular;\n""}\n";

实际运行程序，在水面高度较低时，二者表现基本无异，而当水面高度较高时，二者表现又都不是很理想，毕竟现实生活中水面一高就很难看清水底了，增大相应参数后产生的效果总感觉怪怪的。那么，在程序效果表现较好的范围内，原来的着色器中省去了若干计算的近似就显得非常的划算了。

后记

有些时候，我们不得不用更严谨的态度去深入的研究某项事物，这里也终于是通过查阅资料和自己思考对此项目的核心算法做了我认为比较详细的解释。理论研究与实际应用的侧重点不同，难点也不同，直接比较两者谁更难或是谁更有用并不合适，不过偶尔练习一下自己的理论分析能力和学习能力也是不错的。
一开始我就是想在Qt上实现一个表现较好的水波特效，然后用这种偏应用的关键词去搜索相关文章，倒是有一些代码可供参考，但是其中很少有在原理上做出详细解释的，这里也不得不说CSDN文章的推荐算法相关度真是蛮高的，后来在推荐的文章中看到的曙光，并受到启发使用了类似模拟、仿真、算法等偏向于理论的关键词去搜索，最终找到了想要的资料，获得了期望的答案。

如果想要了解更多的水面模拟知识，获得更加真实的水面效果，可能需要在计算机图形学领域深入的学习，最后再给出一篇相关文章的链接吧。
游戏中的实时水体模拟技术