ヒルベルト曲線について

    • ご無沙汰しております(´・ω・`) なんとなく、前世紀の大数学者ダフィット・ヒルベルトさんの考えた空間充填曲線であるhilbert曲線を三次元においてBlenderというフリーのレンダリングソフトで表現したいなぁという思いから、大学時代のTurboPascal以来というほどに、プログラム言語をいじってその座標点を求めるスプリクトをC言語ベースで書いた。
    • もちろん、自分用なので細かいエラー判定などは省略している。
    • コンパイル環境はVisual Windows for BorlandC++である。C言語はきちんと勉強してないので適当な部分も多いかもしれないので気にしないでください(´・ω・`)

#include
#include
#include
#include

// 三次元空間におけるHilbert曲線の頂点座標{x,y,z}を計算するアルゴリズムである
// Hilbert曲線とは空間充填曲線といわれるフラクタル図形の一種で、位数nのHilbert曲線をH(n)とすると、
// H(n+1)は適当に回転か平行移動した8つ(2の3乗個)のH(n)を最小単位辺(この場合は1)の線分7本でつなげた図形である
// H(n)の極限を考えると、空間が一筆書きで埋めることができる

/* 構造体宣言 */
//Hilbert曲線の頂点座標{x,y,z}の定義
//一般に、H(n)は2^n個の{0,1,2,3,……,(2^n)-1}の組で8^n=(2^3)^n=2^(3n)=(2^n)^3個の頂点がある(※ ^は累乗)
struct point{
int x,y,z;
};

/* プロトタイプ宣言 */
int power(int m, int n);
//mのn累乗(m,n>=0)

void f(int n, struct point *p, struct point *q, int m);

//Hilbert曲線の位数を一段階上げる写像f(n)の定義
// f(1) f(2) f(3)
// H(1) → H(2) → H(3) → H(4) → …
//     f(n-1) f(n) f(n+1)
//  … → H(n-1) → H(n) → H(n+1) → H(n+2) → …
// fi(n)(p)=q(i=1,2,…,8 n=1,2,3…) m:要素組数


