マーチングキューブズの旅(3):各エントリの相対配置について

暇つぶしとしてのマーチングキューブズ自作研究も3回目を迎えたよ。
今回は、各グループにおける、各コンフィギュレーションエントリと、そのプライマリコンフィギュレーションとの相対配置関係を情報化することを考えよう。
マーチングキューブズにおいて予め定義せよとされる14通りの基本テセレーション構成は、即ち、前回までに抽出した、14ある各グループのプライマリコンフィギュレーションについてのみ用意すればよいわけなのだが、当然ながらこれは、逆に言えば、プライマリでない残りの242(=256-14)コンフィギュレーションについては、そのプライマリコンフィギュレーションとの位置関係情報を把握&保持していなくてはならない(!)ということを意味する。

実はここが、マーチングキューブ法を実装するにあたって、実装者が途方にくれてしまいがちな最大の難関なのである。
必要なのは、非プライマリなコンフィギュレーションの各々について、己が属するグループのプライマリコンフィギュレーションとぴったり重ね合わせたときに、8つの頂点のそれぞれが、プライマリコンフィギュレーションにおける、どの頂点に相当するのかという情報だ。より具体的には、

プライマリコンフィギュレーション頂点ID: 7 6 5 4 3 2 1 0 本来の頂点位置
非プライマリコンフィギュレーション頂点ID: a b c d e f g h ←重ね合わせるために操作した後の頂点位置(このa-hを解決し記録しておく)

つまり、上記の、a-h の、即ち0から7の値をとる8つの整数値を、記録しておく必要がある。0から7の値は3bitあればよいので、即ち、1コンフィギュレーションあたり、3×8 = 24byte の情報量があれば十分だ。
さて「十分」と言った。実は、相対配置情報としては1コンフィギュレーションあたり、24バイトも必要ない。よくよく考えてみると、a-hの数値というものは、ある程度、決まりきったものになる。回転・対称・ビット反転という操作の組み合わせは、果たして何通り可能かということを考えてみてほしい。サイコロで考えてみよう。

●回転操作
机の上にサイコロを違った向きで置く置きかたは何通りあるだろうか。3の目が上になる置き方は全部で4通り。なので1の目から6の目までそれぞれを上にして、置く置きかたは24通り。
別の考え方もある。XYZの3次元空間において、ある角を原点にして、サイコロを第一象限に置くとき、その置き方は3通りある。サイコロの角は全部で8つあるので、置き方は全部で3×8=24通り。(24通りの置き方をプログラムで体系的に辿るのは後者の考え方の方が向いているので 実際の実装では後者の考え方を採用した。)
●対称操作
XY面、YZ面、XZ面、3軸分の対象性、つまり鏡が縦横奥方向に3つあるわけなのだが、鏡は2枚置くともとの世界と同じになってしまうので、2枚分は打ち消して、1枚しかないことになり、サイコロと鏡の中のサイコロはどういう置き方をしても重ならないので、回転操作と対称操作とを合わせたサイコロの置き方はこれで48通り。
●ビット反転
サイコロ自体の色が変わるとでも考えて、回転・対称・ビット反転操作とあわせたサイコロの置き方は全部でで96通りとなる。

以上の考察より、プライマリコンフィギュレーションからの可能な操作はたかだか全部で96通りしかないことになる。1コンフィギュレーションあたり、そのプライマリコンフィギュレーションとの位置関係を語るに必要な情報量はたかだか7ビットですむのである。アルゴリズムの一部として最終的に表のようなものを用意するにしても、7ビットですむのなら7ビットで済ませたいところだ。最終的に使いたいa-hの24ビット情報への展開など、コンピュータにやらせればいいのである。
この96通りの操作をどう体系づけてコンピュータに辿らせるか。そのためのアルゴリズムを考えてみよう。以下のようなもので十分だろう。。

  1. 全ての可能な回転操作を行う。
  2. 鏡像(対称)操作(X軸反転1回のみ)を行う。
  3. 1を行う。
  4. ビット反転操作をする
  5. 1、2を行う。

