ヒルベルト曲線順

ヒルベルト曲線は空間充填曲線のひとつだ。このヒルベルト曲線でより高次の空間をたどるとキャッシュ効率がよくなるらしい。

さっそく３次元のモデルの面をヒルベルト曲線順で並び替えてみる。
以下はstanford bunnyのそれぞれの面の順番をHSVのH（色相）に対応させて表示させたものだ。
デフォルトの順番は次のようになる。

かなりめちゃくちゃだ。位置が近ければ近い色に表示されるのが望ましい。
コレを３次元空間をヒルベルト曲線順でたどってみると次のようになる。

だいぶ位置が近い順になった。ただ境目（順番が飛ぶ場所）が結構多い。そこでヒルベルト曲線の階層を１階層にとどめると次のようになる。

なんだか僕にはこちらの方が自然に思える・・・
ヒルベルト曲線順がどういう効果をもたらすか２次元画像をヒルベルト曲線順で色付けして確認してみよう。

どうやらヒルベルト曲線順にすると

• 大きな境目は少なくなる
• 境目の数は多くなる

これは大きく順番が飛ぶ場所を無くす変わりに順番が飛ぶ場所を増やして分散させるといった効果があるようだ。

でキャッシュ効率なんだけど、どうもいまひとつ効果の違いを見出す実験が発見できなかった。自作レイトレーサとOpenGLレンダラでやってみたんだけれども。
そんなわけでここには結果は載せない。
ただ、キャッシュ効率なんてのはプラットフォームしだいで全然変わってくるものだから、プラットフォームが変わるたびに実験して計測してみるといいかもしれない。

以下にソースを示す。
次は面をヒルベルト曲線順に並べ替えるルーチン

struct mface{
	size_t i0;
	size_t i1;
	size_t i2;
};

struct face_id{
  size_t i;//
  vector3 pos;
};

typedef std::vector<face_id>::iterator iterator;

struct sorter : public std::binary_function<face_id,face_id,bool>{
  sorter(int DIR){
    plane_ = DIR>>1;
    inv_ = (DIR&1);
  }
  bool operator()(const face_id& a, const face_id& b)const{
    if(inv_){
      return a.pos[plane_] < b.pos[plane_];
    }else{
      return a.pos[plane_] > b.pos[plane_];
    }
  }
  int plane_;
  bool inv_;
};

void hilbert_sort(iterator a, iterator b, int L, int R, int D, int U, int B, int F);

void hilvert_next(iterator a, iterator b, int L, int R, int D, int U, int B, int F, int Bit)
{
#if 1
//ここを無効にすると１階層のみになる
	switch(Bit){
		case 0: hilbert_sort(a,b, D,U,B,F,L,R);break;//0,0,0
		case 1: hilbert_sort(a,b, B,F,L,R,D,U);break;//0,0,0
		case 3: hilbert_sort(a,b, B,F,L,R,D,U);break;//0,0,0
		case 2: hilbert_sort(a,b, R,L,U,D,B,F);break;//0,0,0
		case 6: hilbert_sort(a,b, R,L,U,D,B,F);break;//0,0,0
		case 7: hilbert_sort(a,b, F,B,L,R,U,D);break;//0,0,0
		case 5: hilbert_sort(a,b, F,B,L,R,U,D);break;//0,0,0
		case 4: hilbert_sort(a,b, D,U,F,B,R,L);break;//0,0,0
	}
#endif

}

void hilvert_x(iterator a, iterator b, int L, int R, int D, int U, int B, int F, int Bit)
{
  if(Bit == 2||Bit == 4){
    std::sort(a,b, sorter(L) );
    iterator m = a+ ((b-a)>>1);
    hilvert_next(a,m,L,R,D,U,B,F,Bit|1);
    hilvert_next(m,b,L,R,D,U,B,F,Bit);
  }else{
    std::sort(a,b, sorter(R) );
    iterator m = a+ ((b-a)>>1);
    hilvert_next(a,m,L,R,D,U,B,F,Bit);
    hilvert_next(m,b,L,R,D,U,B,F,Bit|1);
  }  
}

