領域探索 (kd 木)

説明

点の集合が与えられる．領域が与えられたとき，その領域に含まれる点をすべて報告せよ．という問題を領域探索問題という．これは点位置決定問題の双対問題である．幾何学的な問題だけでなく「年齢 20〜30 歳，年収 500〜1000 万の人をすべて探せ」などのデータベース問い合わせもこの問題に帰着する．

この問題を解くために使える構造が kd 木である．これは一次元における二分探索木の拡張であって，探索木の各段階で x における比較，y における比較を切り替えるものである．点が一様に散らばっている場合，二分探索木と同様に平均計算量 O(log n) となる．

以下では問い合わせ領域は直交領域，すなわち x 軸と y 軸に平行な辺をもつ長方形に限っている．しかしどんな形状の領域でも，bounding box をとるなどしてそれが半平面のどちらに位置するかを高速に判定できれば問い合わせ可能である．

この方法は容易に高次元に拡張できる．しかし次元が高くなるほど「高次元空間中の一様分布」は発生しづらくなるので，フラクショナルカスケーディングなどの技法が要求されることもある．

計算量

• 挿入: O(log n)
• 検索: O(log n + k)，ただし k は範囲内の点の総数．

ただしどちらも点が一様に散らばっている場合の期待計算量である．

使い方

void insert(const point &p)

kd 木に新しい点 p を挿入する．

void search(const point &ld, const point &ru, OUT out)

左下点 ld, 右上点 ru で張られる四角形内に入る点を out に出力する．

TODO

• どこかで verify．

ソースコード

struct kdtree {
  struct node {
    point p;
    node *l, *r;
    node(const point &p)
      : p(p), l(NULL), r(NULL) { }
  } *root;
  kdtree() : root(NULL) { }

#define compare(d, p, q) (d ? real(p) < real(q) : imag(p) < imag(q))
  void insert(const point &p) {
    root = insert(root, 0, p);
  }
  node *insert(node *t, int d, const point &p) {
    if (t == NULL) return new node(p);
    if (compare(d,p,t->p)) t->l = insert(t->l, !d, p);
    else                   t->r = insert(t->r, !d, p);
    return t;
  }
  template <class OUT>
  void search(const point &ld, const point &ru, OUT out) {
    search(root, 0, ld, ru, out);
  }
  template <class OUT>
  void search(node *t, int d, const point &ld, const point &ru, OUT out) {
    if (t == NULL) return;
    const point &p = t->p;
    if (real(ld) <= real(p) && real(p) <= real(ru) &&
        imag(ld) <= imag(p) && imag(p) <= imag(ru)) *out++ = p;
    if (!compare(d,p,ld)) search(t->l, !d, ld, ru, out);
    if (!compare(d,ru,p)) search(t->r, !d, ld, ru, out);
  }
};

Verified