09 Sep 2009

Marching Cubes 作業メモ

研究で3次元の物体データの表面をメッシュ化したい、かつそれを可視化したいという要望がでてきた。具体的には、距離カメラでとった距離データをメッシュ化したい。距離データはカメラ座標系でのxyz座標として与えられる。

ちなみに使っている距離カメラはこれ、

Mesa Imaging AG - SwissRanger SR4000 - miniature 3D time-of-flight range camera. 3D Camera

Marching Cubes

Marching Cubesとは、ボリュームデータの表面をレンダリングする手法の一種。アルゴリズムがシンプルで理解しやすく、そのわりに結果が良いらしい。アルゴリズムの解説はここがわかりやすかった。

メタボールを作る2

要は、データ全体を格子状に区切り、格子内の状態の全パターンについて、その結果（どうレンダリングするか）を準備しておく。あとはこのルールに従ってデータ全体を調べればOKということらしい。

インプットが濃度の分布であれば良いので、3次元物体だけでなく、色々なデータの可視化に使えそう。またレンダリングの際のしきい値をかえることで、濃度値が同じ面が抽出されるので、何かの濃度の分布を等高線みたいに表示させることもたぶんできそう。

今回の研究の用途では、距離カメラでデータをとっているので、濃度は1か0のバイナリになる。

Polygonising a scalar field

marching cubeで検索したら2番目にヒットしたこのページ。書かれたのは15年前だけど、よく見られているようだし（Deliciousでもたくさんブックマークされている）、サンプルプログラムもあるので、ここを参考にすることに。

一番上にあるmarchingsource.cppというファイル。なんとなくコンパイルしてみたら、何となく動いた。しかし、OpenGLが全くわからないので、コードがよくわからない。特にどこがインプットのデータなのかがいまいちわからなかった（どうも自動で生成してるみたいだったけど）。OpenGLを勉強するのもめんどうだったので、いったんパス。

次に見たのはその下にあるrchandra.zipという実装。こちらは上のmarchingsource.cppを参考に作った CIsoSurface という名前のクラスで、入力の生成や表示系はなく、アルゴリズム部だけだったので、わりかし見やすい。

サーフェイスを生成するメインの処理を行うのは、GenerateSurface()という関数。引数は8つで、

ptScalarField: ボリュームデータの配列
tIsoLevel: レンダリングの際の閾値
nCellsX, nCellsY, nCellsZ: ボリュームデータのX, Y, Zそれぞれの要素数
fCellLengthX, fCellLengthY, fCellLengthZ: ボリュームデータの各格子間の距離。サーフェイスの生成には直接は使われない。

ptScalarFieldは3次元配列（XYZ要素に対応）にして、その中の各要素はその点の濃度（今回は物体があるかないかなのでバイナリ）。nCell*はptScalarFieldのXYZのそれぞれの大きさにすれば良い。fCellLength*はあとでカメラの仕様を調べてそれを入力しよう。

計算結果は、主に m_ppt3dVertices, m_piTriangleIndices, m_pvec3dNormalsの3つ。それぞれ、

m_ppt3dVertices: 生成されたサーフェイスを構成する頂点のリスト
m_piTriangleIndices: サーフェイスを構成する三角形の点を、上の m_ppt3dVertices のインデックスで表したもの。[vertex1 of triangle1], [vertex2 of triangle1], [vertex3 of triangle1], [vertex1 of triangle2], ... という順番で並んでいる。
m_pvec3dNormals: たぶん各頂点の法線ベクトル。

今回STL形式で結果を出力しようと思っているので、法線ベクトルは頂点のものじゃなくて三角形のものが欲しい。なので、ちょっとコードに手を入れることに。とは言っても、 m_pvec3dNormals の計算過程で三角形の法線ベクトルは計算されているので、新しくベクタを作ってそれを格納するだけ。

STL形式

ところで、STL形式というのは、3次元の形状を表現するデータ形式の一種。C++のSTLとは関係ない。こちらを参考にした。

STLファイルフォーマット

アスキーとバイナリの2種類があり、各三角形の頂点と法線ベクトルを、フォーマットにのっとって列挙しただけというわかりやすいもの。

STL形式のアスキーフォーマットで出力

上記の CIsoSurface クラスに、計算結果をSTLのアスキーフォーマットで出力する printSTLAscii() という関数を追加した。とはいっても、ただ配列の中身を形式に沿ってプリントするだけなので簡単。

ここまでの修正のパッチはこんな感じ。

diff -u /Users/kosei/Downloads/rchandra/CIsoSurface.cpp ../rchandra/CIsoSurface.cpp
--- /Users/kosei/Downloads/rchandra/CIsoSurface.cpp	2007-08-28 21:10:17.000000000 +0900
+++ ../rchandra/CIsoSurface.cpp	2009-09-08 23:15:44.000000000 +0900
@@ -6,10 +6,14 @@
 //
 // Description: This is the implementation file for the CIsoSurface class.
 
-#include "stdafx.h"
 #include 
+#include 
+#include 
+#include 
 #include "CIsoSurface.h"
 
