Raycasterを実装する｜ブログ｜スタジオスプーン株式会社

この記事で取り扱うRaycasterは、オブジェクトとの当たり判定を計算する機能です。

Webサイト上で3Dシーンを扱う場合に、「オブジェクトをマウスでクリックしたりホバーしたときに何かインタラクションさせたい」ということはよくあります。

Three.jsではRaycasterというオブジェクトとの当たり判定を計算するクラスが用意されていますが、ネイティブのWebGPU APIを扱う場合、このRaycasterを実装する必要があります。

実装サンプルのモデルは、なぜ傾いているのか

3DCGにはシーングラフという概念があります。
実装サンプルでも、Scene - Group - Modelというシーングラフを作成して、Groupに回転や移動をさせています。
シーングラフはModelのワールドマトリクスを解決するための概念で、このModelのワールドマトリクスが適切に更新され、後述するAABBに影響を及ぼさないかを確認するためです。

手順

以下の手順で実装をしていきます。

画面上を動くポインターからRayを生成する
オブジェクトの頂点座標に対してAABB*を求める
RayがAABBと交差するか判定する
Rayがオブジェクトと交差するか判定する

*AABB: Axis-Aligned Bounding Box

Rayの生成

この工程は必ずしも必要ではないですが、Raycasterの主な用途として「マウスポインターがオブジェクトに当たったかどうかを判定する」場合が多いため、画面座標からRayを生成します。

最終的に生成されるRayは、ワールド座標系のOrigin（Rayの原点）と、ワールド座標系のDirection（Rayの方向）になります。

まず、mousemoveやpointermoveイベントから、マウスポインターの位置を-1 ~ 1の範囲で取得します。y方向はWebGPUの座標系に合わせるためflipさせます。

const x = (e.clientX / canvas.clientWidth) * 2 - 1
const y = (1 - e.clientY / canvas.clientHeight) * 2 - 1

次に、cameraの射影行列 $P$ 、ビュー行列 $V$ からOriginとDirectionを求めていきます。
ワールド座標系の位置 $R_{world}$ から、正規デバイス座標系の位置 $R_{ndc}$ は以下のように求まります。

R_{clip} = PVR_{world} \\[1em] R_{world} = V^{-1}P^{-1}R_{clip} \\[1em]

また、 $R_{ndc}$ は、 $R_{clip}$ を透視除算することで求まります。透視除算は名前の通り、射影行列 $P$ を作用させた後のクリップ座標の $w$ 成分で割ることを意味しています。

R_{ndc} = \frac{R_{clip}}{w_{clip}}

このことから、 $R_{world}$ は以下のように求まります。

R_{view}^{\prime} \coloneqq P^{-1}R_{ndc} \\[1em] R_{view} = \frac{R_{view}^{\prime}}{w_{view}^{\prime}} \\[1em] R_{world} = V^{-1}R_{view}

詳しくみる

上記の式の $R_{view}^{\prime}$ から、 $R_{view}$ を求める手順を追います。
$R_{view}^{\prime}$ は、以下のように書き表せます。

\begin{align*} R_{view}^{\prime} &= P^{-1}R_{ndc} \\[1em] &= P^{-1} \frac{R_{clip}}{w_{clip}} \\[1em] &= \frac{P^{-1}PR_{view}}{w_{clip}} \\[1em] &= \frac{R_{view}}{w_{clip}} \tag{1} \end{align*}

ここで、 $R_{view}$ を以下のように置くと、

R_{view} = \begin{bmatrix} x \\ y \\ z \\ 1 \\ \end{bmatrix}

$R_{view}^{\prime}$ は、式 $(1)$ より以下のようになります。

R_{view}^{\prime} = \frac{1}{{w_{clip}}} \begin{bmatrix} x \\ y \\ z \\ 1 \\ \end{bmatrix} = \begin{bmatrix} x/w_{clip} \\ y/w_{clip} \\ z/w_{clip} \\ 1/w_{clip} \\ \end{bmatrix}

$w$ 成分である $w_{view}^{\prime}$ は、以下のようになります。

w_{view}^{\prime} = \frac{1}{w_{clip}}

式 $(1)$ から、

\begin{align*} R_{view}^{\prime} &= \frac{R_{view}}{w_{clip}} \\[1em] R_{view} &= R_{view}^{\prime} × w_{clip} \\[1em] &= \frac{R_{view}^{\prime}}{w_{view}^{\prime}} \end{align*}

つまり、 $R_{view}$ は、 $P^{-1}R_{ndc}$ を計算したあとの $R_{view}^{\prime}$ について、その $w$ 成分で除算した値と等しくなります。

