いろいろな画素形状での離散ボロノイ図

このページでは、2000年度に行われた、いろいろな形状の画素における離散ボロノイ図の研究についてまとめる。

目次

付録

離散ボロノイ図の研究｜村島研究室

2次元直交座標系（平面座標系）

平面上の位置は、2つの実数値をつかって一意に表すことが出来ます。2次元直交座標系は直交する2つの軸、x軸とy軸をつかって点を表現します。

離散ボロノイ図においても、母点座標は実数値としてこの座標系で与えられます。また、格子間の距離を測る場合もこの座標系が用いられます。

戻る

セルインデックス

離散ボロノイ図を作成する空間を格子分割した各々の小領域をセルと呼びます。正方格子、六角格子、三角格子で平面を格子分割すると、平面はセルによって隙間なく埋め尽くされます。2次元平面を埋め尽くすセルは2つの整数値で番号付けすると便利です。この番号をセルインデックスと呼びます。

戻る

セルの代表点

セルは大きさを持ちますが、セル間の距離を計算する場合はセルを代表する点を使った方が便利です。ここではセルの重心をもってその代表点としています。格子点とは異なるので注意。

戻る

セルの大きさ

正方格子、六角格子、三角格子上に作成した離散ボロノイ図を比較することを考えると、平面を公平な精度で分割することが望まれます。たとえば、各セルの一辺の長さを1にした場合、六角格子のセルは大きく、三角格子のセルは小さくなります。

ここでは、同一面積の領域を同数のセルに分割することが平等な精度であるという考えに基づき、各セルの面積が1となるようにしました。

戻る

格子座標系

セルインデックス i,j を目盛りに持つような座標系を格子座標系と呼びます。格子座標系の取り方は一般に一通りではありません。たとえば、正方格子で分割した空間において、j 軸を45度傾けて取ることも可能です。通常は、目的に応じて都合のよい座標系を用いることになります。

戻る

正方格子座標系

正方格子の場合は簡単です。正方格子では正方形のセルが水平・垂直に並んでいるため、格子座標系も直交座標系で表現可能です。

正方格子座標系は、上図のようにインデックス 0,0 のセルの代表点が平面座標系の原点に重なり、インデックス i の増加する方向がx軸方向、インデックス j の増加する方向がy軸方向となるように取ります。

ひとつのセルの面積が1となるようにすると、インデックス 1,0 のセルの代表点が点 (1,0) に、インデックス 0,1 のセルの代表点が点 (0,1) に重なることになります。

正方格子座標のセルインデックスを i,j とし、そのセルの代表点の平面座標を (rx,ry) とすると、

rx = i
ry = j

が成り立ちます。また、セルインデックス i,j から (rx,ry) を得るコードを次に示します。

/*
 * 正方格子の代表点の座標の取得
 */
void getXY(int i,int j,double *x,double *y){
  *x=i;
  *y=j;
}

また、平面上の点 (x,y) が与えられたとき、その点を含むセルのインデックス i,j は次の式で与えられます。

i = floor ( x + 0.5 )
j = floor ( y + 0.5 )

また、(x,y) からセルインデックス i,j を得るコードを以下に示します。

/*
 * 正方格子インデックスの取得
 */
void getIJ(double x,double y,int *i,int *j){
  *i=(int)floor(x+0.5);
  *j=(int)floor(y+0.5);
}

ただし、セルの境界線上に (x,y) があった場合、右または上のセルに所属するものとします。つまり、インデックス i,j のセルの占める領域は、

i - 0.5 ≦ x ＜ i + 0.5
j - 0.5 ≦ y ＜ j + 0.5

です。（下図）

戻る

六角格子座標系

六角格子の場合、1つのセルに隣接するセルは6個あることから、座標軸は直交しない方が自然になります。いろいろな座標軸の取り方が考えられますが、ここでは最も対称性がよいと思われる以下の座標系を採用しました。

六角格子座標系は、上図のようにインデックス 0,0 のセルの代表点が平面座標の原点に重なり、インデックス i の増加する方向がx軸方向であり、インデックス j の増加する方向は i 軸から反時計まわりに60度の方向です。したがって、y軸と j 軸とは一致しません。

また、セルの面積を1に規格化すると、セルの大きさは下図のようになります。

ここで、六角格子に固有の定数として、

HEXUX = 12 -1/4
HEXUY = 108 -1/4

を定義します。

インデックス i,j のセルの代表点の平面座標を (rx,ry) とすると、

rx = ( i * 2 + j ) * HEXUX
ry = j * 3 * HEXUY

が成り立ちます。また、セルインデックス i,j から (rx,ry) を得るコードを次に示します。

/*
 * 六角格子の代表点の座標の取得
 *  HEXUX = pow( 12,-0.25)
 *  HEXUY = pow(108,-0.25)
 */
void getXY(int i,int j,double *x,double *y){
  *x=(i*2+j)*HEXUX;
  *y=j*3*HEXUY;
}

平面内の点 (x,y) が与えられたとき、その点を含むセルのインデックス i,j は次のコードで求められます。

/*
 * 六角格子インデックスの取得
 *  HEXUX = pow( 12,-0.25)
 *  HEXUY = pow(108,-0.25)
 */
void getIJ(double x,double y,int *i,int *j){
  int offx=int(floor(x/HEXUX/4));
  int offy=int(floor(y/HEXUY/6));
  double ux=x/HEXUX-offx*4;
  double uy=y/HEXUY-offy*6;
  int a=int(ux);
  int b=int(uy);
  switch(b){
  case 0: *j=0; *i=(a==0)?0:(a==3)?2:1; break;
  case 1:
    if((a==1 || a==3)?(uy-b >= ux-a):(uy-b >= 1-ux+a)){
      *i=a/2;
      *j=1;
    }
    else{
      *i=(a==0)?0:(a==3)?2:1;
      *j=0;
    }
    break;
  case 2:
  case 3: *j=1; *i=a/2; break;
  case 4:
    if((a==0 || a==2)?(uy-b >= ux-a):(uy-b >= 1-ux+a)){
      *i=(a==0)?-1:(a==3)?1:0;
      *j=2;
    }
    else{
      *i=a/2;
      *j=1;
    }
    break;
  case 5: *j=2; *i=(a==0)?-1:(a==3)?1:0; break;
  }
  *j+=offy*2;
  *i+=offx*2-offy;
}

簡単に説明します。このアルゴリズムは、平面を (4*HEXUX) x (6*HEXUY) の大きさのブロックに分割し、そのブロック単位に与えられた点 (x,y) を平行移動して下図の赤い基本ブロックの中の点 (ux,uy) に持ってきます。

さらに、このブロック内を 4 x 6 の小領域に分割し、平行移動した点 (ux,uy) がそのどの領域に属するかを a,b に求めます。a,b が決まれば、(ux,uy) がどのセルに属するかは何通りかの判断を行えば簡単に求めることが出来ます。

こうして基本ブロック内でどのセルに属しているかを求めた後、最後に平行移動による変化分を補正して、正しいインデックスを求めています。

六角格子の場合も正方格子と同じく、セルの境界線上に点がある場合は右側・上側が優先になります。

上のコードでは、計算を簡単にするために ux,uy を求めるときに、水平・垂直をそれぞれ 4*HEXUX と 6*HEXUY で割って、1x1 の大きさの正方形の内部に正規化しています。

