koga study

VOXEL分割

改訂情報

　ソースプログラムのGetPlace()関数で，近傍VOXELの探索の高速化に伴い改訂． May, 21, 2002

最初に

　アトラクタの解析は，一般に膨大な処理時間を必要とするものである．私は，そう高速に処理できるマシンをもっているわけでもないし，日頃はVisual C++でぼちぼちとプログラムを走らせているのが実状である．昔は全点探索など気にもとめていなかったし，そういうものだと思っていたのだが，ふと，高専時代にレイ・トレーシングの研究で使っていたVOXELという概念をつかうと効率がよいことがわかってきた．しかし，VOXELを高次元空間に拡張すると膨大なメモリが必要となるという欠点から，「5次元くらいの解析までしかつかえないかなぁ」と思ってプログラムを組んでいたわけだが，最近組みなおしてみたら，理論的には何次元でもOKなプログラムが完成した．動作テストは15次元くらいまでしかやってないが多分十分だろう（エラーがでるとすればlistsuu（後述）のオーバーフローだろうか？）．この概念をつかうことで， FNNやリアプノフ解析・時系列予測などの近傍点探索プログラムを格段と高速化できる．

コンピュータグラフィクスの分野のVOXEL分割

　ディスプレイに表示された図形は一般に格子点上に配置した小正方形(画素)の集合で表される．この方法と同じように，立体を3次元の格子点上の小立方体の集合で表すのがボクセル表現である． 2次元の小正方形を画素(ピクセル:pixel)とよぶのに対し，3次元の小立方体をボクセル(voxel)という．境界表現やCSGと比較すると，データ量が膨大で精度が劣り，立体の表示・移動・回転などは手間がかかるという問題がある．しかし，データ構造が簡単で，集合演算がきわめて簡単な原理で行える．人工的で規則的なものでなく，自然界にある非常に不規則な形状を表す場合に向いている．

VOXELの拡張とアトラクタへの応用

　一般的な m 次元空間に再構成されたアトラクタ上の点 X _i に最も近い近傍点を探索する処理を考える．この処理は，未知の構造をもつアトラクタに対しては全ての軌道点に対して距離を計算し（低次元までの距離計算で打ち切る場合を含む），最小となる一点を抽出しなくてはならない．
　しかしながら，座標軸を均等に分割し（実際は均等である必然性は無いが），アトラクタに埋め込まれた座標点のグループ分けをおこなうことで，点 X _i の最近傍点は同じグループの点であることが確約されており，他のグループの点について処理をおこなう必要は無くなるだろう．
　私は，コンピュータグラフィクスにおけるモデリング，レイ・トレーシングなどの処理で利用される 3次元空間のVOXEL分割を m 次元の超空間へと拡張し，そこに埋め込まれたアトラクタの VOXELグループ分けを行い，リアプノフ解析などのアトラクタ上の解析に積極的に利用している．この手法を K-VOXELと命名しよう．（KはKOGAのK？）
　簡単化のため m 次元アトラクタ上に仮定されるVOXELは超立方体とする．（もちろん立方体である必然性はない）各超立方体は空間に隙間無く整然と配置され，上下左右前後 1 VOXELの位置には他のVOXELが存在する．図にすると３次元空間では下図のようなルービックキューブのようなものになる．

　全てのVOXELは座標データを表す位置情報を有しており，隣り合うVOXELどうしの距離を 1 VOXELと定義する．さて，先ほど「点 X _i の最近傍点は同じグループの点であることが確約されて・・・」と述べたが，この真意は偽である．点 X _i がVOXELの端点や辺に近い位置に配置されている場合は，最近傍点は 1 VOXEL離れた隣のVOXEL内に存在する可能性があるし，アトラクタの密度が疎であれば， 2 VOXEL離れたところに存在するかもしれない．即ち，点 X _i が入るVOXELのみでなく，そこからある程度離れたVOXELにたいしても探索をおこなう必要性が生じてくる．何VOXELまで探索するかは経験に基づくところであるが，VOXELは p VOXEL以内の距離にどのような VOXELが存在するか知っておく必要がある．

