Display 3D Reaction Diffusion with Marching Cubes

　この記事は CAMPHOR- Advent Calendar 2019 24日目の記事です．23日目の記事は 5ebec の 誰得なAtomのSyntax themeを作った でした．今回の記事はMarching Cubesという手法に関して，研究ノートをまとめる雰囲気で書いています．万人受けする内容にはなっていませんが，Marching Cubesを触る人の知見になれば幸いです．

About Marching Cubes

　8点からなる立方体（ボクセル）を1つの単位として，各点のスカラー値から等値面*1を作成しポリゴンメッシュを生成するアルゴリズムである．LorensenとClineによって1987年にSIGGRAPHで発表された．主な使用用途に，CT，MRIなどで得られる断層画像の明度をスカラー値として，撮影物の立体形状を再構築することが挙げられる．

Motivation

　3Dの反応拡散を可視化したい．反応拡散は，1種または複数種の作用因子を相互作用させて自律的にパターンをつくる数理モデルである．詳しくはWikipediaなど参照．今回は合成因子・抑制因子2つの相互作用を考えている．はじめは2次元の計算を3次元に書き換えて動くか確かめるためだけにボクセルに色をつけて表示したのだが，この表示方法だと内部構造が見えない（図，黒が濃いほど合成因子の濃度が高い）．そこで，合成因子の濃度が高い領域だけをポリゴンメッシュ化してより見やすくなるようにしたかった．

3次元の反応拡散（Processingで実装）

Algorithm of Marching Cubes

　各ボクセルを構成する8点にはそれぞれ0もしくは1の値が付与されているとする．データのスカラー値が連続型だった場合，閾値を設定し，2値化することで解決する．各点のもつ値の組み合わせによって，張られる等値面の位置・形状は変わってくる．単純計算で2⁸=256通りが考えられるが，Marching Cubesでは，等値面がない場合も含めると等値面は下表のように15種類と決まっている．これは，0と1の分布を逆転させても同じ等値面が張られることから256/2=128通りになり，さらに回転対称性を考慮しているためである．表に示されている等値面を見てみると，例えば，全点で0もしくは1ならば等値面は張られない．1点だけ他の7点と値が違う場合，その点を隔てるようにして三角形の等値面が一枚張られる．等値面を構成する点の位置は二値化されている場合辺の中点となるが，スカラー値が連続型の場合，それを重みとして点の位置が辺の中点からずれる．断層画像を使ってポリゴンメッシュを生成する場合は，隣接した2層の間に立方体が挟まっていると考えれば良い．

The originally published 15 cube configurations (Referenced from Wikipedia)

Coding

　「Polygonising a scalar field」でコードの構成がわかりやすく説明されており，それを踏襲した一通りのテンプレートが，同ページに添付されているrchandra.zip内に存在する．既にこのコードを用いてMarching Cubesを試用している人が書いているように，marchingsource.cppはOpenGLでの描画を前提として書かれているため，読みづらいので，rchandraを使用した．

　 　rchandraは，CIsoSurface.cppとVectors.cpp，それに対応する.hファイルから構成される．重要なのはCIsoSurfaceの方で，CIsoSurfaceクラスを定義し，メンバ関数のGenerateSurface()によって等値面を生成する．GenerateSurface()の呼び出しに必要な引数は以下の4つ．

• ptScalarField：スカラー値をもつ点群の配列（xyzに対応した三次元配列）
• tIsoLevel：等値面の閾値
• nCellsX, nCellsY, nCellsZ：xyz方向の立方体要素の数
• fCellLengthX, fCellLengthY, fCellLengthZ：点群が構成する各格子間のxyz方向における距離

　そして，主に以下の出力結果が得られる．

• m_bValidSurface：等値面を生成できたか返す真偽値
• m_ppt3dVertices：三角形メッシュを構成する点のインデックスを格納する配列．点のxyz座標値を保持する要素数3の配列POINT3D型の配列である．m_ppt3dVertices[i][0]，m_ppt3dVertices[i][1]，m_ppt3dVertices[i][2]はそれぞれi番目の点のx,y,z座標を指す．m_ppt3dVerticesの要素数はm_nVerticesが保持している．
• m_trivecTriangles：三角形メッシュを構成する3点のインデックスを全メッシュについて記述する配列．3点のインデックスはTRIANGLE構造体として定義されている．インデックスはm_ppt3dVerticesのインデックスと対応している．m_trivecTrianglesの要素数はm_nTrianglesが保持している．
• m_pvec3dNormals：三角形メッシュを構成する点ごとの法線ベクトルを格納する配列．

