計算幾何学とは、幾何学的問題を解くアルゴリズムを研究する分野です。扱うのは点、線、多角形といった幾何学的図形にまつわる問題であり、多岐に渡る応用の基礎となっています。

数値や文字を扱うアルゴリズムは、コンピュータ上で自然に表現して処理することができます。一方で幾何学的な問題となると状況は異なり、点や線といった基本的な要素を扱う問題ですら、複雑な処理が必要になります。

幾何学的な問題は、人間が一目見れば一瞬で解けるような問題が多いというのもひとつの特徴です。例えば2つの線分が交差しているかどうかは、人間が一目見れば一瞬で分かるはずです。しかしこれをコンピュータを用いて判定しようとすると、それなりに複雑な計算が必要とされます。

幾何学的な問題の難しさは、考慮しなければならない例外の多さに由来します。例えば上で例にあげた線分の交差であれば、

  • 2 つの線分が点を共有している場合
  • 2 つの線分が平行な場合
  • 2 つの線分が完全に重なっている場合

といったところがあげられるでしょうか。複雑な問題になれば、より多くの例外を考慮しなければならないでしょう。

図形の定義

今回扱う図形は平面上の幾何学的図形とします。今回は点、線分、多角形を以下のように定義します:

// 点
struct point {
  int x, y;
};

// 線分
struct line_segment {
  point p1, p2;
};

// 多角形
typedef point polygon[N_MAX];

点は X, Y 座標の組とします。座標は整数としています。実数とする場合もありますが、大半の場合は整数で十分ですし、速度的な優位性も大きく、また計算がいくらか単純になります。

線分は 2 点の組とします。2 点を結ぶ線分として表現します。

多角形は点の配列とします。配列の先頭から順に点をめぐる図形とします。多角形の表現にはいくつかバリエーションが考えられますが、今回はこのような点の配列とします。

多角形は点や線分と比べるとバリエーションに富む図形です。形によっては、専用のデータ構造を設けたほうが、その図形の特徴を利用することで、よりよいアルゴリズムを導くことができる場合もあります。長方形や正方形はその代表的な例と言えるでしょう。

線分の交差

まず、2 つの線分が交差するかどうかを判定する問題を扱います。以下のようなパターンが考えられます:

fig1.png

冒頭で述べた通り、人間が目で見れば、答えは明らかです。

1 番目は文句なしに交差しています。

2 番目は一方の線分上に、もう一方の線分の端点が存在しています。線分の一部として端点を含むものとして、この場合は交差していると判定します。

3 番目と 4 番目は交差していません。3 番目は、一方の線分をそのまま引き伸ばすと、もう一方の線分と交差します。が、線分としてはやはり交差しません。4 番目では 2 つの線分は平行であり、線分をどこまで引き伸ばしても交差しません。

第一の解(課題あり)

さて、このような 2 つの線分が交差しているかどうかの判定方法を考えてみます。交差判定の関数のシグネチャは次のようなものとします:

bool intersect(line_segment s1, line_segment s2);

2 つの線分 s1, s2 が交差していれば TRUE を、交差していなければ FALSE を返します。

直接的な判定方法としては、各線分を含む直線同士の交点を求め、その交点が各線分上に存在するかどうか判定する方法があります。平面上の直線は ax + by + c = 0 の式で表せるので、直線は次のようなデータ構造で表現できます:

// 直線
struct line {
  int a, b, c;
};

線分から直線に変換するには、線分の端点の両方を含む直線を求めればよく、先の式に p1(x1, y1), p2(x2, y2) を代入すると、次のような連立方程式が得られます:

a*x1 + b*y1 + c = 0
a*x2 + b*y2 + c = 0

次の連立方程式を解けば、

a = y2 - y1
b = x1 - x2
c = y1*(x2 - x1) - x1*(y2 - y1)

が得られます。これをプログラムに落とせば、

line to_line(line_segment s) {
  int dx = s.p2.x - s.p1.x;
  int dy = s.p2.y - s.p1.y;
  line l = { dy, -dx, y1*dx - x1*dy };
  return l;
}

と書けます。

次に、2 直線の交点を求める方法を考えます。2 直線の交点を x, y とすると、

