はじめに

この記事は、自分が忘れたとき用のメモです。Qiita に残すべきではないかもしれません。コメントで指摘があれば削除するかもしれません。

GLSL は縦ベクトル

今の僕は GLSL を「プログラミング言語」と認識しています。つまり C 言語や Python と同列です。GLSL では、ベクトルや行列がプリミティブ型として定義されています。たとえば vec3 v; はベクトルを表し、mat3 m; は行列を表します。GLSL のベクトルは縦ベクトルなので、vec3 r = m * v; は許されていますが vec3 r = v * m; は許されていません。

「GLSL は縦ベクトル」という特徴は JavaScript にまで影響します。モデル変換行列を M、ビュー変換行列を V、プロテクション変換行列を P としたとき、これらの行列をひとつにまとめる場合は P * V * M の順序でなければなりません。M * V * P としてしまうと、これに座標 v を掛けるときに M * V * P * v となり、プロテクション変換から適用されてしまうからです。適用する順序は M -> V -> P でなければなりません。

JavaScript でのベクトルや行列の扱い方を教えて

JavaScript でベクトルと行列を扱う場合は、ライブラリを使うのが楽です。いくつか選択肢がありますが、glMatrix を使えばとりあえずは安心です。理由は、GitHub のスター数が 3000 以上だからです。しかし、演算の際に多少癖のある書き方をしなければなりません。具体的には、行列同士の積の結果は、関数の戻り値として得るのではなく、引数として渡した変数の中に格納されます。このような方式にしているのはおそらく理由があります。たとえばメモリ節約のためだとか。しかし、個人で小さなプログラムを開発するだけであれば、この方式は無駄に冗長です。glMatrix をラップして扱いやすくするのがいいです。

// glMatrix をそのまま使用
const one = glMatrix.mat4.fromValues(...);
const two = glMatrix.mat4.fromValues(...);
const out = glMatrix.mat4.create();
glMatrix.mat4.multiply(out, one, two);

// glMatrix をラップした自作のクラスを使用
const one = Mat4(...);
const two = Mat4(...);
const out = one.multiply(two);

Mat4 というのは関数で、Matrix4 クラスのインスタンスを返します。なぜ new Matrix4(...) ではなく Mat4(...) としているかというと、微妙に長いからです。

とにかく三角形を描画したい！

index.html と index.js と gl-matrix.js の3つのファイルを用意します。gl-matrix.js はこちらにある可能性が高いです。もしリンクが切れている場合は、glMatrix 公式サイトのどこかにあります。このリンクが切れることはおそらくないでしょう。index.html と index.js は次の内容にします。

index.html

<!DOCTYPE html>
<html>
  <head><title></title></head>
  <body>
    <canvas id="canvas" width="640" height="480"></canvas>
    <script src="gl-matrix.js"></script>
    <script src="index.js"></script>
  </body>
</html>

index.js

const mat4 = glMatrix.mat4;

const vss = `
attribute vec3 position;
uniform mat4 mvpMatrix;
void main(void) {
  gl_Position = mvpMatrix * vec4(position, 1.0);
}
`;
const fss = `
void main(void) {
  gl_FragColor = vec4(1.0, 1.0, 1.0, 1.0);
}
`;

window.addEventListener("DOMContentLoaded", () => {
  const gl = document.getElementById("canvas").getContext("webgl");
  const prg = gl.createProgram();
  {
    const vs = gl.createShader(gl.VERTEX_SHADER);
    const fs = gl.createShader(gl.FRAGMENT_SHADER);
    gl.shaderSource(vs, vss);
    gl.shaderSource(fs, fss);
    gl.compileShader(vs);
    gl.compileShader(fs);
    gl.attachShader(prg, vs);
    gl.attachShader(prg, fs);
    gl.linkProgram(prg), gl.useProgram(prg);
  }
  const position = [
    -1, 0, 0,
     1, 0, 0,
     0, 1, 0,
  ];
  {
    const vbo = gl.createBuffer();
    gl.bindBuffer(gl.ARRAY_BUFFER, vbo);
    gl.bufferData(gl.ARRAY_BUFFER, new Float32Array(position), gl.STATIC_DRAW);
    const loc = gl.getAttribLocation(prg, "position");
    gl.enableVertexAttribArray(loc);
    gl.vertexAttribPointer(loc, 3, gl.FLOAT, false, 0, 0);
  }
  {
    const mMatrix = mat4.create();
    const vMatrix = mat4.create();
    const pMatrix = mat4.create();
    const mvpMatrix = mat4.create();

    mat4.fromTranslation(mMatrix, [0, 0, 0]);
    mat4.lookAt(vMatrix, [0, 0, 10], [0, 0, 0], [0, 1, 0]);
    mat4.perspective(pMatrix, 45, gl.canvas.width / gl.canvas.height, 0.1, 100);

    mat4.multiply(mvpMatrix, mMatrix, mvpMatrix);
    mat4.multiply(mvpMatrix, vMatrix, mvpMatrix);
    mat4.multiply(mvpMatrix, pMatrix, mvpMatrix);

    const loc = gl.getUniformLocation(prg, "mvpMatrix");
    gl.uniformMatrix4fv(loc, false, mvpMatrix);
  }
  gl.clearColor(0, 0, 0, 1);
  gl.clear(gl.COLOR_BUFFER_BIT);
  gl.drawArrays(gl.TRIANGLES, 0, 3);
  gl.flush();
});

