参考実装の動画：
固定端バージョン → YouTube https://youtu.be/XxfNnbH-fhI
自由端バージョン → YouTube https://youtu.be/ghV9T7o9abc

Qiita も YouTube の動画を埋め込められるようにして欲しい…。

はじめに

　この文章は NVIDIA からダウンロードできる Vertex Texture Fetch Water (PDF ファイル)に関する簡易的な解説文書です（水面の形状部分のみ解説）。対象とする PDF は 2004 年に公開されたホワイトペーパーですが、計算量が少なく、理論も比較的簡単な方法にて水面形状を生成するものであり、2019 年の今現在読んでみても結構面白いのではないかと思います。

　このホワイトペーパーは、ちょっと不思議な書き方をしていて、数式の機械的な展開は非常に丁寧に行っているにも関わらず、なぜそのようになるのか、といった部分の説明は結構省略されています。親切なんだか不親切なんだか、よく分からない文章です。自分が読んでみて、ここはもう少し説明が必要では？という部分を中心に書いてみました。

　参考として、このホワイトペーパーを参考にして 1 次元の水面（？）シミュレータを実装してみました。こちらは実際に動作するサンプルプログラムを掲載しています（動作の様子は YouTube https://youtu.be/XxfNnbH-fhI にて確認してみて下さい）。

ちなみに、重要な部分だけ知りたいという人は、式 5 の 2 行目から 3 行目への展開が間違っていて、正しくは、

p(t_2)=(1+\alpha)\cdot p(t_1)-\alpha\cdot p(t_0)+\frac{1}{2}a(t_1)h^2

となります（$\alpha\cdot p(t_0)$ の前の符号がプラスではなくマイナスでないと、正しくシミュレーションできません）。

Vertex Texture Fetch Water の座標系

　読んでいけば分かりますが、y 軸が水面の高さを表す軸です。x と z は水面の位置です。なので、もし時刻 $t$ における波の高さを表す関数が $W_t$ と表されるならば、$y=W_t(x,z)$ という形で時刻 $t$ における水面の形状は表現されます。

式 1 の意味

　式 1 として掲載されているこの式ですが、

\frac{\partial^2y}{\partial t^2}=c^2\left(
\frac{\partial^2y}{\partial x^2}+\frac{\partial^2y}{\partial z^2}
\right)

これは、周辺の高さから求められる加速度です。位置 $y$ を時間 $t$ で微分すると速度となり、速度を $t$ で微分すると加速度になります。なので、$y$ を $t$ で 2 階偏微分したものは加速度に他なりません。

この式は後ほどまた出てくるので、頭の片隅にでも覚えておく必要があります。

式 2 の意味

　式 2 が意味するところは、プログラム的には直接は何も言っていません。ここで言いたいのは、上で述べたことと逆で、加速度を 2 階積分すると位置になる、ということです。この伏線は次の式 3 につながっていきます。

式 3 における関数 a 関連

　式 3 では、加速度を表す関数 $a$ が導入され、式 3 中で $\frac{1}{2}a(t_1)(t_2-t_1)^2$ という項を形成しています。この部分が式 2 を紹介した伏線に関連する部分です。

　高校で物理を習った方は、重力定数 $g$ に対して、$\frac{1}{2}gt^2$ という公式を覚えているかもしれません。これは加速度 $g$ の下で時間 $t$ の間に移動した移動量を表すものでした。初速 0 で加速度 $g$ の下で $t$ 秒間運動すると、$t$ で $g$ を積分して速度の関数 $gt$ を得ます。これをさらに $t$ で積分して、位置の関数 $\frac{1}{2}gt^2$ を得る、という訳です。

さて、先程の式にもどって考えてみると、全く同じ形であることが分かります。式 1 によって、周囲の状況から決まる加速度を $a$ という関数で表した場合、時刻 $t_1$ から $t_2$ までの変位量（＝位置の変化）は $a$ に関する部分は $\frac{1}{2}a(t_1)(t_2-t_1)^2$ となるわです。

式 3 の構成は、この周囲によって決まる加速度のみならず、今着目している点の速度 $v(t_1)$ から変位量についても計算しています。なので、

新しい y の位置＝今の y の位置 + 現在の速度からの変位量(=v(t_1)(t_2-t_1)) + 周囲の状況により決定する加速度からの変位量（=\frac{1}{2}a(t_1)(t_2-t_1)^2）

と読むことができます。

式 5