ワールド座標系のNear面とFar面をこの計算で求めます。
そして、Origin $O$ とDirection $D$ は以下のように求まります。

\begin{align*} O &= R_{near} \\[1em] D &= normalize(R_{far} - R_{near}) \end{align*}

コードをみる

unproject(ndc: [number, number, number]): [number, number, number] {
  const p = vec4.create(...ndc, 1)
 
  const view = vec4.transformMat4(p, this.projectionMatrixInverse)
  // prettier-ignore
  view.set([
    view[0] / view[3],
    view[1] / view[3],
    view[2] / view[3],
    1,
  ])
 
  const world = vec4.transformMat4(view, this.viewMatrixInverse)
  return [world[0], world[1], world[2]]
}

Camera.ts（新しいタブで開く）

const worldNear = this.camera.unproject([x, y, 0])
const worldFar = this.camera.unproject([x, y, 1])
 
this.setOrigin(worldNear)
this.setDirection(vec3.sub(worldFar, worldNear)) // setDirection内部でnormalizeしている

PointerRay.ts（新しいタブで開く）

AABBを求める

AABBとは、Axis-Aligned Bounding Boxのことで、「座標軸に平行なオブジェクトを包絡するBox」という意味です。

AABBは、オブジェクトの頂点座標に対して求めていきます。これは、モデルがローカル行列をもとに移動や回転、スケールをしたとしても、頂点座標自体に変更がないためAABBの再計算する必要をなくすためです。

AABBは、オブジェクトの全頂点の位置座標について、軸ごとにmin, maxをとります。

for (let i = 0; i < this.positions.length / 3; i++) {
  const x = this.positions[i * 3 + 0]
  const y = this.positions[i * 3 + 1]
  const z = this.positions[i * 3 + 2]
  vec3.min(this.boundingBox.min, [x, y, z], this.boundingBox.min)
  vec3.max(this.boundingBox.max, [x, y, z], this.boundingBox.max)
}

Geometry.ts（新しいタブで開く）

RayとAABBの交差を求める

Rayとオブジェクトの交差を求める前に、AABBとの交差を求めます。
オブジェクトの交差判定では、オブジェクトを構成する多量の三角形に対して計算を行うため、負荷軽減のためオブジェクトにかすらないRayに対してはループ計算をしないようにします。

求めたAABBは頂点座標系で、Rayはワールド座標系だったので、対象のオブジェクトのワールドマトリクスを使って、Rayを頂点座標系に変換します。

ワールド座標系は、頂点座標系にワールドマトリクスを作用させると求まりました。

R_{world} = M_{world}R_{local}

そのため、 $R_{local}$ は $M_{world}^{-1}$ を使って以下のように求まります。

R_{local} = M_{world}^{-1}R_{world}

// local origin
const O = vec3.transformMat4(this.ray.worldOrigin, target.worldMatrixInverse)
// local direcition
const D = vec3.transformMat4Upper3x3(this.ray.worldDirection, target.worldMatrixInverse)

w成分に注意

Directionはベクトルで平行移動要素が必要ないため、行列を掛ける時の $w$ 成分は0とします。
wgpu-matrixでは、vec3.transformMat4を使うと $w$ 成分を1として計算を行い、vec3.transformMat4Upper3x3を使うと $w$ 成分を0として計算を行います。

RayのDirection上のある位置 $P$ が、AABBのそれぞれの成分の範囲内に収まっていれば、交差していると言えます。

P = O + t\bm{D} \\[1em] min \le P \le max \\[1em] \frac{min - O}{\bm{D}} \le t \le \frac{max - O}{\bm{D}}

つまり、以下を満たす $t$ が同時に存在すれば、交差しているということです。

\frac{x_{min} - O_{x}}{D_{x}} \le t \le \frac{x_{max} - O_{x}}{D_{x}} \\[1em] \frac{y_{min} - O_{y}}{D_{y}} \le t \le \frac{y_{max} - O_{y}}{D_{y}} \\[1em] \frac{z_{min} - O_{z}}{D_{z}} \le t \le \frac{z_{max} - O_{z}}{D_{z}} \\[1em]

各成分の「最小値の最大」と、「最大値の最小」をとり、その大小関係を見ればよいです。

例えば、以下のように置いて、

t_{min,x} \coloneqq \frac{x_{min} - O_{x}}{D_{x}}

max(t_{min,x}, t_{min,y}, t_{min,z}) \le t \le min(t_{max,x}, t_{max,y}, t_{max,z})

コードをみる