/* メイン */
int main(void)
{
  struct point p[8] ={{0,0,0},{0,0,1},{1,0,1},{1,0,0},{1,1,0},{1,1,1},{0,1,1},{0,1,0}};
  //スタートするグレイコード(x,y,z)を表す位数1のHilbert曲線H(1)
//グレイコードは隣り合う座標の差が1であるような数の組。Hilbert曲線の頂点座標はそれらで構成されている

int ORDER; //Hilbert曲線H(n)の位数(1以上。上限は7ぐらいがよいが現在は6まで)
int WRAPPED; //表示する頂点座標の折り返し数
int i,j,k,l,m,n,idx;
struct point *q,*r,*s,*t,*u;

/* 出力ファイルの設定 */
FILE *file;
file = fopen("Hilbert test.txt", "w");


/* hilbert曲線の位数と表示折り返しの入力 */
printf("三次元hilbert曲線の頂点座標をHilbert test.txt(CSVファイル)に出力します");
printf("\n");
printf("求めるhilbert曲線の位数を入力してください(1以上6以下の整数)>");
scanf("%d" , &ORDER);
printf("\n");
printf("何点ごとに座標を折り返し表示するか入力してください(8以上16ぐらいの4の倍数推奨)>");
scanf("%d" , &WRAPPED);
printf("\n");

/* メモリの確保 */
q = (struct point*)malloc(sizeof(struct point)*(8*8));
r = (struct point*)malloc(sizeof(struct point)*(8*8*8));
s = (struct point*)malloc(sizeof(struct point)*(8*8*8*8));
t = (struct point*)malloc(sizeof(struct point)*(8*8*8*8*8));
u = (struct point*)malloc(sizeof(struct point)*(8*8*8*8*8*8));

/* 座標の表示とCSVファイル出力 */
if (ORDER == 1) // ORDER = 1なら H(1)は2^1個の{0,1}の組で8(=8^1=2^3)点
{ for (i = 0; i < 8; i++){
idx = i;
printf( " %04d:(%2d,%2d,%2d)", idx+1 , p[idx].x, p[idx].y, p[idx].z);
if( 0 == (idx +1) % WRAPPED) printf("\n");
fprintf(file,"%d,%d,%d\n", p[idx].x, p[idx].y, p[idx].z);
}
}
else if (ORDER == 2) // ORDER = 2なら H(2)は2^2個の{0,1,2,3}の組で64(=8^2=(2^3)^2=2^6=4^3)点
{ f(1, p, q, 8);
for (j = 0 ; j < 8; j++){
for (i = 0 ; i < 8; i++){
idx = j*8+i;
printf(" %04d:(%2d,%2d,%2d)", (idx +1) , q[idx].x, q[idx].y, q[idx].z);
if( 0 == (idx +1) % WRAPPED) printf("\n");
fprintf(file,"%d,%d,%d\n", q[idx].x, q[idx].y, q[idx].z);
}
}
}
else if (ORDER == 3) // ORDER = 3なら H(3)は2^3個の{0,1,2,……,7}の組で512(=8^3=(2^3)^3=2^9=8^3)点
{ f(1, p, q, 8);
f(2, q, r, 8*8);
for (k = 0 ; k < 8; k++){
for (j = 0 ; j < 8; j++){
for (i = 0 ; i < 8; i++){
idx = k*64+j*8+i;
printf(" %04d:(%2d,%2d,%2d)", (idx +1) , r[idx].x, r[idx].y, r[idx].z);
if( 0 == (idx +1) % WRAPPED) printf("\n");
fprintf(file,"%d,%d,%d\n", r[idx].x, r[idx].y, r[idx].z);
}
}
}
}
else if (ORDER == 4) // ORDER = 4なら H(4)は2^4個の{0,1,2,3,……,15}の組で4096(=8^4=(2^3)^4=2^12=16^3)点
{ f(1, p, q, 8);
f(2, q, r, 8*8);
f(3, r, s, 8*8*8);free(q);
for (l = 0 ; l < 8; l++){
for (k = 0 ; k < 8; k++){
for (j = 0 ; j < 8; j++){
for (i = 0 ; i < 8; i++){
idx = l*512+k*64+j*8+i;
printf(" %06d:(%2d,%2d,%2d)", (idx +1) , s[idx].x, s[idx].y, s[idx].z);
if( 0 == (idx +1) % WRAPPED) printf("\n");
fprintf(file,"%d,%d,%d\n", s[idx].x, s[idx].y, s[idx].z);
}
}
}
}
}
else if (ORDER == 5) // ORDER = 5なら H(5)は2^5個の{0,1,2,3,……,31}の組で3万2768(=8^5=(2^3)^5=2^15=32^3)点
{ f(1, p, q, 8);
f(2, q, r, 8*8);
f(3, r, s, 8*8*8);free(q);
f(4, s, t, 8*8*8*8);free(r);
for (m = 0 ; m< 8; m++){
for (l = 0 ; l < 8; l++){
for (k = 0 ; k < 8; k++){
for (j = 0 ; j < 8; j++){
for (i = 0 ; i < 8; i++){
idx = m*4096+l*512+k*64+j*8+i;
printf(" %06d:(%2d,%2d,%2d)", (idx +1) , t[idx].x, t[idx].y, t[idx].z);
if( 0 == (idx +1) % WRAPPED) printf("\n");
fprintf(file,"%d,%d,%d\n", t[idx].x, t[idx].y, t[idx].z);
}
}
}
}
}
}
else if (ORDER == 6) // ORDER = 6なら H(6)は2^6個の{0,1,2,3,……,63}の組で26万2144(=8^6=(2^3)^6=2^18=64^3)点
{ f(1, p, q, 8);
f(2, q, r, 8*8);
f(3, r, s, 8*8*8);free(q);
f(4, s, t, 8*8*8*8);free(r);
f(5, t, u, 8*8*8*8*8);free(s);
for (n = 0 ; n< 8; n++){
for (m = 0 ; m< 8; m++){
for (l = 0 ; l < 8; l++){
for (k = 0 ; k < 8; k++){
for (j = 0 ; j < 8; j++){
for (i = 0 ; i < 8; i++){
idx = n*32768+m*4096+l*512+k*64+j*8+i;
printf(" %06d:(%2d,%2d,%2d)", (idx +1) , u[idx].x, u[idx].y, u[idx].z);
if( 0 == (idx +1) % WRAPPED) printf("\n");
fprintf(file,"%d,%d,%d\n", u[idx].x, u[idx].y, u[idx].z);
}
}
}
}
}
}
}
free(q);
free(r);
free(s);
free(t);
free(u);


/* 座標の数の表示 */
printf("\n");
printf("位数%dのhilbert曲線の頂点座標:%d個\n", ORDER, power(8,ORDER));
printf("\n");
printf("終了するにはなにかキーを押してください");
getch();
fclose(file);
return 0;
}


