WebGL2でGPGPU Trails

WebGL2でGPGPU Trailsを実装してみたので解説記事を残しておこうと思います。そもそもGPGPU Trailsという言葉が一般的なのかよくわからないですが、ライン状のオブジェクトの位置と速度をGPU上で計算して大量にうねうね動かすやつです。

GitHubにソースコード全体を置いておいたので、必要があればそちらを確認してください。
https://github.com/aadebdeb/WebGL2_GPGPU_Trails

---

まずは単純な線のTrailsから解説していきます。
https://aadebdeb.github.io/WebGL2_GPGPU_Trails/line
![line.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/a8bbf73c-b1bd-f32b-98d9-6c9071ed01a4.png)


はじめに、trailの位置と速度を格納するテクスチャを作成しておきます。今回はテクスチャのx要素が何番目のtrailかを表し、y要素が何番目の頂点かを表すようにしています。例えば、テクスチャの(x, y) = (100, 200)は100番目のtrailの200番目の頂点の情報を格納しています。

```js

function createTexture(gl, sizeX, sizeY, internalFormat, format, type) {
  const texture = gl.createTexture();
  gl.bindTexture(gl.TEXTURE_2D, texture);
  gl.texImage2D(gl.TEXTURE_2D, 0, internalFormat, sizeX, sizeY, 0, format, type, null);
  gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_MAG_FILTER, gl.NEAREST);
  gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_MIN_FILTER, gl.NEAREST);
  gl.bindTexture(gl.TEXTURE_2D, null);
  return texture;
}

function createTrailFbObj(gl, textureSize) {
  const framebuffer = gl.createFramebuffer();
  gl.bindFramebuffer(gl.FRAMEBUFFER, framebuffer);
  // trailの位置を格納するテクスチャの作成
  const positionTexture = createTexture(gl, textureSize.x, textureSize.y, gl.RGBA32F, gl.RGBA, gl.FLOAT);
  gl.bindTexture(gl.TEXTURE_2D, positionTexture);
  gl.framebufferTexture2D(gl.FRAMEBUFFER, gl.COLOR_ATTACHMENT0, gl.TEXTURE_2D, positionTexture, 0);
  // trailの速度を格納するテクスチャの作成
  const velocityTexture = createTexture(gl, textureSize.x, textureSize.y, gl.RGBA32F, gl.RGBA, gl.FLOAT);
  gl.bindTexture(gl.TEXTURE_2D, velocityTexture);
  gl.framebufferTexture2D(gl.FRAMEBUFFER, gl.COLOR_ATTACHMENT1, gl.TEXTURE_2D, velocityTexture, 0);
  // Multiple Render Targetsを利用できるようにする
  gl.drawBuffers([gl.COLOR_ATTACHMENT0, gl.COLOR_ATTACHMENT1]);
  gl.bindFramebuffer(gl.FRAMEBUFFER, null);
  gl.bindTexture(gl.TEXTURE_2D, null);
  return {
    framebuffer: framebuffer,
    positionTexture: positionTexture,
    velocityTexture: velocityTexture,
  };
}

//　テクスチャの(x, y)座標はx番目のtrailのy番目の頂点の情報を表す
const trailTextureSize = new Vector2(trailNum, vertexNum);
// 一つ前の位置や速度をもとに次の位置や速度を更新するため、フレームバッファーを2つ作成してスワップできるようにしておく
let trailFbObjR = createTrailFbObj(gl, trailTextureSize);
let trailFbObjW = createTrailFbObj(gl, trailTextureSize);
const swapTrailFbObj = () => {
  const tmp = trailFbObjR;
  trailFbObjR = trailFbObjW;
  trailFbObjW = tmp;
};
```

作成した位置テクスチャと速度テクスチャを以下のフラグメントシェーダーで適当に初期化しておきます。初期状態は各trailごとにすべての頂点で同じ位置になるようにしています。