a1*x + b1*y + c1 = 0
a2*x + b2*y + c2 = 0

という連立方程式が成立します。これを解くと、

x = (b1*c2 - b2*c1) / (a1*b2 - a2*b1)
y = (a2*c1 - a1*c2) / (a1*b2 - a2*b1)

が得られます。どちらの分母は共通です。この分母が 0 のとき、2 直線は平行であり、交点は存在しません。また、2 直線が完全に重なる場合も分母が 0 になります。今回はこれも平行と同じように扱うことにします。これを踏まえてプログラムに落とせば、

bool intersection_point(line l1, line l2, point* p) {
  int d = (l1.a * l2.b) - (l2.a * l1.b);
  if (d == 0)
    return FALSE;

  int x = (l1.b * l2.c) / d;
  int y = (l2.a * l1.c) / d;
  p->x = x;
  p->y = y;
  return TRUE;
}

といった具合に書けます。

次に任意の点が線分上に存在するかどうかの判定を考えます。線分の端点を A, B とおき、線分上に存在するかどうか判定する点を C とすると、ベクトル AB (以下 a) と ベクトル AC (以下 b) が完全に重なり、かつ a が b 以上の大きさであることを調べればよさそうです。

まず、a と b が重なるかどうかは、内積の定義より求められます。

a・b = |a||b|cosθ

a, b が完全に重なるとき、2 つのベクトルがなす角θは 0 なので、 cosθ = 1 となります。よって、

a・b = |a||b|

かどうかを判定すれば、ふたつのベクトル a, b が重なるかどうかを判定できます。

ベクトルの大きさは |a|, |b| で表現できるので、

|a| >= |b|

かどうかを判定すれば、a が b 以上かどうかを判定できます。これらをあわせると、

int dot(int v1[2], int v2[2]) {
  return v1[0]*v2[0] + v1[1]*v2[1];
}

#define _ABS(x) ((x) < 0 ? (x) * -1 : (x))

bool point_on_segment(point p, line s) {
  int a[2] = { s.p2.x - s.p1.x, s.p2.y - s.p1.y };
  int b[2] = { p.x - s.p1.x, p.y - s.p1.y };
  double as = _ABS(sqrt(a[0]*a[0] + a[1]*a[1]));
  double bs = _ABS(sqrt(b[0]*b[0] + b[1]*b[1]));
  double n = dot(a, b) - as*bs;
  return _ABS(n) < 0.000001 && as >= bs;
}

最後にこれまで作成してきた関数を結集して、もともとの問題であった 2 線分が交差するかどうかの判定を実装すると、

bool intersect(line_segment s1, line_segment s2) {
  line l1 = to_line(s1);
  line l2 = to_line(s2);
  point p;
  return intersection_point(l1, l2, &p)
      && point_on_segment(p, s1)
      && point_on_segment(p, s2);
}

これでようやく所望のプログラムを得ることができました。当初の発想としては直感的な方法を採ったつもりが、このように複雑な実装が必要になってしまいました。こういった点が幾何学的アルゴリズムの難しい点であると言えます。

ちなみにこのプログラム、 intersection_point 関数で求める 2 直線の交点を求めた結果を整数にしているので、正確な交点を求めることができていません。つまり、少なくともこの方法では、点の X, Y 座標を整数で持つことはできません。

また、 intersection_point 関数に完全に一致する 2 直線を渡した場合、正確な交点を求めることができないため、FALSE を返しています。本来であれば、同関数から FALSE が返ってきた場合でも、平行する 2 線分が交わるかどうかの例外処理を組む必要があるでしょう。

第二の解

ところで、線分交差の判定は、より簡単な方法を採ることができます。発想を少し変えて、一方の線分の両端点と、もう一方の端点の片方との位置関係を調べることで、判定できます。

一方の線分の両端点と、もう一方の線分の片方の端点を順にめぐるときと、もう片方の端点を順にめぐるときの回転方向が異なっているかどうかを両方の線分について調べれば、線分が交差しているかどうか判定することができます。

fig2.png

ここで必要になる関数は、与えられた 3 点を順にめぐったとき、反時計回りになるか時計回りになるかを判定する関数です。Counter ClockWise より、 ccw と命名します。