戻る

三角格子座標系

三角格子は隣接するセルが3つしかないので、やはり直交座標系では表しにくくなります。直交座標系や軸がジグザクになる座標系も考えたのですが、対称性という観点からは次のような「穴のあいた」座標系がシンプルでよいと思われます。

この座標系において、赤で書かれたインデックスは対応するセルが存在しないことを意味しています。たとえば、インデックス 2,0 に対応するセルはありません。

これは、i = 2 のセルはすべて存在しないという意味ではありません。2,1 のセルは存在します。

見て分かるとおり、三角格子のセルには左向きのものと右向きのものとがあります。ここでは原点に当たるインデックス 0,0 のセルは左向きと決めています。原点のセルは存在した方と何かと便利ですし、1,0 や 0,1 もあった方がいいかと言うのかその理由です。このため、-1,0 と 0,-1 のセルがありません。

インデックス i,j が与えられたとき、そのセルが存在しているかどうかは次のようにしてわかります。同時に、存在しているなら右向きか左向きかもわかります。

/*
 * 三角格子のセル i,j のタイプ判定
 *  0 : 左向き
 *  1 : 右向き
 *  2 : 存在しない
 */
int getCellType(int i,int j){
  i=i%3;
  j=j%3;
  if(i<0) i+=3;
  if(j<0) j+=3;
  return (i+j)%3;
}

このアルゴリズムは、セルインデックスを、

0 ≦ i ＜ 3
0 ≦ j ＜ 3

の平行四辺形の領域を単位に、この領域に入るように平行移動した後、そのインデックスがどこにあるかで判定を行っています。

上の図のように、赤で囲まれた範囲内には9個のインデックスがありますが、黒で示した3つはセルが存在しません。この位置は、i+j が2になります。緑で示した点は右向きのセルであり、この位置では i+j は1か4ですから、どちらも3で割ると1余ります。青の点は左向きであり、この位置は i+j は0か3で、3で割った余りは0です。

ただし、負の数を割った余りは負の数となるため、余りが負になっている場合は補正しています。

セルの面積を1に規格化すると、セルの大きさは下図のようになります。

ここで、三角格子に固有の定数として、

TRIUX = 3 1/4 / 3
TRIUY = 3 -1/4

を定義します。

インデックス i,j のセルの代表点の平面座標を (rx,ry) とすると、

rx = ( i * 2 - j ) * TRIUX
ry = j * TRIUY

が成り立ちます。

/*
 * 三角格子の代表点の座標の取得
 *  TRIUX = pow(3,0.25)/3
 *  TRIUY = pow(3,-0.25)
 */
void getXY(int i,int j,double *x,double *y){
  *x=(i*2-j)*TRIUX;
  *y=j*TRIUY;
}

平面内の点 (x,y) が与えられたとき、その点を含むセルのインデックス i,j は次のコードで求められます。

/*
 * 三角格子のインデックスの取得
 *  TRIUX = pow(3,0.25)/3
 *  TRIUY = pow(3,-0.25)
 */
void getIJ(double x,double y,int *i,int *j){
  int offx=int(floor((x+TRIUX*2)/TRIUX/6));
  int offy=int(floor(y/TRIUY/2));
  double ux=(x+TRIUX*2)/TRIUX/6-offx;
  double uy=y/TRIUY/2-offy;
  int a=ux<0.5;
  int b=uy<=ux;
  int c=uy>=1-ux;
  if(b && c){ *i=2; *j=1; }
  else if(!b && !c){ *i=0; *j=1; }
  else if(b && !c){ *i=(a)?0:1; *j=0; }
  else{ *i=(a)?1:2; *j=2; }
  *j+=offy*2;
  *i+=offx*3+offy;
}

これも簡単に説明します。このアルゴリズムでは、平面全体を (6*TRIUX) x (2*TRIUY) の大きさのブロックに分割し、そのブロック単位に与えられた点 (x,y) を平行移動して下図の赤い基本ブロックの中の点 (ux,uy) に持ってきます。

さらに (ux,uy) の関係で図の右の a,b,c の条件を求め、各条件の組み合わせでどの位置にあるかを判定します。組み合わせは以下のとおりです。

a b c 位置
- o o 2,1
- x x 0,1
o o x 0,0
x o x 1,0
o x o 1,2
x x o 2,2

最後に offx と offy で補正して元の位置に戻しています。このコードでも、(ux,uy) を求める際に基本ブロックのサイズが 1x1 になるように正規化を行っています。

戻る

境界矩形

ボロノイ図は空間全体に広がり境界を持ちませんが、離散ボロノイ図の場合、メモリの制約から何らかの境界で大きさを制限する必要があり、適当な大きさの矩形を用いて大きさを制限しています。この矩形を境界矩形と呼びます。

境界矩形 (bx1,by1)-(bx2,by2)
(bx1,by1) ：左下の平面座標
(bx2,by2) ：右上の平面座標

境界矩形を利用して、どのセルを処理対象とするかが決まります。格子の形状によってインデックスの並びが異なりますが、次のような基本インデックス base_i,base_j と幅 width 高さ height を用いて処理対象となるセルを決定することが出来ます。

基本インデックス base_i,base_j
(bx1,by1) を含むセル（基本セル）のインデックス
幅 width
(bx2,by1) を含むセルのインデックスを i,j としたとき i - base_i + 1
高さ height
(bx1,by2) を含むセルのインデックスを i,j としたとき j - base_j + 1

基本セルから始めて、すべてのセルを1回ずつ取り出すことをセルの反復と呼びます。

正方格子の場合、境界矩形とセルの関係は次のようになります。

六角格子の場合は次のようになります。

上の図を見ると分かるように、六角格子のセルには、上下左右にその一部が境界矩形の内部にありながら処理対象セルにならないものが存在します。したがって、矩形内にランダムに母点を分布させるような場合、母点に対応するセルが必ず処理対象となるようにするためには、母点の分布する矩形領域を以下のようにしなければなりません。

(bx1+HEXUX,by1+HEXUY) - (bx2-HEXUX,by2-HEXUY)

三角格子の場合は次のようになります。

または、

上の2つの図の違いは、基本セルが左向きか右向きかという点です。三角格子を用いる場合、インデックスに穴があいているため、幅 width の計算に注意する必要があります。基本セルが左向きの場合に、境界の右端のセルが右向きであった場合、幅を1つ増やさなければなりません。そうしないと、そのひとつ上の行の右端のセルが処理対象とならないからです。具体的には以下のコードが必要となります。

if(!base_right && getCellType(base_i+width-1,base_j)==1)
  width++;

この場合は、上下の境界矩形内に処理対象とならないセルが存在します。したがって、矩形内にランダムに母点を分布させるような場合、母点に対応するセルが必ず処理対象となるようにするためには、母点の分布する矩形領域を以下のようにしなければなりません。

(bx1,by1+TRIUY) - (bx2,by2-TRIUY)

戻る

セルの矩形領域内判定

セル i,j がこの離散ボロノイ図の境界矩形内のセルかどうかを判定する方法について考えます。この判定は、離散ボロノイ図の外側にあるセルにアクセスしないようにするために必要です。

ここでは、離散ボロノイ図全体の中にあるかどうかだけでなく、セル bi,bj から w,h で示される範囲に限定して、インデックス i,j のその領域内にあるかどうかを判定するメソッドを定義します。