index.html をブラウザで表示すると、画面に三角形が描画されます。

GLSL を別ファイルに記述したいんだけど

頑張ってください。async/await と fetch() を使えば綺麗に書けますよ。頂点シェーダの拡張子は .vert 、フラグメントシェーダの拡張子は .frag にするといいです。VSCode でシンタックスハイライトが付きます。セキュリティ周りでエラーが出る場合は、python3 -m http.server 8000 で Web サーバを立ち上げればいいです。

頂点シェーダとフラグメントシェーダの処理単位の違い

頂点シェーダは頂点の数だけ、フラグメントシェーダは描画されるピクセルの数だけ実行されます。そして、頂点シェーダからフラグメントシェーダに何らかの情報を渡すことができます。疑問を感じませんか？1対1の関係ではないのに、なぜ情報を渡すことができるのでしょう。結論としては、うまく補間されているのです。この補間は開発者の知らないところで行われています。興味があれば調べるといいです。「ラスタライズ処理」といいます。このあたりの話は、こちらのページで図を交えて丁寧に解説されています。

球体を表示したい

けっこう大変です。IQが要ります。頭に3次元を思い浮かべて頂点の割り振り方を考える必要があります。ノートに図示してみるのもいいです。ただ、球体はどこから見ても円になるので、意外と図示するのが難しいです。地球の経線と緯線のような線を書けばいいです。その手間すら省きたいのであれば、次のコードをコピペすればいいです。…と思ってコードを探したのですが、見つかりませんでした。どうやら削除したようです。ということで、頑張って実装してください。実装したときの達成感は半端ないですよ。

影を表現したい

実は今の僕はよく理解できていません。頑張って勉強してください。あわよくば、以下のコードに含まれる法線ベクトル周りのバグを取り除いてください。Polygon#rotate() にバグがあるらしい（サイコロを描画したときの影の付き方が微妙におかしい）のですが、今の僕にはわかりませんでした。

以下のコードの簡単な説明をします。前提として、glMatrix がなければ動作しません。

Polygon というクラスは、画面に描画する物体を表します。決して三角形「だけ」を表すわけではありません。この名付け方に微妙に納得がいってませんが、Object だと汎用的すぎるため、とりあえず Polygon にしています。もっといい名前があれば変更してください。

コンストラクタとして2つの引数を取ります。まず data について。data は頂点の情報で、配列です。要素は position と color と normal の3つのキーを持つオブジェクトです。現時点ではこの3つで十分です（法線ベクトル周りのバグを取り除くために、更にいくつかのキーを追加する必要があるかもしれません）。position は頂点の座標です。これは必須中の必須です。color は頂点の色です。これも必須です。normal は法線ベクトルです。これも必須です。キーが欠けていると、予期しないエラーが発生するかもしれません。次に index について。index は頂点番号の配列です。要素数は3の倍数です。まあ詳しく説明せずともわかると思います（正直なところ説明が面倒）。

Model#prepare() を呼べば、そのモデルを描画する準備が整います。準備が整うと、あとは gl.drawElements() を呼ぶだけで画面にモデルが描画されます。しかし、今の僕は悩んでいます。gl.drawElements(gl.TRIANGLES, model.index.length, gl.UNSIGNED_SHORT, 0); というように、モデルの index の長さを参照しなければいけないんですよね。これはあまりよくないです。ということで、最適な設計は「Model 自身に描画メソッドを実装すること」だと思っています。model.draw(gl); のような感じですね。ここまでくれば、Model#prepare() も必要なくなるように思います。まあどのような実装にするかはあなたに任せます。

const vec3 = glMatrix.vec3;
const mat4 = glMatrix.mat4;

class Polygon {
  constructor(data = [], index = []) {
    this.data = data;
    this.index = index;
  }
  add(poly) {
    const index = this.index.concat(
      poly.index.map(i => i + this.data.length)
    );
    const data = this.data.concat(poly.data);
    return new Polygon(data, index);
  }
  rotate(axis, angle) {
    return new Polygon(this.data.map(d => {
      const m = mat4.create();
      const p = vec3.create();
      const n = vec3.create();
      mat4.fromRotation(m, Math.PI / 180 * angle, axis);
      vec3.transformMat4(p, d.position, m);
      vec3.transformMat4(n, d.normal, m);
      return { ...d, position: Array.from(p), normal: Array.from(n) };
    }), this.index);
  }
  translate(x, y, z) {
    return new Polygon(this.data.map(d => {
      const m = mat4.create();
      const p = vec3.create();
      mat4.fromTranslation(m, [x, y, z]);
      vec3.transformMat4(p, d.position, m);
      return { ...d, position: Array.from(p) };
    }), this.index);
  }
  primitive() {
    return {
      index: this.index,
      position: this.data.map(d => d.position).flat(),
      color: this.data.map(d => d.color).flat(),
      normal: this.data.map(d => d.normal).flat(),
    };
  }
  model(gl) {
    return new Model(gl, this);
  }
}