int ccw(point p0, point p1, point p2) {
  int dx1 = p1.x - p0.x;
  int dy1 = p1.y - p0.y;
  int dx2 = p2.x - p0.x;
  int dy2 = p2.y - p0.y;
  if (dx1*dy2 > dx2*dy1) return +1;
  if (dx1*dy2 < dx2*dy1) return -1;
  if ((dx1*dx2 < 0) || (dy1*dy2 < 0)) return -1;
  if ((dx1*dx1+dy1*dy1) < (dx2*dx2+dy2*dy2)) return +1;
  return 0;
}

関数の戻り値は 3 値とし、反時計回りのときは 1, 時計回りのときは -1 とします。 例外として 3 点が一直線上に存在する場合がありますが、これについては後述します。 0 についても後述します。

p0 から p1 (線分A)と、 p0 から p2 (線分 B) の両線分の傾きを求め、その大小を比較することによって、反時計回りになるか時計回りになるかを判定できます。線分 A の傾きは dy1/dx1 で、線分 B の傾きは dy2/dx2 で求めることができます。線分 B の傾きが線分 A の傾きよりも大きいとき 3 点は反時計回りに、線分 B の傾きが線分 A の傾きより小さいとき 3 点は時計回りに回ります。つまり、

if (dy2/dx2 > dy1/dx1) return +1;  // 反時計回り
if (dy2/dx2 < dy1/dx1) return -1;  // 時計回り

というプログラムが得られます。このとき、dx1, dx2 は、線分が垂直の場合に 0 になり、その例外を考慮しなければなりません。これは両辺に dx1*dx2 をかけることで、対処できます。よって、

if (dx1*dy2 > dx2*dy1) return +1;  // 反時計回り
if (dx1*dy2 < dx2*dy1) return -1;  // 時計回り

という判定文が得られます。

傾きが等しいとき、3 点は反時計回りでも時計回りでもなく、直線上に存在します。 ccw 関数の基本的な使われ方として、線分上の両端点を第一、第二引数とし、線分とは別の点を第三引数とし、三点の位置関係を調べるというものがあります。ここではそのような使い方を前提として 3 値のとり方を考えます。p1 から p2 を結ぶ線分を線分 C とし、線分 A と線分 C が重なるかどうかを境界として、重なるときは 0 以下の数値, 重ならないときは正の数値を返すようにします。

まず p0p1p2 の間にあるとき -1 を返すことにします。

if ((dx1*dx2 < 0) || (dy1*dy2 < 0)) return -1;

p0p1p2 の間にあれば、線分 A と線分 C は重なります。このとき dx1, dx2 または dy1, dy2 の符号が異なるはずですので、それを利用して判定します。

次に p2 が線分 A の延長線上にあるとき +1 を返すことにします。

if ((dx1*dx1+dy1*dy1) < (dx2*dx2+dy2*dy2)) return +1;

三平方の定理を用いて、線分 A の長さと線分 B の長さを比較しています。線分 B のほうが長ければ、 p2 は線分 A の延長線上にあります。このとき線分 A と線分 B は重なりません。

このような ccw 関数を用いて線分交差の判定を実装すると、

bool intersect(line_segment s1, line_segment s2) {
  return (ccw(s1.p1, s1.p2, s2.p1) * ccw(s1.p1, s1.p2, s2.p2) <= 0)
      && (ccw(s2.p1, s2.p2, s1.p1) * ccw(s2.p1, s2.p2, s1.p2) <= 0);
}

と書けます。

こちらの解も計算量は決して少なくありませんが、必要になるプログラムの量は、最初に導いたものよりずいぶん短くなっています。 ccw 関数のような、よく練られた基本的な関数を用いることで、例外の影響を低減することができます。

単純な閉路

任意の数の点集合において、すべての点を通り、かつ途中で交差せずに出発点に戻るような図形を求める方法を考えます。このような図形は 単純な閉路 と呼ばれます。

fig3.png

ここでは各点を最短距離で通るというような解を求めるといった高度な問題は扱わず、単純にすべての点を交差せずにまわることだけを考えます。

単純な閉路を求める問題は、次のようにして解くことができます:

  1. 点の集合から、Y 座標がもっとも小さい座標を基準点とする
  2. 基準点と各点に引いた直線と、基準点から水平方向に引いた直線とのなす角度を計算する
  3. 2 で求めた角度が小さい順に並べる