　式 4 はベレの式の紹介で、「ほら、似てるでしょ？」的な意味合いなので、実装作業に関しては意味がありません。しかし、式 5 はそうではありません。プログラムの中心的な役割を果たす、非常に重要な式です。このホワイトペーパーでは、困ったことに、この非常に重要な式 5 に 間違いがあります。

　式 5 の 2 行目から 3 行目への展開が間違っていて、正しくは、

p(t_2)=(1+\alpha)\cdot p(t_1)-\alpha\cdot p(t_0)+\frac{1}{2}a(t_1)h^2

となります。$\alpha\cdot p(t_0)$ の前の符号が プラスではなくマイナスです。正しい式でないと、正しくシミュレーションできず悩むこととなります（私も少し悩みました）。

基本的には、この正しい式 5 を使えばシミュレータは実装できるのですが、周辺の状況から決定される加速度 $a$ の値を確定していかなければなりません。ホワイトペーパーでも、その話に入っていきます。

式 6

　周辺の状況による加速度は式 1 で求めることができるわけですが、そもそも、水面の形を形成する関数も分からないので、微分値を求めることができません。そこで、このホワイトペーパーでは、3 次関数として近似するものとして話をすすめていこう、という説明です。

　ここは非常に重要な部分なので、補足説明しますが、3 次関数を使って水面を近似するわけではありません。3 次関数を使って近似したとするとという考え方であるので、その辺りを勘違いすると訳が分からなくなるかと思います。

式 7

　式 7 は、このホワイトペーパーにおいて最も訳の分からない部分のハイライトとなっています。でも、ここで言っていることは重要であり、また、実はシンプルです。

　このホワイトペーパーでは $X$,$Z$ 平面に広がる水面を取り扱ってきているのに、いきなり 1 変数関数 $f$ が出てきて面食らったかと思います。これは $X$ 軸方向や $Z$ 軸方向など、ある方向に沿って 3 次関数 $f$ で近似したとすると、という文脈で説明されているからです。

言ってしまえば、式 1 における $\frac{\partial^2 y}{\partial x^2}$ や$\frac{\partial^2 y}{\partial z^2}$ をそれぞれ計算する方法を考えましょう、ということです。

ある点 $(x_0,z_0)$ における x 軸方向に関する水面の形状がある関数 $f$ を用いて、$y=f(x)$ として近似されているならば、 $\frac{\partial^2 y}{\partial x^2}$ は $f''(x_0)$ として近似されることとなります。

それを踏まえて、じゃあ、その $f''(x_0)$ ってどうやって計算すりゃ良いのよ？という話になります。

ゲームでは水面をグリッド（格子）で構成されたポリゴンで表現されることが多い状況です。そのため、今着目している点 $(x_0,z_0)$ における 2 重偏微分値である $f''(x_0)$ などを、これらグリッド上の値を使って、低コストで計算できると非常に嬉しい。それが可能であることをここでは説明しています。

今、着目している点 $x_0$ を原点として座標軸を取り直し、3 次関数 $f(x)=c_0x^3+c_1x^2+c_2x+c_3$ で近似されたとすると、$f''(0)$ はもちろん、

f''(0)=6c_0x+2c_1=2c_1

となります。

もちろん、この $c_1$ という係数は未知の定数です。これをどうやって、未知の定数のまま処理しつつ、具体的な値として取り出す方法を考えなければなりません。

この説明では、今は $x$ 軸方向の近似関数として $f$ を考えていました。「近似できているものとする」訳ですから、隣のグリッドの位置を $(x_0-1,z_0)$ および $(x_0+1,z_0)$ とするならば、これらのグリッドにおける今現在の $y$ の値も当然この関数 $f$ によって近似されているはずです。

その前提の下に、計算を行ってみると、

(x_0-1,z_0) での y の値 = これを f_{-1} とする = f(-1)=-c_0+c_1-c_2+c_3

同様にして、

(x_0+1,z_0) での y の値 = これを f_{+1} とする = f(+1)=c_0+c_1+c_2+c_3

を得ます。

着目している点の値も既知ですが、計算してみると、

(x_0,z_0) での y の値 = これを f_0 とする = f(0)=c_3

となります。

これらの材料を上手く使って、具体的な近似式（つまり、具体的な $c_0,c_1,c_2,c_3$ ）を計算することなく、$f''(0)=2c_1$ を計算することを考えます。

$f_{-1}+f_{+1}$ より、$2c_1+2c_3$ を得ます（実際に足し算するだけです）。

そして、今 $f_0=c_3$ だったので、結局、

f''(0)=f_{-1}+f_{+1}-2f_0=2c_1