class Model {
  constructor(gl, polygon) {
    const { position, color, normal, index } = polygon.primitive();

    this.position = gl.createBuffer();
    gl.bindBuffer(gl.ARRAY_BUFFER, this.position);
    gl.bufferData(gl.ARRAY_BUFFER, new Float32Array(position), gl.STATIC_DRAW);
    this.color = gl.createBuffer();
    gl.bindBuffer(gl.ARRAY_BUFFER, this.color);
    gl.bufferData(gl.ARRAY_BUFFER, new Float32Array(color), gl.STATIC_DRAW);
    this.normal = gl.createBuffer();
    gl.bindBuffer(gl.ARRAY_BUFFER, this.normal);
    gl.bufferData(gl.ARRAY_BUFFER, new Float32Array(normal), gl.STATIC_DRAW);

    this.index = index;
  }
  prepare(gl, prg) {
    const list = [
      ["position", this.position],
      ["color", this.color],
      ["normal", this.normal],
    ];
    list.forEach(([name, vbo]) => {
      gl.bindBuffer(gl.ARRAY_BUFFER, vbo);
      const loc = gl.getAttribLocation(prg, name);
      gl.enableVertexAttribArray(loc);
      gl.vertexAttribPointer(loc, 3, gl.FLOAT, false, 0, 0);
    });

    const index = gl.createBuffer();
    gl.bindBuffer(gl.ELEMENT_ARRAY_BUFFER, index);
    gl.bufferData(gl.ELEMENT_ARRAY_BUFFER, new Int16Array(this.index), gl.STATIC_DRAW);
  }
}

なぜ GPU は数千コアとかあるのに CPU とそれほど値段が変わらないの？

わからないです。頑張って調べてください。

ゴルフゲーム作りたいんだけど、衝突判定どうすればいいの？

ベクトルの内積と外積を使えばいいです。この2つでたいてい事足ります。内積と外積の感覚を取り戻したいのであれば、こちらの記事を一通り読めばいいです。ものすごくわかりやすく解説されています。ベクトルの内積と外積も glMatrix で計算しますが、行列のときと同じくそのままでは扱いづらいです。うまくラップしましょう。最初は速度やメモリのことは考えなくてもいいです。とりあえず動くものを作ることを目指しましょう。実装するとなると結構難しいですよ。速度やメモリよりも集中すべき問題があります。それは浮動小数点数の誤差や、フレーム特有の離散的であるがゆえに発生する問題などです。これらの問題を解決するのは難しいです。まあ頑張ってください。僕は途中で投げ出しました。

もう少し詳しく書き残しておきます。

「ゴルフボールが地面に衝突する」とはどういうことでしょうか。地面はポリゴンの組み合わせで作られており、ゴルフボールもまた同じです。ということで、ポリゴン同士の衝突判定ってどうするの？という問いに答える必要があります。しかし、この問題を解くのはたぶん難しいです。よって、「ポリゴン同士の衝突判定」の問題に帰結すべきではないです。それよりも、ゴルフボールを球体と考えて「ポリゴンと球体の衝突判定」の問題に帰結したほうがずっと解きやすいです。このことを念頭に置いて考察を続けます。

飛んでいるゴルフボールは離散的です。スローカメラでボールの飛んでいる様子を映したとき、あるフレームで (1, 2, 3) にあったボールが、次のフレームで突然 (2, 3, 4) に移動したりします。これによりひとつの懸念点が出てきます。「ボールが地面をすり抜けてしまうかもしれない」という懸念です。この問題の解決策は、問題を「移動前と移動後の座標を結んだ線分が地面に接しているか」と言い換えればいいです。離散的だったものを連続的に考えるのです。これにより、ボールが地面をすり抜けることはありません。

眠たいので続きはまた明日書きます。残しておきたいコードがあるのでここに残しておきます。

const vec3 = glMatrix.vec3;

class Vector3 {
  constructor(x, y, z) {
    if (x.length != null) {
      [this.x, this.y, this.z] = x;
    } else {
      [this.x, this.y, this.z] = [x, y, z];
    }
  }
  add(that) {
    const out = vec3.create();
    vec3.add(out, this.primitive(), that.primitive());
    return new Vector3(out);
  }
  sub(that) {
    const out = vec3.create();
    vec3.sub(out, this.primitive(), that.primitive());
    return new Vector3(out);
  }
  dot(that) {
    return vec3.dot(this.primitive(), that.primitive());
  }
  cross(that) {
    const out = vec3.create();
    vec3.cross(out, this.primitive(), that.primitive());
    return new Vector3(out);
  }
  scale(n) {
    const out = vec3.create();
    vec3.scale(out, this.primitive(), n);
    return new Vector3(out);
  }
  length() {
    return vec3.length(this.primitive());
  }
  normalize() {
    const out = vec3.create();
    vec3.normalize(out, this.primitive());
    return new Vector3(out);
  }
  quat(q) {
    const out = vec3.create();
    vec3.transformQuat(out, this.primitive(), q);
    return new Vector3(out);
  }
  primitive() {
    return [this.x, this.y, this.z];
  }
}

