CUDAでMC法の変更点

追加された行はこの色です。
削除された行はこの色です。
CUDAでMC法へ行く。
*CUDAでMarching Cube法 [#g5e1f9e1]
Marching Cube法((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.))
は流体シミュレーションなどで得られた陰関数曲面から
三角形メッシュを作成する手法である．手順は以下．

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

より詳しい手順やソースコードは以下のサイトを参考に．
-http://local.wasp.uwa.edu.au/~pbourke/geometry/polygonise/

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

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

**ファイル構成 [#d9e62cf1]
-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:http://gpgpu.org/developer/cudpp]]ライブラリ(CUDA Data Parallel Primitives library)を用いている(cudpp.h)．

**元データ [#k0dea681]
メッシュ化する陰関数値はこのサンプルではデバイス関数として，もしくは，サンプルボリュームとして与えられる．
これらは define.h の
 #define SAMPLE_VOLUME 1
で切り替えられる．

**処理の流れ [#m856ec77]
CUDAによるメッシュ化は，marchingCubes.cpp内の computeIsosurface関数で行われる．
手順としては，
+グリッドごとにCube 8頂点の関数値を調べ，Cubeごとのメッシュ頂点数を求める．
この処理を行うカーネル関数は classifyVoxel で，カーネルスレッド数はグリッド数(numVoxels)である．
カーネル関数内の処理手順は，
++グリッド8頂点の関数値を field[8] に格納
++8頂点の関数値が閾値より大きいか，小さいかを各ビットの値とした整数値cubeindexを求める．
++テーブル numVertsTable (テクスチャ numVertsTex) をcubeindexで参照し，頂点数を求め，uint型配列 d_voxelVerts に格納．
++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)である．
カーネル関数内の処理は，
++各グリッドの8頂点座標 v[8] を参照
++1と同様に8頂点の関数値 field[8] を求め，テーブル参照用の整数値cubeindexを算出
++v と field の値から12辺について陰関数値の線型補間によりメッシュ頂点の位置を計算，vertlist[12]に格納
++グリッド内のメッシュ数は頂点数/3になるので，cubeindexから頂点数 numVerts を求め，メッシュごとに以下の処理を繰り返す．
+++各メッシュがどの頂点から構成されるかの情報をテーブル triTable (テクスチャ triTex)から取得(edge)
+++vertlist[edge]を参照することで3頂点座標を取得(v[3])
+++3頂点座標より外積で法線を算出(calcNormal関数)
+++結果を配列 pos, normにそれぞれ格納
+cudaGLUnmapBufferObject()によりVBOのマップを外す．

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

**Prefix sum (Scan)とは [#k4d4640f]
大きさ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:http://en.wikipedia.org/wiki/Prefix_sum]](Wikipadia)

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

-[[Parallel Prefix Sum (Scan) with CUDA:http://developer.download.nvidia.com/compute/cuda/1_1/Website/projects/scan/doc/scan.pdf]](注 PDF)
CUDAでMC法 の変更点