　さて，ここまでの話は単純だが，実際に m = 10 次元のアトラクタで，各軸の分割数を q = 16 とすると，必要なVOXEL数はいくつになるだろうか？
　正解は，16 ¹⁰ ≒ 10 ¹² である．
　これでは，VOXELを構成するためのメモリが 4.4 [TB] ほど必要になってしまう．しかし，よくよく考えてみるとアトラクタに埋め込む点の数 N はせいぜい 100 万点が限度だろうから，全ての点がそれぞれ 1 つのVOXEL に入ってもVOXELの数は 10 ⁶ であり，のこりは空のVOXELだということになる．経験的に，カオス力学系の場合，アトラクタを構成するVOXELの最大数は N に依存しない．（想定した空間の外に飛び出すことがないため）．ただし，一様乱数などを埋め込んだ場合は軌道点が空間内に一様に分布するため， 10 ¹² 個のVOXELをフルに利用することになる．

K-VOXELのコーディング

　コーディング作業は以下の流れでおこなう．説明がわかりにくいと思われますので，そのときはプログラムリストを参照してください．

　１．アトラクタの N 個の点がどのVOXELに入るかしらべ，それぞれのVOXELの番地を求めておく．

　２．全ての番地をクイックソートし，必要なVOXELの個数を求める．

　３．以下のデータ構造のK-VOXELを作成する

	{
		int* 番地情報; ※
		int  VOXELに入る１つの点への指標; †
		int* このVOXELから p VOXELの距離にあるVOXELへの指標;
	}

　※　番地情報は m 個の座標データとするのが最も単純である．しかし，後々の計算量を省くため，私は 5次元ごとに一まとめにしている．
　†　VOXELに入る１つの点への指標（インデックス）は，このVOXELがアクセスされたときに，返信するアドレス値のことであり，とび先のアドレスには次のアドレスが格納されている(後述 ix[] )．

　４．アトラクタの番地にリンクする飛び先アドレスとVOXEL番号の指標を作る

	{
		int ix[N];  //飛び先，初期値はすべてINT_MAXにする
		int idx[N]; //自分の点が入るVOXELのアドレス
	}

　これは， X _i が選択された場合，voxel[i]の番地のVOXEL情報が読み込まれ， VOXELから呼び出された場合は，距離計算などの処理のあと，ix[i]のさすアドレスの点に飛ぶことを表す．

　５．全てのVOXELどうしの距離を計算し， p VOXEL以内の距離にあるVOXELは指標により連結する．

プログラムリスト

C++言語によるCKVoxelというクラスを作っています．ファイルは，KVoxel.h，KVoxel.cppに分かれます．
上記の文字に対応していない場合があります．コメントを参照してください．

KVoxel.hの内部構造


// KVoxel.h: interface for the CKVoxel class.
//
//////////////////////////////////////////////////////////////////////

#if !defined(AFX_VOXEL_H__968ECF45_7769_11D4_BA93_CE05854FCE21__INCLUDED_)
#define AFX_VOXEL_H__968ECF45_7769_11D4_BA93_CE05854FCE21__INCLUDED_

#if _MSC_VER > 1000
#pragma once
#endif // _MSC_VER > 1000
class CKVoxel  
{
public:
	int NMAX;				// アトラクタの軌道点の総数
	int ARM;				// 近傍ボクセルとみなす距離の最大値
	int DIVISION;				// 各軸の分割数
	_int64 Maxsuu;				// VOXELの数

	//配列情報
	// アトラクタ軌道点が保持する情報
	int* idx;				// アトラクタの軌道点が入るVOXELの番号
	int* ix;				// 次に呼び出す軌道点へのアドレス
	// VOXELが保持する情報
	int** placeID;				// VOXELの近傍ボクセルのVOXEL番号
	int* voxel;				// VOXELに入る軌道点の先頭番地

	// 処理によってアクセスしてもらう軌道点のリスト
	int* ret_list;

	int past_get_place;			// 過去にGetPlace関数で呼び出されたVOXEL番号

	//関数
	void SavePlace1();			// 1次元の時系列データを遅れ時間で埋め込む
	void SavePlace2();			// m次元の配列データをそのまま埋め込む
	void GetPlace(int t, int arm = 1, bool all = false);		// 外部処理からVOXELを呼び出し適切な数の軌道点の番地をゲットする