を得ます。つまり、具体的な近似値を一切計算することなく、周辺からの加速度の元となる $f$ の 2 階微分値が「両隣の値の和-2×注目地点の値」で具体的に求めることができますよ、ということです。

この不思議な式 7 の意味するところはこういう意味です。

式 8

　式 1 では $X$ 軸方向、$Z$ 軸方向、それぞれの方向の 2 階偏微分を足し合わせたものが周辺から定まる加速度であることに他ならないので、それを式 7 の結果を用いて計算しているだけです。

$y_{-1,0}$ とかありますが、$y_{k,l}$ とすると、着目点を $(x_0,z_0)$ とするならば、$y_{k,l}=点(x_0+k,z_0+l)における y の値$ という意味です。

今回は、水面の形状を生成する部分に着目した解説文書ですので、解説部分はここで終わりです。

参考実装

　平面ではなく、例えば $X$ 軸方向のみを考慮した水平面（水平線？）の生成プログラムを以下に示しておきます（全 44 行）。マウスで水面の高さに影響できる、インタラクティブなものです。コードは Processing にて記述しております。

固定端バージョン

　まずは、実装がシンプルな固定端バージョンから示します（より水の挙動に近い自由端バージョンはこの後示します）：

float[] gPos=new float[51];
float[] gPosNew=new float[51];
float[] gPosOld=new float[51];

void setup() {
  size(500,100);
}

void draw() {
  if(mousePressed==true){
    gPos[mouseX/10]=(50-mouseY)/1.0f;
  }
  update();
  background(255,255,255);
  strokeWeight(1);
  stroke(0,0,0);
  for(int i=0; i<50; i++) {
    line(i*10,-gPos[i]+50,(i+1)*10,-gPos[i+1]+50);
  }
  fill(0,64,255);
  noStroke();
  for(int i=0; i<50; i++) {
    beginShape();
      vertex(i*10,-gPos[i]+50+1);
      vertex(i*10,100);
      vertex((i+1)*10,100);
      vertex((i+1)*10,-gPos[i+1]+50+1);
    endShape(CLOSE);
  }
}

void update() {
  float h=1.0f;
  float d=0.99f;
  float c2=0.8f;
  for(int i=1; i<50; i++) {
    float a=c2*(gPos[i-1]+gPos[i+1]-2*gPos[i]);
    gPosNew[i]=constrain((1+d)*gPos[i]-d*gPosOld[i]+0.5f*a*h*h,-40,+40);
  }
  for(int i=0; i<51; i++) {
    gPosOld[i]=gPos[i];
    gPos[i]=gPosNew[i];
  }
}

水面なので、本当は両端は自由端であるべきですが、実装が面倒だったので固定端にしております。
両端をつなげたトーラスにするのが楽かなーとか思いつつ実装してません。興味のある人はやってみると良いかと思います。

自由端バージョン（2019.02.12 追記）

　せっかくなので自由端バージョンも作成してみました。トーラス処理（左右の端をつなげる処理）のため、少しばかりコードが増えています。あと、水量を変化させないような処理も追加されています。

float[] gPos=new float[51];
float[] gPosNew=new float[51];
float[] gPosOld=new float[51];

void setup() {
  size(500,100);
}

void draw() {
  if(mousePressed==true){
    int t=mouseX/10;
    float d=(50-mouseY)-gPos[t];
    gPos[mouseX/10]+=d;
    gPos[(t+5)%51]-=d;  // adjust the amount of water.
  }
  update();
  background(255,255,255);
  strokeWeight(1);
  stroke(0,0,0);
  for(int i=0; i<50; i++) {
    line(i*10,-gPos[i]+50,(i+1)*10,-gPos[i+1]+50);
  }
  fill(0,128,255);
  noStroke();
  for(int i=0; i<50; i++) {
    beginShape();
      vertex(i*10,-gPos[i]+50);
      vertex(i*10,100);
      vertex((i+1)*10,100);
      vertex((i+1)*10,-gPos[i+1]+50);
    endShape(CLOSE);
  }
}

void update() {
  float h=1.0f;
  float d=0.99f;
  float c2=0.8f;
  for(int i=0; i<51; i++) {
    float a=c2*(gPos[(i-1+51)%51]+gPos[(i+1)%51]-2*gPos[i]);
    gPosNew[i]=(1+d)*gPos[i]-d*gPosOld[i]+0.5f*a*h*h;
  }
  for(int i=0; i<51; i++) {
    gPosOld[i]=gPos[i];
    gPos[i]=gPosNew[i];
  }
}