function Vec3(x, y, z) {
  return new Vector3(x, y, z);
}

次の関数は、平面ABCと直線DEの交点を求める関数です。こちらのページを参考にしました。

function kouten(A, B, C, D, E) {
  const F = E.sub(D);

  const N = B.sub(A).cross(C.sub(A));
  const a = N.x;
  const b = N.y;
  const c = N.z;
  const d = - a * A.x - b * A.y - c * A.z;

  const t = - (a * D.x + b * D.y + c * D.z + d) / (a * F.x + b * F.y + c * F.z);

  const x = D.x + t * F.x;
  const y = D.y + t * F.y;
  const z = D.z + t * F.z;

  return new Vector3(x, y, z);
}

次の関数は、点Dから平面ABCへ垂線を下ろしたときの交点を求める関数です。

function suisen(A, B, C, D) {
  const N = C.sub(A).cross(B.sub(A)).normalize();
  return D.add(N.scale(A.sub(D).dot(N)));
}

本当は次のように記述したいんですよね。でも叶わぬ願いです。

function suisen(A, B, C, D) {
  const N = ((C - A) * (B - A)).normalize();
  return D + N.scale((A - D) . N);
}

次の関数は、△ABCに点Qが含まれているかどうかを判定する関数です。浮動小数点数の扱い方がわからないため、判定の部分は未実装です。

function contain(A, B, C, Q) {
  const a = B.sub(A).cross(Q.sub(A));
  const b = C.sub(B).cross(Q.sub(B));
  const c = A.sub(C).cross(Q.sub(C));
  // TODO
}

明日になりました。上のTODOにどんなコードを書けばいいかがわかりました。「2つのベクトルが同じような方向を向いているか」は、2つのベクトルの内積を取ることでわかります。内積が正であれば同じ向きを向いています。これを実装に落とし込むと次のようになります。

function contain(A, B, C, Q) {
  const a = B.sub(A).cross(Q.sub(A));
  const b = C.sub(B).cross(Q.sub(B));
  const c = A.sub(C).cross(Q.sub(C));
  return a.dot(b) >= 0 && a.dot(c) >= 0;
}

こんなにスマートに解決できるんですね。ベクトルが大好きになりました。3Dゲームを作るとき、ベクトルの感覚を掴んでいればかなり楽しく開発を進められそうです。

仕事中に色々考えていたのですが、まず問題が不明確でした。ただ漠然と「ゴルフゲームを実装したい」と思っていただけで、それを具体的な問題に言い換えることをしませんでした。ですので、考えることがわからずに悩んでいただけでした。この状態はよくないので、問題を明確化しようと思います。

僕がまず実装したいのは、ゴルフボールを飛ばして、地面を何度かバウンドしながら転がり続けて、最終的に止まるという機能です。しかしこれでもまだ問題が不明確です。もっと問題を明確にします。

ゴルフボールを球体と考えます。その球体は重さ $m$ と半径 $r$ を持ちます。最初、ゴルフボールには初速 $\overrightarrow{v_0} $ が与えられます。ボールには常に重力加速度 $\vec{g}$ がかかります。それと同様に、風の力 $\vec{W}$ がかかります。この $2$ つの力を受けながらボールは飛び、いずれ地面に着地します。地面は $N$ 個の三角形で構成されており、$i$ 番目の三角形の頂点を $A_i, B_i, C_i (1 \le i \le N)$ とおきます。ボールは離散的に飛ぶものとします。初期位置は $H$ とし、そのときの時間を $0$ とします。

ここでひとつのルールを決めておきます。それは、「速度を求めてからボールを移動させる」というルールです。「ボールを移動させてから速度を求める」との2択ですが、どちらでも特に変わらない気がするので、てきとーに一方を選びます。

現在の時間を $t$ とします。このときのボールの速度を $\vec{v_t}$ 、位置を $\vec{h_t}$ とします。まずは時間 $t+1$ のときの速度を求めてみます。ボールには、重力加速度 $\vec{g}$ と風の力 $\vec{F}$ がかかります。風の力を加速度に変換すると $\vec{F} / m$ となります。よって、$t + 1$ のときの速度は $\vec{v_t} + \vec{g} + \vec{F} / m$ です。これに時間 $1$ を掛けると位置の変化量になります。よって、$t + 1$ のときの位置は $\vec{h_t} + \vec{v_t} + \vec{g} + \vec{F} / m$ です。

\begin{cases}
\vec{v_{t+1}} = \vec{v_t} + \vec{g} + \frac{\vec{F}}{m} \\
\vec{h_{t+1}} = \vec{h_t} + \vec{v_t} + \vec{g} + \frac{\vec{F}}{m}
\end{cases}