正方格子の場合は、j 軸の傾きがないので単純です。

/*
 * 矩形領域内判定
 */
int inside(int i,int j,int bi,int bj,int w,int h){
  return bi<=i && i<bi+w && bj<=j && j<bj+h;
}
int inside(int i,int j){
  return inside(i,j,base_i,base_j,width,height);
}

六角格子では、j 軸の傾きによる補正が必要です。

/*
 * 矩形領域内判定
 */
int inside(int i,int j,int bi,int bj,int w,int h){
  int bi1=bi-(j-bj)/2;
  return bi1<=i && i<bi1+w && bj<=j && j<bj+h;
}
int inside(int i,int j){
  return inside(i,j,base_i,base_j,width,height);
}

三角格子の場合、補正の向きが逆になり、さらに基本セルの向きも関係します。

/*
 * 矩形領域内判定
 */
int inside(int i,int j,int bi,int bj,int w,int h,int br){
  int bi1=bi+(j-bj+br)/2;
  return bi1<=i && i<bi1+w && bj<=j && j<bj+h;
}
int inside(int i,int j){
  return inside(i,j,base_i,base_j,width,height,base_right);
}

ここで、base_right は基本セルが右向きかどうかを表すフラグで、

base_right = getCellType(base_i,base_j)

です。

戻る

セル配列

セルの情報は配列として記憶します。C 言語の配列は基本的に1次元配列ですから、セルインデックス i,j を配列のインデックス a に変換する、またはその逆が必要となります。

ここでは、セルインデックスの値の妥当性に対するチェックは行っていません。

セルの状態を配列として準備する場合、正方格子と六角格子では使わないインデックスはないので width x height の大きさの配列を準備する必要があります。ここでは、以下のようにセル配列が確保されているとして説明します。

int cell_size=width*height;
unsigned int *cell=new unsigned int[cell_size];

/*
 * セルへのアクセス
 * OB : 領域外をあらわすコード
 */
void setCell(int a,unsigned int c){
  if(a<0 || a>=cell_size) return;
  cell[a]=c;
}
unsigned int getCell(int a){
  if(a<0 || a>=cell_size) return OB;
  return cell[a];
}

正方格子の場合は、次のコードで相互に変換できます。

/*
 * セルインデックス i,j から配列のインデックス a へ
 */
int c2a(int i,int j){
  if(!inside(i,jt)) return -1;
  int ai=i-base_i;
  int aj=j-base_j;
  return ai+aj*width;
}

/*
 * 配列のインデックス a からセルインデックス i,j へ
 */
int a2c(int a,int *i,int *j){
  if(a<0 || a>=cell_size) return 0;
  *i=a%width+base_i;
  *j=a/width+base_j;
  return 1;
}

c2a() では、与えられたインデックス i,j が範囲外にある場合 -1 を返します。

六角格子の場合は、インデックス j 軸が傾いているので、行が上がってくると i の補正が必要となります。j が2増えるごとに i が1ずつ減少しますから、その分を補正すると、

/*
 * セルインデックス i,j から配列のインデックス a へ
 */
int c2a(int i,int j){
  if(!inside(i,j)) return -1;
  int aj=j-base_j;
  int ai=i-base_i+aj/2;
  return ai+aj*width;
}

/*
 * 配列のインデックス a からセルインデックス i,j へ
 */
int a2c(int a,int *i,int *j){
  if(a<0 || a>=cell_size) return 0;
  int ai=a%width;
  int aj=a/width;
  *i=ai+base_i-aj/2;
  *j=aj+base_j;
  return 1;
}

となります。

三角格子の場合、格子座標に穴があいているため、幅が width で高さが height のであっても、実際には width x height よりも少ない数のセルしか存在しません。行数は height ですが、各行のセル数は基本セルの向きによって次のように変化します。

基本セル width
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 …
左向き 1 2 - 3 4 - 5 6 - 7 8 - 9 10 - 11 12 - …
右向き 1 - 2 3 - 4 5 - 6 7 - 8 9 - 10 11 - 12 …

ハイフン ( - ) は穴の部分のインデックスで、対応するセルが存在しないことを意味します。この表から、

1行のセル数 ta_width = width - int ( width / 3 )

であることが分かります。したがって、三角格子のセル配列は、ta_width x height の大きさあれば十分であることになります。以下では、次のようにセル配列が確保されているとして説明します。

int cell_size=ta_width*height;
Cell cell=new Cell[cell_size];

/*
 * セルへのアクセス
 * OB : 領域外をあらわすコード
 */
void setCell(int a,int c){
  if(a<0 || a>=cell_size) return;
  cell[a]=c;
}
int getCell(int a){
  if(a<0 || a>=cell_size) return OB;
  return cell[a];
}

次に、基準を 0 に揃えた場合の水平方向のインデックス i と配列のインデックス n の関係は、以下のようになります。

この図から、左端のセルが右向きか左向きかによって i と n の関係は異なっていますが、ピンクで示した穴のインデックスは使用されることはないので、i から n を求める場合はどちらも、

n = i - int ( ( i + 1 ) / 3 )

という式で表現できることが分かります。しかし、逆に n から i を求める場合、左端のセルが左向きか右向きかが影響します。いま、r は左端のセルが右向きならば 1 、左向きなら 0 であるものとします。このとき、n から i を得る式は、

i = n + int ( n / 2 ) + ( ( n % 2 ) ? r : 0 )

となります。

各行の左端のセルの向きは、基本セルの向きと行番号によって決まります。j - base_j は最下行からの行数を表し、これが偶数ならばその行の左端のセルの向きは基本セルの向きと同じ、そうでないなら逆になります。したがって、基本セルの向きを base_right で表すとすると、インデックス j の行の左端のセルの向き r は、

r = ( ( j - base_j ) % 2 ) ? ( 1 - base_right ) : ( base_right )

で表されます。

三角格子のインデックス j 軸は傾いているため、行が上がっていくと i 方向の補正が必要になりますが、これにも基本セルの向きが関係します。基本セルが左向きの場合、j が2増えるごとに 1 づつ増加します。基本セルが右向きの場合、最初の1回だけ、j が1増えると i も1増え、その後は j が2増えるごとに 1 づつ増加します。これを式で表すと、j に対する i の変化量は、

int ( ( j - base_j + base_right ) / 2 )

で表されます。よってその分を補正すると、

/*
 * セルインデックス i,j から配列のインデックス a へ
 */
int c2a(int i,int j){
  if(!inside(i,j)) return -1;
  int aj=j-base_j;
  int ai=i-base_i-(aj+base_right)/2;
  ai-=(ai+1)/3;
  return ai+aj*ta_width;
}

/*
 * 配列のインデックス a からセルインデックス i,j へ
 */
void a2c(int a,int *i,int *j){
  if(a<0 || a>=cell_size) return 0;
  int ai=a%ta_width;
  int aj=a/ta_width;
  int r=((aj-base_j)%2)?(1-base_right):(base_right);
  ai=ai+ai/2+((ai%2)?r:0);
  *i=ai+base_i+(aj+base_right)/2;
  *j=aj+base_j;
  return 1;
}

となります。

戻る

セルの反復