上記アルゴリズムで辿れる各配置に0から96までの通し番号をつけることもできるだろう。
それでは、上述の242ケースの非プライマリコンフィギュレーションの各々について、プライマリからの相対配置情報を調べ上げる方法について考えよう。
「グループ」の定義から、任意のプライマリコンフィギュレーションについて、上記1)〜5)の操作を行うことで、グループの全てのコンフィギュレーションを辿ることができるはずである。このとき、あるコンフィギュレーションにたどり着いた時点での配置情報(上述のa-h)を、そのコンフィギュレーション番号とセットで記録しておく。
(下図は、実際の実装のイメージである。例えばZ軸まわり回転なら、プライマリにおける頂点ID値と頂点の2値パターンを同時に回転している。その結果コンフィギュレーション#11となった。右側の図で、7番から0番の各頂点ID値が本来の頂点位置でいうところの、どこにおるかという情報を、a-h情報として記録しておく。下の実装では実際にはh->aの向きで出力させている)

言うまでもないが、前回までの考察により、プライマリコンフィギュレーションが全て抽出されているからこそ、このようなアプローチが可能となるのである。
以上の処理を行うプログラムをVisualStudio2008SP1を使い、c++0xで記述してみたのが、以下のコードだ。

#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 );
		idh[8]=(idh[8]==0)?1:0;
	}
	static void symValY( int& val, array9& idh ){ 
		trval_sub( val, &tr::symIdY, idh ); 
		idh[8]=(idh[8]==0)?1:0;
	}
	static void symValZ( int& val, array9& idh ){ 
		trval_sub( val, &tr::symIdZ, idh ); 
		idh[8]=(idh[8]==0)?1:0;
	}
	static void complement( int& val, array9& idh ){ 
		val = (~val)&0xff; 
		idh[8]=(idh[8]==0)?1:0; 
	}
private:
	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 ); }
	
	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 ){
		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;
		}
		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;
	}
};

//typedef unordered_map< int, tr::array9 > ConfigurationMap;
typedef map< int, tr::array9 > ConfigurationMap;


void turn3way( int& cfg, tr::array9& ids, ConfigurationMap* mp  )
{
	for( int i=0; i<3; ++i ){	
		ConfigurationMap::iterator it = mp->find( cfg );
		if( it==mp->end() ){
			mp->operator[]( cfg ) = ids;
		}
		tr::rotValX( cfg, ids );
		tr::rotValZ( cfg, ids );
	}
	return;
}

void try_4( int& cfg, tr::array9& ids, ConfigurationMap* mp  )
{
	for( int i=0; i<4; ++i ){
		turn3way( cfg, ids, mp );
		tr::rotValY( cfg, ids );
	}// i
	return;
}

void try_front_and_back( int& cfg, tr::array9& ids, ConfigurationMap* mp  )
{
	try_4( cfg, ids, mp );
	tr::rotValX( cfg, ids );
	tr::rotValX( cfg, ids );
	try_4( cfg, ids, mp );
	tr::rotValX( cfg, ids );
	tr::rotValX( cfg, ids );
	return;
}

void try_rot_and_sym( int& cfg, tr::array9 ids, ConfigurationMap* mp  )
{
	try_front_and_back( cfg, ids, mp );
	tr::symValX( cfg, ids );
	try_front_and_back( cfg, ids, mp );
	tr::symValX( cfg, ids );
	return;
}


void try_whole_op( int cfg, tr::array9 ids, ConfigurationMap* mp  )
{
	try_rot_and_sym( cfg, ids, mp );
	tr::complement( cfg, ids );
	try_rot_and_sym( cfg, ids, mp );	
	return;
}