```glsl
#version 300 es

precision highp float;

#define PI 3.14159265359

layout (location = 0) out vec3 o_position;
layout (location = 1) out vec3 o_velocity;

float random(float x){
  return fract(sin(x * 12.9898) * 43758.5453);
}

vec3 randomInSphere(float v) {
  float z = random(v * 0.42 + 213.23) * 2.0 - 1.0;
  float phi = random(v * 0.19 + 313.98) * PI * 2.0;
  float r = random(v * 0.35 + 192.75);

  float a = pow(r, 1.0 / 3.0);
  float b = sqrt(1.0 - z * z);

  return vec3(a * b * cos(phi), a * b * sin(phi), a * z);
}

void main(void) {
  // trailの番号をランダムシードとして与えることで同じtrailではすべての頂点位置を同じにしている
  o_position = 1.0 * randomInSphere(gl_FragCoord.x);
  o_velocity = vec3(0.0);
}
```

ここまでで初期化処理が完了したので、ここからはレンダリングループ内で行います。

以下のフラグメントシェーダーで各trailの位置と速度を更新します。各trailの先頭の頂点はカールノイズのベクトル場から加速度を求めて位置と速度を更新しています。各trailの先頭の頂点以外は一つ前方にある頂点の情報をそのままコピーしておきます。

```glsl
#version 300 es

precision highp float;

layout (location = 0) out vec3 o_position;
layout (location = 1) out vec3 o_velocity;

uniform sampler2D u_positionTexture;
uniform sampler2D u_velocityTexture;
uniform float u_time;
uniform float u_deltaTime;
uniform float u_maxSpeed;
uniform float u_maxForce;
uniform float u_boundaryRadius;
uniform float u_noiseScale;

...

// 3次元のカールノイズ
vec3 curlnoise(vec3 x, float time) {
  // 実装は省略
}

vec3 limit(vec3 v, float max) {
  if (length(v) > max) {
    return normalize(v) * max;
  }
  return v;
}

void main(void) {
  ivec2 coord = ivec2(gl_FragCoord.xy);

  vec3 nextPosition, nextVelocity;
  if (coord.y == 0) { // 各trailの先頭の頂点は一つ前の位置と速度をもとに新しい位置と速度を計算
    vec3 velocity = texelFetch(u_velocityTexture, coord, 0).xyz;
    vec3 position = texelFetch(u_positionTexture, coord, 0).xyz;

    // 現在位置のカールノイズの値から加速度を求める
    vec3 acceleration = curlnoise(position * u_noiseScale, u_time * 0.2);
　
    // 遠くに飛んでいかないように一定範囲外ならば中心に戻るように加速度を加える
    acceleration = u_maxForce * mix(acceleration, -normalize(position), smoothstep(u_boundaryRadius, u_boundaryRadius * 1.05, length(position)));

    nextVelocity = limit(velocity + u_deltaTime * acceleration, u_maxSpeed);
    nextPosition = position + u_deltaTime * nextVelocity;
  } else { // 先頭以外の頂点は一つ先の頂点の情報を移しておく
    nextPosition = texelFetch(u_positionTexture, ivec2(coord.x, coord.y - 1), 0).xyz;
    nextVelocity = texelFetch(u_velocityTexture, ivec2(coord.x, coord.y - 1), 0).xyz;
  }
  o_position = nextPosition;
  o_velocity = nextVelocity;
}
```

trailをGeometry Instancingでレンダリングします。レンダリングに必要な頂点の情報は位置テクスチャに含まれているので`gl.bindBuffer`や`gl.bindVertexArray`で渡す情報はありません。頂点の数が`vertexNum`個のラインを`trailNum`個レンダリングできるように`gl.drawArraysInstanced`を呼び出します。

```js

gl.drawArraysInstanced(gl.LINE_STRIP, 0, vertexNum, trailNum);
```

シェーダー側では`gl_InstancedID`にtrailの番号が、`gl_VertexID`に頂点の番号が含まれているので、それをもとに位置テクスチャから頂点位置を取得します。

```glsl
#version 300 es

uniform sampler2D u_positionTexture;
uniform mat4 u_vpMatrix;

void main(void) {
  vec3 position = texelFetch(u_positionTexture, ivec2(gl_InstanceID, gl_VertexID), 0).xyz;
  gl_Position = u_vpMatrix * vec4(position, 1.0);
}
```

---

単純なラインだとボリューム感がないので次に平面メッシュをtrailにしてみます。

https://aadebdeb.github.io/WebGL2_GPGPU_Trails/plane
![plane.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/85f74368-79c4-eca7-2198-a11b3a0b4565.png)