正方格子の場合、base_i,base_j から始めて、水平方向に width 回の反復を垂直方向に height 回繰り返せばよいことになるので、以下のコードですべてのセルにアクセスすることが出来ます。

for(int j=base_j;j<base_j+height;j++){
  for(int i=base_i;i<base_i+width;i++){

    // ... 処理

  }
}

六角格子の場合、base_i,base_j から横に width 個のセルが並んでいますが、j 軸が傾いているため、2回に1回は i の開始位置を左に1つずらしてやる必要があります。したがって、反復のためのコードは次のようになります。

int bi=base_i,c=0;
for(int j=base_j;j<base_j+height;j++){
  for(int i=bi;i<bi+width;i++){

    // ... 処理

  }
  if((++c)%2==0) bi--;
}

三角格子の場合、基本セルの向きによって j が i に与える影響が異なります。基本セルが右向きの場合、1回早く i の補正が始まります。したがって反復のためのコードは次のようになります。

int bi=base_i,c=base_right;
for(int j=base_j;j<base_j+height;j++){
  for(int i=bi;i<bi+width;i++){

    // ... 処理

  }
  if((++c)%2==0) bi++;
}

なお、この場合の i,j はすべてのセルインデックスを回りますから、中にはセルが存在しないインデックスも現れます。i,j を処理に使用する際は注意が必要です。

戻る

セル間の距離

離散ボロノイ図を距離を比較して作成する場合、インデックス i0,j0 のセルに対して、i1,j1 および i2,j2 のどちらがより近いかを判定する必要があります。離散ボロノイ図の定義に従えば、2つ以上の母点から等距離にあるセルは、母点の順序番号の若いものに所属しなければなりません。

セルの代表点間の距離を実数演算をつかって求めた場合、完全に等距離にあるが位置関係は異なるような2つのセル間の距離は、数値誤差によってどちらが大きくなるかはっきりしません。したがって、数値誤差の影響のない整数演算によって距離を判定する必要があります。

正方格子の場合、セルインデックスと平面座標が一致するため、簡単です。

この場合、セル間の距離を求めるコードと、判定用のコードは次のようになります。

/*
 * セル間の距離
 */
double dist(int i1,int j1,int i2,int j2){
  return sqrt(dist2(i1,j1,i2,j2));
}
int dist2(int i1,int j1,int i2,int j2){
  int dx=i2-i1;
  int dy=j2-j1;
  return dx*dx+dy*dy
}

/*
 * 3つの格子の距離の比較
 * i1,j1 の方が i2,j2 よりも i0,j0 に近いかどうか
 */
int nearer(int i0,int j0,int i1,int j1,int i2,int j2){
  return dist2(i0,j0,i1,j1)<dist2(i0,j0,i2,j2);
}

六角格子の場合、少し厄介です。下の図のように、j 軸が傾いているため、i 方向の差は単にインデックス i の差とはなりません。

セル間の距離の計算は、セルの面積を 1 に規格化しないほうが簡単になります。また、セル間の距離を比較することが目的の場合、単位は特に関係しません。そこで、六角形の1辺の長さを 2 として、

このようなサイズを仮定すると計算が楽です。この単位系で計算した距離は、必要に応じて HEXUX,HEXUY を用いて面積 1 に規格化された単位系に変換できます。

上の図の長さを用いると、i1,j1 と i2,j2 との平面座標での差 dx,dy は次のようになります。

dx = ( ( i2 - j1 ) * 2 + j2 - j1 ) * sqrt(3)
dy = ( j2 - j1 ) * 3

これより、i1,j1 と i2,j2 との間の距離 d の2乗は、

d2 = dx2 + dy2 = 3 * (( i2 - i1 ) * 2 + j2 - j1 )2 + 9 * ( j2 - j1 )2

で計算でき、整数演算により2点の距離の比較を行うことが出来ます。

以上より、六角格子の場合のセル間の距離と、どちらが近いかの判定コードは次のようになります。

/*
 * 規格化されないセル間の距離、距離の2乗
 */
double dist(int i1,int j1,int i2,int j2){
  return sqrt(dist2(i1,j1,i2,j2));
}
int dist2(int i1,int j1,int i2,int j2){
  int dx=(i2-i1)*2+j2-j1;
  int dy=(j2-j1)*3;
  return 3*dx*dx+9*dy*dy;
}

/*
 * 3つの格子の距離の比較
 * i0,j0 に i1,j1 の方が i2,j2 よりも近いかどうか
 */
int nearer(int i0,int j0,int i1,int j1,int i2,int j2){
  return dist2(i0,j0,i1,j1)<dist2(i0,j0,i2,j2);
}

三角格子の場合は、下の図のようになります。

三角格子の場合も、セルの面積を 1 に規格化せず、

上の長さを使って距離を計算します。セル i1,j1 と i2,j2 との間の平面座標での差 dx,dy は、

dx = ( ( i2 - i1 ) * 2 - j2 + j1 ) * sqrt(3)
dy = ( j2 - j1 ) * 3

となります。これによりセル間の距離 d の2乗は、

d2 = dx2 + dy2 = 3 * ( ( i2 - i1 ) * 2 - j2 + j1 )2 + 9 * ( j2 - j1 )2

で計算でき、整数演算により2点の距離の比較を行うことが出来ます。

以上より、三角格子の場合のセル間の距離と、どちらが近いかの判定コードは次のようになります。

/*
 * セル間の距離、距離の2乗、距離の4乗
 * 距離の4乗は、実際には 27*dist^4 であり、正確な値。
 */
double dist(int i1,int j1,int i2,int j2){
  return sqrt(dist2(i1,j1,i2,j2));
}
int dist2(int i1,int j1,int i2,int j2){
  int di=(i2-i1)*2-j2+j1;
  int dj=j2-j1;
  return 3*di*di+9*dj*dj;
}

/*
 * 3つの格子の距離の比較
 * i0,j0 に i1,j1 の方が i2,j2 よりも近いかどうか
 * ただし、3つのセルは存在していること
 */
int nearer(int i0,int j0,int i1,int j1,int i2,int j2){
  return dist2(i0,j0,i1,j1)<dist2(i0,j0,i2,j2);
}

戻る

離散ボロノイ図表示スクリプト drawDVD.ps

ダウンロード： drawDVD.ps [7.0KB] （← 別途 DVD データファイルが必要です）

Adobe Postscript で記述した、正方格子、六角格子、三角格子での離散ボロノイ図をページ出力するためのスクリプトです。離散ボロノイ図のデータファイルを読み込み、図として出力します。以下の特徴があります。

読み込み可能な離散ボロノイ図データのフォーマットは次のとおりです。なお、データファイル名は "drawDVD.ps" の DVD-FILE-NAME を直接いじってください。

例 説明 
(HEX)

0 0
10 10
3
true
true
true
red_i red_j
1 0 0
0 1 0
0 0 1
0 0 0 0 1 1 1 ...	
0 0 0 1 1 1 1 ...
...
0 2 2 2 2 2 1 ...
2 2
4 6
3 3
 
離散ボロノイ図の種類
  (SQR) (HEX) (TRI) のいずれか
図の左下のインデックス
図の大きさ
母点数
格子線を表示するか
母点を表示するか
赤点を表示するか
赤点のインデックス（表示する場合のみ）
色0 (R G B)
色1
色2
左下から横優先で並べた
離散点の母点番号データ
...
...
母点位置0
母点位置1
母点位置2