	void CountListMax(int & j, int & k, int & length);	// 近傍VOXELの最大数をカウントする
	
	void CKVoxelMain();			// ボクセルの登録処理のメイン

	CKVoxel(int m_dim, int nmax,            double** x, int div = 16);
	CKVoxel(int m_dim, int nmax, int m_tau, double*  x, int div = 16);
	virtual ~CKVoxel();
private:
	double* x1;				// 1次元の時系列データ
	double** x2;				// m次元の配列データ
	int tau;				// 遅れ時間
	int dim;				// 埋め込み次元
};

#endif // !defined(AFX_VOXEL_H__968ECF45_7769_11D4_BA93_CE05854FCE21__INCLUDED_)

KVoxel.cppの内部構造


// KVoxel.cpp: implementation of the CKVoxel class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "KVoxel.h"
#include "Sort.h"
#include <math.h>

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

// 多次元配列データを埋め込む場合のクラスの定義
//   m_dim : 埋め込み次次元
//   nmax  : アトラクタのデータ点の総数
//   **x   : データポインタ
//   di    : 各座標軸の分割数
CKVoxel::CKVoxel(int m_dim, int nmax, double **x, int div)
{
	DIVISION = div;							//分割数
	x1 = NULL, x2 = x;							//アトラクタのポインタ
	NMAX = nmax;								//点数
	dim = m_dim;								//埋め込み次元

	//アトラクタの配置
	placeID = new int*[nmax];						//座標データを格納するための配列の定義
	for(int i = 0;i < nmax; i++){
		placeID[i]=new int[m_dim / 5 + 1];			//５次元を一まとめにする
	}
	SavePlace2();							//埋め込み処理の呼び出し

	//VOXELの定義
	CKVoxelMain();							//VOXELのメイン処理へ移行
}

// 1次元時系列データを埋め込む場合のクラスの定義
//   m_dim : 埋め込み次次元
//   nmax  : アトラクタのデータ点の総数
//   m_tau : 埋め込みの遅れ時間
//   *x    : データポインタ(★注意★　必ずx[0-(m_dim-1)*m_tau]にデータが存在すること)
//   di    : 各座標軸の分割数
CKVoxel::CKVoxel(int m_dim, int nmax, int m_tau, double *x, int div)
{
	DIVISION = div;							//分割数
	x2 = NULL, x1 = x;							//アトラクタのポインタ
	tau = m_tau;								//遅れ時間
	NMAX = nmax;								//点数
	dim = m_dim;								//埋め込み次元

	//アトラクタの配置
	placeID = new int*[nmax];						//座標データを格納するための配列の定義
	for(int i = 0; i < nmax; i++){
		placeID[i]=new int[m_dim / 5 + 1]; //５次元を一まとめにする
	}
	SavePlace1();							//埋め込み処理の呼び出し

	//VOXELの定義
	CKVoxelMain();							//VOXELのメイン処理へ移行
}

