WebGLでマーチングキューブ法

マーチングキューブ法はボクセルデータからメッシュを作成する手法です。空間を格子上に分割することで作成された各直方体(セル)について、それを構成する8つの頂点の値に応じて三角形ポリゴンを生成することでメッシュを作成します。
マーチングキューブ法 - Wikipedia

WebGLで実装するにあたり、以下のサイトを参考にしました。このサイトの方法では、1つのセルに対して8つの頂点がとりうる値のパターンは256あるので、そのパターンすべてについて三角形ポリゴンの作成方法をルックップテーブルとしてあらかじめ格納しておきます。各セルについてルックアップテーブルを参照しながら三角形ポリゴンを必要な数だけ生成しています。
Polygonising a scalar field (Marching Cubes)

マーチングキューブ法は対象にするボクセルデータにより必要となる三角形ポリゴンの数が変わります。リアルタイムにメッシュを生成する場合はジオメトリシェーダーを使ってGPU側でポリゴンを作成するのが一般的だと思います。しかし、WebGLにはジオメトリシェーダーがないため、シェーダー側でポリゴンを作成する方法がありません。マーチングキューブ法では各セルで生成される三角形ポリゴンの最大数は5つです。そのため今回のWebGLの実装では、各セルについて5つの三角形を描画するようにドローコールを行い、不要な三角形はシェーダー側でピクセルを破棄して描画しないようにしています。

ソースコード全文はGithubに置いておきました。
aadebdeb/WebGL_MarchingCubes: Real time Marching Cubes in WebGL

まずはじめに、ルックアップテーブルを格納するテクスチャを作成しておきます。シェーダー側にルックアップテーブルを直接定義すると、シェーダーの定数制限に引っかかるのでテクスチャで渡しています。参考にしているサイトにはedgeTableとtriTableの2種類がありますが、edgeTableは今回の実装では使用しなので、triTableに対応したテクスチャのみを作成しています。ArrayBufferViewからテクスチャを作成する方法については私が以前に書いた記事を参考にしてください。

const TRI_TABLE = new Uint32Array([
  -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
  0, 8, 3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
  // ... 長いので省略
  0, 3, 8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
  -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1
]);

function createTriTexture(gl) {
  const triTexture = gl.createTexture();
  gl.bindTexture(gl.TEXTURE_2D, triTexture);
  gl.texImage2D(gl.TEXTURE_2D, 0, gl.R32UI, 4096, 1, 0, gl.RED_INTEGER, gl.UNSIGNED_INT, TRI_TABLE);
  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 triTexture;
}

マーチングキューブ法を実行するためのパラメータとドローコールは以下のようになります。marchingSpaceがマーチングキューブ法を実行する空間全体の大きさ、cellNumがセル数、cellSizeがセルの大きさを表しています。cellNum.x*cellNum.y*cellNum.zが総セル数であり、各セルについて三角形が最大5つ描画される(3*5=15)ので、頂点数vertexNumはcellNum.x*cellNum.y*cellNum.z*15になります。

const marchingSpace = new Vector3(10.0, 10.0, 10.0);
const cellNum = new Vector3(50, 50, 50);
const cellSize = new Vector3(marchingSpace.x / cellNum.x, marchingSpace.y / cellNum.y, marchingSpace.z / cellNum.z);
const vertexNum = cellNum.x * cellNum.y * cellNum.z * 15;
...
gl.drawArrays(gl.TRIANGLES, 0, vertexNum);

マーチングキューブ法を実行する頂点シェーダーは以下のようになります。　アルゴリズム自体は参考にしているサイトと同じです。main関数では、まずgl_VertexIDをもとに現在どのセルのどの三角形の頂点を対象にしているかを決定します。セルの8つの頂点の値から、ルックアップテーブルを参照するためのインデックスを作成します。ルックアップテーブルから取得した値により現在対象にしている三角形の頂点がどの辺にあるかがわかるので、辺の2つの端点の位置と値から頂点位置を決定しています。ルックアップテーブルの値が-1のときは三角形を描画しないということなので、gl_Positionには適当な値を入れておき、フラグメントシェーダーでピクセルを破棄するためのフラグを立てておきます。

#version 300 es

precision highp usampler2D;

out vec3 v_normal; // ワールド座標系における法線
flat out int v_discard; // 1のときピクセルを破棄

uniform usampler2D u_triTexture; // triTableの値を格納するテクスチャ
uniform mat4 u_mvpMatrix;
uniform mat4 u_normalMatrix;
uniform ivec3 u_cellNum; // セル数
uniform vec3 u_cellSize; // セルの大きさ
uniform float u_time;

// 符号付き距離関数を返す関数
// サンプルで実際に使用しているものを載せるのは冗長なため、ここでは単純な球の距離関数を使用している
float getDistance(vec3 p) {
  return length(p) - 2.0;
}

// 法線を返す関数
vec3 getNormal(vec3 p) {
  float e = 0.01;
  return normalize(vec3(
    getDistance(p + vec3(e, 0.0, 0.0)) - getDistance(p - vec3(e, 0.0, 0.0)),
    getDistance(p + vec3(0.0, e, 0.0)) - getDistance(p - vec3(0.0, e, 0.0)),
    getDistance(p + vec3(0.0, 0.0, e)) - getDistance(p - vec3(0.0, 0.0, e))
  ));
}