void hilvert_y(iterator a, iterator b, int L, int R, int D, int U, int B, int F, int Bit)
{
	if(Bit == 4){
		std::sort(a,b, sorter(D) );
		iterator m = a+ ((b-a)>>1);
		hilvert_x(a,m,L,R,D,U,B,F,Bit|2);
		hilvert_x(m,b,L,R,D,U,B,F,Bit);
	}else{
		std::sort(a,b, sorter(U) );
		iterator m = a+ ((b-a)>>1);
		hilvert_x(a,m,L,R,D,U,B,F,Bit);
		hilvert_x(m,b,L,R,D,U,B,F,Bit|2);
	}
}

void hilvert_z(iterator a, iterator b, int L, int R, int D, int U, int B, int F, int Bit)
{
	  std::sort(a,b, sorter(F) );
	  iterator m = a+ ((b-a)>>1);
	  hilvert_y(a,m,L,R,D,U,B,F,Bit);
	  hilvert_y(m,b,L,R,D,U,B,F,Bit|4);
}

void hilbert_sort(iterator a, iterator b, int L, int R, int D, int U, int B, int F)
{
	if(b-a<=1)return;	
	hilvert_z(a,b,L,R,D,U,B,F,0);
}

void hilbert_sort(iterator a, iterator b)
{
	hilbert_sort(a,b,0,1,2,3,4,5);
}


//面をヒルベルト曲線順に並べ替える関数
//mfが面の配列
//mvが頂点配列
void hilbert_sort_face(std::vector<mface>& mf, std::vector<vector3>& mv)
{
	if(mf.empty())return ;
	size_t fsz = mf.size();


	std::vector<face_id> tmp(fsz);
	for(size_t i = 0;i<fsz;i++){
		tmp[i].i = i;
		tmp[i].pos = real(1.0/3)*(mv[mf[i].i0]+mv[mf[i].i1]+mv[mf[i].i2]);
	}
	hilbert_sort(tmp.begin(), tmp.end());

	std::vector<mface> other(fsz);
	for(size_t i = 0;i<fsz;i++){
		other[i] = mf[tmp[i].i];
	}

	mf.swap(other);
}

以下は２次元のヒルベルト曲線順を求めるルーチン

class hilbert2D{
public:
	hilbert2D(int order)
	{
		int wid = 1<<order;
		o = order;
		x = 0;
		y = 0;
		i = 0;	
		w = wid;
		v.resize(w*w);
	}

	void step(int dir){
		switch(dir&3){
			case 0:x=x+1;break;
			case 1:y=y+1;break;
			case 2:x=x-1;break;
			case 3:y=y-1;break;
		}

		v[w*y+x]=i;
		i++;
	}
	int operator[](int id)const{return v[id];}
	size_t size()const{return v.size();}

	void hilbert(int dir, int rot, int order){
		if(order==0)return;
		dir = dir + rot;
		hilbert( dir, -rot, order-1);
		step(dir);
		dir = dir - rot;
		hilbert( dir,  rot, order-1);
		step(dir);
		hilbert( dir,  rot, order-1);
		dir = dir - rot;
		step(dir);
		hilbert( dir, -rot, order-1);
	}

	void hilbert(){
		x = -1;
		y = 0;
		i=0;
		step(0);
		hilbert(0,1,o);
	}

private:
	int o;
	int w;
	std::vector<int> v;  
	int x;
	int y;
	int i;
};
//使用例
//hilbert2D h2d(order);
//h2d.hilbert();
//int w = 1<<order;
//for(int y=0;y<w;y++){
//	for(int x=0;x<w;x++){
//		int i = w*y+x;
//		std::cout << h2d[i] <<", ";//x,yは何番目？
//	}
//	std::cout << std::endl;
//}