まずは平面メッシュをJavaScript側で作成しておきます。通常のメッシュの場合は各頂点についてオブジェクト座標系の位置と法線を格納する配列があれば十分ですが、今回はそれ以外にメッシュの各頂点をどのtrailの頂点と紐づけるかという情報を保持する配列も作成しておきます。レンダリング時にはこの情報をもとに平行移動行列を作成し、ワールド座標系の位置に変換します。

```js

function createPlaneMesh(planeWidth, trailVertexNum) {
  const vertexNum = 2 * trailVertexNum;
  const indexNum = 6 * (trailVertexNum - 1);

  const indices = new Int16Array(indexNum); 
  const trailVertices = new Float32Array(vertexNum); // 紐づくtrailの頂点の番号を格納する配列
  const positions = new Float32Array(3 * vertexNum); // 位置を格納する配列
  const normals = new Float32Array(3 * vertexNum); // 法線を格納する配列

  const halfWidth = planeWidth * 0.5;

  let posCount = 0;
  let normalCount = 0;
  for (let i = 0; i < trailVertexNum; i++) {
    // trailの頂点位置をもとにワールド座標系に変換するのでオブジェクト座標系の位置と法線はすべて同じでいい
    posCount = addVertex(positions, posCount, 0.0, halfWidth, 0.0); 
    normalCount = addVertex(normals, normalCount, 1.0, 0.0, 0.0);
    posCount = addVertex(positions,  posCount, 0.0, -halfWidth, 0.0);
    normalCount = addVertex(normals, normalCount, 1.0, 0.0, 0.0);
    trailVertices[2 * i] = trailVertices[2 * i + 1] = i;
  }

  let indexCount = 0;
  for (let i = 0; i < trailVertexNum - 1; i++) {
    indexCount = addQuad(indices, indexCount,
      2 * i, 2 * i + 1, 2 * i + 2, 2 * i + 3
    );
  }

  return {
    indices: indices,
    trailVertices: trailVertices,
    positions: positions,
    normals: normals,
  };
}
```

メッシュの場合は向きがあるために、位置と速度のほかにどちらを上向きにするかという情報を格納するテクスチャも作成しておきます。

```js
function createTrailFbObj(gl, textureSize) {
  const framebuffer = gl.createFramebuffer();
  gl.bindFramebuffer(gl.FRAMEBUFFER, framebuffer);
  ...
  const upTexture = createTexture(gl, textureSize.x, textureSize.y, gl.RGBA32F, gl.RGBA, gl.FLOAT);
  gl.bindTexture(gl.TEXTURE_2D, upTexture);
  gl.framebufferTexture2D(gl.FRAMEBUFFER, gl.COLOR_ATTACHMENT2, gl.TEXTURE_2D, upTexture, 0);
  gl.drawBuffers([gl.COLOR_ATTACHMENT0, gl.COLOR_ATTACHMENT1, gl.COLOR_ATTACHMENT2]);
  ...
  return {
    framebuffer: framebuffer,
    positionTexture: positionTexture,
    velocityTexture: velocityTexture,
    upTexture, upTexture
  };
}
```

テクスチャの初期化時には、上方向も適当に初期化しておきます。

```glsl
#version 300 es

precision highp float;

#define PI 3.14159265359

layout (location = 0) out vec3 o_position;
layout (location = 1) out vec3 o_velocity;
layout (location = 2) out vec3 o_up;

...

void main(void) {
  o_position = 1.0 * randomInSphere(gl_FragCoord.x);
  o_velocity = vec3(0.0);
  o_up = randomOnSphere(gl_FragCoord.x * 0.21 + random(3424.34));
}
```

trailの更新はラインのときと同じように行います。平面メッシュでは上方向も新たに求めます。

```glsl
#version 300 es

precision highp float;

layout (location = 0) out vec3 o_position;
layout (location = 1) out vec3 o_velocity;
layout (location = 2) out vec3 o_up;

uniform sampler2D u_positionTexture;
uniform sampler2D u_velocityTexture;
uniform sampler2D u_upTexture;
...

void main(void) {
  ivec2 coord = ivec2(gl_FragCoord.xy);

  vec3 nextPosition, nextVelocity, nextUp;
  if (coord.y == 0) {
    ...
    // 現在の速度と一つ前の上方向をもとに新たな上方向を求める
    vec3 front = normalize(nextVelocity);
    vec3 right = cross(front, up);
    nextUp = normalize(cross(right, front));
  } else {
    ...
    nextUp = texelFetch(u_upTexture, ivec2(coord.x, coord.y - 1), 0).xyz;
  }
  o_position = nextPosition;
  o_velocity = nextVelocity;
  o_up = nextUp;
}

```

