比較，kd-treeとoctree

Last updated at 2025-03-13Posted at 2025-03-12

はじめに

点間の距離計算が頻繁に求められる事例においては，それを効率良くできるように点群の空間配置を表現することが重要です．常道として用いられるのは $k$d木（kd-tree） と 八分木（octree） です．どちらもよく知られたものですが，結局どちらを使えば良いのか？は当然気になるところです．それに答えるために，Zeoでの実装を紹介しつつ，最近傍点探索と$r$近傍点群探索の二つの評価結果を比較してみようというのが本記事の趣旨です．

$k$d木＝任意の$k$次元の木ですが，実質的にその構造が意味を持つのはデカルト空間においてくらいだし，$k$があまり大きくなると計算効率も落ちるので，実用上は$k=3$と思ってほとんど差し支えありません．また，octreeは元来デカルト空間に特化した表現です．したがって，本記事で言う点群とはデカルト空間における点群のことですので，ご注意下さい．

Zeoのkd-tree実装

データ型の定義と点の追加処理

kd-treeのデータ型は，次のように定義しています．

typedef struct{
  int id;          /*!< identifier of a tree node */
  zAxis split;     /*!< split axis index */
  zVec3D point;    /*!< spliting vertex */
  zAABox3D region; /*!< region to cover */
} zVec3DTreeData;

typedef struct __zVec3DTree{
  struct __zVec3DTree *parent;
  struct __zVec3DTree *child[2];
  uint size;
  zVec3DTreeData data;
} zVec3DTree;

zVec3DTreeDataのメンバ変数splitが，領域の分離軸IDになります．
zAxisは列挙型で，0が$x$軸，1が$y$軸，2が$z$軸にそれぞれ対応します．
pointには3次元の点の座標が入り，この点を通るように分離軸が定められます．
regionは当該ノードがカバーする直方体領域を規定する構造体で，その型zAABox3D（Axis-Aligned Box）は次のように定義されます．

typedef struct{
  zVec3D min; /*!< minimum coordinates */
  zVec3D max; /*!< maximum coordinates */
} zAABox3D;

各ノードが分離軸の座標値さえ持っていれば，根元側の値を伝播することでカバー領域を知ることは出来るので，全ての値を持たせるのは冗長ではあるのですが，領域判定を単純にするために導入しました．
最初のメンバ変数idはおまけで，後述するように同一の木の中で何番目に生成されたノードかを保存するようにしています．kd-treeと配列の相互変換を実装するための措置ですが，無くても支障はありません．

zVec3DTreeは，実際には上記のコードを直接書いているわけではなく，ZEDAのzTreeClass()マクロで生成しています．これはノードのデータ型と同時に，初期化関数zVec3DTreeInit()，破棄関数zVec3DTreeDestroy()，ノード生成関数zVec3DTreeNodeAlloc()まで自動生成するものです．それらの関数の中身については割愛し，肝心の木に点を追加する処理を次に示します．

static zVec3DTree *_zVec3DTreeCreateLeaf(zAxis split, const zVec3D *point, int id)
{
  zVec3DTree *leaf;

  if( !( leaf = zAlloc( zVec3DTree, 1 ) ) ){
    ZALLOCERROR();
    return NULL;
  }
  leaf->size = 1;
  leaf->data.id = id;
  leaf->data.split = split;
  zVec3DCopy( point, &leaf->data.point );
  leaf->child[0] = leaf->child[1] = NULL;
  return leaf;
}

static int _zVec3DTreeChooseBranch(const zVec3DTree *node, const zVec3D *point)
{
  return point->e[(int)node->data.split] >= node->data.point.e[(int)node->data.split] ? 0 : 1;
}

static zVec3DTree *_zVec3DTreeAddPoint(zVec3DTree *node, const zVec3D *point, int id)
{
  int b;
  zVec3DTree *leaf;

  node->size++;
  if( node->child[( b = _zVec3DTreeChooseBranch( node, point ) )] )
    return _zVec3DTreeAddPoint( node->child[b], point, id );
  if( !( leaf = _zVec3DTreeCreateLeaf( ( node->data.split + 1 ) % 3, point, id ) ) )
    return NULL;
  node->child[b] = leaf;
  zAABox3DCopy( &node->data.region, &leaf->data.region );
  if( b == 0 )
    leaf->data.region.min.e[(int)node->data.split] = node->data.point.e[(int)node->data.split];
  else /* b == 1 */
    leaf->data.region.max.e[(int)node->data.split] = node->data.point.e[(int)node->data.split];
  return leaf;
}