最初は (SQR) , (HEX) , (TRI) のいずれかでなければなりません。これはそれぞれ、正方格子、六角格子、三角格子の離散ボロノイ図を意味します。カッコも必要です。

次の2つの値は、表示する離散ボロノイ図の左下の格子のインデックスです。格子のインデックスについては、下の各格子座標系の説明を読んでください。

次の2つの値は、表示する離散ボロノイ図の大きさです。これもインデックス座標で指定します。

次の値は母点数です。後に続く色と母点位置は、この母点数と同じだけある必要があります。

次の2つは true か false を書きます。これらはそれぞれ、格子のグリッド線を表示するかどうか、母点の位置に黒丸を表示するかどうかを指定します。

その次は赤点を表示するかどうかを true か false で指定します。赤点とは指定したインデックスの位置に表示する赤い丸印で、特定のインデックスがどこであるかを正確に知るために使用します。これを true にした場合、続いてインデックス red_i red_j も指定しなければなりません。

その次に、各母点の領域を塗りつぶす色を RGB で指定します。各成分は [0,1] の範囲の実数値で指定し、0 が最も暗く、1 が最も明るくなります。1 0 0 は赤、1 1 1 は白を表します。

続いて、各格子の母点番号が並びます。この母点番号は、離散ボロノイ図の左下のインデックスから、i 方向優先でたどった順番に並んでいなければなりません。各格子のインデックスの順番については、下の各格子座標系の説明を読んでください。

最後に母点の位置を母点数分だけ並べます。母点の位置はその母点の格子のインデックスで指定します。このデータは母点を表示しない場合は不要です。

各データの間は1つ以上の空白文字（空白、タブ、改行）で区切ります。改行はデータの区切りならばどこにあっても構いません。

なお、離散ボロノイ図の範囲外の領域は色なしで表示されます。

戻る

離散ボロノイ図の定義

離散ボロノイ図とは、与えられた N 個の母点 p のセルインデックス kip,kjp ( p = 0,1, ... N-1 ) をもとにして、すべてのセルについてどの母点が最も近いかで色分けして出来る図です。このとき、近さを計る物差しとしてセル間の距離を使用します。2つ以上の母点から等距離にあるセルは、番号 p が小さい母点に所属すると決めます。

母点 p の離散ボロノイ領域 Vp は、その母点に所属しているセル i,j の集合として次のように定義されます。

この式の中で、d(i1,j1,i2,j2) はセル i1,j1 と i2,j2 との距離を意味します。

実際に離散ボロノイ図のデータを計算機を使って作成する場合は、母点ごとに色分けする代わりに各セルに母点の番号 p を格納していきます。番号は整数値ですから、プログラムでのセルは整数を格納できるような配列として表現されます。セル配列の項でも述べましたが、通常は次のように宣言します。

unsigned int *cell=new unsigned int[cell_size];

セル i,j を母点 p に所属させるには、

setCell(c2a(i,j),p);

とします。また、セル i,j が所属している母点の番号を得るには、

getCell(c2a(i,j));

とします。

戻る

全探索法による作成

全探索法とは、離散ボロノイ図のすべてのセルについて、すべての母点との距離を計算し所属を決定します。最も定義に忠実で理解しやすいですが、母点数やセル数が増えると作成に要する時間が急速に増大します。

ここでは詳しい作成法や作成コードは省略します。

戻る

波面法による作成

すべての母点から、同時に同心円状に領域を広げていくことで離散ボロノイ図を作成する方法が波面法です。セル数が一定ならば、母点数によらずほぼ一定の時間で離散ボロノイ図を作成できます。（ただし、1母点あたりのセル数が 100 を超えたあたりから一定時間ではなくなります）。

波面法により離散ボロノイ図作成するには、順序表という表が必要です。この表はセルの形状に応じて準備する必要がありますが、いったん作っておくと、離散ボロノイ図を作成するのにいちいち作り直す必要はありません。順序表の作成方法については順序表の作成を見てください。

格子ごとの順序表が準備できると、波面法アルゴリズムの大きな流れは格子によらず同一になります。以下にその流れを示します。

  1. 母点を母点番号の小さい順に並べた母点リストを準備する。
  2. 順序表を準備する。
  3. すべての母点のカウンタを0にする。
  4. すべてのセルを未所属状態にする。
  5. 順序表の先頭の要素 e を読む。
  6. 母点リストの先頭の母点 kp を読む。
  7. e1 を e にする。
  8. e1 の相対位置のセルおよびその対称位置のセルが未所属であれば、そのセルを母点 kp に所属させる。
  9. 1つでもセルを所属させることが出来れば、母点 kp のカウンタを0にする。そうでなければ、カウンタを1つ増やす。
  10. e1 の次の要素が等距離であれば、e1 を１つ進めて 8. へ。
  11. kp のカウンタが e1 の停止判定数に達したら、kp を母点リストから外す。
  12. リストに kp の次の母点があれば、kp を1つ進めて 7. へ。
  13. 母点リストに母点があるか、または e の次の順序表要素があれば、e を1つ進めて 6. へ。
  14. 終了。

このアルゴリズムの中で使われている対称位置だけが、セルの形によって異なります。

正方格子の場合、対称位置は一般に8箇所あります。

順序表の中に用意されているセルは、第1象限の下半分（黄色部分）にある青のセルだけですから、これを対称移動して全方向の同じ位置関係のセルを所属させます。上の図は一般の場合ですが、i,j が特殊な値のときは、この8箇所よりも少なくなります。

母点の位置が i,j の時の、対称位置は以下のようになります。

条件 対称位置
di = dj = 0 ( i , j )
di = dj != 0
( i+di , j+di ) ( i-di , j-di )
( i-di , j+di ) ( i+di , j-di )
dj = 0
( i+di , j     ) ( i-di , j     )
( i     , j+dj ) ( i     , j-di )
その他
( i+di , j+dj ) ( i+di , j-dj )
( i+dj , j+di ) ( i+dj , j-di )
( i-dj , j+di ) ( i-dj , j-di )
( i-di , j+dj ) ( i-di , j-dj )

六角格子の場合、対称位置は一般に12個所あり、次の図のようになります。

順序表に含まれているのは黄色の部分のに含まれる青いセル di,dj のみであり、これを対称移動させて残りの11個のセルの位置を算出します。di,dj の値によっては、対称位置は12箇所よりも少なくなります。

条件 対称位置
di = dj = 0 ( i , j )
di = dj != 0
( i + di , j + dj ) ( i - di , j + di + dj ) ( i - di - dj , j + dj )
( i - di , j - dj ) ( i + di , j - di - dj ) ( i + di + dj , j - dj )
dj = 0
( i + di , j ) ( i , j + di ) ( i - di , j + di )
( i - di , j ) ( i , j - di ) ( i + di , j - di )
その他
( i + di , j + dj ) ( i + dj , j + di ) ( i - dj , j + di + dj )
( i - di , j + di + dj ) ( i - di - dj , j + di ) ( i - di - dj , j + dj )
( i - di , j - dj ) ( i - dj , j - di ) ( i + dj , j - di - dj )
( i + di , j - di - dj ) ( i + di + dj , j - di ) ( i + di + dj , j - dj )

三角格子の場合、対称点は一般に6箇所あり、次の図のようになります。