dx を 2 点の水平方向の距離、dy を 2 点の垂直方向の距離とするとき、基準点と各点とを結んだ線分がなす角度は、逆正接関数 (atan) を用いて次のように求めることができます。C の標準関数だと、 atan2 関数がより使いやすいので、今回はこちらを利用します:

atan2(dy, dx)

基準点は点集合の中で Y 座標がもっとも小さい点なので、 atan2 の戻り値は 0 <= N <= pi になります。あとはこれらを単純に組み合わせて、

bool less_y(const point& p1, const point& p2) {
  return p1.y != p2.y ? p1.y < p2.y : p1.x < p2.x;
}

struct less_deg {
  point p0;

  less_deg(const point& p) : p0(p) { }

  bool operator()(const point& p1, const point& p2) const {
    int dx1 = p1.x - p0.x;
    int dy1 = p1.y - p0.y;
    int dx2 = p2.x - p0.x;
    int dy2 = p2.y - p0.y;
    double a1 = atan2(dy1, dx1);
    double a2 = atan2(dy2, dx2);
    if (a1 != a2) return a1 < a2;
    if (dy1 != dy2) return a1 < M_PI/2 ? dy1 < dy2 : dy1 > dy2;
    return dx1 < dx2;
  }
};

void closed_path(polygon ps, int N) {
  if (N <= 2)
    return;

  point p0 = *std::min_element(ps, ps + N, less_y);
  std::sort(ps, ps + N, less_deg(p0));
}

と書けます。

less_y 関数は Y 座標が最小の点を見つけるための比較関数です。Y 座標が同一である場合の例外として、その場合には X 座標の大小比較を行います。これにより、Y 座標がもっとも小さく、かつその中でももっとも X 座標が小さい点を見つけることができます。

less_deg 構造体は、基準点と各点に引いた直線と水平方向に引いた直線とがなす角度を用いて、2 点を比較する関数オブジェクトです。

この関数では 2 つの例外を扱う必要があります。

各点のなす角度が等しい場合、なす角度が 90 度未満の場合には Y 座標が小さい点を、90 度以上の場合は Y 座標が大きい点を先にたどるようにします。基準点は Y 座標がもっとも小さい点なので、90 度を境として行きから帰りが切り替わります。行きでは下から順に、帰りでは上から順にたどるようにします。

Y 座標が等しい場合には、X 座標がより小さい点を先にたどるようにします。これは基準点と Y 座標が同じ点の場合を指します。基準点は Y 座標がもっとも小さいもののうち、X 座標がもっとも小さいものを選択しているので、小さい順にたどっていけば交差がないことを保証できます。

点の多角形による包含

次に考える問題は、任意の点が多角形内に含まれているかどうかを判定するという問題です。点と多角形が与えられ、点が多角形内に含まれていれば TRUE を、含まれていなければ FALSE を返す関数を考えます。

まず、任意の数の頂点を持つ、一般的な多角形について考えます。今回は、与えられた点から任意の方向へ十分に長い線分(テスト線分)をのばし、その線分と多角形との交点の数を調べるという方法を採ってみます。交点の数が奇数であれば点は多角形の内側に、偶数であれば外側にあると判定できます。

fig4.png

赤い点が包含を判定する点で、そこからのびている赤い点線がテスト線分です。

ここでもやはり考慮すべき例外があります。テスト線分が多角形の辺を含む場合です。

図中の 3 において、例外を考慮せずにテスト線分が多角形の点のちょうど上を通る場合も交点の数を 1 とすると、交点の数は 4、つまり偶数となります。しかし図を見れば点が多角形に含まれていることは明らかです。また、図中の 4 はその逆が成り立ちます。やはり例外を考慮する必要性があることが分かります。

この例外を考慮して、テスト線分が多角形の点のちょうど上を通るとき、交点の数を 1 とするか 0 とするかを場合分けする必要があります。