zVec3DTree *zVec3DTreeAddPoint(zVec3DTree *tree, const zVec3D *point)
{
  if( tree->data.split == zAxisInvalid ){
    tree->size = 1;
    tree->data.id = 0;
    tree->data.split = zX;
    zVec3DCopy( point, &tree->data.point );
    return tree;
  }
  return _zVec3DTreeAddPoint( tree, point, tree->size );
}

処理の本体は，再帰的に領域を絞り込んで与えられた点を登録する関数_zVec3DTreeAddPoint()です．特に難しいことはやっていませんが，葉ノードを新たに生成する際に，親ノードのsplitの値を1増やして3で割った余りをsplitの値としているところがポイントでしょうか．こうすることで，根から葉に向かって分離軸が$x$→$y$→$z$→$x$→$\cdots$ と順番に変わっていきます．

※分割された直方体領域の最も長い辺に沿った軸で分割するという流儀もありますが，$-\infty$ / $\infty$を領域境界にとることを許容しているため，処理がややこしくなるのを避けて最も簡単な規則を採用しました．

葉ノードの領域は，いったん親ノードの領域をコピーした後，分離軸に沿った最小値or最大値（どちらにするかは，点がどちらの領域に属するかで変わります）を，その点の分離軸成分に置き換えるようにしています．点が存在しない領域に対応するノードは生成されません．

ノードのメンバ変数sizeには，自分自身を含めた部分木のノード数が入りますので，idには，そのノードが木全体で何番目に作られたものであるかを示す番号が入ることになります．繰り返しになりますが，これらは無くても木構造の操作に支障はありません．

r近傍点探索

$r$近傍点探索は，与えられた点からの距離が$r$未満である点を集めた部分点群を作る処理です．

static bool _zVec3DTreeVicinityTest(const zVec3DTree *node, const zVec3D *p, double radius_sqr, zVec3DData *vicinity)
{
  if( zVec3DSqrDist( &node->data.point, p ) < radius_sqr ){
    if( !zVec3DDataAdd( vicinity, &node->data.point ) ) return false;
  }
  return true;
}

static zVec3DData *_zVec3DTreeVicinityOpp(zVec3DTree *node, const zVec3D *point, double radius_sqr, zVec3DData *vicinity)
{
  if( !_zVec3DTreeVicinityTest( node, point, radius_sqr, vicinity ) ) return NULL;
  if( node->child[0] && _zVec3DTreeIsOverlap( node->child[0], point, radius_sqr ) )
    if( !_zVec3DTreeVicinityOpp( node->child[0], point, radius_sqr, vicinity ) ) return NULL;
  if( node->child[1] && _zVec3DTreeIsOverlap( node->child[1], point, radius_sqr ) )
    if( !_zVec3DTreeVicinityOpp( node->child[1], point, radius_sqr, vicinity ) ) return NULL;
  return vicinity;
}

static zVec3DData *_zVec3DTreeVicinity(const zVec3DTree *tree, const zVec3D *point, double radius_sqr, zVec3DData *vicinity)
{
  int b;
  zVec3DTree *ob; /* opposite branch */

  if( tree->child[( b = _zVec3DTreeChooseBranch( tree, point ) )] )
    if( !_zVec3DTreeVicinity( tree->child[b], point, radius_sqr, vicinity ) ) return NULL;
  if( !_zVec3DTreeVicinityTest( tree, point, radius_sqr, vicinity ) ) return NULL;
  ob = tree->child[1-b];
  if( ob && _zVec3DTreeIsOverlap( ob, point, radius_sqr ) )
    if( !_zVec3DTreeVicinityOpp( ob, point, radius_sqr, vicinity ) ) return NULL;
  return vicinity;
}

zVec3DData *zVec3DTreeVicinity(const zVec3DTree *tree, const zVec3D *point, double radius, zVec3DData *vicinity)
{
  zVec3DDataInitAddrList( vicinity );
  return _zVec3DTreeVicinity( tree, point, _zSqr(radius), vicinity );
}