これで、地面を考慮しない場合のボールの速度と位置を求めることができました。次は地面との衝突を考えていきます。「地面と衝突している」というのは、地面とボールの中心との距離が $r$ 以下であるということです。しかし、ここで離散的であるが故の問題が浮上します。$t$ 時点では地面より $r$ 以上上の位置にボールがあって、$t + 1$ 時点で $r$ 以上下の位置にボールがある場合があるのです。この場合でも、地面と衝突したとみなせます。これをどのように考えればいいでしょうか。簡単です。「地面より上」をプラスで表し、「地面より下」をマイナスで表すのです。条件はそのまま「地面とボールの中心との距離が $r$ 以下」です。これで衝突判定ができました。

数式を使って明確にします。地面といっても $N$ 個の三角形があります。今回は、$i$ 番目の三角形との衝突を考えます（これができれば、あとは $N$ 個の三角形で同じことをすればいいだけです。…なんだか計算量的に嫌な予感がしますが、とりあえず今は考えないことにします）。三角形 $A_iB_iC_i$ とボールとの距離はどのようにして求められるでしょうか。これは、ベクトルの内積と外積を用いることで求めることができます。まずは外積を使って、三角形の法線ベクトルを求めます。法線ベクトルとは、三角形の面に対して垂直なベクトルのことです。法線ベクトルは $\vec{AB} \times \vec{AC}$ で求まります。次に、正規化した法線ベクトルと $\vec{Ah_{t+1}}$ との内積を取ります。実はこれが、三角形の面とボールの中心との距離です。

まとめると、三角形の面と衝突した場合は次の不等式を満たします。

normalize(\vec{AB} \times \vec{AC}) \cdot \vec{Ah_{t+1}} \le r

「三角形」ではなく「三角形の面」と言っている部分がミソです。つまりこの条件だけでは、三角形と衝突しているかどうかがわからないのです。どうすればわかるでしょうか。「三角形にボールが衝突している」というのは、衝突したときのボールと面の接点が三角形に含まれることです。この「三角形の中に点が含まれるかどうか」というのは、外積を使えば簡単にわかります。よって、まずはボールと面の接点を求めることだけを考えればいいです。

ボールと面との接点を求めるために、まずはボールの軌道と面の交点を求めることにします。ボールの軌道の直線の方程式は $\vec{h_t} + m \cdot \vec{h_th_{t+1}}$（$m$ は媒介変数）です。この直線と三角形 $A_iB_iC_i$ の面との交点は、こちらの方法で求めることができます。これで交点を求めることができたので、ここから接点を求めていきます。ここで、先ほどの交点を $P$ とおき、求めたい接点を $Q$ とおき、座標 $\vec{h_t}$ から面 $A_iB_iC_i$ に垂線を下ろしたときの交点を $S$ とおき、ボールと面が接したときのボールの中心を $O$ とおき、面 $A_iB_iC_i$ の法線ベクトルを $\vec{n}$ とおきます。すると、$\triangle{POQ}$ と $\triangle{Ph_tS}$ は相似となります。そして、$|\vec{OQ}|$ と $|\vec{h_tS}|$ と $\vec{S}$ は次のように求めることができます。

\begin{cases}
|\vec{OQ}| = r \\
|\vec{h_tS}| = \vec{Ph_t} \cdot \vec{n} \\
\vec{S} = \vec{h_t} - |\vec{h_tS}| \cdot \vec{n}
\end{cases}

よって、$Q$ を次の式で求めることができます。

Q = \vec{P} + |\vec{OQ}| / |\vec{h_tS}| \cdot \vec{PS}

少し長かったですが、これでようやくボールと面との接点を求めることができました。あとは、この接点が三角形に含まれているかどうかです。この記事のずっと上のほうに書いた気がしますが、改めて書きます。三角形内に点が含まれているかどうかは、外積を使えば簡単に求まります。具体的には、$\vec{A_iB_i} \times \vec{A_iQ}$ と $\vec{B_iC_i} \times \vec{B_iQ}$ と $\vec{C_iA_i} \times \vec{C_iQ}$ がすべて同じ方向を向いていればいいです。これらの外積を次の文字に置き換えたとき、

\vec{a} = \vec{A_iB_i} \times \vec{A_iQ} \\
\vec{b} = \vec{B_iC_i} \times \vec{B_iQ} \\
\vec{c} = \vec{C_iA_i} \times \vec{C_iQ}

「すべて同じ方向を向いているかどうか」は、$\vec{a} \cdot \vec{b} \ge 0$ かつ $\vec{a} \cdot \vec{c} \ge 0$ を満たすかどうかでわかります。

ここまで計算して、やっと半分といったところです。しかし、ここまでついてこれたのであれば、ゴールまではもうすぐです。あと考えるべき問題は「衝突したときにボールの位置はどうなるか」だけだからです。続きは風呂に入りながら考えます。

風呂から出たあとに、上の文章を読み返しました。思っていたよりもしっかりと理論立てて考えられていました。嬉しかったです。この調子で続きも考えていきます。