　等値面を生成するための最小のコードは以下のようになる．これをmain.cppとしてg++ main.cpp CIsoSurface.cpp Vectors.cppを叩くと，等値面が生成される．C++の実行環境は人によって違うと思うので，コマンドはそれに合わせて適宜変更してほしい．

　注意点として，nCellsX, nCellsY, nCellsZは（xyzの各軸上に並ぶ点の数-1）にする．先程のコードだとpnum-1にあたる．先人のブログだと，この数は点の総数と同値になっていたが，それで実行すると等値面が正しく生成されなかった．

　rchandra内のコードのみだと等値面メッシュのデータが出力されないので，.vtk，.objファイル形式でデータを出力する関数をCIsoSurfaceクラス内に作成した．.vtk，.objファイル形式についての説明は割愛．三角形メッシュを構成する点の座標とインデックスがあれば作成できるので，m_ppt3dVerticesとm_trivecTrianglesを各ファイル形式に合わせてプリントするだけでよい．Paraview，MeshLabなどでメッシュを表示する場合は，法線ベクトルが.vtk，.objファイルに記載されていなくても，インデックスの順番から面の法線の方向を判断して描画してくれるので，今は法線ベクトルの出力をしていない．参考として，.vtkファイル出力用のコードを以下に載せる．

Display Isosurfaces with Paraview

　合成因子/抑制因子の相互作用による反応拡散を計算し，合成因子の濃度をスカラー値としてメッシュを生成した．等値面の閾値は適当に設定している．定常状態になると，迷宮模様の立体構造が現れた（図）．Paraviewでは，Reload Files→Find New Filesで新規生成されたデータファイルを読み込めるので，起動し直さなくても最新の計算結果を描画することが可能である．初期の色彩設定ではメッシュの裏表の明暗が明瞭でないため，形状がわかりにくい．面の法線方向はメッシュの裏表で統一されているので，これを利用して，今回は表面の色を黒に近づけた．方法は，Paraview内のパレットのアイコンをクリックし，Edit Current Paletteを選択する．その中の Color used when solid coloring surfaces and faces. を黒にすると，以下のように空洞部分が黒くなって形状が比較的わかりやすくなる．表裏で白黒を反転させたい場合は，表面は白のまま，メッシュを選択した状態で Properties 下に現れるBackface color を黒にする．アドバイスをくれた @iftianjiar に感謝．

Paraviewで表示した結果

Paraviewでメッシュの表面を黒くした場合

Marching Tetrahedron

　Marching TetrahedronはMarching Cubesの応用である．立方体メッシュを単位として等値面を張ると穴があく場合があるので，それを防ぐために立方体要素をさらに6つの四面体に分割し，各四面体の中でメッシュを張る．Reaction Diffusionのシミュレーションで穴が目立つ場合は導入を考えていたが，表示結果を見たら問題なさそうだったので，着手していない．

References

• Lorensen, W. E.; Cline, Harvey E. (1987). "Marching cubes: A high resolution 3d surface construction algorithm". ACM Computer Graphics. 21 (4): 163–169. doi:10.1145/37402.37422.
• Reaction Diffusion
  • https://en.wikipedia.org/wiki/Reaction%E2%80%93diffusion_system
  • https://ja.wikipedia.org/wiki/%E5%8F%8D%E5%BF%9C%E6%8B%A1%E6%95%A3%E7%B3%BB
• Polygonising a scalar field
  • http://paulbourke.net/geometry/polygonise/
• https://en.wikipedia.org/wiki/Marching_cubes
• http://please-sleep.cou929.nu/20090908.html
• https://en.wikipedia.org/wiki/Marching_tetrahedra

　明日はshiba6vの担当です．例年CAMPHOR- Advent Calendarの最終日には，代表による1年の総括が綴られています．歴代の総括を読むと，CAMPHOR-が年ごとに違うカラーを出しながらも団体の維持・発展に努めてきたことがわかりました．今年はどのような総括になるのでしょうか？お楽しみに！

　CAMPHOR-に興味を持った方はこちらの記事をどうぞ．それでは良いクリスマスを．

*1:同じ値をもつ点の集まりでできる面のこと．等高線を面に拡張したものだと思ってもらえばよい．今回は筆者が設定した閾値と同じ値をもつ等値面のみを表示するようにしている．