CUDAでMC法

CUDAでMarching Cube法†

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

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

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

http://local.wasp.uwa.edu.au/~pbourke/geometry/polygonise/

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関数で行われる．手順としては，

グリッドごとに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を格納する．
(空のボクセルを以降の処理でスキップする場合のみ)グリッド占有の状態 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 個のスレッドで実行する．
手順2と同様にして，グリッドの頂点数 d_voxelVerts についてもScanし(d_voxelVertsScan)，最後の要素を参照することで，総頂点数 totalVerts を得る．
結果の頂点位置，法線をOpenGLのVBOに格納するために，cudaGLMapBufferObject()によりVBOをデバイスメモリ上にマップする
(空でない)グリッドごとにメッシュ頂点座標，頂点の並びを計算する．この処理を行うカーネル関数は陰関数の参照に関数を用いる場合は 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にそれぞれ格納
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はデータパラレルプログラミングで良く用いられる．より詳しくは，

Prefix sum(Wikipadia)

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

Parallel Prefix Sum (Scan) with CUDA(注 PDF)

Menu