本文会对一个基于波动方程来模拟水面实现的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使用的算法的确是基于波动方程的,但是在上述文章的基础上进一步做了近似和简化,下面会给出一部分推导过程,请结合上方的文章学习。
关键的不同之处是对于水波衰减的处理。上方文章中波动方程的无阻尼形式是下方推导的基础。
那么首先来看看二维波动方程的无阻尼形式:
∂ 2 z ∂ t 2 = c 2 ⟮ ∂ 2 z ∂ x 2 + ∂ 2 z ∂ y 2 ⟯ (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∂