まずは $\vec{h_t}$ と $\vec{O}$ と $\vec{h_{t+1}}$ の関係から考えます。これらは同じ直線上にあるので、$\vec{O} = \vec{h_t} + m \cdot \vec{h_{t+1}}$（$m$ は媒介変数）と表せます。そして $0 \le m \le 1$ を満たします。これについては、それぞれのベクトルの意味を理解していればすぐにわかります。実はこの $m$ は、三角形にボールが衝突したときの時間を表しています。たとえば $m = 0.7$ であれば $\vec{O} = \vec{h_{t+0.7}}$ です。まず知りたいことは衝突時の時間であるため、この $m$ を求めます。次に $1 - m$ を求めます。…ここで脳が活発に動きすぎたため、頭の中で話が完結してしまいました。順序立てて説明するのが若干面倒なので、ここからは勢いで書きます。

初期位置 $\vec{O}$ 、初速 $\vec{v_0}$ で、$1 - m$ 後の速度と位置はどこか？を求めればいいです。今回は、速度は初速から変わらないものとします。$\vec{O}$ は既に求まっています。$1 - m$ も計算しました。残りは速度です。ここが今回の肝です。三角形にボールが衝突したときの速度の変化を求めるのに、三角形の面の向きが重要になってきます。そして、これは $\vec{n}$ として先ほど求めました。地面の弾性が $100%$（地面がスーパーボールで出来ていると考える）で摩擦係数が $0$ のとき、衝突後の速度は、四次元数を使うことで面白いように簡単に求まります。$\vec{n}$ を軸として $\vec{v_t}$ を $180^\circ$ 回転させてやればいいのです（読み返して気付いたのですが、回転させる前に、速度ベクトルを逆方向に向ける必要がありますね）。あとはこうして求まった速度ベクトルに、反発係数や動摩擦係数を適用すればいいです。ここで、速度ベクトルを2つの成分に分解します。面に垂直なベクトルと、面に平行なベクトルです。そして、面に垂直なベクトルのほうに反発係数を適用し、面に平行なベクトルのほうに動摩擦係数を適用し、最後にその2つを足してやります。それで $\vec{v_0}$ が求まりました。速度ベクトルを2つの成分に分解するにはどうすればよいでしょうか。簡単です。$\vec{n}$ と、回転前の速度ベクトルの内積を取り、それに $\vec{n}$ を掛ければいいです。こうして面に垂直なベクトルが求まれば、あとは引き算で面に平行なベクトルが求まります。終わりました。

簡単でしたね。いや、難しかったか。というか考えるのが楽しすぎます。ベクトルの感覚を掴むとこうも楽しくなるのですね。

考察を終えたので、これから実装します。

形式ばった説明

実装の前に、形式ばった説明に言い換えます。

問題

半径 $r$ 、重さ $m$ の球体のボールがあります。時間 $t$ にボールが $h_t$ の位置にあり、速度 $\vec{v_t}$ で動いています。動いているボールには、重力加速度 $\vec{g}$ と風の力 $\vec{F}$ が常にかかります。単位時間後のボールの位置 $h_{t+1}$ と速度 $\vec{v_{t+1}}$ を求めてください。ただし、ボールの下には $N$ 個の三角形 $A_iB_iC_i$ が隙間なく詰められており、単位時間後にボールと衝突する可能性があります。衝突した場合、反発係数 $e$ 、動摩擦係数 $\mu$ で速度が減衰するものとします。あと、ボールは単位時間ごとに離散的に動くものとします。そして、$h_{t+1} = h_{t} + \vec{v_{t+1}}$ を満たすものとします。

回答

三角形との衝突は考えないことにします。すると、単位時間後の速度は $\vec{v_{t+1}} = \vec{v_t} + (\vec{g} + \vec{F} / m) \times 1 = \vec{v_t} + \vec{g} + \vec{F} / m$ になります。そして、位置をベクトルで表すと、単位時間後の位置は $\vec{h_{t+1}} = \vec{h_t} + \vec{v_t} + \vec{g} + \vec{F} / m$ になります。

\begin{cases}
\vec{v_{t+1}} = \vec{v_t} + \vec{g} + \vec{F} / m \\
\vec{h_{t+1}} = \vec{h_t} + \vec{v_t} + \vec{g} + \vec{F} / m
\end{cases}

次は三角形との衝突を考えます。ここで、三角形 $A_iB_iC_i$ の正規化された法線ベクトルを $\vec{n_i} (= normalize(\vec{A_iB_i} \times \vec{A_iC_i}))$ とすると、三角形の面とボールが衝突している場合、まず次の条件を満たします。

\vec{n_i} \cdot \vec{A_ih_{t+1}} \le r

次に、衝突時の三角形の面とボールとの接点 $Q_i$ が、三角形に含まれるかどうかを判定します。三角形の面と直線 $\vec{h_t} + x \cdot \vec{h_th_{t+1}}$ との交点を $P_i$ とすると、$P_i$ はこちらに記載されている方法で求めることができます。さらに、$\vec{h_t}$ から面 $A_iB_iC_i$ に垂線を下ろしたときの交点を $S_i$ 、ボールと面が接したときのボールの中心を $O_i$ とおきます。$|\vec{O_iQ_i}|$ 、$|\vec{h_tS_i}|$ 、$\vec{P_iS_i}$ は次の式で求めることができます。