bool inside(point t, polygon poly, int N) {
  int count = 0, j = 0;
  line_segment lt, lp;
  polygon p;

  std::copy(poly, poly + N, &p[1]);
  p[0] = p[N];
  p[N + 1] = p[1];

  lt.p1 = t;
  lt.p2 = t;
  lt.p2.x = INT_MAX;

  for (int i = 1; i <= N; ++i) {
    if (ccw(lt.p1, lt.p2, p[i]) != 0) { // (1)
      if (i == j + 1) { // (2)
        lp.p1 = p[i];
        lp.p2 = p[j];
        if (intersect(lp, lt))
          ++count;
      }
      else {
        if (ccw(lt.p1, lt.p2, p[i]) * ccw(lt.p1, lt.p2, p[j]) < 0) // (3)
          ++count;
      }
      j = i;
    }
  }

  return count & 1;
}

lt はテスト線分を表現するオブジェクトです。判定する点 t から水平方向へ線分をのばしています。j は多角形の頂点のうち、テスト線分上になかった最後の頂点をさす添え字です。

基本ルートは (1) と (2) の両条件を満たした場合です。つまり、多角形の頂点がテスト線分上にない場合です。多角形の各辺とテスト線分とが交差するかどうかを判定し、交差していれば交点の数 count+1 します。

例外ルートを考えます。(1) の条件を満たさない場合、頂点 p[i] はテスト線分上にあると判定できます。その場合には j を更新せずに次の頂点を調べます。このように j を更新しなかった後で (1) の条件を満たすと、(2) の条件を満たさず、(3) の処理に流れます。(3) では p[i]p[j] についてテスト線分との ccw を取っています。p[i]p[j] はいずれもテスト線分上に存在しない点であり、これらがテスト線分の両側に位置するときのみ、交点の数を +1 します。

この解では、与えられる多角形 poly の頂点が反時計回りに並んでおり、かつ poly[0] が、全頂点のうち Y 座標がもっとも小さく、かつその内でもっとも X 座標がもっとも小さい点をであることを前提としています。

やはり、例外を考慮すると、それなりに複雑なプログラムになってしまいます。

点の長方形による包含

先の問題では、与えられた点が任意の数の頂点を持つ多角形に含まれているかどうかを判定しました。このような一般的な多角形ではなく、長方形のようなより単純で特徴のある多角形であれば、より単純な実装方法が考えられます。

// 長方形
struct rect {
  int x1, y1, x2, y2;
};

bool insiderect(point t, rect r) {
  return p.x >= r.x1
      && p.x <= r.x2
      && p.y >= r.y1
      && p.y <= r.y2;
}

長方形は左下と右上の頂点座標で表現します。x1, y1 は左下の頂点の X 座標と Y 座標を、x2y2 は右上の頂点の X 座標と Y 座標を持つものとします。点の X 座標と Y 座標の両方が、長方形の上下左右の X 座標、Y 座標の範囲内におさまっているかどうかをしています。至って自明なプログラムです。

一般的な多角形の場合と比べると、ずいぶん簡単なプログラムになりました。長方形は一般的な多角形に比べて極めて単純であり扱う例外の数が少ないため、このように単純なプログラムに落としこむことができます。

領域探索

最後に、領域探索について考えます。領域探索とは、点集合の中から、与えられた領域の中に含まれる点をすべて探しだすという問題です。一次元的な要素であれば、二分探索を利用することで簡単かつ効率的に探索することができます。今回扱いたいのは二次元的な要素であり、二分探索をそっくりそのまま使うことはできません。

逐次探索

まず単純に思いつく方法として、点集合に含まれる点それぞれについて、与えられた領域内に含まれるかどうかを判定していく方法です。逐次探索と呼ぶことにします。

class Range {
private:
  struct node {
    point p;
    node* next;
  };
  node* z;
  node* head;

public:
  Range();
  void insert(point p);
  int search(rect range);
};

Range::Range() {
  z = new node;
  head = z;
}

void Range::insert(point p) {
  node* n = new node;
  n->p = p;
  n->next = head;
  head = n;
}

int Range::search(rect range) {
  int count = 0;
  for (node* n = head; n != z; n = n->next) {
    if (insiderect(n->p, range))
      ++count;
  }
  return count;
}

もう少し工夫して、X 座標で二分探索を実施し、次に X 座標の領域探索に引っかかった点集合について Y 座標で二分探索を実施すれば、包含判定の実施回数は減らせそうです。ここまでやれば、小さいサイズの問題に対しては十分でしょう。

