計算幾何 〜 2次元凸包 理論と実践 〜

in English

はじめに

　前回すこしふれた2次元凸包アルゴリズムについて説明し、C++による実装例を示す。

アルゴリズム

　今回紹介するアルゴリズムは逐次添加法と呼ばれ、$O(n\log{n})$の計算量を持つ¹。

入力：$A=\{P_1,P_2,\cdots,P_n\}$（添字を固有番号と呼ぶことにする）
出力：凸包 $C(A)$

1. $A$に属する点を $x$ 座標の小さい順にならべる。この点列を $A^{\prime}=\{P_{(1)},P_{(2)},\cdots,P_{(n)}\}$ と書き、添字を整列番号と呼ぶことにする。$x$ 座標が同じ値を持つ場合は固有番号の大きいものから並べる。整列番号から固有番号へ変換するテーブル $T_{{\rm sorted}\rightarrow{\rm original}}$ を保持しておく。
2. $A^{\prime}$ の最初の3点から三角形を作り、これを $C(\{P_{(1)},P_{(2)},P_{(3)}\})$ とする。時計回りにたどるときの点列 $V_{\rm clockwise}$ と反時計回りにたどるときの点列 $V_{\rm counterclockwise}$ を保持する。たとえば、時計回りにまわったとき整列番号 $a$ の次に整列番号 $b$ の点がくるなら$V_{\rm clockwise}(a)=b$となるようにすればよい。
3. 整列番号 $(3)$ から $(n)$ までの点を順に追加する。各点の追加ごとにアルゴリズム $X$（後述）を実行する。
4. 出力値 $C(\{P_{(1)},P_{(2)},\cdots,P_{(n)}\})$ を得る。$V_{\rm clockwise}$ あるいは $V_{\rm counterclockwise}$ をたどれば、凸包を得ることができる。

　先のページで述べた記号摂動法を適用する際、各座標に以下のように摂動を加える。 $$\begin{eqnarray} (x_i,y_i)&\rightarrow&(x_i+\epsilon^{M^i},y_i+\epsilon^{M^{n+i}}) \nonumber \\ &\equiv&(x_i^{\ast},y_i^{\ast}) \end{eqnarray}$$ ここで、$n$ は点の総数、$M$ は十分大きな量、$\epsilon$ は十分小さな量である。添字 $i$ として固有番号を採用することにより、たとえ $x$ 座標値が同じであっても、仮想的に $x$ 座標値の小さな順に点を添加していくことができる。たとえば、$x$座標値が同じ3点を考えてみる（下図）。括弧で囲んだ数字は整列番号、囲んでいない番号は固有番号である。$x$ 座標が同じ値を持つ場合は固有番号の大きいものから並べたので、各点の摂動量（青線）と摂動後の点（赤丸）は図のようになる。整列番号の小さい点は摂動量も小さくなることが分る。

　次にアルゴリズム $X$ について述べる。

入力：$C(\{P_{(1)},P_{(2)},\cdots,P_{(k)}\})$、$P_{(k+1)}$、$V_{\rm clockwise}$、$V_{\rm counterclockwise}$、$T_{{\rm sorted}\rightarrow{\rm original}}$
出力：$C(\{P_{(1)},P_{(2)},\cdots,P_{(k+1)}\})$

1. $P_{(s)}\leftarrow P_{(k)}$
2. 3点 $P_{(k+1)},P_{(s)},P_{(t)}$ を考える。ここで、$t=V_{\rm counterclockwise}(s)$である。
3. 3点をこの順にたどったとき右に曲がるなら $s\leftarrow t$ として2へ、左に曲がるなら4へ（下図参照）。
4. 左に曲がった整列番号$(u)$を保持する。
5. $P_{(s)}\leftarrow P_{(k)}$
6. 3点 $P_{(k+1)},P_{(s)},P_{(t)}$ を考える。ここで、$t=V_{\rm clockwise}(s)$である。
7. 3点をこの順にたどったとき左に曲がるなら $s\leftarrow t$ として6へ、右に曲がるな8へ。
8. 右に曲がった整列番号$(l)$を保持する。
9. $V_{\rm counterclockwise}(l)=k+1$、$V_{\rm counterclockwise}(k+1)=u$ とする。
10. $V_{\rm clockwise}(u)=k+1$、$V_{\rm clockwise}(k+1)=l$ とする。

　右に曲がる、左に曲がるという判定は先のページで取り上げた判定式 $F$ を使えば良い。

1. $F(P_{(i)},P_{(j)},P_{(k)})$ を考える。
2. $T_{{\rm sorted}\rightarrow{\rm original}}$ を使って、整列番号を固有番号に変換する。 $$\begin{eqnarray} a&=&T_{{\rm sorted}\rightarrow{\rm original}}(i) \nonumber \\ b&=&T_{{\rm sorted}\rightarrow{\rm original}}(j) \nonumber \\ c&=&T_{{\rm sorted}\rightarrow{\rm original}}(k) \nonumber \end{eqnarray}$$
3. 固有番号のリスト $L=(a,b,c)$ を大きい順に並べ替え、リスト $L^{\prime}=(a^{\prime},b^{\prime},c^{\prime})$ を作成する。並べ替えの際、奇置換なら $p=-1$、偶置換なら $p=1$ とする。
4. $F(P_{a^{\prime}},P_{b^{\prime}},P_{c^{\prime}})$ の正負を判定する。その結果に$p$の値を乗じたものが $F(P_{a},P_{b},P_{c})$ の符号となる。正なら左に折れ、負なら右に折れる。

C++による実装

ソースはここ。main関数は以下の通りである。

1. 16行目から39行目：乱数による点の生成を行っている。
2. 43行目：テンプレートクラスConvexHull2dのオブジェクトを作成している。ここではテンプレート引数としてintを指定したが、boost::multiprecision::cpp_intを指定することもできる。
3. 44行目：アルゴリズムを実行している。
4. 46行目以降：最初に生成した点と凸包を描画している。
5. 57行目：凸包を構成する点を取り出している。

デモ

開発環境

1. Mac OS X 10.8.4
2. プロセッサ：3.06 GHz Intel Core 2 Duo
3. メモリ：4GB
4. Xcode4.6.2 with Apple LLVM 4.2（C++ Language Dialect → C++11, C++ Standard Library → libc++）
5. boost-1.53.0（Apple LLVM 4.2でコンパイルしたもの。こちら）
6. opencv-2.4.5（Apple LLVM 4.2でコンパイルしたもの。こちら）

参考文献

1. 計算幾何プログラミング, 杉原厚吉, 岩波書店