private insideAABB(target: Mesh, O: Vec3, D: Vec3) {
  const aabb = target.geometry.getBoundingBox()
 
  let near = -Infinity
  let far = Infinity
 
  for (let i = 0; i < 3; i++) {
    if (D[i] === 0) {
      // この軸以外の面上にあるRay（例えば、Dx=0の場合、yz平面上のRayとなる）
      if (O[i] < aabb.min[i] || O[i] > aabb.max[i]) {
        // 例えば、Dx=0の場合、Oxの範囲を見るだけでAABBに含まれるかどうかを判別できる
        return false
      }
      continue
    }
 
    const t1 = (aabb.min[i] - O[i]) / D[i]
    const t2 = (aabb.max[i] - O[i]) / D[i]
 
    near = Math.max(near, Math.min(t1, t2))
    far = Math.min(far, Math.max(t1, t2))
 
    if (near > far) return false
  }
 
  return far >= Math.max(near, 0)
}

Raycaster.ts（新しいタブで開く）

Rayとオブジェクトの交差を求める

オブジェクトは無数の三角形ポリゴンから構成されるため、Rayと三角形ひとつひとつに対して交差判定を行います。

三角形ABC上の点 $P$ は、Rayと三角形の頂点座標で、それぞれ以下のように表せます。

\begin{align*} P &= O + t\bm{D} \\[1em] P &= A + u\bm{AB} + v\bm{AC} \end{align*}

以下の連立方程式を解いて、未知数 $u, v, t$ を求めていきます。

O + t\bm{D} = A + u\bm{AB} + v\bm{AC} \\[1em] u\bm{AB} + v\bm{AC} - t\bm{D} = O - A \\[1em] u\bm{AB} + v\bm{AC} - t\bm{D} = \bm{AO} \\[1em]

ベクトル $\bm{AB}, \bm{AC}, \bm{-D}, \bm{AO}$ をそれぞれ $b, c, d, o$ と置いてマトリクスの形にすると、 $u, v, t$ は以下のように求まります。

\begin{bmatrix} b_{x} & c_{x} & d_{x} \\ b_{y} & c_{y} & d_{y} \\ b_{z} & c_{z} & d_{z} \\ \end{bmatrix} \begin{bmatrix} u \\ v \\ t \\ \end{bmatrix} = \begin{bmatrix} o_{x} \\ o_{y} \\ o_{z} \\ \end{bmatrix} \\[1em] \begin{bmatrix} u \\ v \\ t \\ \end{bmatrix} = \begin{bmatrix} b_{x} & c_{x} & d_{x} \\ b_{y} & c_{y} & d_{y} \\ b_{z} & c_{z} & d_{z} \\ \end{bmatrix}^{-1} \begin{bmatrix} o_{x} \\ o_{y} \\ o_{z} \\ \end{bmatrix}

逆行列の計算は割愛します。

このとき、以下の条件を満たせば交差点が存在することになります。

$t \ge 0$ ：交差点がRayの前方にある
$u \ge 0, v \ge 0, u + v \le 1$ ：交差点が三角形内部にある

コードをみる

for (let i = 0; i < target.geometry.indices.length / 3; i++) {
  const iA = target.geometry.indices[i * 3 + 0]
  const iB = target.geometry.indices[i * 3 + 1]
  const iC = target.geometry.indices[i * 3 + 2]
 
  this.setPosition(target, iA, A)
  this.setPosition(target, iB, B)
  this.setPosition(target, iC, C)
 
  // O + tD = A + u(B - A) + v(C - A)
  // u(B - A) + v(C - A) - tD = O - A
  vec3.sub(B, A, AB)
  vec3.sub(C, A, AC)
  vec3.sub(O, A, AO)
 
  mat3.set(
    AB[0], AB[1], AB[2],
    AC[0], AC[1], AC[2],
    -D[0], -D[1], -D[2],
    M,
  )
 
  // Rayが面に対して平行な場合は、交点を持たず、逆行列が存在しない（行列式|M| = 0 となる）
 
  if (Math.abs(mat3.determinant(M)) > 1e-8) {
    mat3.inverse(M, invM)
    // uvt = invM x AO
    vec3.transformMat3(AO, invM, uvt)
    const [u, v, t] = [...uvt]
 
    // 条件
    // t: 交差点がrayの前方にある
    // u, v: 交差点が三角形内部にある
    if (t >= 0 && u >= 0 && v >= 0 && u + v <= 1) {
      // 交点をもつ
      const localIntersectPos = vec3.add(O, vec3.mulScalar(D, t))
      const worldIntersectPos = vec3.transformMat4(localIntersectPos, target.worldMatrix)
      // ...
    }
  }
}

Raycaster.ts（新しいタブで開く）

Raycasterを実装する

手順

Rayの生成

AABBを求める

RayとAABBの交差を求める

Rayとオブジェクトの交差を求める

セクションへのリンクRELATED ARTICLES