// v0, v1の値をもとにp0, p1を補間した値を返す
vec3 interpolate(vec3 p0, vec3 p1, float v0, float v1) {
  return mix(p0, p1, -v0 / (v1 - v0));
}

void main(void) {
  int cellId = gl_VertexID / 15; // セルのID
  int vertexId = gl_VertexID % 15; // セル内での頂点のID

  ivec3 cellIdx = ivec3(
    cellId % u_cellNum.x, (cellId % (u_cellNum.x * u_cellNum.y)) / u_cellNum.x, cellId / (u_cellNum.x * u_cellNum.y));
  vec3 cellCorner = (0.5 * vec3(u_cellNum) - vec3(cellIdx)) * u_cellSize;

  // 現在のセルの各頂点位置を求める
  vec3 c0 = cellCorner;
  vec3 c1 = cellCorner + u_cellSize * vec3(1.0, 0.0, 0.0);
  vec3 c2 = cellCorner + u_cellSize * vec3(1.0, 1.0, 0.0);
  vec3 c3 = cellCorner + u_cellSize * vec3(0.0, 1.0, 0.0);
  vec3 c4 = cellCorner + u_cellSize * vec3(0.0, 0.0, 1.0);
  vec3 c5 = cellCorner + u_cellSize * vec3(1.0, 0.0, 1.0);
  vec3 c6 = cellCorner + u_cellSize * vec3(1.0, 1.0, 1.0);
  vec3 c7 = cellCorner + u_cellSize * vec3(0.0, 1.0, 1.0);

  // 現在のセルの各頂点の値を求める
  float v0 = getDistance(c0);
  float v1 = getDistance(c1);
  float v2 = getDistance(c2);
  float v3 = getDistance(c3);
  float v4 = getDistance(c4);
  float v5 = getDistance(c5);
  float v6 = getDistance(c6);
  float v7 = getDistance(c7);

  // セルの各頂点の値からルックアップテーブルを参照するためのインデックスを求める
  int cubeIdx = 0;
  if (v0 < 0.0) cubeIdx |= 1;
  if (v1 < 0.0) cubeIdx |= 2;
  if (v2 < 0.0) cubeIdx |= 4;
  if (v3 < 0.0) cubeIdx |= 8;
  if (v4 < 0.0) cubeIdx |= 16;
  if (v5 < 0.0) cubeIdx |= 32;
  if (v6 < 0.0) cubeIdx |= 64;
  if (v7 < 0.0) cubeIdx |= 128;

  int tri = int(texelFetch(u_triTexture, ivec2(cubeIdx * 16 + vertexId, 0), 0).x);
  if (tri == -1) {
    // -1のとき描画する三角形はないので、ピクセルを破棄するようにする
    gl_Position = vec4(vec3(0.0), 1.0);
    v_discard = 1;
  } else {
    // 描画する三角形があるときは、頂点位置を求める
    vec3 position;
    if (tri == 0) {
      position = interpolate(c0, c1, v0, v1);
    } else if (tri == 1) {
      position = interpolate(c1, c2, v1, v2);
    } else if (tri == 2) {
      position = interpolate(c2, c3, v2, v3);
    } else if (tri == 3) {
      position = interpolate(c3, c0, v3, v0);
    } else if (tri == 4) {
      position = interpolate(c4, c5, v4, v5);
    } else if (tri == 5) {
      position = interpolate(c5, c6, v5, v6);
    } else if (tri == 6) {
      position = interpolate(c6, c7, v6, v7);
    } else if (tri == 7) {
      position = interpolate(c7, c4, v7, v4);
    } else if (tri == 8) {
      position = interpolate(c0, c4, v0, v4);
    } else if (tri == 9) {
      position = interpolate(c1, c5, v1, v5);
    } else if (tri == 10) {
      position = interpolate(c2, c6, v2, v6);
    } else if (tri == 11) {
      position = interpolate(c3, c7, v3, v7);
    }
    gl_Position = u_mvpMatrix * vec4(position, 1.0);
    v_normal = (u_normalMatrix * vec4(getNormal(position), 0.0)).xyz;
    v_discard = 0;
  }
}

フラグメントシェーダーは以下のようになります。v_discard == 1のときはdiscardでピクセルを破棄し、それ以外の場合は普通にシェーディングしています。

#version 300 es

precision highp float;

in vec3 v_normal;
flat in int v_discard;

out vec4 o_color;

const vec3 LIGHT_DIR = normalize(vec3(0.5, 0.5, 1.0));

void main(void) {
  if (v_discard == 1) {
    discard;
  } else {
    vec3 normal = normalize(v_normal);
    vec3 color = vec3(0.3, 0.3, 1.0) * max(0.0, dot(LIGHT_DIR, normal));
    o_color = vec4(color, 1.0);
  }
}

実際に実行してみると以下のようになりました。球やトーラスなど曲面で構成されるオブジェクトの場合はきれいにレンダリングされますが、直方体の角部分では法線が正しく求まらないのかアーティファクトが目立つ結果となっています。セルの数を増やせば目立たなくなりますが、その分負荷も大きくなります。