int main( int argc, char* argv[] )
{
	// 各プライマリコンフィギュレーションへの操作を通して、
	// 同じグループのコンフィギュレーションを全て導き出してみよう。
	cout <<endl<<endl<<endl;

	int val=30;
	array< int, 14 > groups={ 0, 1, 3, 6, 7, 15, 22, 23, 24, 25, 27, 30, 60, 105 };
	vector< ConfigurationMap > mpv;
	int n_config=0;

	for each( const int& cfg in groups ){
		ConfigurationMap mp;
		tr::array9 ids={0,1,2,3,4,5,6,7,0};// 最後のはビット反転の有無
		try_whole_op( cfg, ids, &mp );
		cout<< cfg <<": "<<endl;
		int cnt=0;
		for each( const ConfigurationMap::value_type& pr in mp ){
			cout<< "  "<<pr.first <<"  (";
			for( int i=0; i<8; ++i ) 
				cout<< pr.second[i];
			cout<< ":"<<pr.second[8];
			cout<<")";
			++cnt;
			if( cnt%3==0 ) cout<<endl;
		}
		cout<< endl;
		n_config += mp.size();
	}// i
	cout<< endl<< n_config <<endl;
	return 0;
}

そして出力は以下のようになった。

group #0:
  0  (01234567:0)  255  (01234567:1)
group #1:
  1  (01234567:0)  2  (15370426:0)  4  (23670145:0)
  8  (31752064:0)  16  (46570213:0)  32  (54761032:0)
  64  (62734051:0)  127  (73516240:1)  128  (73516240:0)
  191  (62734051:1)  223  (54761032:1)  239  (46570213:1)
  247  (31752064:1)  251  (23670145:1)  253  (15370426:1)
  254  (01234567:1)
group #3:
  3  (01234567:0)  5  (02461357:0)  10  (31752064:0)
  12  (23670145:0)  17  (04152637:0)  34  (15370426:0)
  48  (54761032:0)  63  (67452301:1)  68  (62734051:0)
  80  (46570213:0)  95  (57134602:1)  119  (37261504:1)
  136  (37261504:0)  160  (57134602:0)  175  (46570213:1)
  187  (62734051:1)  192  (67452301:0)  207  (54761032:1)
  221  (15370426:1)  238  (04152637:1)  243  (23670145:1)
  245  (31752064:1)  250  (02461357:1)  252  (01234567:1)

group #6:
  6  (01234567:0)  9  (20316475:0)  18  (04152637:0)
  20  (02461357:0)  33  (10543276:0)  40  (15370426:0)
  65  (40625173:0)  72  (23670145:0)  96  (46570213:0)
  111  (54761032:1)  123  (62734051:1)  125  (31752064:1)
  130  (31752064:0)  132  (62734051:0)  144  (54761032:0)
  159  (46570213:1)  183  (23670145:1)  190  (40625173:1)
  215  (15370426:1)  222  (10543276:1)  235  (02461357:1)
  237  (04152637:1)  246  (20316475:1)  249  (01234567:1)

group #7:
  7  (01234567:0)  11  (13025746:0)  13  (20316475:0)
  14  (32107654:0)  19  (04152637:0)  21  (02461357:0)
  31  (75643120:1)  35  (10543276:0)  42  (15370426:0)
  47  (67452301:1)  49  (45016723:0)  50  (51407362:0)
  55  (76325410:1)  59  (62734051:1)  69  (26043715:0)
  76  (23670145:0)  79  (54761032:1)  81  (40625173:0)
  84  (64207531:0)  87  (73516240:1)  93  (57134602:1)
  112  (46570213:0)  115  (37261504:1)  117  (31752064:1)
  138  (31752064:0)  140  (37261504:0)  143  (46570213:1)
  162  (57134602:0)  168  (73516240:0)  171  (64207531:1)
  174  (40625173:1)  176  (54761032:0)  179  (23670145:1)
  186  (26043715:1)  196  (62734051:0)  200  (76325410:0)
  205  (51407362:1)  206  (45016723:1)  208  (67452301:0)
  213  (15370426:1)  220  (10543276:1)  224  (75643120:0)
  234  (02461357:1)  236  (04152637:1)  241  (32107654:1)
  242  (20316475:1)  244  (13025746:1)  248  (01234567:1)

group #15:
  15  (01234567:0)  51  (04152637:0)  85  (02461357:0)
  170  (15370426:0)  204  (23670145:0)  240  (46570213:0)

