特にドロー系の画像編集ソフトなどでは、指定された2次元領域に関係しているオブジェクトを素早く特定して処理を行うために空間インデックスといわれる類のデータ構造をしばしば用いている。その内の一つである領域四分木というデータ構造とアルゴリズムを見ていき、Pythonで素朴に実装してみる。具体的なチンパンコード付き。

データ構造とアルゴリズム

四分木とは、次のような問題を解くことのできるアルゴリズムの1つであり、場合にもよるがQ-treeだとかQuad-treeとも呼ばれる。

脳死アルゴリズムによってこの問題を解くと計算ステップのオーダはO(MN)である。二分木といった木構造を用いると計算ステップのオーダはO((N+M)log N)となる。

二分木とは

四分木にいくまえにその簡単な場合である二分木についてみていく。二分木とは、大小関係を持ったデータをある意味でソートしながら記憶し、ある点に最も近いデータ、あるいはある区間内にあるデータの列挙といったクエリに対して効率良く答えることが出来る。ここでいう大小関係とは、数直線上の点であればそのまま数値の大小関係を用いればよいし、文字列であれば辞書に出てくる順番の前後関係を大小関係と見なせばよい。

具体的にどうするかというと、次の図のようなデータ点の1つ1つが頂点に対応する木構造を構築する。

一番上から順番に自分との大小関係を比較していき、自分のいる側に移動していって最終的にそれ以上エッジが無くなったらその先に自分を追加すればよい。全ての頂点は高々2つのエッジを持つという特徴がある。クエリに対しては直接的に二分探索を使って応答することが出来る。

四分木とは

二分木は大小関係を持ったデータに対して適用できるとあったが平面上の点というデータに対しては使えるだろうか？とにかく大小関係を与えればよい、という点では例えば「原点との距離」を用いれば可能だろうが先に上げた例のようなクエリに対してはふさわしくないだろう。あるいはx座標のみ参照して大小関係を与えるという方法もあるがやや微妙。そこで四分木というデータ構造がある。やることは簡単で、xの大小関係とyの大小関係のそれぞれ2つによってエッジを生成していく。そのほかのやり方は二分木の場合と同じ。

領域四分木

単純な四分木でも先の問題に答えることは出来そうだが、まだやや問題がある。二分木とその素朴な拡張による四分木のどちらに対しても言えることとして、生成される木構造がデータの入力順に依存しているということがある。二分木の場合の極端な例として、値の小さいものから順番にデータを木構造に追加していくと右側エッジの先にしかノードが生えず、結果として単なるリストと等価なものが出来上がってしまう。これではlogNにはならないので二分木はデータ追加の過程でルートノードを付けかえる木構造のバランシング処理を行うことが普通だ。四分木の場合にも同じ問題は起こり、コーディングの面とバランス判定基準の調整の面の両方でバランシング処理はより面倒になる。

四分木の場合には、木構造をデータ入力の順番に依存させない領域四分木という方法がある（領域二分木というものが無いわけではないが聞いたことが無い、なぜか？）。領域四分木では全データのバウンディングボックスを再帰的に4分割していくことによって木構造を構成していく。分割のイメージはページ上部のものを参考にされたい。

さっきまでのデータ構造とは異なり、それぞれの頂点は領域を表しデータは木構造の「葉」としてグラフの末端に格納される。もし例の問題のように各座標の上下限がわからなくても一度データを舐めればオーダO(N)でバウンディングボックスの算出は可能である。明らかな特徴として、領域四分木ではデータを入力する順番を変えても得られる木構造は同じになる。

PythonによるQ-treeのコーディング

最初に言っておかなければならないのは、これは実験であり車輪の再開発である。どんな言語でもあると思うがPythonにもpyqtreeというQ-treeの実装を与えるパッケージはあるのですぐに何かをやりたい場合はそれを使うのが賢明だろう。

Pythonプログラム

import numpy as np

QmaxD = 50 # maximum depth
QmaxN = 5  # maximum items per boxel
Qtree = {}
Qbbox = [-1, 1, -1, 1] #xmin, xmax, ymin, ymax

def Qreset(bbox, QmaxN=5, QmaxD=50):
  Qtree = {}
  Qbbox = list(bbox)

def Qinsert(point):
  qtree, qbbox = Qtree, list(Qbbox)

  for _ in range(QmaxD):
    key, x0, y0 = 0, (qbbox[0]+qbbox[1])/2, (qbbox[2]+qbbox[3])/2
    if point[0] < x0:
      qbbox[1] = x0
      key += 1
    else:
      qbbox[0] = x0
    if point[1] < y0:
      qbbox[3] = y0
      key += 2
    else:
      qbbox[2] = y0

    if not key in qtree:
      qtree[key] = []
      break
    if type(qtree[key]) is list:
      break

    qtree = qtree[key]

  if len(qtree[key]) < QmaxN:
    qtree[key].append(point)
  else:
    subtree = {}
    x0, y0 = (qbbox[0]+qbbox[1])/2, (qbbox[2]+qbbox[3])/2
    for point in qtree[key]+[point]:
      subkey = 0
      if point[0] < x0: subkey += 1
      if point[1] < y0: subkey += 2
      if subkey in subtree:
        subtree[subkey].append(point)
      else:
        subtree[subkey] = [point]
    qtree[key] = subtree