この場合も、di,dj の条件によって次のように6箇所よりも少ない対称位置となります。

条件 対称位置
di = dj = 0 ( i , j )
di = dj != 0
( i + di , j + dj ) ( i - di , j ) ( i , j - dj )
dj = 0
( i + di , j ) ( i , j + di ) ( i - di , j - di )
その他
( i + di , j + dj ) ( i + dj , j + di ) ( i - dj , j + di - dj )
( i - di , j - di + dj ) ( i - di + dj , j - di ) ( i + di - dj , j - dj )

ただし、この対称位置は母点のセルが左向きの場合です。母点が右向きのセルだった場合は、対称位置が左右反対になります。三角格子のセルインデックス座標系において左右を反転することは、j = 2 i という軸に対して反転させることになります。（下図）

この図から、左向きのセルを中心に考える場合の相対位置を di,dj とすると、右向きのセルが中心に来たときの同じ位置関係のインデックス di ',dj ' は、

di ' = - di + dj
dj ' = dj

となります。なぜなら、まず直線 j = 2 i は i 軸に垂直ですから、線対称移動しても j の値は変化しません。次に、j = j0 のときの j = 2 i 上の i の値 i0 は、

i0 = j0 / 2

ですから、これに対して i の値を対称移動するためには、

  - ( i - i0 ) + i0
= - i + 2 * i0
= - i + j0

を計算すればよいことになるからです。

したがって、三角格子の離散ボロノイ図を波面法で作成する場合、母点が左向きならば di,dj をそのまま用い、右向きならば di ' = - di + dj の変換を行い、di ',dj 使えばよいことになります。

戻る

順序表の作成

正方格子、六角格子、三角格子上の離散ボロノイ図を波面法を用いて作成する場合、それぞれの格子状で作成された順序表が必要になります。ここでは、各格子上の順序表について説明します。

順序表とは、母点を中心にして波面を広げていく際に参照する、原点からの距離が小さい順にセルの相対位置を並べた表です。この表には、着目するセルの相対位置 di,dj のほかに、波面の拡大を停止させるための停止判定数 stop も格納されています。正方格子の順序表の最初の一部分を以下に示します。

k 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 …
dik 0 1 1 2 2 2 3 3 3 4 4 3 4 5 4 5 5 4 5 6 6 …
djk 0 0 1 0 1 2 0 1 2 0 1 3 2 0 3 1 2 4 3 0 1 …
stopk -1 -1 -1 -1 -1 -1 -1 7 8 9 9 10 10 11 12 12 12 13 13 14 14 …

実際に表に格納されるは下の3行です。

波面法により母点の領域を拡大するとき、実際には360度全方向に領域を拡大していく必要がありますが、表の中にはそのうちの一部分の座標だけが含まれています。これは、その他の領域の座標は、表の中の座標を線対称・点対称移動させることによって簡単に得ることが出来るからです。

停止判定数とは、母点の領域の拡大を止めるための条件を表しています。波面法によって離散ボロノイ図を作成する場合、まわりがすべて他の母点の領域になってしまった母点は、それ以上領域を拡大することがありませんから、速やかに処理を中止して他の母点の処理に専念する必要があります。このためには、自分の領域のまわりがすべて他の母点の領域になったかどうかを判定する必要があります。

正方格子の場合の順序表の作成の様子を以下に示します。

作成に際して与えられる条件は、セルが存在する範囲の領域幅 width と停止判定に用いる円の円周幅 Ra です。width は母点の領域が広がる可能性のある最大幅よりも大きい値である必要があります。例えば 400x600 の縦長な領域内に離散ボロノイ図を作成する場合、width は 600 以上の値でなければなりません。また Ra は飛び地の発生を考慮して、必要な精度が保たれるような値にする必要があります。最低でも sqrt(2) 以上でないと、正しい離散ボロノイ図が作成されません。正方格子の場合、1段飛び地に対処するためには 3 程度を用います。停止判定の詳細は飛び地についてを参照してください。

width が与えられると、順序表内のセル A,B,C のインデックスが決まり、順序表に必要な要素数 size が求まります。正方格子の場合の順序表のサイズは簡単で、セル B のインデックスを iB,0 とすると、最後の列に iB + 1 個のセルがあり、全部で iB + 1 列あるので、

size = ( 1 + iB + 1 ) * ( iB + 1 ) / 2 = ( iB + 2 ) ( iB + 1) / 2

です。

六角格子の場合、下の図のようになります。

このように、六角格子では width の値によって、順序表の右端のセルの位置が3通りのパターンを持ちます。B0,B1,B2 の違いは、

B0 : iB mod 3 = 0
B1 : iB mod 3 = 1
B2 : iB mod 3 = 2

です。波面法で離散ボロノイ図を作成する際に、領域を広げるためのセルはあるのに順序表の要素がなくなると、所属が決まらないセルが残ってしまいます。これを避けるにために、与えられた幅 width ぎりぎりまで領域が広がったとしても十分に順序表要素が足りるように準備しておかなければなりません。

上の図で、水平方向に B0 のセルまで必要な場合、赤で示したセルまでは順序表に含まれていなければなりません。緑、青のセルはそれぞれ B1,B2 の場合を示しています。これを見ると、B0 と B1,B2 では、C0,C1,C2 を過ぎてからのセルの選ばれ方に微妙な違いがあることが分かります。つまり、B0 では、C0 まで進んだ後、次の列はセル数が1つ減り、その次の列からは順に2つずつ減っていきます。ところが B1,B2 は C1,C2 まで進んだ後、次の列からセル数が2つずつ減っていきます。

順序表の要素数を多少多めに取っても構わないという場合、

iB' = floor( ( iB + 1 ) / 3 ) * 3 + 1

で新しい iB' を計算し、これを使うことで、十分な要素を含んだ順序表を1通りの場合（この場合は緑のケース）で作成可能です。

緑のケースでの順序表のサイズは、D1 のインデックスを iD,0 とすると、D1 までに、

1 + 2 + … + iD + ( iD + 1 ) = ( iD + 2 ) ( iD + 1) / 2

個のセルがあり、さらに D1 の次から B1 までの iB - iD 列では、B1 の側から数えると、

2 + 4 + 6 + …
= 2 * ( 1 + 2 + 3 + … )
= 2 * ( iB - iD ) ( 1 + iB - iD ) / 2
= ( iB - iD ) ( 1 + iB - iD )

個のセルがあります。よって、これらをあわせると、

size = ( iD + 2 ) ( iD + 1) / 2 + ( iB - iD ) ( 1 + iB - iD )

となります。ここで、緑のケースでは、

iD = floor( 2 * iB / 3 ) + 1

です。

三角格子の場合、以下のようになります。

三角格子の場合、上の図のような縦長の三角形領域になりますが、幅 width の範囲のセルしか必要でないので、上部の灰色の部分は順序表に入れる必要はありません。また、三角格子のセルインデックスのうち、赤で示したインデックスはセルが存在しませんので、これも順序表に入れる必要はありません。

六角格子と同じく、三角格子の場合でも幅 width の位置を含むセルが B0,B1 の2通りあり、インデックスを端まで進んだ後の i,j の変え方に少しの違いがあります。しかしこれも、順序表の様子数が多少増えることが構わない場合は、どちらかに統一することで一通りのコードで順序表を作成できます。