group #22:
  22  (01234567:0)  41  (15370426:0)  73  (23670145:0)
  97  (46570213:0)  104  (73516240:0)  107  (62734051:1)
  109  (54761032:1)  121  (31752064:1)  134  (31752064:0)
  146  (54761032:0)  148  (62734051:0)  151  (73516240:1)
  158  (46570213:1)  182  (23670145:1)  214  (15370426:1)
  233  (01234567:1)
group #23:
  23  (01234567:0)  43  (15370426:0)  77  (23670145:0)
  113  (46570213:0)  142  (31752064:0)  178  (54761032:0)
  212  (62734051:0)  232  (73516240:0)
group #24:
  24  (01234567:0)  36  (04152637:0)  66  (02461357:0)
  126  (15370426:1)  129  (15370426:0)  189  (02461357:1)
  219  (04152637:1)  231  (01234567:1)
group #25:
  25  (01234567:0)  26  (10543276:0)  28  (26043715:0)
  37  (04152637:0)  38  (13025746:0)  44  (31752064:0)
  52  (40625173:0)  56  (57134602:0)  61  (73516240:1)
  62  (64207531:1)  67  (02461357:0)  70  (20316475:0)
  74  (37261504:0)  82  (45016723:0)  88  (62734051:0)
  91  (76325410:1)  94  (51407362:1)  98  (54761032:0)
  100  (67452301:0)  103  (75643120:1)  110  (46570213:1)
  118  (32107654:1)  122  (23670145:1)  124  (15370426:1)
  131  (15370426:0)  133  (23670145:0)  137  (32107654:0)
  145  (46570213:0)  152  (75643120:0)  155  (67452301:1)
  157  (54761032:1)  161  (51407362:0)  164  (76325410:0)
  167  (62734051:1)  173  (45016723:1)  181  (37261504:1)
  185  (20316475:1)  188  (02461357:1)  193  (64207531:0)
  194  (73516240:0)  199  (57134602:1)  203  (40625173:1)
  211  (31752064:1)  217  (13025746:1)  218  (04152637:1)
  227  (26043715:1)  229  (10543276:1)  230  (01234567:1)

group #27:
  27  (01234567:0)  29  (20643175:1)  39  (10325476:1)
  46  (31752064:0)  53  (04152637:0)  58  (51734062:1)
  71  (02461357:0)  78  (32761054:1)  83  (40516273:1)
  92  (62734051:0)  114  (54761032:0)  116  (64752031:1)
  139  (13570246:1)  141  (23670145:0)  163  (15370426:0)
  172  (73625140:1)  177  (45670123:1)  184  (57134602:0)
  197  (26370415:1)  202  (37261504:0)  209  (46570213:0)
  216  (76543210:1)  226  (75316420:1)  228  (67452301:0)

group #30:
  30  (01234567:0)  45  (13025746:0)  54  (04152637:0)
  57  (10543276:0)  75  (20316475:0)  86  (02461357:0)
  89  (26043715:0)  99  (45016723:0)  101  (40625173:0)
  106  (73516240:0)  108  (76325410:0)  120  (75643120:0)
  135  (32107654:0)  147  (51407362:0)  149  (64207531:0)
  154  (57134602:0)  156  (62734051:0)  166  (31752064:0)
  169  (15370426:0)  180  (67452301:0)  198  (37261504:0)
  201  (23670145:0)  210  (54761032:0)  225  (46570213:0)

group #60:
  60  (01234567:0)  90  (02461357:0)  102  (04152637:0)
  153  (15370426:0)  165  (46570213:0)  195  (23670145:0)

group #105:
  105  (01234567:0)  150  (15370426:0)

14のプライマリコンフィギュレーションに対して1)〜5)の操作を行うことで、各グループのコンフィギュレーションが全て現れ、結果的に全256コンフィギュレーションとその配置情報とが、関連付けて出力された。
プライマリのみについてテセレーションを持つ以上、鏡像(対称)操作もしくはビット反転操作のほどこされた相対配置を持つコンフィギュレーションについては、後にテセレーションの三角の向き(CW、CCW)を逆にしなくてはならないので、相対配置情報の一番後ろに、:1がついているはずだ。
あとは14プライマリコンフィギュレーションについてのテセレーションについて考えるだけだ。それはまた次回に。