+using namespace std;
+
 template  const unsigned int CIsoSurface::m_edgeTable[256] = {
 	0x0  , 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c,
 	0x80c, 0x905, 0xa0f, 0xb06, 0xc0a, 0xd03, 0xe09, 0xf00,
@@ -681,7 +685,7 @@
 	vecIterator = m_trivecTriangles.begin();
 	m_nTriangles = m_trivecTriangles.size();
 	m_piTriangleIndices = new unsigned int[m_nTriangles*3];
-	for (i = 0; i < m_nTriangles; i++, vecIterator++) {
+	for (unsigned int i = 0; i < m_nTriangles; i++, vecIterator++) {
 		m_piTriangleIndices[i*3] = (*vecIterator).pointID[0];
 		m_piTriangleIndices[i*3+1] = (*vecIterator).pointID[1];
 		m_piTriangleIndices[i*3+2] = (*vecIterator).pointID[2];
@@ -704,8 +708,9 @@
 	}
 
 	// Calculate normals.
-	for (i = 0; i < m_nTriangles; i++) {
-		VECTOR3D vec1, vec2, normal;
+	for (unsigned int i = 0; i < m_nTriangles; i++) {
+		VECTOR3D vec1, vec2;
+		vector  normal(3);
 		unsigned int id0, id1, id2;
 		id0 = m_piTriangleIndices[i*3];
 		id1 = m_piTriangleIndices[i*3+1];
@@ -728,10 +733,11 @@
 		m_pvec3dNormals[id2][0] += normal[0];
 		m_pvec3dNormals[id2][1] += normal[1];
 		m_pvec3dNormals[id2][2] += normal[2];
+		triNormals.push_back(normal);
 	}
 
 	// Normalize normals.
-	for (i = 0; i < m_nNormals; i++) {
+	for (unsigned i = 0; i < m_nNormals; i++) {
 		float length = sqrt(m_pvec3dNormals[i][0]*m_pvec3dNormals[i][0] + m_pvec3dNormals[i][1]*m_pvec3dNormals[i][1] + m_pvec3dNormals[i][2]*m_pvec3dNormals[i][2]);
 		m_pvec3dNormals[i][0] /= length;
 		m_pvec3dNormals[i][1] /= length;
@@ -739,6 +745,43 @@
 	}
 }
 
+template  int CIsoSurface::printSTLAscii()
+{
+  string solid = "CIsoSurface";
+
+  if (!IsSurfaceValid())
+    {
+      cout << "Isosurface is not velid" << endl;
+      return -1;
+    }
+
+  cout << "solid " << solid << endl;;
+
+  for (int i=0; i;
 template class CIsoSurface;
 template class CIsoSurface;
diff -u /Users/kosei/Downloads/rchandra/CIsoSurface.h ../rchandra/CIsoSurface.h
--- /Users/kosei/Downloads/rchandra/CIsoSurface.h	2007-08-28 21:10:17.000000000 +0900
+++ ../rchandra/CIsoSurface.h	2009-09-08 22:02:44.000000000 +0900
@@ -14,6 +14,8 @@
 #include 
 #include "Vectors.h"
 
+using namespace std;
+
 struct POINT3DID {
 	unsigned int newID;
 	float x, y, z;
@@ -48,6 +50,9 @@
 	// valid.
 	int GetVolumeLengths(float& fVolLengthX, float& fVolLengthY, float& fVolLengthZ);
 
+	// Print isosurface as a STL Ascii format
+	int printSTLAscii();
+
 protected:
 	// The number of vertices which make up the isosurface.
 	unsigned int m_nVertices;
@@ -112,6 +117,9 @@
 	// Lookup tables used in the construction of the isosurface.
 	static const unsigned int m_edgeTable[256];
 	static const unsigned int m_triTable[256][16];
+
+	// Normals of each triangles
+	vector  > triNormals;
 };
 #endif // CISOSURFACE_H
 
diff -u /Users/kosei/Downloads/rchandra/Vectors.cpp ../rchandra/Vectors.cpp
--- /Users/kosei/Downloads/rchandra/Vectors.cpp	2007-08-28 21:10:17.000000000 +0900
+++ ../rchandra/Vectors.cpp	2009-09-08 23:13:09.000000000 +0900
@@ -5,7 +5,6 @@
 //
 // Description: This is the implementation file for POINT3DXYZ class.
 
-#include "stdafx.h"
 #include "Vectors.h"
 
 POINT3DXYZ operator+(const POINT3DXYZ& pt3dPoint1, const POINT3DXYZ& pt3dPoint2)

使い方はこんな感じ。

#include "CIsoSurface.h"

int main(int argc, char ** argv)
{
  CIsoSurface <short> *ciso = new CIsoSurface <short> ();

  short volume[3][3][3] = {
    0, 1, 1, 1, 1, 0, 0, 0, 0,
    0, 1, 1, 1, 1, 0, 0, 0, 0,
    0, 1, 1, 1, 1, 0, 0, 0, 0,
  };

  ciso->GenerateSurface(&volume[0][0][0], 1, 3, 3, 3, 1, 1, 1);
  ciso->printSTLAscii();

  return 0;
}