\begin{cases}
|\vec{O_iQ_i}| = r \\
|\vec{h_tS_i}| = \vec{n_i} \cdot \vec{Ph_t} \\
\vec{P_iS_i} = \vec{h_t} + |\vec{h_tS_i}| \vec{n_i} - \vec{P_i}
\end{cases}

条件により $\triangle{P_iO_iQ_i}$ と $\triangle{P_ih_tS_i}$ は相似であるため、$\vec{P_iQ_i}$ は次のように求められます。

\vec{P_iQ_i} = (|\vec{O_iQ_i}| / |\vec{h_tS_i}|) \vec{P_iS_i}

よって、$Q_i$ は次の式で求められます。

Q_i = (|\vec{O_iQ_i}| / |\vec{h_tS_i}|) \vec{P_iS_i} + \vec{P_i}

次は、$Q_i$ が三角形 $A_iB_iC_i$ に含まれるかどうかを判定します。$\vec{a_i}$ 、$\vec{b_i}$、$\vec{c_i}$ を次のように定義すると、

\begin{cases}
\vec{a_i} = \vec{A_iB_i} \times \vec{A_iQ_i}\\
\vec{b_i} = \vec{B_iC_i} \times \vec{B_iQ_i}\\
\vec{c_i} = \vec{C_iA_i} \times \vec{C_iQ_i}
\end{cases}

三角形 $A_iB_iC_i$ に $Q_i$ が含まれている場合は、次の条件を満たします。

\vec{a_i} \cdot \vec{b_i} \ge 0 \land \vec{a_i} \cdot \vec{c_i} \ge 0

まとめると、衝突している場合は次の条件を満たします。

\vec{n_i} \cdot \vec{A_ih_{t+1}} \le r \land \vec{a_i} \cdot \vec{b_i} \ge 0 \land \vec{a_i} \cdot \vec{c_i} \ge 0

$N$ 個の三角形すべてが上の条件を満たさなかった場合、ボールはどの三角形とも衝突しなかったことを意味します。よって、最初に求めた $\vec{v_{t+1}}$ と $\vec{h_{t+1}}$ が答えです。

衝突時の時間を $t + s$ とおくと、$s$ は次の式で求められます。

s = |\vec{h_tO_i}| / |\vec{h_th_{t+1}}|

