CUDAでMarching Cube法

Marching Cube法*1 は流体シミュレーションなどで得られた陰関数曲面から 三角形メッシュを作成する手法である.手順は以下.

  1. 立方体グリッド(Cube)を空間に配置
  2. 各Cubeの辺上で陰関数値が0となる点を線型補間や二分法,ニュートン法などで算出
  3. Cubeに含まれる点の配置と数から三角形メッシュを生成

より詳しい手順やソースコードは以下のサイトを参考に.

CUDAのSDK内に含まれるサンプルにMC法のCUDA実装がある. ここではそれについて解説する.

サンプルの場所

MC法のCUDAサンプルは

(SDKインストールフォルダ)\C\src\marchingCubes

にあります.Visual Studioの場合,(SDKフォルダ)\C\src以下にRelease.slnがあるので, それを開くと全サンプルのプロジェクトを含んだソリューションが開かれる. その中の"marchingCubes"プロジェクトがMC法のサンプルである.

ファイル構成

  • defines.h : marchingCubes_kernel.cuでインクルードされる.各種定義など
  • marchingCubes.cpp : メインファイル
  • marchingCubes_kernel.cu : MC法カーネル,デバイス関数など(カーネルを呼び出す関数(launch_*)も含まれる)
  • tables.h : Cube内のメッシュ構成判定テーブルなど

これらの他に,FBOの管理などを行うrendercheck_gl.h, rendercheck_gl.cpp, GLSLコードのコンパイルなどを行う nvShaderUtils.h などが用いられている. これらのファイルは (SDKフォルダ)\C\common\inc\ にヘッダがある. また,Prefix sum(Scan)を行うために,CUDPPライブラリ(CUDA Data Parallel Primitives library)を用いている(cudpp.h).

元データ

メッシュ化する陰関数値はこのサンプルではデバイス関数として,もしくは,サンプルボリュームとして与えられる. これらは define.h の

#define SAMPLE_VOLUME 1

で切り替えられる.

処理の流れ

CUDAによるメッシュ化は,marchingCubes.cpp内の computeIsosurface関数で行われる. 手順としては,

  1. グリッドごとにCube 8頂点の関数値を調べ,Cubeごとのメッシュ頂点数を求める. この処理を行うカーネル関数は classifyVoxel で,カーネルスレッド数はグリッド数(numVoxels)である. カーネル関数内の処理手順は,
    1. グリッド8頂点の関数値を field[8] に格納
    2. 8頂点の関数値が閾値より大きいか,小さいかを各ビットの値とした整数値cubeindexを求める.
    3. テーブル numVertsTable (テクスチャ numVertsTex) をcubeindexで参照し,頂点数を求め,uint型配列 d_voxelVerts に格納.
    4. uint型配列 d_voxelOccupied にボクセル内に頂点があれば1, なければ0を格納する.
  2. (空のボクセルを以降の処理でスキップする場合のみ)グリッド占有の状態 d_voxelOccupied をPrefix sum (Scan)して (結果は d_voxelOccupiedScan),メッシュを生成するべきグリッドのみがつめて格納されたuint型配列 d_compVoxelArray を生成する.ScanにはcuDPPライブラリのcudppScan関数を用いる. cudppScan関数によりScanされた結果 d_voxelOccupiedScan の最後の値は, exclusive scanを用いているので要素の数をnとすると0からn-2番目までの要素の合計になっている. これに d_voxelOccupied 配列の最後の値(n-1番目の値)を足すことでメッシュを生成すべきグリッド数 activeVoxels を求める. そして,d_voxelOccupiedScanに従い,空でないグリッドについて,
    d_compVoxelArray[ d_voxelOccupiedScan[i] ] = i;
    の処理を行うカーネル関数 compactVoxels を numVoxels 個のスレッドで実行する.
  3. 手順2と同様にして,グリッドの頂点数 d_voxelVerts についてもScanし(d_voxelVertsScan), 最後の要素を参照することで,総頂点数 totalVerts を得る.
  4. 結果の頂点位置,法線をOpenGLのVBOに格納するために,cudaGLMapBufferObject()によりVBOをデバイスメモリ上にマップする
  5. (空でない)グリッドごとにメッシュ頂点座標,頂点の並びを計算する. この処理を行うカーネル関数は陰関数の参照に関数を用いる場合は generateTriangles , サンプルボリュームを用いる場合は,generateTriangles2 で,カーネルスレッド数は手順2を行った場合は, 空でないグリッド数(activeVoxels),そうでない場合は,グリッド数(numVoxels)である. カーネル関数内の処理は,
    1. 各グリッドの8頂点座標 v[8] を参照
    2. 1と同様に8頂点の関数値 field[8] を求め,テーブル参照用の整数値cubeindexを算出
    3. v と field の値から12辺について陰関数値の線型補間によりメッシュ頂点の位置を計算,vertlist[12]に格納
    4. グリッド内のメッシュ数は頂点数/3になるので,cubeindexから頂点数 numVerts を求め,メッシュごとに以下の処理を繰り返す.
      1. 各メッシュがどの頂点から構成されるかの情報をテーブル triTable (テクスチャ triTex)から取得(edge)
      2. vertlist[edge]を参照することで3頂点座標を取得(v[3])
      3. 3頂点座標より外積で法線を算出(calcNormal関数)
      4. 結果を配列 pos, normにそれぞれ格納
  6. cudaGLUnmapBufferObject()によりVBOのマップを外す.

このサンプルでは,メッシュ頂点座標をメッシュごとに並べたものを生成している. これはglDrawArrays で描画する場合は便利ではあるが, メッシュを保存してモデリングしたり,レンダリングする際には頂点座標とその幾何構造(接続情報)で構成されるメッシュの方が用いやすい. これについては別のページで実装を述べる. また,環境にもよるが,cudaGLMapBufferObject, cudaGLUnmapBufferObject関数の処理が重く, あらかじめ確保していたデバイスメモリ上に直接結果を書き込み,cudaMemcpyでデバイスからホストメモリにコピーした方が高速な場合もある.

Prefix sum (Scan)とは

大きさnのある配列aがあったときScanは以下のようになる.

Inclusive scan : [a[0], a[0]+a[1], a[0]+a[1]+a[2], ... , a[0]+a[1]+a[2]+...+a[n-2]+a[n-1]]
Exclusive scan : [0, a[0], a[0]+a[1], a[0]+a[1]+a[2], ... , a[0]+a[1]+a[2]+...+a[n-2]]

Scanはデータパラレルプログラミングで良く用いられる. より詳しくは,

を参照.また,CUDAによる並列Scanについては以下参照.


*1 Lorensen, W. E. and Cline, H. E. "Marching cubes: a high resolution 3D surface construction algorithm", Computer Graphics 21 (Proc. SIGGRAPH87), pp.163-169, 1987.

トップ   編集 凍結 差分 履歴 添付 複製 名前変更 リロード   新規 一覧 検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2011-10-27 (木) 15:09:13