三角格子における順序表の大きさを求めるのはちょっと手間がかかるので、プログラムで数えるのが楽です。次のコードは、インデックス 0,0 から始めて、順序表のすべてのセルを順番にたどっていき、セルの個数をかぞえます。

  // 順序表のサイズをかぞえる
  // width : 幅
  // サイズ→ size
  int iB,jB,iC,jC;
  getIJ(width,0,&iB,&jB);
  getIJ(width,width,&iC,&jC);
  if(iB%2==0) iB++;
  if(getCellType(iC,++jC)==2) iC++;
  int size=0;
  for(int i=0;i<=iC;i++){
    if(i<=iB){
      for(int j=0;j<=i;j++){
        if(getCellType(i,j)!=2) size++;
      }
    }
    else if(i<=jC){
      for(int j=(i-iB)*2;j<=i;j++){
        if(getCellType(i,j)!=2) size++;
      }
    }
    else{
      for(int j=(i-iB)*2;j<=jC;j++){
        if(getCellType(i,j)!=2) size++;
      }
    }
  }

三角格子のセルには左向きと右向きとがあり、上の順序表は左向きのセルの母点を中心にして作成されています。したがって、このまま、母点が右向きのセルであるような場合に適用すると、間違った図を作成してしまいます。

幸いにして、左向きの母点用の順序表は、わずかな計算で右向きの母点用に利用することが出来ますので、わざわざ右向き用の順序表を作成する必要はありません。この件については、波面法のアルゴリズムのところで詳しく説明します。

戻る

飛び地について

飛び地とは、同一母点に所属する領域が連続しないで存在している場合の、母点を含まない側の領域の総称です。離散化されていないボロノイ図では、1つのボロノイ領域は凸多角形であり、ボロノイ領域に飛び地は存在しません。しかし、離散ボロノイ図では離散化によるサンプリングの粗さによってボロノイ図に飛び地が発生する可能性があります。また、重みつきボロノイ図においては、離散化誤差の影響ではなく、離散化の有無にかかわらず本質的に飛び地が発生する可能性があります。

ここでは、重みつきボロノイ図については触れず、離散化誤差によって生じる飛び地について説明します。

平面を格子状に離散化し、母点とセルの代表点との距離で離散ボロノイ図を作成した場合、セルの代表点は隙間を持って並んでいるため、離散化されていないボロノイ領域がこの隙間の間に伸びて広がっている場合、飛び地が発生します。飛び地が発生するケースの例を以下に示します。

赤枠で示したセルが飛び地です。上の図の飛び地は、本当ならば地続きになっているべきものですが、セルの代表点が隙間を持つために不連続な領域になってしまいます。

上で示した飛び地は、どれか1つのセルが青い領域になれば連続した領域となることが出来ます。1つのセルによって地続きになることが出来る飛び地を、1段飛び地と呼びます。

一般に、少なくとも N 個のセルが自分の領域にならないと地続きの領域になれないような飛び地を N 段飛び地と呼ぶことにします。母点の配置によっては、2段、3段飛び地なども存在する可能性があります。離散平面の大きさが有限で、母点がセルの代表点にしかないような場合、発生する飛び地の段数は有限に抑えられると考えられますが、理論的に何段飛び地まで発生しうるかはまだ明らかになっていません。

飛び地自体は、離散ボロノイ図の定義にのっとった正しい領域です。が、問題となる可能性があるのは、「1つのボロノイ領域は連続した凸多角形領域である」という仮定に基づいた処理がなされる場合です。

例えば、波面法によって離散ボロノイ図を作成する場合、ある母点の領域が「完成する」とその母点からの波面の拡大は停止させます。ある母点の領域が完成したかどうかは、その母点を中心とする同心円の円周上のすべての点が他の母点の領域に属しているかどうかで行います。つまり、その母点の領域は、母点から連続であるという仮定に基づいているわけです。

ところが飛び地が存在する場合、いったんは円周上のすべての点が他の母点の所属になった後に、再度自分の領域が現れる場合があることになります。このようのことが起こると、正しい領域である飛び地が発生しないで消えてしまうことになるのです。

波面法ではこの問題を低減させるために、次のような対策をとっています（理論的には、飛び地はかなり先の段数まで発生する可能性があるため、完全に回避するとまでは、今のところいえません）。

波面法では、すべての母点から同時に同心円状に領域を広げていくことで離散ボロノイ図を作成します。このとき、領域拡大を停止させるために、同心円の円周上のすべてのセルが他の母点の所属になったかどうかを判定しています。ここで、円周にいくらかの幅を持たせることによって、特定の段数までの飛び地が消えてしまわないようにすることが出来ます。

具体的な方法を、格子ごとに述べます。

正方格子の場合、1段、2段、3段飛び地が発生する場合で、それぞれ最も遠い位置関係にあるケースは次の図のようになります。

したがって、1段飛び地が消えないようにするためには sqrt(5) 以上の停止判定円の円周幅があればよいことになります。これらの実際の値は、

R1 = pow(5,0.5) = 2.236068
R2 = pow(13,0.5) = 3.6055513
R3 = 5.0

です。

六角格子の場合、3段飛び地までの最も遠い配置は次の図のようになります。

これらの数字の実際の値は、

R1 = pow(6,0.5)*pow(3,-0.25) = 1.86121
R2 = pow(14,0.5)*pow(3,-0.25) = 2.843045
R3 = pow(26,0.5)*pow(3,-0.25) = 3.874417

です。

三角格子の場合、3段飛び地までの最も遠い配置は次の図のようになります。

これらの数字の実際の値は、

R1 = pow(52,0.5)*pow(3,-0.5)*pow(3,-0.25) = 3.16345
R2 = pow(84,0.5)*pow(3,-0.5)*pow(3,-0.25) = 4.02067
R3 = pow(172,0.5)*pow(3,-0.5)*pow(3,-0.25) = 5.753383

です。三角格子で3段飛び地まで対応するには、停止判定円の円周幅 R を 6 にすれば大丈夫ということです。

戻る

隣接情報の抽出

出来上がった離散ボロノイ図から隣接情報を抽出する場合、図のすべてのセルについてそれに隣接するセルを調べます。あるセルの所属母点番号が k1 で、そのセルの隣接セルの母点が k2 のとき、k1 != k2 ならば、母点 k1 と k2 とは隣接していることになります。

母点の隣接関係は行列を用いて表現することが出来ます。母点数が N の場合、次のような N x N 行列を準備します。

n(0,0) n(0,1) … n(0,N-1)
n(1,0) n(1,1) … n(1,N-1)
： ： 　 ：
n(N-1,0) n(N-1,1) … n(N-1,N-1)

n(k1,k2) は母点 k1 と母点 k2 とが隣接しているかどうかを表すフラグです。n(k1,k2) = 0 ならば母点 k1 と k2 は隣接していません。n(k1,k2) = 1 ならば2つの母点は隣接しています。

しかし、この表現には無駄があります。この行列は母点 k1 と k2 とが隣接していることを、2つの要素、n(k1,k2) と n(k2,k1) の両方を使って表しています。実際には片方だけあれば十分ですから、次のような表で足りることが分かります。

1 2 … k2 … N-2 N-1
 0 　 n(0,1) n(0,2) … n(0,k2) … n(0,N-2) n(0,N-1)