言うまでもなく、逐次探索の実行時間は問題のサイズに線形に増加します。問題のサイズが大きい場合には適さない手法と言えます。

バケット法

次に考える方法は、平面を正方形の区画で分割し、探索する領域を少なくするというものです。この正方形の区画をバケットと呼び、この領域探索方法をバケット法と呼ぶことにします。今回はバケットを平面上に縦横同じ数だけ並べることにします。

fig5.png

この図では、青い枠で表した検索領域と交わる 2, 3, 6, 7 のバケットについてのみ、その区画内に含まれる点に対して検索領域内に含まれるかどうかを調べます。

バケット法を実装するにあたって問題になるのは、平面上をいくつのバケットで分割するかという点と、バケットの辺の長さをどれだけにするかという点です。

  • 平面上に存在する点の数を N
  • 1バケットあたりの点の数(の期待値)を M
  • バケットの行数(および列数)を G
  • バケットの辺の長さを S

とおきます。約 N / M 個のバケットで平面を分割すればよいので、

G = ceil(sqrt(N / M))
S = ceil(座標の最大値 / G)

と導くことができます。 座標の最大値 とは、すべての点の X および Y 座標のうち、もっとも大きい数値のことを指します。平面の座標空間が広くなければ、平面の縦幅、横幅のうち長いほうを採用するという方法で選んでもいいかもしれません。 N は問題により与えられる値であり、 G, S は上記の式から導くことができました。 M は自分で適当な値を決める必要があるでしょう。

以下のプログラムでは、G および S があらかじめ定数として定義されていると仮定しています:

class Range {
private:
  struct node {
    point p;
    node* next;
  };
  node* z;
  node* buckets[G][G];

public:
  Range();
  void insert(point p);
  int search(rect range);
};

Range::Range() {
  z = new node;
  for (int i = 0; i < G; ++i)
    for (int j = 0; j < G; ++j)
      buckets[i][j] = z;
}

void Range::insert(point p) {
  node *n = new node;
  n->p = p;
  n->next = buckets[p.x / S][p.y / S];
  buckets[p.x / S][p.y / S] = n;
}

buckets は点のリスト (node オブジェクトへのポインタ) の配列として表現しています。コンストラクタ内で全要素を空のリストとして初期化します。

点の挿入では、点の X 座標と Y 座標をそれぞれバケットの辺の長さ S で割ることで、点が属するべきバケットを割り出しています。

int Range::search(rect range) {
  int count = 0;
  for (int i = range.x1 / S; i <= rect.x2 / S; ++i) {
    for (int j = range.y1 / S; j <= rect.y2 / S; ++j) {
      for (node* n = buckets[i][j]; n != z; n = n->next) {
        if (insiderect(n->p, range))
          ++count;
      }
    }
  }
  return count;
}

点の検索では、検索するべきバケットを検索領域の端点から割り出しています。挿入時と同じように、 S で割ることで割り出しています。検索対象のバケットに含まれるすべての点について insiderect 関数を用いて領域内に含まれるかどうかを判定します。

バケット法は、点が平面上に一様に分散している場合に強さを発揮します。一方で、点が一部のバケットに集中している場合、そのバケットをに引っかかる検索をするとすべての点を調べる必要があるので、あまり役に立ちません。また、 GS が動的かつ頻繁に変わるような問題の場合には buckets の再構築のコストが無視できなくなるでしょう。

2D 木 (kD 木)

最後に 2D 木(2 次元木)による解を考えます。2D 木とは領域探索に適したデータ構造であり、基本的な考え方は二分木と同じであると言えます。またバケット法の弱点のひとつである問題の動的さへの弱さおいても、二分木と同じような性質を持つ 2D 木であれば特に問題はありません。

2D 木と二分木の異なる点は、木を構築する際に用いるキーを、階層によって交互に切り替えるという点です。

  • ルートとその子(第一階層)を分かつのは Y 座標
  • 第一階層と第二階層を分かつのは X 座標
  • 第二階層と第三階層を分かつのは Y 座標
  • ...

といった具合です。

fig6.png