// VOXELの登録処理のメイン
void CKVoxel::CKVoxelMain()
{
	int i, j, k, t, p, r;

	ARM = 2;									// 近傍VOXELとの距離の定義

	//番地データ(placeID)をクイックソートする
	int* index=new int[NMAX];							// ソート前の番地がわかるように指標データの作成
	for(i = 0; i < NMAX; index[i] = i, i++);
	//多次元ソート
	CSort sort;
	sort.QuickSortRiseIP(placeID, index, 0, NMAX-1, 0);			//indexとplaceIDのポインタを同時にソートする（後述）
	//各層をソートする（5次元ごとにまとめてある）
	for(r = 1; r <= dim / 5; r++){
		//局所アドレスを移動しながらソートする
		for(i = 0; i < NMAX-1;   ){
			for(j = i + 1; j < NMAX; j++){
				for(t = r-1; t >= 0 && placeID[i][t] == placeID[j][t]; t--);
				if(t >= 0 && placeID[i][t] != placeID[j][t]) break;
			}
			t = j;								//ソート対象アドレスの決定
			if(i != t - 1){						//ソート処理
				sort.QuickSortRiseIP(placeID, index, i, t-1, r);
			}
			i = t;								//局所アドレスの移動
		}
	}

	//必要なVOXELの個数(Maxsuu)をカウントする
	Maxsuu = 1;
	for(i = 0; i < NMAX-1; i++){
		for(t = r-1; t >= 0 && placeID[i][t] == placeID[i + 1][t]; t--);
		if(t >= 0 && placeID[i][t] != placeID[i + 1][t]) Maxsuu++;
	}

	//VOXELから最初に呼び出される軌道点の番地を宣言
	voxel = new int[Maxsuu];
	for(i = 0; i < Maxsuu; voxel[i++]=INT_MAX);				//(初期値はINT_MAX)

	//次々に呼び出される軌道点のとび先番地を宣言
	ix = new int[NMAX];
	for(i = 0;i < NMAX; ix[i++]=INT_MAX);					//(初期値はINT_MAX)

	//ボクセルの格納
	idx = new int[NMAX];								//各軌道点の格納されるボクセルの番号の宣言
	int** ID = new int*[Maxsuu];						//ボクセルの番地を保存しておくＩＤ（後で使用する）
	j = 0;										//ボクセルの個数のカウンタ
	for(i = 0; i < NMAX;   ){

		//各軸における番地の保存
		ID[j] = new int [dim];
		k = dim - 1;
		for(r = dim / 5; r >= 0; r--){
			p = placeID[i][r];
			for(   ; k >= r * 5; k--){
				ID[j][k] = p % DIVISION, p /= DIVISION;
			}
		}

		//同じVOXELに含まれるグループの登録
		p = -1;									//１つ前のインデックスの値
		for(;;){
			if(p < 0){
				voxel[j] = index[i];					//最初の１点はVOXELに登録
			}else{
				ix[p] = index[i];					//次からは飛び先アドレスとしてixに登録
			}
			p = index[i];							//１つ前のインデックスの値の更新
			idx[ index[i] ] = j;						//軌道点の格納されるVOXELのアドレスの保存

			if(i == NMAX - 1){ i++; break; }

			//同じ座標のVOXELである限りi++しながら登録し続ける
			for(t = dim / 5; t >= 0 && placeID[i][t] == placeID[i + 1][t]; t--);
			i++;
			if(t >= 0) break;
		}
		j++;
	}

	//placeIDの一時削除(後で再び定義して利用する)
	for(i = 0; i < NMAX; i++) delete [] placeID[i];
	delete [] placeID;

	//自分の近傍のVOXELのリストを作る
	//リストの個数のカウント
	int listsuu=0;						//リストの最大数
	int length=0;						//探索中のVOXELどうしの距離
	j=0;								//現在のループ数
	CountListMax(j,listsuu,length);

	//placeIDの再定義(自分の近傍のVOXELの番地を保存する)
	placeID = new int*[Maxsuu];
	for(i = 0;i < Maxsuu; i++){
		placeID[i] = new int [listsuu + 1];				//個数がわからないので一応最大値とする
		index[i] = 0;								//ボクセルに対するリストの個数の初期化
		for(j = 0; j <= listsuu; placeID[i][j++] = INT_MAX);	//初期化，INT_MAXは探索終了コード
	}

	//全てのVOXELどうしの距離を計算し，ARM以内であればリストに登録する
	for(i = 0; i < Maxsuu - 1; i++){
		for(j = i + 1; j < Maxsuu; j++){
			int arm = 0;							//VOXELどうしの距離
			for(r = 0; arm <= ARM && r < dim; r++){
				arm += abs(ID[i][r] - ID[j][r]);
			}
			if(arm <= ARM){						//登録する
				placeID[i][ index[i]++ ]=j;
				placeID[j][ index[j]++ ]=i;
			}
		}
		//いらなくなったIDの削除
		delete [] ID[i];
		
		//メモリの節約のため，不要な部分のリストは削除する．
		int* copy_place=new int[index[i] + 1];
		int* pointer;
		for(j = 0; j <= index[i]; j++) copy_place[j] = placeID[i][j];
		pointer = placeID[i], placeID[i] = copy_place, copy_place = pointer;
		delete [] copy_place;
	}

	//i == Maxsuu - 1に対する処理
	i = Maxsuu - 1;
	//いらなくなったIDの削除
	delete [] ID[i];
	//メモリの節約のため，不要な部分のリストは削除する．
	int* copy_place=new int[index[i] + 1];
	int* pointer;
	for(j = 0; j <= index[i]; j++) copy_place[j] = placeID[i][j];
	pointer = placeID[i], placeID[i] = copy_place, copy_place = pointer;
	delete [] copy_place;

	//その他の要素の削除
	delete [] ID;
	delete [] index;

	//get_placeによって呼び出される軌道点のリストの宣言
	ret_list=new int[NMAX + 1];

	//過去に呼び出されたVOXEL番号の新規登録
	past_get_place = -1;
}

