テセレーションという言葉は、「形状表面を、互いに重ならない複数の三角の集合で隙間なく埋め尽くすときの、個々の三角形の決め方・選び方」を意味する。もっと簡単に「ある形状領域（エリア）を複数の三角形に分割するやり方を決めること」と言ってしまってもいいかも知れない。マーチングキューブズでは、即ち、セル単位で見た各コンフィギュレーションにおける等値面を、複数の三角形で表現するため、その個々の三角の張り方までをも最終的には決めてやらなくてはならない。
ところでテセレーションとは言ったものの、考えなくてはならないヤマが2段階ある。それは
- 各コンフィギュレーションで描くべき形状（ループ）の定義
- 定義したループ形状のテセレーション
となる。
１は、各コンフィギュレーションにおける、等値面形状の輪郭形状を決定することを意味する。各セル単位で見た等値面の輪郭形状は必ず閉じた輪の形状になっているので、これをループと呼んだりもする。２は、１で決めたループ形状内部を、どういう三角形の集まりで表現してやるかという問題である。
想像に難くないかも知れないが、上記の１と２、どちらも結構奥が深い。何故なら、どちらも実は正解が、無数通りとは言わないが、不確定数通りあり得るためである。答えが一意に決まるという意味ではむしろ、前回までにやってきた14通りのグループ分けだとか、相対配置情報の算出とかの方がずっと簡単だ。マーチングキューブズ法はなかなか一筋縄では達成できないアルゴリズムなのである。
●等値面輪郭形状（ループ）の決定
マーチングキューブアルゴリズムの話題としてよく上がるのが、考えている３次元格子の間隔よりも形状が入り組んでいる場合に生じる、不確定性である。つまり一つのセルの中に二つ以上のループが考えられるケースだ。後でわかることだが、そのようなケースは前回までに挙げた14グループの中でも、７グループほど存在するようだ。これらのグループのケースでは、これから求めようとしている等値面輪郭形状の決め方が複数通りある場合があり、故に不確定性として問題視？される。しかしこれはアルゴリズムに本質的なものであり、避けようがない。不確定性が問題となるようなケースでは、サンプリングの格子間隔を細かくするしか結局のところ解決策はない。
等値面を構成する三角の頂点は、格子のそれではなく、然るべき格子エッジ上のどこかに存在する、補間点（エッジポイント）だ。つまりループ形状は、最終的には、エッジポイントの配列として定義されるわけだ。しかしここで問題としたいのはエッジポイントの配列を求めることではなく、エッジポイントが乗っているところのエッジの配列として、ループ形状を求めておきたいのである。つまり、補間点の座標そのものを計算するのではなく、まず補間点はどこのエッジ上で計算すればよいのかということを、各ケース別に求めておきたいのである。（上図）
ここではハーフエッジクアッドメッシュデータ構造としてキューブを再構築し、以下のようなルールに則ってループを決定して行くことにより、14グループの輪郭ループ形状セットを全て自動で一意のもとに決定する。（図も参照）
以下「有効頂点」は閾値よりも低い値を持つ頂点を表し、無効頂点は、閾値以上の値を持つ頂点を表すものとする。
- 有効エッジ（片側に有効頂点、片側に無効頂点を持つエッジ）上に定義された、有効頂点を始点とする任意のハーフエッジをカレントハーフエッジとする
- カレントハーフエッジの乗っているエッジをループに登録
- カレントハーフエッジが定義されているFace上の残りの３つのハーフエッジを、Prevハーフエッジから始めて反時計周りに、探索していく。探索ハーフエッジの乗っているエッジが有効エッジだったなら、それをループに登録
- カレントハーフエッジを上記で登録したエッジ由来の探索ハーフエッジのOppositeを次のカレントハーフエッジにする
- ループが最初のエッジに戻ってくるまで、3を繰り返す。
- ループに登録されていない未処理の有効エッジがなくなるまで、1に戻る。
（なお、３の探索において複数候補がある場合、候補の中にループの先頭のエッジが含まれる場合は、それを優先し、そこでループを完結させ、次なるループの作成にうつる（これがために、個人的にはこれを「せっかちなループ生成の方法」と呼んでいる））
各プライマリコンフィギュレーションにおけるループセットの出力結果は、以下のようになった。
cfg #0 [00000000] (0 loops) cfg #1 [00000001] (1 loops) loop1: 3 edge points: 0-1 4-0 0-2 cfg #3 [00000011] (1 loops) loop1: 4 edge points: 0-2 3-1 1-5 4-0 cfg #6 [00000110] (2 loops) loop1: 3 edge points: 0-1 3-1 1-5 loop2: 3 edge points: 0-2 2-6 2-3 cfg #7 [00000111] (1 loops) loop1: 5 edge points: 4-0 2-6 2-3 3-1 1-5 cfg #15 [00001111] (1 loops) loop1: 4 edge points: 4-0 2-6 7-3 1-5 cfg #22 [00010110] (3 loops) loop1: 3 edge points: 0-1 3-1 1-5 loop2: 3 edge points: 0-2 2-6 2-3 loop3: 3 edge points: 4-0 5-4 6-4 cfg #23 [00010111] (1 loops) loop1: 6 edge points: 3-1 1-5 5-4 6-4 2-6 2-3 cfg #24 [00011000] (2 loops) loop1: 3 edge points: 4-0 5-4 6-4 loop2: 3 edge points: 3-1 2-3 7-3 cfg #25 [00011001] (2 loops) loop1: 4 edge points: 0-1 5-4 6-4 0-2 loop2: 3 edge points: 3-1 2-3 7-3 cfg #27 [00011011] (1 loops) loop1: 6 edge points: 0-2 2-3 7-3 1-5 5-4 6-4 cfg #30 [00011110] (2 loops) loop1: 5 edge points: 0-1 0-2 2-6 7-3 1-5 loop2: 3 edge points: 4-0 5-4 6-4 cfg #60 [00111100] (2 loops) loop1: 4 edge points: 0-2 2-6 7-3 3-1 loop2: 4 edge points: 4-0 1-5 5-7 6-4 cfg #105 [01101001] (4 loops) loop1: 3 edge points: 0-1 4-0 0-2 loop2: 3 edge points: 3-1 2-3 7-3 loop3: 3 edge points: 1-5 5-7 5-4 loop4: 3 edge points: 2-6 6-4 6-7
ここでは、エッジは両側の頂点IDによって表記している。各ループの向きは、有効頂点側が表側となるべく、即ち表側から見たループ頂点（エッジ上の点）の定義順は全て反時計周りになっている。
ループを２つ以上持つグループが半分ほどあることがわかることと思うが、このような完全自動のアプローチをとることにより、それらあいまいなケースについてのループ決定をいかに解決したかということの説明もしやすいと思うのだが、いかがだろうか。
そして以下が上記ループ決定を行うための完全なC++ソース（vc++2008sp1互換）であった。
#include <iostream> #include <set> #include <vector> #include < algorithm > using namespace std; #include <array> #include <map> #include <unordered_map> using namespace std::tr1; // x: 0->001 // y: 0->010 // z: 0->100 class tr { public: enum Axis{ X=0, Y, Z }; typedef array< int, 9 > array9; public: static void rotValX( int& val, array9& idh ){ trval_sub( val, &tr::rotIdX, idh ); } static void rotValY( int& val, array9& idh ){ trval_sub( val, &tr::rotIdY, idh ); } static void rotValZ( int& val, array9& idh ){ trval_sub( val, &tr::rotIdZ, idh ); } static void symValX( int& val, array9& idh ){ trval_sub( val, &tr::symIdX, idh, !0 ); } static void symValY( int& val, array9& idh ){ trval_sub( val, &tr::symIdY, idh, !0 ); } static void symValZ( int& val, array9& idh ){ trval_sub( val, &tr::symIdZ, idh, !0 ); } static void complement( int& val, array9& idh ){ val = (~val)&0xff; idh[8]=(idh[8]==0)?1:0; } static void rotIdX( int& id ){ rotid_sub( id, Y, Z ); } static void rotIdY( int& id ){ rotid_sub( id, X, Z ); } static void rotIdZ( int& id ){ rotid_sub( id, X, Y ); } static void symIdX( int& id ){ sym_planar( id, X ); } static void symIdY( int& id ){ sym_planar( id, Y ); } static void symIdZ( int& id ){ sym_planar( id, Z ); } private: static void bit_mod( int& val, int bitid, bool s ){ val &= ~(1<<bitid); val |= s<<bitid; } static void rotid_sub( int& id, Axis a, Axis b ) { bool u = !(id&(1<<b)); bool v = id&(1<<a); bit_mod( id, a, u ); bit_mod( id, b, v ); return; } static void sym_planar( int& id, Axis a) { bool v = !( id&(1<<a) ); bit_mod( id, a, v ); return; } template< class TrIdFunc > static void trval_sub( int& val, TrIdFunc fn, array9& idh, bool is_flip=0 ){ int ids[8]; for( int i=0; i<8; ++i ){ int id=i; fn( id ); ids[i] = id; } subst( val, ids ); // for( int i=0; i<8; ++i ){ int id=idh[i]; fn( id ); idh[i] = id; } if( is_flip ) idh[8]=(idh[8]==0)?1:0; return; } static void subst( int& val, int ids[8] ){ int tmp = val; val = 0; for( int i=0; i<8; ++i ) val |= ( (tmp&(1<<i))!=0 ) << ids[i]; return; } }; // 14あるプライマリコンフィギュレーションの // テセレーションについて考えよう class Cube { public: class Edge; class HalfEdge; typedef map< pair< int, int >, Cube::Edge* > EdgeMap; typedef array< int, 4 > array4; Cube( void ){} ~Cube( void ){ for each( EdgeMap::value_type ipair in m_edge_map ){ Edge* ie = ipair.second; delete ie->m_he->m_opp; delete ie->m_he; }// ie } class Edge { public: Edge( void ){ m_flag=0; m_used = 0; } bool Flag( void ){ return m_flag; } bool Used( void ){ return m_used; } HalfEdge* m_he; bool m_flag; bool m_used; }; class HalfEdge { public: HalfEdge( int v ){ m_opp=0; m_vertex_id = v; } HalfEdge* m_next; HalfEdge* m_opp; Edge* m_edge; int m_vertex_id; }; public: void AddFace( array4& ids ){ unsigned int arr[5]; copy( ids.begin(), ids.end(), arr ); arr[4] = ids.front(); HalfEdge* he_arr[4]; for( int i=0; i<4; ++i ){ int s=arr[i]; int e=arr[i+1]; HalfEdge*& he = he_arr[i]; he = new HalfEdge( s ); if( s>e ){ swap( s, e ); } EdgeMap::iterator it = m_edge_map.find( make_pair( s, e ) ); if( it==m_edge_map.end() ){ Edge* new_edge = new Edge(); m_edge_map.insert( make_pair( make_pair( s, e ), new_edge ) ); new_edge->m_he = he; he->m_edge = new_edge; }else{ it->second->m_he->m_opp = he; he->m_opp = it->second->m_he; he->m_edge = it->second; } }// i for( int i=0; i<3; ++i ){ he_arr[i]->m_next = he_arr[i+1]; }// i he_arr[3]->m_next = he_arr[0]; } public: EdgeMap m_edge_map; }; template< class Func > void rot_h( Cube::array4& ids, Func fn ) { for( int i=0; i<4; ++i ){ fn( ids[i] ); }// id return; } int main( void ) { // サイコロメッシュ（クアッドメッシュ、６面体）を構築します Cube* cube = new Cube(); { Cube::array4 ids={ 0, 1, 5, 4 }; for( int i=0; i<4; ++i ){ cube->AddFace( ids ); rot_h( ids, &tr::rotIdX ); }// i for( int i=0; i<2; ++i ){ rot_h( ids, &tr::rotIdZ ); cube->AddFace( ids ); rot_h( ids, &tr::rotIdZ ); }// i } // 各プライマリコンフィギュレーションについてテセレーションを // 作成していきます。 array< int, 14 > groups={ 0, 1, 3, 6, 7, 15, 22, 23, 24, 25, 27, 30, 60, 105 }; for each( const int& cfg in groups ){ // 有効値を持つ頂点と持たない頂点をそれぞれ一つずつ両端に持つ // エッジのフラグを立てます。 int cnt=0; for each( Cube::EdgeMap::value_type ipair in cube->m_edge_map ){ Cube::Edge* ie = ipair.second; bool s_flag = (1<<ie->m_he->m_vertex_id) & cfg; bool e_flag = (1<<ie->m_he->m_next->m_vertex_id) & cfg; if( s_flag ^ e_flag ){ ie->m_flag = !0; ++cnt; }else{ ie->m_flag = 0; } }// ipair // フラグのついたエッジについて、順次、隣接する、やはりフラグのついたエッジを辿ることにより、 // テセレーションすべきループを作成していきます。 vector< Cube::Edge* > ev; for each( Cube::EdgeMap::value_type ipair in cube->m_edge_map ){ Cube::Edge* ie = ipair.second; if( ie->m_flag ){ ev.push_back( ie ); ie->m_used = 0; } }// ipair // ループを作成 vector< vector< Cube::Edge* > > loop_vec; while( !ev.empty() ){ vector< Cube::Edge* > loop; Cube::HalfEdge* cur=ev.front()->m_he; { // 始頂点側に有効ビット頂点が来るようにします。 if( cfg & (1<<cur->m_vertex_id) ){ }else{ cur = cur->m_opp; } } cur->m_edge->m_used=!0;// 使用済みにする loop.push_back( cur->m_edge ); while( !0 ){ // 現在の loop を閉じることを何よりも優先する if( cur->m_next->m_next->m_next->m_edge == loop.front() ){ break; }else if( cur->m_next->m_next->m_edge == loop.front() ){ break; }else if( cur->m_next->m_edge == loop.front() ){ break; } // used がなければ if( cur->m_next->m_next->m_next->m_edge->m_flag && !cur->m_next->m_next->m_next->m_edge->m_used ){ cur = cur->m_next->m_next->m_next->m_opp; }else if( cur->m_next->m_next->m_edge->m_flag && !cur->m_next->m_next->m_edge->m_used ){ cur = cur->m_next->m_next->m_opp; }else if( cur->m_next->m_edge->m_flag && !cur->m_next->m_edge->m_used ){ cur = cur->m_next->m_opp; }else{ break; } cur->m_edge->m_used=!0;// 使用済みにする loop.push_back( cur->m_edge ); }// while; loop_vec.push_back( loop ); ev.erase( remove_if( ev.begin(), ev.end(), mem_fun( &Cube::Edge::Used) ), ev.end() );// 使用済みは削除 }// while cout << "cfg #"<<cfg<<" ["; { for( int i=0; i<8; ++i ) cout<< ( 0!= ((1<<(7-i))&cfg) ); } cout<<"] ("; cout << loop_vec.size()<<" loops)"<<endl; int loop_cnt=0; for each( const vector< Cube::Edge* >& loop in loop_vec ){ cout << " loop"<<loop_cnt+1<<": "<< loop.size() << " edge points: "; for each( Cube::Edge* const & ie in loop ){ cout<< ie->m_he->m_vertex_id << "-"<<ie->m_he->m_next->m_vertex_id<<" "; }// ie ++loop_cnt; cout<< endl;; }// loop cout<< endl;; }// cfg return 0; }
trクラスは第一回より利用してきているものの改良でID値の回転などのために用いている。Cubeクラスは、今回考えたハーフエッジクアッドメッシュをサポートするクラスであり、キューブモデルの構築は、main関数の冒頭にて行っている。main関数の残りはループ形状抽出コアと、結果の出力処理だ。
この場合のハーフエッジデータ構造は、実はエッジとハーフエッジしか必要ないので、Cubeクラスの実装もそのようにしている。Faceの概念は必要ないし、Vertexの概念も今回は0から7の頂点ID、即ち整数値のみで十分だ。