`createPlaneMesh`メソッドで作成した平面メッシュをGeometry Instancingで`trailNum`分だけレンダリングします。

```js
gl.bindVertexArray(meshVao);
gl.drawElementsInstanced(gl.TRIANGLES, mesh.indices.length, gl.UNSIGNED_SHORT, 0, trailNum);
```

レンダリング時の頂点シェーダーでは、各メッシュの頂点について速度の向きと上方向から回転行列を、紐づいたtrailの頂点位置から平行移動行列を求めてモデル行列を作成しています。そのモデル行列でオブジェクト座標の位置と法線を変換することでワールド座標の位置と法線を求めています。

```glsl
#version 300 es

layout (location = 0) in vec3 i_position;
layout (location = 1) in vec3 i_normal;
layout (location = 2) in float i_trailVertex; // 紐づけるtrialの頂点の番号

out vec3 v_normal;

uniform sampler2D u_positionTexture;
uniform sampler2D u_velocityTexture;
uniform sampler2D u_upTexture;
uniform mat4 u_vpMatrix;

mat4 getLookMat(vec3 front, vec3 up) {
  vec3 z = -normalize(front);
  vec3 y = up;
  vec3 x = cross(y, z);

  return mat4(
    x.x, x.y, x.z, 0.0,
    y.x, y.y, y.z, 0.0,
    z.x, z.y, z.z, 0.0,
    0.0, 0.0, 0.0, 1.0
  );
}

void main(void) {
  vec3 instancePosition = texelFetch(u_positionTexture, ivec2(gl_InstanceID, int(i_trailVertex)), 0).xyz;
  vec3 velocity = texelFetch(u_velocityTexture, ivec2(gl_InstanceID, int(i_trailVertex)), 0).xyz;
  vec3 up = texelFetch(u_upTexture, ivec2(gl_InstanceID, int(i_trailVertex)), 0).xyz;

  // 現在の速度の方向と上方向から回転行列を求める
  mat4 lookMat = getLookMat(normalize(velocity), up);

  // 紐づいたtrailの頂点位置と回転行列からモデル行列を求める
  mat4 modelMatrix = mat4(
    1.0, 0.0, 0.0, 0.0,
    0.0, 1.0, 0.0, 0.0,
    0.0, 0.0, 1.0, 0.0,
    instancePosition.x, instancePosition.y, instancePosition.z, 1.0
  ) * lookMat;

  mat4 mvpMatrix = u_vpMatrix * modelMatrix;
  v_normal = (lookMat * vec4(i_normal, 0.0)).xyz;
  gl_Position = mvpMatrix * vec4(i_position, 1.0);
}
```

---

メッシュを平面から円柱にしてみます。円柱になるとクオリティも一気にあがりますが、ポリゴン数も増えるので負荷も一気にあがります。

https://aadebdeb.github.io/WebGL2_GPGPU_Trails/cylinder
![cylinder.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/f08aa182-0623-055b-dfb3-872ba21f9ef9.png)

円柱のGPGPU trailsと平面のGPGPU Trailsではメッシュの作り方しか基本的には違いがありません。なのでメッシュ作成部分のJavaScriptだけをここには載せておきます。平面の場合と同じようにレンダリング時の頂点シェーダーで紐づいたtrailの頂点位置をもとにモデル行列を作るので、オブジェクト座標系ではXY平面に頂点を配置しています。