//自分の近傍のVOXELのリストの最大個数を再帰ループによりカウントする
//　j : ループの深さ（次元）
//　k : リストの個数
//　length : VOXELどうしの距離
void CKVoxel::CountListMax(int & j,int & k,int & length){
	int i;
	j++;
	for(i = -ARM; i <= ARM; i++){
		length += abs(i);
		if(length <= ARM){
			if(j == dim){
				k++;
			}else{
				CountListMax(j, k, length);
			}
		}
		length -= abs(i);
	}
	j--;
}

// デコンストラクタ
CKVoxel::~CKVoxel()
{
	int i;
	for(i = 0; i < Maxsuu; i++){
		delete [] placeID[i];
	}
	delete [] placeID;
	delete [] ret_list;
	delete [] idx;
	delete [] ix;
	delete [] voxel;
}

// 1次元時系列データの埋め込み
void CKVoxel::SavePlace1()
{
	int i, j, p, k, q;
	double dmin, dmax, length;

	//データの最大・最小値の計算
	dmin = dmax = x1[0];
	for(i = -(dim-1)*tau;i < NMAX; i++){
		if(x1[i] < dmin) dmin = x1[i];
		else
		if(x1[i] > dmax) dmax = x1[i];
	}

	//アトラクタの直径から分割距離の逆数を計算
	length = DIVISION / (dmax * 1.01 - dmin);

	//各軌道点からVOXELの座標を計算する
	for(i = 0; i < NMAX; i++){
		p = i;
		j = 0;
		for(k = 1; k <= dim / 5 + 1; k++){				//5次元ごとに一まとめにする
			q = 0;
			for(   ; j < k * 5 && j < dim; j++, p-=tau)
				q = q * DIVISION + (int)((x1[p] - dmin) * length);
			placeID[i][k - 1] = q;					//IDの登録
		}
	}
}

// 多次元配列データの埋め込み
void CKVoxel::SavePlace2()
{
	int i, j, k, q;
	double dmin, dmax, length;

	//データの最大・最小値の計算
	dmin = dmax = x2[0][0];
	for(j=0; j < dim; j++)
		for(i=0; i < NMAX; i++){
			if(x2[j][i] < dmin) dmin = x2[j][i];
			else 
			if(x2[j][i] > dmax) dmax = x2[j][i];
		}

	//アトラクタの直径から分割距離の逆数を計算
	length = DIVISION / (dmax * 1.01 - dmin);

	//各軌道点からVOXELの座標を計算する
	for(i = 0; i < NMAX; i++){
		j = 0;
		for(k = 1; k <= dim / 5 + 1; k++){				//5次元ごとに一まとめにする
			q = 0;
			for(   ; j < k * 5 && j < dim; j++)
				q = q * DIVISION + (int)((x2[j][i] - dmin) * length);
			placeID[i][k - 1] = q;					//IDの登録
		}
	}
}