def Qquery(bbox):
  result = []
  Qstack = [(Qtree, list(Qbbox), 0)]
  pc = 0
  while len(Qstack) > 0:
    pc += 1
    qtree, qbbox, key = Qstack[-1]
    if type(qtree) is list:
      for point in qtree:
        if bbox[0] <= point[0] and point[0] <= bbox[1] and bbox[2] <= point[1] and point[1] <= bbox[3]:
          result.append(point)
      Qstack.pop()
      continue

    if key > 3:
      Qstack.pop()
      continue

    x0, y0, qbbox = (qbbox[0]+qbbox[1])/2, (qbbox[2]+qbbox[3])/2, list(qbbox)
    if key % 2 == 1:
      qbbox[1] = x0
    else:
      qbbox[0] = x0
    if key // 2 == 1:
      qbbox[3] = y0
    else:
      qbbox[2] = y0

    Qstack[-1] = (Qstack[-1][0], Qstack[-1][1], key+1)
    if key in qtree:
      if qbbox[0] <= bbox[1] and bbox[0] <= qbbox[1] and qbbox[2] <= bbox[3] and bbox[2] <= qbbox[3]:
        Qstack.append((qtree[key], qbbox, 0))
  return result

はい。結構長くなってしまっている。読むのは何かしらのトレーニングにはなるだろうが苦痛なのでやめたほうがいいだろう。これをmyqtree.pyとしていろいろビジュアライズしていこう。

プログラムの解説

解説といってもコードを追うようなことはしないが、Pythonのdictを使いまくっている。パフォーマンスを追求することは目的ではないのが体感クエリに関しては相当早く応答する（dictの特性からどんなにPythonの中で変なことしてても速いのは当たり前だが）。ちなみに予想どおりdictをアロケートしまくるインサートは遅い。

領域四分木では一つの頂点（領域）にぶら下がる葉（データ）の数に上限を設け、その上限を超えたときにその領域を分割する。このプログラムではその上限をQmaxN=5としている。ここで１つだけ問題があり、まったく同一の点がQmaxN個以上入力されると、実装にもよるが「ボクセル内のデータ数がQmaxNを超えてしまう」か「無限に領域を分割するループに入る」のいずれかが起こってしまう。後者の対策としてQmaxDとして深さの上限を設けているが50だと全然ヒットしないので放置してある。ついでにこの実装だとおそらく前者が起こるので気にしなくていいはず。

それにしてもこの手の車輪をC/C++で再開発していたころに比べて大分やりやすくなったもんだな。

疑似データによる実行結果

バウンディングボックスは[-1,+1]^2の範囲としてその中に一様乱数を用いて2000個のランダム点を散らす。3つのクエリに対する応答を可視化した結果次のようになる。これに用いたコードは次のようになる。

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as ptch
from myqtree import Qreset, Qinsert, Qquery

points = []
for i in range(2000):
  point = np.random.rand(2)*2-1
  points.append(point)
  Qinsert(point)

bboxes = [
  [0.25, 0.4, 0.25, 0.4],
  [-0.5, -0.3, 0.2, 0.9],
  [-0.9, -0.4, -0.5, -0.4]
]

fig, ax = plt.subplots(figsize=(5,5))
ax.set_xlim([-1,1])
ax.set_ylim([-1,1])
ax.set_aspect('equal')

P = np.array(points).T
ax.scatter(P[0], P[1], s=3)

for bbox in bboxes:
  result = Qquery(bbox)
  if len(result) > 0:
    R = np.array(result).T
    ax.scatter(R[0], R[1], s=6)
    ax.add_patch(ptch.Rectangle((bbox[0], bbox[2]), bbox[1]-bbox[0], bbox[3]-bbox[2], fill=False, lw=1))
plt.show()
plt.close(fig)

これに対する結果は以下。

プロットした点の大きさの関係ではみ出ているように見えるやつもあるがクエリで与えられたボックスの中にあるものをもれなく列挙出来ている。

計算ステップの計測実験

ここまでであれば高々2000点だし脳死アルゴリズムとの違いは大してわからない。そこで、クエリ数M=3（先の例で使ったもの）を固定してデータ数Nを大きくしていったときに、計算量の目安としてループカウンタがどこまで大きくなるかを計測した。ついでに予想される計算量の曲線もプロットしていく。素朴に考えるとインサートにNlogN、クエリには3logNのオーダで曲線は伸びていくはずである。

これが結果。insertのループカウンタについてはおおむね予想通りの回数が出ている。ところがqueryについては大幅に大きなステップ数（の増加オーダ）がでています。なぜだろうと考えると、このqueryのループカウンタは四分木のエッジを辿るところにも1ループ使っているので本来ならば2^Nのオーダが出るはずの深さ優先のイテレーションです。深さ優先の利点は記憶域のオーダが高々木構造の深さ倍で済むところで、さらに今回の場合はエッジの先にある領域が完全にバウンディングボックスからはみ出ている場合にはそのエッジを切り落とすことができます（枝刈り）。この切り落とした回数をcutとしてプロットしてみると枝刈りによって相当な恩恵を受けていることがわかります。実際、2^NのオーダはプロットしていませんがNlogNよりも急峻なはずですからqueryの計算ステップは2^Nと比べて大幅に削減されていることがわかります。

まとめ

以上でQ-treeのデータ構造とアルゴリズム、さらにPythonのチンパンコードによる軽い実験をおわる。最後のところはM=O(N)としてみるとN^2と(N+M)logN=NlogNの比較になるので領域四分木が速いだろうことはわかるだろう。ついでにチンパンコードから得た所感はクエリに対してはべらぼうに早いのでM=大だともっと強くなるだろうね。対数目盛でわかりにくいがN=5,000,000までPython3.8.1+Windows10の環境で動いているのでPythonはメモリ1G喰っても文句言わずに動いて偉いですな。