渲染三维场景时经常会遇到需要渲染各种水体的情况，比如湖泊、河流、海洋等，不仅需要水体表面要有接近真实的随时间而变化的波动，还要有令人信服的颜色、反射光、透明度等细节。实时渲染水面的方法有很多，从简单的若干正弦波叠加，到《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个点中计算线性内插即可。
两种波形。第一个波峰较尖，用来绘制细小的水波，第二个波峰较宽，用来绘制波长较长的水波。
|
1 2 3 4 5 6 7 8 9 10 |
static const GLfloat gerstner_pt_a[22] = { 0.0,0.0, 41.8,1.4, 77.5,5.2, 107.6,10.9, 132.4,17.7, 152.3,25.0, 167.9,32.4, 179.8,39.2, 188.6,44.8, 195.0,48.5, 200.0,50.0 }; static const GLfloat gerstner_pt_b[22] = { 0.0,0.0, 27.7,1.4, 52.9,5.2, 75.9,10.8, 97.2,17.6, 116.8,25.0, 135.1,32.4, 152.4,39.2, 168.8,44.8, 184.6,48.5, 200.0,50.0 }; |
线性内插函数：
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 |
static float gerstnerZ(float w_length, float w_height, float x_in, const GLfloat gerstner[22]) { x_in = x_in * 400.0 / w_length; while(x_in < 0.0) x_in += 400.0; while(x_in > 400.0) x_in -= 400.0; if(x_in > 200.0) x_in = 400.0 - x_in; int i = 0; float yScale = w_height/50.0; while(i<18 && (x_in<gerstner[i] || x_in>=gerstner[i+2])) i+=2; if(x_in == gerstner[i]) return gerstner[i+1] * yScale; if(x_in > gerstner[i]) return ((gerstner[i+3]-gerstner[i+1]) * (x_in-gerstner[i]) / (gerstner[i+2]-gerstner[i]) + gerstner[i+3]) * yScale; } |
用一个结构体来保存时间和每个波的波长、振幅、方向、频率和起始坐标：
|
1 2 3 4 5 6 7 8 |
static struct { GLfloat time; GLfloat wave_length[WAVE_COUNT], wave_height[WAVE_COUNT], wave_dir[WAVE_COUNT], wave_speed[WAVE_COUNT], wave_start[WAVE_COUNT*2]; } values; |
计算每个网格点的z值：
|
1 2 3 4 5 6 |
wave = 0.0; for(int w=0; w<WAVE_COUNT; w++){ d = (pt_strip[index] - values.wave_start[w*2] + (pt_strip[index+1] - values.wave_start[w*2+1]) * tan(values.wave_dir[w])) * cos(values.wave_dir[w]); wave += values.wave_height[w] - gerstnerZ(values.wave_length[w], values.wave_height[w], d + values.wave_speed[w] * values.time, gerstner_pt_a); } pt_strip[index+2] = START_Z + wave; |
法线的计算
为了计算光照，需要知道每个顶点的法线方向。法线是根据网格上相邻4点的坐标计算的。
由于每个点与四个网格面相邻，每个面有独立的法线，所以应该把4个面的法线的均值作为这个顶点的法线。所谓的“均值”就是把4个规范化的法线相加在除以4。CG中涉及法线的计算一般都需要在计算前规范化法线向量（使其长度为1），但因为每个网格面在x和y方向上的尺寸相同，所以法线不需要规范化就可以相加，只对加和后的法线做一次规范化即可。这样就减少了四次规范化和一次除法。
规范化函数：
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 |
static int normalizeF(float in[], float out[], int count) { int t=0; float l = 0.0; if(count <= 0.0){ printf("normalizeF(): Number of dimensions should be larger than zero.\n"); return 1; } while(t<count && in[t]<0.0000001 && in[t]>-0.0000001){ t++; } if(t == count){ printf("normalizeF(): The input vector is too small.\n"); return 1; } for(t=0; t<count; t++) l += in[t] * in[t]; if(l < 0.0000001){ l = 0.0; for(t=0; t<count; t++) in[t] *= 10000.0; for(t=0; t<count; t++) l += in[t] * in[t]; } l = sqrt(l); for(t=0; t<count; t++) out[t] /= l; return 0; } |
这个函数里需要注意的是，由于使用浮点数进行计算，在小数点后第8位开始容易产生误差，如果输入向量每个分量都很小，则很难保证计算结果的正确性；同时，如果输入向量的长度很小，则需要先把输入向量扩大10000倍，再进行计算，以减小误差对计算结果的影响。这个函数还支持原地计算（即输入和输出为同一个向量），适用范围是很广的。
法线计算形式如下：
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 |
int p0 = index-STRIP_LENGTH*3, p1 = index+3, p2 = index+STRIP_LENGTH*3, p3 = index-3; float xa, ya, za, xb, yb, zb; if(i > 0){ if(j > 0){ xa = pt_strip[p0] - pt_strip[index], ya = pt_strip[p0+1] - pt_strip[index+1], za = pt_strip[p0+2] - pt_strip[index+2]; xb = pt_strip[p3] - pt_strip[index], yb = pt_strip[p3+1] - pt_strip[index+1], zb = pt_strip[p3+2] - pt_strip[index+2]; pt_normal[index] += ya*zb-yb*za; pt_normal[index+1] += xb*za-xa*zb; pt_normal[index+2] += xa*yb-xb*ya; } if(j < STRIP_LENGTH-1){ xa = pt_strip[p1] - pt_strip[index], ya = pt_strip[p1+1] - pt_strip[index+1], za = pt_strip[p1+2] - pt_strip[index+2]; xb = pt_strip[p0] - pt_strip[index], yb = pt_strip[p0+1] - pt_strip[index+1], zb = pt_strip[p0+2] - pt_strip[index+2]; pt_normal[index] += ya*zb-yb*za; pt_normal[index+1] += xb*za-xa*zb; pt_normal[index+2] += xa*yb-xb*ya; } } if(i < STRIP_COUNT-1){ if(j > 0){ xa = pt_strip[p3] - pt_strip[index], ya = pt_strip[p3+1] - pt_strip[index+1], za = pt_strip[p3+2] - pt_strip[index+2]; xb = pt_strip[p2] - pt_strip[index], yb = pt_strip[p2+1] - pt_strip[index+1], zb = pt_strip[p2+2] - pt_strip[index+2]; pt_normal[index] += ya*zb-yb*za; pt_normal[index+1] += xb*za-xa*zb; pt_normal[index+2] += xa*yb-xb*ya; } if(j < STRIP_LENGTH-1){ xa = pt_strip[p2] - pt_strip[index], ya = pt_strip[p2+1] - pt_strip[index+1], za = pt_strip[p2+2] - pt_strip[index+2]; xb = pt_strip[p1] - pt_strip[index], yb = pt_strip[p1+1] - pt_strip[index+1], zb = pt_strip[p1+2] - pt_strip[index+2]; pt_normal[index] += ya*zb-yb*za; pt_normal[index+1] += xb*za-xa*zb; pt_normal[index+2] += xa*yb-xb*ya; } } if(normalizeF(&pt_normal[index], &pt_normal[index], 3)) printf("%d\t%d\n", index/3/STRIP_LENGTH, (index/3)%STRIP_LENGTH); index += 3; |
网格的绘制
最初编写代码时，我考虑过在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（储存在显存中）。变换代码如下：
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 |
for(int c=0; c<(STRIP_COUNT-1); c++) { for(int l=0; l<2*STRIP_LENGTH; l++) { if(l%2 == 1){ pt = c*STRIP_LENGTH + l/2; }else{ pt = c*STRIP_LENGTH + l/2 + STRIP_LENGTH; } index = STRIP_LENGTH*2*c+l; for(int i=0; i<3; i++){ vertex_data[index*3+i] = pt_strip[pt*3+i]; normal_data[index*3+i] = pt_normal[pt*3+i]; } } } |
对法线的处理也是如此。产生和绑定VAO的代码如下：
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 |
names.attributes.position = glGetAttribLocation(names.program, "position"); glGenBuffers(1, &names.vertex_buffer); names.attributes.normal = glGetAttribLocation(names.program, "normal"); glGenBuffers(1, &names.normal_buffer); glBindBuffer(GL_ARRAY_BUFFER, names.vertex_buffer); glBufferData(GL_ARRAY_BUFFER, sizeof(vertex_data), vertex_data, GL_STATIC_DRAW); glVertexAttribPointer(names.attributes.position, 3, GL_FLOAT, GL_FALSE, sizeof(GLfloat)*3, (void*)0); glEnableVertexAttribArray(names.attributes.position); glBindBuffer(GL_ARRAY_BUFFER, names.normal_buffer); glBufferData(GL_ARRAY_BUFFER, sizeof(normal_data), normal_data, GL_STATIC_DRAW); glVertexAttribPointer(names.attributes.normal, 3, GL_FLOAT, GL_FALSE, sizeof(GLfloat)*3, (void*)0); glEnableVertexAttribArray(names.attributes.normal); |
绘制网格的代码如下：
|
1 2 |
for(int c=0; c<(STRIP_COUNT-1); c++) glDrawArrays(GL_TRIANGLE_STRIP, STRIP_LENGTH*2*c, STRIP_LENGTH*2); |
|
1 2 3 4 5 6 7 8 |
glm::mat4 Projection = glm::perspective(45.0f, (float)(SCREEN_WIDTH/SCREEN_HEIGHT), 1.0f, 100.f); glm::mat4 viewTransMat = glm::translate(glm::mat4(1.0f), glm::vec3(0.0f, 0.0f, -2.5f)); glm::mat4 viewRotateMat = glm::rotate(viewTransMat, -45.0f, glm::vec3(1.0f, 0.0f, 0.0f)); glm::mat4 ModelViewMat = glm::scale(viewRotateMat, glm::vec3(1.0f, 1.0f, 1.0f)); glm::mat3 NormalMat = glm::transpose(glm::inverse(glm::mat3(ModelViewMat))); glUniformMatrix4fv(glGetUniformLocation(names.program, "modelViewMat"), 1, GL_FALSE, glm::value_ptr(ModelViewMat)); glUniformMatrix4fv(glGetUniformLocation(names.program, "perspProjMat"), 1, GL_FALSE, glm::value_ptr(Projection)); glUniformMatrix3fv(glGetUniformLocation(names.program, "normalMat"), 1, GL_FALSE, glm::value_ptr(NormalMat)); |
为了使glm::perspective()的参数符合自己的习惯，我把glm/gtc/matrix_transform.inl的第254到264行改成如下内容：
|
1 2 3 4 5 6 7 8 9 10 11 |
valType const rad = glm::radians(fovy); #endif valType tanHalfFovy = tan(rad); detail::tmat4x4<valType, defaultp> Result(valType(0)); Result[0][0] = valType(1) / (tanHalfFovy); Result[1][1] = (aspect) / (tanHalfFovy); Result[2][2] = - (zFar + zNear) / (zFar - zNear); Result[2][3] = - (valType(2) * zFar * zNear) / (zFar - zNear); Result[3][2] = - valType(1); |
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 |
#version 130 attribute vec3 position; attribute vec3 normal; uniform mat4 modelViewMat; uniform mat4 perspProjMat; uniform mat3 normalMat; uniform float time; varying vec2 texture_coord; varying vec3 normalVect; varying vec3 lightVect; varying vec3 eyeVect; varying vec3 halfWayVect; varying vec3 reflectVect; void main() { gl_Position = perspProjMat * modelViewMat * vec4(position, 1.0); float tex_x = (position.x + time/20.0) / 8.0 + 0.5; float tex_y = 0.5 - (position.y + time/25.0) / 5.0; texture_coord = vec2(tex_x, tex_y); vec3 eyePos = vec3(0.0, 0.0, 5.0); vec3 lightPos = vec3(1.0, 3.0, 0.0); vec3 ptVertex = vec3(modelViewMat * vec4(position, 1.0)); eyeVect = normalize(eyePos - ptVertex); lightVect = normalize(lightPos - ptVertex); halfWayVect = eyeVect + lightVect; normalVect = normalMat * normal; reflectVect = 1.0 * eyeVect - 2.0 * dot(-1.0*eyeVect, normalVect) * normalVect; } |
绘制的网格效果如下：
添加材质和光照
我使用了一张512*512大小的贴图，格式为tga，图像可以在文章末尾所给的链接获取。
读取和绑定贴图：
|
1 2 3 4 5 6 7 |
names.diffuse_texture = initTexture("water-texture-2.tga"); names.uniforms.diffuse_texture = glGetUniformLocation(names.program, "textures[0]"); glUniform1i(names.uniforms.diffuse_texture, 0); names.normal_texture = initTexture("water-texture-2-normal.tga"); names.uniforms.normal_texture = glGetUniformLocation(names.program, "textures[1]"); glUniform1i(names.uniforms.normal_texture, 1); |
initTexture()函数用来读取贴图文件和生成贴图对象，这里就不详细介绍了，具体细节可以参考文章末尾给出的代码链接。为了试验法线贴图，我使用了两张贴图，在每次绘制时要交替启用两张贴图，以保证fragment shader能同时读到两张贴图：
|
1 2 3 4 5 |
glActiveTexture(GL_TEXTURE1); glBindTexture(GL_TEXTURE_2D, names.normal_texture); glActiveTexture(GL_TEXTURE0); glBindTexture(GL_TEXTURE_2D, names.diffuse_texture); |
添加材质后的效果：
经过比较，我觉得Ward的各向异性光照效果比较好，其原理可以参考：D3DBook:(Lighting) Ward – GPWiki。
fragment shader代码如下：
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 |
#version 130 uniform sampler2D textures[2]; uniform vec4 materAmbient, materSpecular; uniform vec4 lightDiffuse, lightAmbient, lightSpecular; uniform vec4 envirAmbient; varying vec2 texture_coord; varying vec3 normalVect; varying vec3 lightVect; varying vec3 eyeVect; varying vec3 halfWayVect; varying vec3 reflectVect; void main() { vec4 diffuse, ambient, globalAmt; vec4 specular; vec3 eyeDir, lightDir, normalDir, halfWayDir, reflectDir; float NdotL, NdotH, NdotR, S, temp, delta; float alpha = 0.4; eyeDir = normalize(eyeVect); lightDir = normalize(lightVect); normalDir = normalize(normalVect); halfWayDir = normalize(halfWayVect); reflectDir = normalize(reflectVect); NdotL = max(dot(normalDir, lightDir), 0.0); NdotH = max(dot(normalDir, halfWayDir), 0.0); NdotR = max(dot(normalDir, reflectDir), 0.0); delta = acos(NdotH); temp = -1.0 * tan(delta) * tan(delta) / alpha / alpha; S = pow(2.71828, temp) / 4.0 / 3.14159 / alpha / alpha / pow(NdotL*NdotR, 0.5); diffuse = texture2D(textures[0], texture_coord) * lightDiffuse; globalAmt = envirAmbient * materAmbient; ambient = envirAmbient * lightAmbient; specular = materSpecular * lightSpecular; gl_FragColor = NdotL * (diffuse + specular * S) + globalAmt; } |
添加光照后的效果：
全部代码和贴图可以在这里查看：
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）。
顶点法向量的计算方法很巧妙啊！
顺便问一下，GPU Gems http://http.developer.nvidia.com/GPUGems/gpugems_ch01.html 上的公式 (12) 是怎么求的？
我按照他说的 N = B x T, 却怎么也消不去求和的那些项
你好，也就是说你这个只是改变了顶点的高度值，并没有改变顶点在水平面的坐标？
是的，在计算高度的过程中已经考虑了坐标点在水平方向的平移。所以这种很均匀的网格只适合用来模拟比较平缓的水面。
博主，能说说gerstnerZ和calcuWave方法具体实现原理吗？我现在想用分形噪声和海浪谱来模拟海面，您有没有做过这块的研究，跪求指导和建议！
因为我的代码中没有实时计算Gestner波形，只是预先保存了两个参数不同的半波（坐标范围(0,0)-(200,50)），然后用gerstnerZ来缩放和内插，适应不同的波长和幅度。而calcuWave里三个循环块的作用分别是在每个方格点叠加所有波的Z值；为每个方格点计算法向；和把密集保存在数组中的点重新排列成一条条的三角形带，便于OpenGL直接调用绘制多边形。
由于目前CG并不是我的研究重点，关于水面模拟的其他技术也不是十分了解，所以恐怕很难帮得上你：）
您好，请问参数d和x_in分别表示什么？怎么由插值函数控制波动的？
d就是公式z=g(x,y)下面图中的d，d的方向代表单个波的传播方向。模板波是二维的（x_in和z），实际的波是三维的（x，y和z）。首先根据三维空间水平网格上x，y坐标值推算出每个波各自传播方向上的d（也就是图中的d），作为用于求该波高度的横坐标。d+波速*时间就是当前时刻的x_in，通过对x_in插值求出z，这个z就是该波在三维的网格空间中的高度。随着时间逐渐增加，x_in也逐渐变化，波动就产生了。
百度网盘下载下来，配置好了，编译不过去是不是有啥要修改？
编译不通过的错误信息是怎样的呢？
您好，求波高的公式里面的d是如何计算的，根据您的代码我看的不是很懂，希望您帮忙解答一下，不胜感激
wave = 0.0;for(intw=0; w<WAVE_COUNT; w++){d = (pt_strip[index] - values.wave_start[w*2] + (pt_strip[index+1] - values.wave_start[w*2+1]) *tan(values.wave_dir[w])) *cos(values.wave_dir[w]);wave += values.wave_height[w] - gerstnerZ(values.wave_length[w], values.wave_height[w], d + values.wave_speed[w] * values.time, gerstner_pt_a);}pt_strip[index+2] = START_Z + wave;d其实就是根据z=g(x,y)那个公式计算的，只是把sin拆成cos*tan，然后提出cos。这个公式是为了把任意方向波面的x坐标转换到水平方向的公共坐标系内。
你好，我现在有波高，波向，具体点的经纬度，请问该如何模拟近岸海浪，该片区域大概1w个点
我不太确定用Gerstner波模拟近岸海浪能否达到你的要求，如果你觉得形状可以接受的话，已知波高，波向和点的经纬度，就可以直接根据坐标来绘制了。如果点的经纬坐标不规则，关键要找到一个方法把散点用多边形面表示出来，最基本的就属Delaunay三角剖分了吧。
我们的游戏里面要用自己做一个水面，看了这篇文章好几次了，深有学习。
非常高兴我的文章能够对你有所帮助！
相当有帮助啊！可以跟你用qq交流请教下不？我的是：43008117.
前辈，可以帮我看看我的程序应该怎么用gerstner波么？我不太看懂这个波和您的处理方法。
void Initmesh()
{
//创建海浪二维网格
float Lx = 0.0f;
float Lz;
for (int x = 0; x < SIZEX; x++)
{
Lz = 0.0;
//海浪网格
for (int z = 0; z < (SIZEZ/2 + STARTZ); z++)
{
//网格中心点为原点，分布在XZ平面，Y值均设为0
mesh[x][z][0] = ((float) SIZEX*0.5) – Lx;
mesh[x][z][1] = 0.0f;
mesh[x][z][2] = (float) STARTZ – Lz;
Lz ++;
}
Lx ++;
}
}
网格是XZ平面上SIZEX（150）*SIZEZ（150）的 四个角的坐标分别是（-75x,20z,0y)(75x,20z,0y)(-75x,-75z,0y)(75x,-75z,0y) 每个点的坐标都按XZY 的顺序 存在mesh中
现在每个点Y都是0 要怎么处理才能得到Y的gerstner的值
如果是沿着x轴绘制的话，可以这样计算y坐标值：
mesh[x][z][1] = A*cos( (double)(Lx – w*t) );
再计算偏移后的x坐标值：
mesh[x][z][0] = Lx – Q*Di*A*sin( (double)(Di*Lx – w*t) );
其中Q,Di,A,w,t的含义与公式里的对应。需要注意用均匀的矩形网格坐标计算x和y，计算得到的x就不再均匀了。
谢谢前辈。
还有
static const GLfloat wave_para[6][6] = {
{ 1.6, 0.12, 0.9, 0.06, 0.0, 0.0 },
{ 1.3, 0.1, 1.14, 0.09, 0.0, 0.0 },
{ 0.2, 0.01, 0.8, 0.08, 0.0, 0.0 },
{ 0.18, 0.008, 1.05, 0.1, 0.0, 0.0 },
{ 0.23, 0.005, 1.15, 0.09, 0.0, 0.0 },
{ 0.12, 0.003, 0.97, 0.14, 0.0, 0.0 }
};
static const GLfloat gerstner_pt_a[22] = {
0.0,0.0, 41.8,1.4, 77.5,5.2, 107.6,10.9,
132.4,17.7, 152.3,25.0, 167.9,32.4, 179.8,39.2,
188.6,44.8, 195.0,48.5, 200.0,50.0
};
static const GLfloat gerstner_pt_b[22] = {
0.0,0.0, 27.7,1.4, 52.9,5.2, 75.9,10.8,
97.2,17.6, 116.8,25.0, 135.1,32.4, 152.4,39.2,
168.8,44.8, 184.6,48.5, 200.0,50.0
};
您这三块数字有什么说法没？1.6,0.12,0.9,0.6 从哪里来的
坐标范围为什么是（0,0）——（200,50）
麻烦前辈说一下，谢谢
不要这么客气，，我只是个学生
因为我是用6个不同的波合成，wave_para用来存放每个波的参数，依次为波长、振幅、方向、速度和起始坐标。每个参数是多次尝试确定的，只是为了整体效果更好。
这6个波实际上只有两种波形，分别存放在gerstner_pt_a和gerstner_pt_b中，相当于两个模板，是提前计算好的。(200,50)是任选的，(0,0)到(200,50)刚好是半个波形。
绘制新波形时选出一个模板，横向和纵向缩放适应波长和振幅。渲染时网格点的坐标由模板插值得到，不需要用三角函数实时计算Gerstner波，这样也是为了节省计算量。
谢谢前辈，我正在做本科毕设，刚接触图形编程不久，您的代码真实帮了大忙，用到我的程序里，很好的实现gerstner的变换，谢谢您。
HI,你好，我是在开发海体的时候，搜索Gerstner波的时候搜索到你的文章的，非常感兴趣，而且也看到了你的一些素描的东西，更加佩服了，希望可以交个朋友.
当然，能结交志同道合的朋友一直是我的荣幸！
您好，想请问一下，二维Gerstner波表示为y=f(x)的话，这个f(x)方程应该怎么列呢……
Gerstner波一般都用参数方程的形式（x=x(p),y=y(p)）来表示。三角函数的反函数分段表示本身就比较繁琐，如果想用一行公式来表示，也只能以[k*PI,(k+1)*PI)和[(k+1)*PI,(k+2)*PI)分段（当然还要考虑参数D,w,t）从方程y=y(p)解出p=p(y)，再代入x=x(p)，得到x=x(p(y))=q(y)的形式。q(y)的形式一定很复杂，所以不确定能否求得它的反函数以表示为y=f(x)的形式。
楼主真心好人一枚，对于我们这些新手来说真是很好的资源，非常非常感谢！！
抱歉回复的很晚。谢谢你的支持！
博主真心好人一枚
抱歉回复的很晚。谢谢支持！
博主你好，最近我也在研究opengl模拟水面，发现网上资料很少，我又看不懂《GPU GEMS》,直到发现了你这篇博文，不过我将你的github上的源码下下来，发现源码不是很完整，于是我添加了头文件，链接库的声明，也修改了matrix_transform.inl，配好了环境，终于没有报错了，不过运行出来还是黑屏。所以想问博主能给我发个整个工程的打包文件吗？谢谢！
你好，我几天前发过去的邮件被谷歌服务器退回了，给你一个网盘的链接吧：
http://pan.baidu.com/s/1eQnBv2M我这也是，黑的
下载链接已经更新：http://pan.baidu.com/s/1dDFgmB3
而且在GitHub目录里的代码一般是最新的：https://github.com/johnhany/OpenGLProjects/tree/master/GerstnerWave