一番下のzVec3DTreeVicinity()が直接呼び出すべき関数で，vicinityをzVec3Dポインタのリストとして初期化してから，与えられた点pointからの距離がradius未満の点を全て探し出し，vicinityに登録します．

まずはpointが属する領域のノードを再帰的に探し，それぞれのノードに保存されている点とpointとの距離がradius以下であれば，vicinityに追加します．また，反対側の領域のノードについても，をの直方体領域とpointを中心とする半径radiusの球が重なっているならば確認します．最近傍点探索の時と考え方はほぼ同じであることが，お分かり頂けると思います．

Zeoのoctree実装

データ型の定義と点の追加処理

octreeのデータ型については，別記事「分解能可変octreeを作ったよ」で紹介しています．空間分割数を可変にする趣旨で作ったもので，いわゆるMorton数を使った高速領域探索はしていませんので，ご注意下さい．octreeのノードはオクタント(octant)，その子ノードはサブオクタント(suboctant)ともそれぞれ呼ばれます．

r近傍点探索

octreeを使った$r$近傍点探索処理を示します．

static zVec3DData *_zVec3DOctantVicinity(zVec3DOctant *octant, const zVec3D *point, double radius_sqr, zVec3DData *vicinity)
{
  int i;
  zVec3DListCell *cp;

  if( zAABox3DSqrDistFromPoint( &octant->region, point ) >= radius_sqr ) return vicinity;
  if( !zListIsEmpty( &octant->points ) ){
    zListForEach( &octant->points, cp )
      if( zVec3DSqrDist( &cp->data, point ) < radius_sqr )
        if( !zVec3DDataAdd( vicinity, &cp->data ) ) return NULL;
    return vicinity;
  }
  for( i=0; i<8; i++ ){
    if( !octant->suboctant[i] ) continue;
    if( !_zVec3DOctantVicinity( octant->suboctant[i], point, radius_sqr, vicinity ) ) return NULL;
  }
  return vicinity;
}

zVec3DData *zVec3DOctreeVicinity(zVec3DOctree *octree, const zVec3D *point, double radius, zVec3DData *vicinity)
{
  zVec3DDataInitAddrList( vicinity );
  return _zVec3DOctantVicinity( &octree->root, point, _zSqr(radius), vicinity );
}

こちらはよりシンプルで，

与えられた点pointから着目しているノード（オクタント）octant領域までの距離がradian以上ならば何もしない
octantが葉ノードならば，その中に含まれる全ての点のうちpointからの距離がradian未満のものを全てvicinityに登録する
葉ノードで無いならば，高々7個のsuboctantについて再帰的に探索する

という処理になっています．直方体領域と球が重なっているか判定する処理には，kd-treeと同じくzAABox3DSqrDistFromPoint()を使っています．

kd-treeとoctreeの性能比較

r近傍点探索

$r$近傍点探索についても調べました．点群生成条件と試行回数は先程と同じ，octreeの空間分解能も同じく$1\times1\times1$とし，$r$は$1$としました．
結果を示します．

念のため，これも全ての解が一致していることは確認しています．今度もoctreeの方がkd-treeよりも総じて遅いのですが，差は小さくなっており，kd-treeの方が時間がかかっている例もあります．

さらに，点群を直方体領域で一様に分布させるのではなく，中心座標$(5,5,5)$，半径$4$の球面上で一様に分布するように生成してみました．点の数は10000に減らしています．他の条件は全て同じです．このときの結果を示します．

この条件では，kd-treeの方が総じてoctreeよりも遅くなると分かります．何箇所か時間が突出しているのは，$r$近傍点が存在したケースです．つまり，ほとんどの場合では$r$近傍点が無かったということになります．このとき，計算時間の多くは近傍点リストに新たな点を登録するためのメモリ操作に費やしていることが分かりました．

考察というほどでもないもの

上記の結果から，次のことが言えそうです．

最近傍点探索はkd-treeの方が有利
octreeによる最近傍点探索の効率は，空間分解能の設定に大きく依存
一般的なケースでの$r$近傍点探索は，kd-treeとoctreeで大きな性能差は無い
$r$近傍点探索において，点群の分布に偏りがある場合はoctreeの方が有利な可能性がある

点群配置と計算効率の関係については，もっと条件を変えて調べないとなんとも言えないところがありますが，使用時の目安にはなるかと思います．

You get articles that match your needs
You can efficiently read back useful information
You can use dark theme

What you can do with signing up