/* mのn累乗(m,n>=0) */
int power(int m, int n)

{
return ( (n == 0) ? 1: m*power(m,n-1) );
}

/* Hilbert曲線の位数を一段階上げる写像f(n) */
void f(int n, struct point *p, struct point *q, int m)
{
int i;
struct point *p_w;//作業用の座標

   /* メモリの確保 */
p_w = (struct point*)malloc(sizeof(struct point)*(m*8));

for (i = 0 ; i < m ; i++){
if (n == 1)//f1,8-yzx f2,3,6,7-zyx f4,5-xzy:H(1)→H(2)
{ p_w[i].y = p[i].y;
p_w[i].z = p[i].z;
p_w[i].x = p[i].x;}
else if (n == 2)//f1,8-zxy f2,3,6,7-xzy f4,5-yxz:H(2)→H(3)
{ p_w[i].y = p[i].z;
p_w[i].z = p[i].x;
p_w[i].x = p[i].y;}
else if (n == 3)//f1,8-zyx f2,3,6,7-yzx f4,5-xyz:H(3)→H(4)
{ p_w[i].y = p[i].z;
p_w[i].z = p[i].y;
p_w[i].x = p[i].x;}
else if (n == 4)//f1,8-zyx f2,3,6,7-yzx f4,5-xyz:H(4)→H(5)
{ p_w[i].y = p[i].z;
p_w[i].z = p[i].y;
p_w[i].x = p[i].x;}
else if (n == 5)//f1,8-zxy f2,3,6,7-xyz f4,5-yxz:H(5)→H(6)
{ p_w[i].y = p[i].z;
p_w[i].z = p[i].x;
p_w[i].x = p[i].y;}

    //f1(n)
    q[i].x = p_w[i].y;
q[i].y = p_w[i].z;
q[i].z = p_w[i].x;
//f2(n)
q[i+m].x = p_w[i].z + power(2,n);
q[i+m].y = p_w[i].y;
q[i+m].z = p_w[i].x;
//f3(n)
q[i+2*m].x = p_w[i].z + power(2,n);
q[i+2*m].y = p_w[i].y + power(2,n);
q[i+2*m].z = p_w[i].x;
//f4(n)
q[i+3*m].x = -p_w[i].x + power(2,n) -1;
q[i+3*m].y = -p_w[i].z + power(2,n+1) -1;
q[i+3*m].z = p_w[i].y;
//f5(n)
q[i+4*m].x = -p_w[i].x + power(2,n) -1;
q[i+4*m].y = -p_w[i].z + power(2,n+1) -1;
q[i+4*m].z = p_w[i].y + power(2,n);
//f6(n)
q[i+5*m].x = p_w[i].z + power(2,n);
q[i+5*m].y = -p_w[i].y + power(2,n+1) -1;
q[i+5*m].z = -p_w[i].x + power(2,n+1) -1;
//f7(n)
q[i+6*m].x = p_w[i].z + power(2,n);
q[i+6*m].y = -p_w[i].y + power(2,n) -1;
q[i+6*m].z = -p_w[i].x + power(2,n+1) -1;
//f8(n)
q[i+7*m].x = -p_w[i].y + power(2,n) -1;
q[i+7*m].y = p_w[i].z;
q[i+7*m].z = -p_w[i].x + power(2,n+1) -1;
}
free(p_w);

}