1 　 　 n(1,2) … n(1,k2) … n(1,N-2) n(1,N-1)
： 　 　 　 … ： … ： ：
k1 　 　 　 　 n(k1,k2) … n(k1,N-2) n(k1,N-1)
： 　 　 　 　 　 … ： ：
N-3 　 　 　 　 　 　 n(N-3,N-2) n(N-3,N-1)
N-2 　 　 　 　 　 　 　 n(N-2,N-1)
　 　 　 　 　 　 　 　

赤い部分は使わない部分です。したがって、この表の大きさは N * (N-1) / 2 です。また、n(k1,k2) の k1 と k2 が満たすべき条件として、k1 < k2 が必要です。したがって、例えば母点 7 と母点 3 が隣接しているかどうかは、7 と 3 を入れ替えて n(3,7) をチェックすることで判断します。

この表を作るのに要する時間は離散ボロノイ図の大きさに比例しますから、O(S) です。母点数には依存しません。

まずこの表にアクセスするためのメソッドを定義しておきます。

C 言語の配列は1次元配列ですから、2次元行列の要素にアクセスするためには、2次元の添え字 k1,k2 を1次元の配列インデックス a に変換する必要があります。しかし、上の表のように各行の要素数が異なっている場合、ちょっとだけめんどうです。

上の図のように、n(0,1) から始めて縦方向優先で左から右へ順番に配列を割り当てるとすると、赤で示した n(k1,k2) の位置は、

k2 * ( k2 - 1 ) / 2 + k1

番目になることが分かります（緑に部分に k2 * ( k2 - 1 ) / 2 個の要素があり、青の部分に k1 この要素があるので）。したがって、隣接情報にアクセスするコードは以下のようになります。

/*
 * 隣接情報の初期化・設定・取り出し
 * int N     : 母点数
 * char *nei : 大きさ N*(N-1)/2 の配列
 */
void resetNeighbors(){
  int s=N*(N-1)/2;
  for(int i=0;i<s;i++) nei[i]=0;
}
void setNeighbor(int k1,int k2){
  if(k1==k2 || k1<0 || k1>=N || k2<0 || k2>=N) return;
  if(k1>k2){ int d=k1; k1=k2; k2=d; }
  nei[k2*(k2-1)/2+k1]=1;
}
int getNeighbor(int k1,int k2){
  if(k1==k2 || k1<0 || k1>=N || k2<0 || k2>=N) return 0;
  if(k1>k2){ int d=k1; k1=k2; k2=d; }
  return nei[k2*(k2-1)/2+k1];
}

ここではグローバルに母点数 N と隣接情報配列 nei が宣言されているとしました。getNeighbor() は、同じ母点どうし、または範囲外の母点番号が与えられたときは、0 （隣接なし）を返します。

以上を踏まえて、各格子ごとに隣接情報を取り出すアルゴリズムを考えます。

正方格子の場合、上下左右の4箇所を調べる必要があります。上下左右のセルは、着目しているセルのインデックスを i,j とすると、

の4箇所ですから、隣接情報配列を作成するコードは次のようになります。

/*
 * 隣接情報配列の作成
 */
void makeNeighborInfo(){
  resetNeighbors();
  for(int j=base_j;j<base_j+height;j++){
    for(int i=base_i;i<base_i+width;i++){
      KPID id=getCell(c2a(i,j));
      KPID r=getCell(c2a(i+1,j));
      KPID l=getCell(c2a(i-1,j));
      KPID u=getCell(c2a(i,j-1));
      KPID d=getCell(c2a(i,j+1));
      setNeighbor(id,r);        // id と r は隣接している
      setNeighbor(id,l);        // id と l は隣接している
      setNeighbor(id,u);        // id と u は隣接している
      setNeighbor(id,d);        // id と d は隣接している
    }
  }
}

六角格子の場合、隣接しているセルは6箇所あります。着目しているセルのインデックスを i,j とすると、

の6箇所ですから、隣接情報配列を作成するコードは次のようになります。

/*
 * 隣接情報配列の作成
 */
void makeNeighborInfo(){
  resetNeighbors();
  int bi=base_i,c=0;
  for(int j=base_j;j<base_j+height;j++){
    for(int i=bi;i<bi+width;i++){
      KPID id=getCell(c2a(i,j));
      KPID n1=getCell(c2a(i+1,j));
      KPID n2=getCell(c2a(i+1,j-1));
      KPID n3=getCell(c2a(i,j-1));
      KPID n4=getCell(c2a(i-1,j));
      KPID n5=getCell(c2a(i-1,j+1));
      KPID n6=getCell(c2a(i,j+1));
      setNeighbor(id,n1);        // id と n1 は隣接している
      setNeighbor(id,n2);        // id と n2 は隣接している
      setNeighbor(id,n3);        // id と n3 は隣接している
      setNeighbor(id,n4);        // id と n4 は隣接している
      setNeighbor(id,n5);        // id と n5 は隣接している
      setNeighbor(id,n6);        // id と n6 は隣接している
    }
    if((++c)%2==0) bi--;
  }
}

六角格子の場合は、着目しているセルの向きによって2通りの位置関係があります。

したがって、隣接情報配列を作成するコードは次のようになります。

/*
 * 隣接情報配列の作成
 */
void makeNeighborInfo(){
  resetNeighbors();
  int bi=base_i,c=base_right;
  for(int j=base_j;j<base_j+height;j++){
    for(int i=bi;i<bi+width;i++){
      KPID id=getCell(c2a(i,j));
      KPID n1,n2,n3;
      switch(getCellType(i,j)){
      case 0:
        n1=getCell(c2a(i+1,j));
        n2=getCell(c2a(i-1,j-1));
        n3=getCell(c2a(i,j+1));
        break;
      case 1:
        n1=getCell(c2a(i-1,j));
        n2=getCell(c2a(i+1,j+1));
        n3=getCell(c2a(i,j+1));
        break;
      case 2:
        continue;
      }
      setNeighbor(id,n1);        // id と n1 は隣接している
      setNeighbor(id,n2);        // id と n2 は隣接している
      setNeighbor(id,n3);        // id と n3 は隣接している
    }
    if((++c)%2==0) bi++;
  }
}

セルの向きは getCellType() メソッドで得られますので、これが 0 ならば上図の左の形、1 ならば右の形となります。2 ならばそのセルは存在しませんので、何もしないで次の反復に移ります。

戻る

更新履歴

2001.1.16 公開
2001.1.17 三角格子のインデックス変換 a2c(int a,int *i,int *j) の間違いを修正
2001.1.17 drawDVD.ps へリンクを張った
2001.1.25 六角格子で、母点の存在範囲の上下の補正が抜けていたのを修正
2001.1.25 セルの矩形領域内判定を追加
2001.1.25 三角格子の場合の、幅 width の修正に base_right の考慮が抜けていたのを追加
2001.1.25 drawDVD.ps のバグを修正
2001.1.25 セル間の距離の比較を追加
2001.1.26 三角格子のセル間の距離比較が間違っていたのを訂正
2001.1.27 隣接情報の抽出を追加
2001.1.31 drawDVD.ps に赤点表示機能を追加
2001.1.31 順序表の作成を追加
2001.2.01 波面法による作成を追加
2001.2.05 三角格子のセル間距離の定義が間違っていたのを修正

戻る