ゲルストナー波(Gerstner Wave)

ゲルストナー波(Gerstner Wave)

※2019/06/16 各シェーダを修正しました。

 海洋の荒い波は単純なsin波とは異なり、波の山は鋭く、谷はより広く平になります。このような波はゲルストナー波を用いることで表現することができます。以下にGPU Gems:Chapter 1. Effective Water Simulation from Physical Modelsに掲載されているゲルストナー波の式を示します。3次元における表面の座標位置\(P(x,y)\)は以下の通りです。

\( \boldsymbol{P}(x,y,t)= \left[ \begin{align} &x+\sum(Q_{i}A_{i}\times\boldsymbol{D_{i}}.x\times cos(w_{i}\boldsymbol{D_{i}}\cdot (x,y)+\varphi_{i}t)),\\ &y+\sum(Q_{i}A_{i}\times\boldsymbol{D_{i}}.y\times cos(w_{i}\boldsymbol{D_{i}}\cdot (x,y)+\varphi_{i}t)),\\ &\sum(A_{i}sin(w_{i}\boldsymbol{D_{i}}\cdot(x,y)+\varphi_{i}t)) \end{align} \right] \)

\(Q_{i}\)は波の鋭さ、\(A_{i}\)は振幅、 \(\boldsymbol{D}\)は波の進行方向、\(w\)は波数、\(\varphi\)は角周波数です。また、従法線ベクトル\(\boldsymbol B\)及び接ベクトル\( \boldsymbol T\)はそれぞれ以下の式で表せます。

\(\boldsymbol{B}= \left[ \begin{align} &1-\sum(Q_{i}\times \boldsymbol{D_{i}}.x^{2}\times WA \times S),\\ &-\sum(Q_{i}\times \boldsymbol{D_{i}}.x\times \boldsymbol{D_{i}}.y\times WA \times S),\\ &\sum(\boldsymbol{D_{i}}.x\times WA \times C) \end{align} \right] \)

\(\boldsymbol{T}= \left[ \begin{align} &-\sum(Q_{i}\times \boldsymbol{D_{i}}.x\times \boldsymbol{D_{i}}.y\times WA \times S),\\ &1-\sum(Q_{i}\times \boldsymbol{D_{i}}.y^{2}\times WA \times S),\\ &\sum(\boldsymbol{D_{i}}.y\times WA \times C) \end{align} \right] \)

ここで、\(WA\)、\(S\)及び\(C\)は以下の通りです。

\( \begin{align} &WA=w_{i}\times A_{i}\\ &S=sin(w_{i}\boldsymbol{D_{i}}\cdot \boldsymbol{P}+\varphi_{i}t))\\ &C=cos(w_{i}\boldsymbol{D_{i}}\cdot \boldsymbol{P}+\varphi_{i}t)) \end{align} \)

これらの式を用いてシェーダにより海洋の波を表現します。

shder

 頂点シェーダ内で座標位置、従法線ベクトル及び接ベクトルを計算しています。求めた座標位置により平面の各頂点を移動しています。さらに、従法線ベクトルと接ベクトルの外積を計算し、法線を求めています。フラブメントシェーダ内ではライティングの計算及び当ブログに掲載したグリッドシェーダによってグリッド状の線を表示する処理を行っています。

GPU Gemsに掲載されている式とシェーダではy軸とz軸が逆となります。そのため、座標位置、従法線ベクトル及び接ベクトルのy成分とz成分を入れ替えて計算しています。また、法線ベクトル\(\boldsymbol{N}\)は、以下の式に示すように従法線ベクトル \(\boldsymbol{B}\)と接ベクトル\(\boldsymbol{N}\)の外積を計算することで求まります。

$$ \boldsymbol{N}=\boldsymbol{B}\times \boldsymbol{T}= \begin{pmatrix} b_{1}\\ b_{2}\\ b_{3} \end{pmatrix} \times \begin{pmatrix} t_{1}\\ t_{2}\\ t_{3} \end{pmatrix} =\begin{pmatrix} b_{2}t_{3}-b_{3}t_{2}\\ b_{3}t_{1}-b_{1}t_{3}\\ b_{1}t_{2}-b_{2}t_{1} \end{pmatrix} $$

シェーダではy成分とz成分が逆となるため法線ベクトル\(\boldsymbol{N’}\)は以下の式となります。

$$ \boldsymbol{N’}= \begin{pmatrix} b_{2}t_{3}-b_{3}t_{2}\\ b_{1}t_{2}-b_{2}t_{1}\\ b_{3}t_{1}-b_{1}t_{3} \end{pmatrix} $$

この法線ベクトル \(\boldsymbol{N’}\)をy成分とz成分を入れ替えた従法線ベクトル\(\boldsymbol{B’}\)と接ベクトル \(\boldsymbol{T’}\) を用いて計算します。式は以下の通りです。

$$ \boldsymbol{N’}=\boldsymbol{T’}\times \boldsymbol{B’}= \begin{pmatrix} t_{1}\\ t_{3}\\ t_{2} \end{pmatrix} \times \begin{pmatrix} b_{1}\\ b_{3}\\ b_{2} \end{pmatrix} =\begin{pmatrix} b_{2}t_{3}-b_{3}t_{2}\\ b_{1}t_{2}-b_{2}t_{1}\\ b_{3}t_{1}-b_{1}t_{3} \end{pmatrix} $$

このように、従法線ベクトルと接ベクトルを逆にして外積を計算することでy成分とz成分を入れ替えた法線ベクトルを計算することができます。

実行結果

 上記シェーダの実行結果は以下の通りです。sin波より鋭い山となっていることがわかります。

以下のように数値を変更すると

となります。メッシュのエッジ部分が明るくなっており、きれいにライティングができていないことがわかります。

ベクトル計算を変更

 ライティングをきれいに表示するため、各種ベクトルを頂点シェーダではなくフラグメントシェーダで計算するように変更しました。

shader

 従法線ベクトルと接ベクトルをフラグメントシェーダ内で計算し、法線ベクトルを求めています。

実行結果

 上記シェーダの実行結果は以下の通りです。エッジ部分が明るく表示されずに、滑らかにライティングができていることがわかります。

ワールド座標で計算

 以上のシェーダはローカル座標で座標位置の計算を行っていました。ローカル座標で計算した場合、複数の平面を並べ、カメラからの距離に応じてメッシュを変更する、つまり、LODを用いることができません。そこで、ワールド座標で計算するようにシェーダを変更しました。

shader

 ゲルストナー波をワールド座標で計算するシェーダは以下の通りです。

実行結果

 上記シェーダの実行結果は以下の通りです。ワールド座標で計算しているため、平面を動かしてもそれに伴って波が移動することはありません。

参考サイト

GPU Gems:Chapter 1. Effective Water Simulation from Physical Models