```js

function addVertex(vertices, vi, x, y, z) {
  vertices[vi++] = x;
  vertices[vi++] = y;
  vertices[vi++] = z;
  return vi;
};  

function addTriangle(indices, i, v0, v1, v2) {
  indices[i++] = v0;
  indices[i++] = v1;
  indices[i++] = v2;
  return i;
};

function addQuad(indices, i, v00, v10, v01, v11) {
  indices[i] = v00;
  indices[i + 1] = indices[i + 5] = v10;
  indices[i + 2] = indices[i + 4] = v01;
  indices[i + 3] = v11;
  return i + 6;
};  

function createCylinderMesh(radius, angleSegment, trailVertexNum) {
  const vertexNum = (1 + angleSegment) * 2 + trailVertexNum * angleSegment;
  const indexNum = 3 * (angleSegment * 2 + (trailVertexNum - 1) * angleSegment * 2);

  const indices = new Int16Array(indexNum);
  const trailVertices = new Float32Array(vertexNum);
  const positions = new Float32Array(3 * vertexNum);
  const normals = new Float32Array(3 * vertexNum);

  const angleStep = 2.0 * Math.PI / angleSegment;

  let posCount = 0;
  let normalCount = 0;
  let trailVertexCount = 0;

  posCount = addVertex(positions, posCount, 0.0, 0.0, 0.0);
  normalCount = addVertex(normals, normalCount, 0.0, 0.0, -1.0);
  trailVertices[trailVertexCount++] = 0;
  for (let ti = 0; ti < trailVertexNum + 2; ti++) {
    for (let ai = 0; ai < angleSegment; ai++) {
      const angle = ai * angleStep + Math.PI * 0.5;
      const position = new Vector3(radius * Math.cos(angle), radius * Math.sin(angle), 0.0);
      posCount = addVertex(positions, posCount, position.x, position.y, position.z);
      if (ti === 0) { // 前のcap部分
        normalCount = addVertex(normals, normalCount, 0.0, 0.0, -1.0);
        trailVertices[trailVertexCount++] = 0;
      } else if (ti === trailVertexNum + 1) { // 胴体部分
        normalCount = addVertex(normals, normalCount, 0.0, 0.0, 1.0);
        trailVertices[trailVertexCount++] = trailVertexNum - 1;
      } else { // 後ろのcap
        const normal = Vector3.norm(position);
        normalCount = addVertex(normals, normalCount, normal.x, normal.y, 0.0);
        trailVertices[trailVertexCount++] = ti - 1;
      }
    }
  }
  posCount = addVertex(positions, posCount, 0.0, 0.0, 0.0);
  normalCount = addVertex(normals, normalCount, 0.0, 0.0, 1.0);
  trailVertices[trailVertexCount++] = trailVertexNum - 1;

  let indexCount = 0;
  for (let ai = 0; ai < angleSegment; ai++) { // 前のcap部分のインデックス
    const aj = ai !== angleSegment - 1 ? ai + 1 : 0;
    indexCount = addTriangle(indices, indexCount, 0, ai + 1, aj + 1);
  };
  let vertexOffset = angleSegment + 1;
  for (let ti = 0; ti < trailVertexNum - 1; ti++) { // 胴体部分のインデックス
    for (let ai = 0; ai < angleSegment; ai++) {
      const aj = ai !== angleSegment - 1 ? ai + 1 : 0;
      const tj = ti + 1;
      const v00 = ai + ti * angleSegment + vertexOffset;
      const v10 = aj + ti * angleSegment + vertexOffset;
      const v01 = ai + tj * angleSegment + vertexOffset;
      const v11 = aj + tj * angleSegment + vertexOffset;
      indexCount = addQuad(indices, indexCount, v00, v01, v10, v11);
    }
  }
  vertexOffset += angleSegment * trailVertexNum;
  for (let ai = 0; ai < angleSegment; ai++) { // 後ろのcap部分のインデックス
    const aj = ai !== angleSegment - 1 ? ai + 1 : 0;
    indexCount = addTriangle(indices, indexCount, vertexNum - 1, aj + vertexOffset, ai + vertexOffset);
  }

  return {
    indices: indices,
    trailVertices: trailVertices,
    positions: positions,
    normals: normals
  };
}
```

---

というわけでGPGPU Trailsを試してみました。円柱までやるとシェーディングができてかっこいい感じになります。
![cylinder2.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61569/aa0aa7f9-a5e4-c38f-0f70-57bcb60c2fc3.gif)

trailの位置と速度の更新自体はそこまで複雑ではないのでGPGPUパーティクルをやったことがあれば理解できるのではないでしょうか。本当に大変なのは、JavaScript側でメッシュを作るところな気がしています。