2D 木によって平面が分割されていく様子を図にしています。

  1. 点 A を木に追加します。第一階層を Y 座標で分かつため、水平方向に直線を引きます
  2. 点 B を木に追加します。第二階層を X 座標で分かつため、垂直方向に直線を引きます。このとき点 A の Y 座標より上側については気にしません
  3. 点 C を木に追加します。第三階層を Y 座標で分かつため、水平方向に直線を引きます。このとき点 B の X 座標より右側については気にしません
  4. 点 D を木に追加します。第二階層を X 座標で分かつため、垂直方向に直線を引きます。このとき点 A の Y 座標より下側については気にしません
  5. 点 E を木に追加します。第三階層を Y 座標で分かつため、水平方向に直線を引きます。このとき点 D の X 座標より左側については気にしません

こうしたプロセスを得て構築された 2D 木は、次のようになります:

fig7.png

点 A より上に位置している点が A ノードの右側に、点 A より下に位置している点が A ノードの左側にぶら下がっていることが分かります。点 A より上に位置している点のうち、点 D より右側に位置している点は D ノードの右側に、点 D より左側に位置している点は D ノードの左側にぶら下がっていることが分かります。

2D 木のプログラムは至ってふつうの二分木とほとんど同じで、階層ごとに X 座標と Y 座標を切り替えるようにするだけです。

class Range {
private:
  struct node {
    point p;
    node *l, *r;
  };
  node *z, *head;

public:
  Range();
  void insert(point p);
  int search(rect range);

private:
  int search_r(node* t, rect range, int d);
};

Range::Range() {
  z = new node;
  head = new node;
  head->p = { INT_MIN, INT_MIN };
  head->l = z;
  head->r = z;
}

void Range::insert(point p) {
  node *f, *t;
  int d, left;
  for (d = 0; t = head; t != z; d = !d) {
    left = d ? (p.x < t->p.x) : (p.y < t->p.y);
    f = t;
    t = left ? t->l : t->r;
  }

  t = new node;
  t->p = p;
  t->l = z;
  t->r = z;
  if (left)
    f->l = t;
  else
    f->r = t;
}

点の挿入では、挿入するべきノードを木構造を根から葉へ下りながら探します。探す際のキーを深さに応じて X と Y で切り替えています。

int Range::search(rect range) {
  return search_r(head->r, range, 1);
}

int Range::search_r(node* t, rect range, int d) {
  if (t == z)
    return 0;

  int tx1 = range.x1 < t->p.x;
  int tx2 = range.x2 <= t->p.x;
  int ty1 = range.y1 < t->p.y;
  int ty2 = range.y2 <= t->p.y;
  int t1 = d ? tx1 : ty1;
  int t2 = d ? tx2 : ty2;

  int count = 0;
  if (t1)
    count += search_r(t->l, range, !d);
  if (insiderect(t->p, range))
    ++count;
  if (t2)
    count += search_r(t->r, range, !d);

  return count;
}

点の探索でも、木構造を根から葉へ下りながら探します。挿入時と同じように、次に探索するノードを決定するキーを、探索中のノードの深さに応じて X と Y で切り替えています。

挿入のときと異なる点は、左の子ノードを探すか、右の子ノードを探すか、または両方を探すかという場合分けが必要であるという点です。検索領域が現在のノードの左側(下側)におさまっている場合、左の子ノードだけを探索すればよいということになります。逆に現在のノードの右側(上側)におさまっている場合、右の子ノードだけを探索すればよいということになります。検索領域が左右(上下)にまたがっている場合、両方の子ノードを探索する必要があります。

2D 木はバケット法に比べて

  • 点の偏りに強い
  • (場合によっては)使用するメモリ空間が少ない

という特徴があります。

2D 木はその名の通り 2D、つまり二次元空間を検索するためのデータ構造ですが、上記のアルゴリズムをそのまま拡張して、三次元、四次元と、扱う次元数を増やしていくことができます。より一般に、k 次元の空間を検索するためのデータ構造を kD 木と呼びます。実装も上記のアルゴリズムをほとんどそのまま採用でき、X と Y を交互に切り替えていたところを、階層をおうごとに 1, 2, 3, ... k とキーを切り替えていけばよいのです。

参考文献