次に速度ベクトルを求めることを考えます。反発係数 $1$ 、動摩擦係数 $0$ としたとき、衝突後の速度ベクトル $\vec{v_{t+1}'}$ は、$-\vec{v_{t+1}}$ を、$\vec{n_i}$ を軸として $180^\circ$ 回転させることで求められます。次に反発係数と動摩擦係数を考えますが、その前に、$\vec{v_{t+1}'}$ を2つに分解します。そのうち、三角形 $A_iB_iC_i$ に垂直な方向の速度ベクトル $\vec{w}$ は $(\vec{n_i} \cdot \vec{v_{t+1}'}) \vec{n_i}$ で求められるため、三角形 $A_iB_iC_i$ に水平な速度ベクトル $\vec{u}$ は $\vec{v_{t+1}'} - \vec{w}$ で求められます。

\begin{cases}
\vec{w} = (\vec{n_i} \cdot \vec{v_{t+1}'}) \vec{n_i} \\
\vec{u} = \vec{v_{t+1}'} - \vec{w}
\end{cases}

摩擦力は $\mu m(\vec{g} \cdot -\vec{n_i})$ であるため、…ここで理論に穴があるっぽいことに気付いたので少し方向転換（こちらの方法でもたぶん自然なボールの動きになる）。三角形 $A_iB_iC_i$ に水平な方向にかかる反発係数っぽいものを $d$ とおくと、最終的な速度ベクトルは次のようになります。

\vec{v_{t+1}} = e\vec{w} + d\vec{u}

よって、$\vec{h_{t+1}}$ は次のようにして求められます。

\vec{h_{t+1}} = \vec{O_i} + (1 - s)\vec{v_{t+1}}

実装

function* yields(v0, h0, tris) {
  const g = Vec3(0, -9.8 / 60, 0);
  const F = Vec3(0, 0, 0);
  const e = 0.5;
  const d = 0.99;
  const r = 0.04267;
  const m = 0.04593;

  let v = v0;
  let h = h0;
  
  while (true) {
    let vout = v.add(g).add(F.scale(1 / m));
    let hout = h.add(vout);

    for (const tri of tris) {
      const A = tri[0];
      const B = tri[1];
      const C = tri[2];
      
      const n = B.sub(A).cross(C.sub(A)).normalize();
      if (n.dot(hout.sub(A)) > r) {
        continue;
      }

      const P = kouten(A, B, C, h, hout);
      const _OQ = r;
      const _hS = n.dot(h.sub(P));
      const PS = h.add(n.scale(_hS)).sub(P);
      const PQ = PS.scale(_OQ / _hS);
      const Q = PQ.add(P);

      const a = B.sub(A).cross(Q.sub(A));
      const b = C.sub(B).cross(Q.sub(B));
      const c = A.sub(C).cross(Q.sub(C));
      if (a.dot(b) < 0 || a.dot(c) < 0) {
        continue;
      }

      const O = h.sub(P).scale(_OQ / _hS).add(P);
      const _s = O.sub(h).length() / hout.sub(h).length();
      console.assert(0 <= _s && _s <= 1);

      const q = glMatrix.quat.create();
      glMatrix.quat.setAxisAngle(q, n.primitive(), Math.PI);
      const v_ = vout.scale(-1).quat(q);
      const w = n.scale(n.dot(v_));
      const u = v_.sub(w);

      vout = w.scale(e).add(u.scale(d));
      hout = O.add(vout.scale(1 - _s));
      break;
    }

    yield [vout, hout];

    v = vout;
    h = hout;
  }
}

使用例

const L = 1234;
const y = yields(Vec3(0, 10, 10), Vec3(0, 0, 0), [
  [
    Vec3(-L, -10,  L),
    Vec3( L, -10,  L),
    Vec3(-L, -10, -L)
  ],
  [
    Vec3( L, -10,  L),
    Vec3( L, -10, -L),
    Vec3(-L, -10, -L)
  ],
]);
for (let i = 0; i < 300; i++) {
  const [v, h] = y.next().value;
  console.log(v.primitive(), h.primitive());
}

実行結果

時空の狭間に消えるボール…。

原因

法線ベクトルが逆方向になっていました。以下のように修正することで直りました。

// 修正前
const PS = h.add(n.scale(_hS)).sub(P);
// 修正後
const PS = h.add(n.scale(-_hS)).sub(P);

問題点

上記の実装には不満があります。まず、ボールの回転が考慮されていません。次に、摩擦力を「三角形$A_iB_iC_i$に水平な方向にかかる反発係数っぽいもの」に置き換えて誤魔化したのですが、それだとボールの動きが少し不自然になります。どう不自然かというと、地面が少しでも傾いていれば、ボールは静止せずに動き続けるという点です。通常では、多少の傾きがあっても摩擦力により静止します。この2つの問題を解決するために、もう一度理論を構築していこうと思います。

と思ったのですが、剛体の回転が思った以上に難しいので、とりあえず今は置いておきます。まず空気抵抗と動摩擦力を実装しようと思います。といいながら、もう既に空気抵抗は実装済みです。あ、あと転がっている間は風の影響を受けたくないのと、その代わりに動摩擦力をかけるようにしようと思います。いきなりコードに手をつけると訳が分からなくなる恐れがあるので、まずは文章を書きながら設計を考えていこうと思います。

少しだけ楽させてもらいます。空気抵抗は、最初に速度を求めるときのみ考慮し、それ以外のときは考慮しないことにします。つまりは、地面に衝突した後の空気抵抗を考慮しないということです。…転がっている間は空気抵抗の影響をどうするか。本当に動摩擦力だけにするか。考え方としては、転がっている間は風の影響を無視すると決めたので、空気抵抗も風のようなものなので同じく無視する方針のほうが統一感がある気がします。ということで、転がっている間は風だけでなく空気抵抗も無視します。よって「転がっているかどうか」で条件分岐します。

四次元数使ってみたけど、通常の行列とできること変わらなくない？

おそらく変わりません。四次元数は「任意軸の回転を気軽に扱える」らしいですが、通常の行列でも glMatrix を使えば mat4.fromRotation() で気軽に扱えます。更に、回転軸と角度から四次元数を作成する quat.setAxisAngle() の回転軸は、正規化されたベクトルである必要があります。mat4.fromRotation() はそうではありません。回転軸が (1, 1, 0) とかでも許されるわけです。

まあでも、この2つの関数の実装を見ると、四次元数が簡単に作れることが分かります。mat4.fromRotation() は30行以上ありますが、quat.setAxisAngle() は7行です。

export function setAxisAngle(out, axis, rad) {
  rad = rad * 0.5;
  let s = Math.sin(rad);
  out[0] = s * axis[0];
  out[1] = s * axis[1];
  out[2] = s * axis[2];
  out[3] = Math.cos(rad);
  return out;
}

ただ、quat.setAxisAngle() を呼び出す回数はそれほど多くありませんし、3次元座標に四次元数を適用（この言い方が合っているかどうかは微妙）するときに使用する vec3.transformQuat() のコード量もそれほど少なくないです。ぱっと見、3次元座標に行列を適用するときのコード量よりも多いです。計算にかかる時間もあまり変わらない気がします。

今のところ、回転に四次元数を使うのは余計な混乱を生むので使わなくていいです。ただ、16個の数がたった4個の数だけで表せるというのは頭の片隅に入れておけばいいです（元は4つの数（回転軸を表す3つの数と角度を表す1つの数）だったものが別の4つの数で表せて、更にその4つの数は色々と扱いやすいというのは、どこか美しさを感じる）。

おわりに

とりあえず書き残したいことは書けたので、これでようやく Haskell や競プロに集中できます。