//時刻tの点を中心にした点群をゲットする
//　t : 現在の時刻
// arm : ARMの近傍点の探索を更に arm 回続ける場合
// all : VOXELを使わず，全ての点を返し，全点探索とする場合
void CKVoxel::GetPlace(int t, int arm, bool all)
{
	// VOXELを使わず，全ての点を返し，全点探索とする場合
	if(all){
		if(past_get_place == -1){		// 新規設定で全点探索
			ret_list[NMAX] = INT_MAX;				// 終了コード
			for(int j = NMAX; j; ) --j,ret_list[j] = j;
		}
		return;
	}

	//一回前と同じ場所を中心に検索している場合は，同じものを使ってもらう
	int j;						// tに対応するVOXEL番号
	if((j=idx[t])past_get_place) return;		// 一回前と同じ場所を中心に検索している場合
	
	int i, p, k, q, r, lastpos;
	past_get_place=j;					// 今回のVOXEL番号を登録しておく

	int* list=new int[Maxsuu];				//自分以外のVOXELを探索する時のVOXELのリスト
	for(i=Maxsuu; i; list[--i] = -1);			//探索終了キーを全て記載する
	lastpos=past_get_place;					//近隣であると登録されたVOXELの最終番号
	
	//自分の近傍VOXELを登録
	for(k=0; (q=placeID[past_get_place][k]) < INT_MAX; list[lastpos]=q, lastpos=q, k++);

	//更に腕を伸ばして，ほかの点も探しに行く場合（通常は必要ありません）
	if(arm >= 2){
		register int loop, top, bottom;
		for(bottom = past_get_place, loop = 2; loop <= arm; loop++){
			if((top = list[bottom]) < 0) break;
			for(bottom = lastpos, j = top;  ; j = list[j]){
				//未登録情報の検索
				for(k = 0; (q = placeID[j][k]) < INT_MAX; k++)
					if(list[q] < 0 && lastpos != q)		//まだ登録していなかったら
						list[lastpos] = q, lastpos = q;
				if(j == bottom) break;
			}
		}
	}

	//時刻tに対応するVOXELに登録されている軌道点を全てゲットする
	//近傍のVOXELに登録されている軌道点を全てゲットする
	for(j = 0, k = past_get_place; k >= 0; k = list[k])
		for(p = voxel[k], ret_list[j++] = p; ix[p] < INT_MAX; p = ix[p], ret_list[j++] = p);

	//リストの削除
	delete [] list;

	//探索終了コード
	ret_list[j]=INT_MAX;
}

言語は一応C++です．（C言語も混ざっていますが）．
バグがありましたら連絡してください．お願いします．
GetPlace()関数は更に高速検索が可能になるように改良しました． May, 21, 2002

実例

#include "KVoxel.h"					//宣言

ある関数{
	int DIMENSION;
	int m_Nmax;
	int tau;
	double* m_data = new double[m_Nmax];
	double* data = &m_data[(DIMENSION - 1) * tau];

	CKVoxel* kvoxel;

	//例えばこんな感じでVOXELを呼び出す
	kvoxel = new CKVoxel(DIMENSION, m_Nmax - (DIMENSION - 1) * tau, tau, data);

	int t, j, n;
	for(t = 0; t < m_Nmax - (DIMENSION - 1) * tau; t++){
		kvoxel->GetPlace(t,1,false);			//VOXELの呼び出し
		for(j = 0; (n = kvoxel->ret_list[j]) < INT_MAX; j++){
		
			//data[t]とdata[n]との距離計算などの処理

		}
	}

	delete voxel;
}

クイックソート関数

前述のクイックソートの関数は次のように定義される．

// ポインタとindexを一度にソートする関数
//　viewはdata[][view]をみながらソートするという意味
void CSort::QuickSortRiseIP(int **data, int *index, int bottom, int top, int view)
{
	int i, j, k;
	int x;
	int* f;

	x = data[(bottom + top) >> 1][view], i = bottom, j = top;
	for(;;){
		while(data[i][view] < x) i++;
		while(x < data[j][view]) j--;

		if(i >= j) break;

		//入れ替え処理
		f = data[i], data[i] = data[j], data[j] = f;
		k = index[i], index[i] = index[j], index[j] = k;

		i++, j--;
	}
	if(bottom < i - 1) QuickSortRiseIP(data, index, bottom, i-1, view);
	if(j + 1 < top)    QuickSortRiseIP(data, index,  j + 1, top, view);
}

戻る

2002/5
by Hiroyuki Koga