渲染三維場景時經常會遇到需要渲染各種水體的情況，比如湖泊、河流、海洋等，不僅需要水體表面要有接近真實的隨時間而變化的波動，還要有令人信服的顏色、反射光、透明度等細節。實時渲染水面的方法有很多，從簡單的若干正弦波疊加，到《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個點中計算線性內插即可。
我們採用兩種波形。第一個波峰較尖，用來繪製細小的水波，第二個波峰較寬，用來繪製波長較長的水波。
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
};
線性內插函數：
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+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;
}
用一個結構體來保存時間和每個波的波長、振幅、方向、頻率和起始座標：
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值：
wave = 0.0; for(int w=0; w法線的計算
為了計算光照，需要知道每個頂點的法線方向。法線是根據網格上相鄰4點的座標計算的。
![水面的簡單渲染 – Gerstner波]()
由於每個點與四個網格面相鄰，每個面有獨立的法線，所以應該把4個面的法線的均值作為這個頂點的法線。所謂的“均值”就是把4個規範化的法線相加在除以4。CG中涉及法線的計算一般都需要在計算前規範化法線向量（使其長度為1），但因為每個網格面在x和y方向上的尺寸相同，所以法線不需要規範化就可以相加，只對加和後的法線做一次規範化即可。這樣就減少了四次規範化和一次除法。
規範化函數：
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-0.0000001){ t++; } if(t == count){ printf("normalizeF(): The input vector is too small.\n"); return 1; } for(t=0; t這個函數里需要注意的是，由於使用浮點數進行計算，在小數點後第8位開始容易產生誤差，如果輸入向量每個分量都很小，則很難保證計算結果的正確性；同時，如果輸入向量的長度很小，則需要先把輸入向量擴大10000倍，再進行計算，以減小誤差對計算結果的影響。這個函數還支持原地計算（即輸入和輸出為同一個向量），適用範圍是很廣的。
法線計算形式如下：
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（儲存在顯存中）。變換代碼如下：
for(int c=0; c對法線的處理也是如此。產生和綁定VAO的代碼如下：
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);繪製網格的代碼如下：
for(int c=0; c使用GLM計算模型視圖矩陣、透視投影矩陣和法線變換矩陣：
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行改成如下內容：valType const rad = glm::radians(fovy); #endif valType tanHalfFovy = tan(rad); detail::tmat4x4 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);vertex shader代碼如下：
#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; }繪製的網格效果如下：
![水面的簡單渲染 – Gerstner波]()
添加材質和光照
我使用了一張512*512大小的貼圖，格式為tga。
讀取和綁定貼圖：
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能同時讀到兩張貼圖：
glActiveTexture(GL_TEXTURE1); glBindTexture(GL_TEXTURE_2D, names.normal_texture); glActiveTexture(GL_TEXTURE0); glBindTexture(GL_TEXTURE_2D, names.diffuse_texture);添加材質後的效果：
![水面的簡單渲染 – Gerstner波]()
經過比較，我覺得Ward的各向異性光照效果比較好。fragment shader代碼如下：
#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; }添加光照後的效果：
