LoginSignup
1
0
@zhidao(Zhidao)

最小包含球を求めるエレガントなアルゴリズム

Last updated at Posted at 2024-01-29

はじめに

$n$個の点から成る点群$\mathcal{P}=\left\{\boldsymbol{p}_{i}|i=1\sim N\right\}$を内部および表面上に全て含む球のうち，最小半径のものを 最小包含球（Bounding Ball, Smallest Enclosing Ball） と呼びます．最も単純な包含立体の一つであり，領域判定や干渉判定などのいろいろな場面で有用なものです．これを求めるアルゴリズムは幾つかありますが，本記事ではその中でもエレガントな，次の論文で提案されているものを解説します．

  • Emo Welzl, Smallest enclosing disks (balls and ellipsoids), in New Results and New Trends in Computer Science, Lecture Notes in Computer Science, Springer-Verlag, pp. 359--370, 1991.

原著論文は著者自身が公開しているものが読めます．

以下では，点群$\mathcal{P}$の最小包含球を$\mathcal{BB}(\mathcal{P})$と表すことにします．

点の数N≦4のときの最小包含球

点の数$N=1$のときは，$\mathcal{BB}(\mathcal{P})$は$\mathcal{P}$に含まれるただ一つの点を中心とする半径$0$の球となります．最小包含球が実用上の意味を持つのは$N\geq 2$の場合で，このとき次の命題$\star$が成り立ちます．

命題$\star$）$N\geq 2$のとき，$\mathcal{BB}(\mathcal{P})$は次の3つの条件のうちいずれかを必ず満たす．

  1. $\mathcal{P}$に含まれるある2点を直径の両端に持つ．
  2. $\mathcal{P}$に含まれるある3点から成る三角形の外接円を最大断面とする．
  3. $\mathcal{P}$に含まれるある4点から成る四面体の外接球である．

1つめの条件を満たす球は${}_{N}C_{2}$個，2つめの条件を満たす球は${}_{N}C_{3}$個，3つめの条件を満たす球は${}_{N}C_{4}$個あります．$N\leq 4$のときは高々${}_{4}C_{2}+{}_{4}C_{3}+{}_{4}C_{4}=10$通りの候補を実際に求め，その中から最小半径のものを選べば良いわけです．

具体的にそれぞれの球を求めてみます．データ構造としては，球は中心座標$\boldsymbol{p}_{\mathrm{C}}$と半径$r$の組で表されます．

2点を直径の両端点とする球

これは単純に，中心座標$\boldsymbol{p}_{\mathrm{C}}=\frac{\displaystyle\boldsymbol{p}_{i}+\boldsymbol{p}_{j}}{\displaystyle 2}$，半径$r=\frac{\displaystyle \|\boldsymbol{p}_{i}-\boldsymbol{p}_{j}\|}{\displaystyle 2}$と求まります．

3点を頂点とする三角形の外接円を最大断面とする球

image.png

上図のように，三角形の頂点を$\mathrm{P}_{i}(\boldsymbol{p}_{i})$，$\mathrm{P}_{j}(\boldsymbol{p}_{j})$，$\mathrm{P}_{k}(\boldsymbol{p}_{k})$，中心を$\mathrm{P}_{\mathrm{C}}(\boldsymbol{p}_{\mathrm{C}})$とそれぞれおきます．$s_{i}=\triangle\mathrm{P}_{\mathrm{C}}\mathrm{P}_{j}\mathrm{P}_{k}$，$s_{j}=\triangle\mathrm{P}_{\mathrm{C}}\mathrm{P}_{k}\mathrm{P}_{i}$，$s_{k}=\triangle\mathrm{P}_{\mathrm{C}}\mathrm{P}_{i}\mathrm{P}_{j}$とそれぞれおいたとき，

\boldsymbol{p}_{\mathrm{C}}=\frac{s_{i}\boldsymbol{p}_{i}+s_{j}\boldsymbol{p}_{j}+s_{k}\boldsymbol{p}_{k}}{s_{i}+s_{j}+s_{k}}

であることが知られています．
円周角の定理より$\angle\mathrm{P}_{j}\mathrm{P}_{\mathrm{C}}\mathrm{P}_{k}=2\angle\mathrm{P}_{j}\mathrm{P}_{i}\mathrm{P}_{k}$，$\angle\mathrm{P}_{k}\mathrm{P}_{\mathrm{C}}\mathrm{P}_{i}=2\angle\mathrm{P}_{k}\mathrm{P}_{j}\mathrm{P}_{i}$，$\angle\mathrm{P}_{i}\mathrm{P}_{\mathrm{C}}\mathrm{P}_{j}=2\angle\mathrm{P}_{i}\mathrm{P}_{k}\mathrm{P}_{j}$ですので，

\displaylines{
s_{i}=\frac{1}{2}r^{2}\sin 2\angle\mathrm{P}_{j}\mathrm{P}_{i}\mathrm{P}_{k}
\\
s_{j}=\frac{1}{2}r^{2}\sin 2\angle\mathrm{P}_{k}\mathrm{P}_{j}\mathrm{P}_{i}
\\
s_{k}=\frac{1}{2}r^{2}\sin 2\angle\mathrm{P}_{i}\mathrm{P}_{k}\mathrm{P}_{j}
}

さらに$l_{i}=\|\boldsymbol{p}_{j}-\boldsymbol{p}_{k}\|$，$l_{j}=\|\boldsymbol{p}_{k}-\boldsymbol{p}_{i}\|$，$l_{k}=\|\boldsymbol{p}_{i}-\boldsymbol{p}_{j}\|$とおくと，正弦定理

\displaylines{
\sin\angle\mathrm{P}_{j}\mathrm{P}_{i}\mathrm{P}_{k}=l_{i}/2r
\\
\sin\angle\mathrm{P}_{k}\mathrm{P}_{j}\mathrm{P}_{i}=l_{j}/2r
\\
\sin\angle\mathrm{P}_{i}\mathrm{P}_{k}\mathrm{P}_{j}=l_{k}/2r
}

および余弦定理

\displaylines{
\cos\angle\mathrm{P}_{j}\mathrm{P}_{i}\mathrm{P}_{k}=\frac{l_{j}^{2}+l_{k}^{2}-l_{i}^{2}}{2l_{j}l_{k}}
\\
\cos\angle\mathrm{P}_{k}\mathrm{P}_{j}\mathrm{P}_{i}=\frac{l_{k}^{2}+l_{i}^{2}-l_{j}^{2}}{2l_{k}l_{i}}
\\
\cos\angle\mathrm{P}_{i}\mathrm{P}_{k}\mathrm{P}_{j}=\frac{l_{i}^{2}+l_{j}^{2}-l_{k}^{2}}{2l_{i}l_{j}}
}

より，倍角公式を使えば

\displaylines{
s_{i}=\frac{1}{2}r^{2}\cdot\frac{l_{i}}{2r}\cdot\frac{l_{j}^{2}+l_{k}^{2}-l_{i}^{2}}{2l_{j}l_{k}}
\\
s_{j}=\frac{1}{2}r^{2}\cdot\frac{l_{j}}{2r}\cdot\frac{l_{k}^{2}+l_{i}^{2}-l_{j}^{2}}{2l_{k}l_{i}}
\\
s_{k}=\frac{1}{2}r^{2}\cdot\frac{l_{k}}{2r}\cdot\frac{l_{i}^{2}+l_{j}^{2}-l_{k}^{2}}{2l_{i}l_{j}}
}

となります．中心座標を求める際には，これらの比だけが意味を持つので，

\displaylines{
s_{i}=l_{i}^{2}(l_{j}^{2}+l_{k}^{2}-l_{i}^{2})
\\
s_{j}=l_{j}^{2}(l_{k}^{2}+l_{i}^{2}-l_{j}^{2})
\\
s_{k}=l_{k}^{2}(l_{i}^{2}+l_{j}^{2}-l_{k}^{2})
}

としても問題ありません．これらを上の式に代入すれば，$\boldsymbol{p}_{\mathrm{C}}$を得ます．
半径は，$r=\|\boldsymbol{p}_{i}-\boldsymbol{p}_{\mathrm{C}}\|$と求まります．

4点を頂点とする四面体の外接球

頂点を$\mathrm{P}_{i}(\boldsymbol{p}_{i})$，$\mathrm{P}_{j}(\boldsymbol{p}_{j})$，$\mathrm{P}_{k}(\boldsymbol{p}_{k})$，$\mathrm{P}_{l}(\boldsymbol{p}_{l})$，中心を$\mathrm{P}_{\mathrm{C}}(\boldsymbol{p}_{\mathrm{C}})$とそれぞれおきます．
外接球の中心は，辺$\overline{\mathrm{P}_{i}\mathrm{P}_{l}}$の二等分面，辺$\overline{\mathrm{P}_{j}\mathrm{P}_{l}}$の二等分面，辺$\overline{\mathrm{P}_{k}\mathrm{P}_{l}}$の二等分面の交点として求まります．このことから，$\boldsymbol{e}_{i}=\boldsymbol{p}_{i}-\boldsymbol{p}_{l}$，$\boldsymbol{e}_{j}=\boldsymbol{p}_{j}-\boldsymbol{p}_{l}$，$\boldsymbol{e}_{k}=\boldsymbol{p}_{k}-\boldsymbol{p}_{l}$，$\boldsymbol{e}_{\mathrm{C}}=\boldsymbol{p}_{\mathrm{C}}-\boldsymbol{p}_{l}$とそれぞれおくと，次が成り立ちます．

\displaylines{
\boldsymbol{e}_{i}^{\mathrm{T}}\boldsymbol{e}_{\mathrm{C}}=\frac{1}{2}\|\boldsymbol{e}_{i}\|^{2}
\\
\boldsymbol{e}_{j}^{\mathrm{T}}\boldsymbol{e}_{\mathrm{C}}=\frac{1}{2}\|\boldsymbol{e}_{j}\|^{2}
\\
\boldsymbol{e}_{k}^{\mathrm{T}}\boldsymbol{e}_{\mathrm{C}}=\frac{1}{2}\|\boldsymbol{e}_{k}\|^{2}
}

これらは，次のようにまとめられます．

\begin{bmatrix} \boldsymbol{e}_{i}^{\mathrm{T}} \\ \boldsymbol{e}_{j}^{\mathrm{T}} \\ \boldsymbol{e}_{k}^{\mathrm{T}} \end{bmatrix}
\boldsymbol{e}_{\mathrm{C}}=
\frac{1}{2}
\begin{bmatrix} \|\boldsymbol{e}_{i}\|^{2} \\ \|\boldsymbol{e}_{j}\|^{2} \\ \|\boldsymbol{e}_{k}\|^{2} \end{bmatrix}

さらにこれは，$\boldsymbol{e}_{\mathrm{C}}$について次のように直接的に解けます．

\displaylines{
\boldsymbol{e}_{\mathrm{C}}=
\frac{1}{2}
\begin{bmatrix} \boldsymbol{e}_{i}^{\mathrm{T}} \\ \boldsymbol{e}_{j}^{\mathrm{T}} \\ \boldsymbol{e}_{k}^{\mathrm{T}} \end{bmatrix}^{-1}
\begin{bmatrix} \|\boldsymbol{e}_{i}\|^{2} \\ \|\boldsymbol{e}_{j}\|^{2} \\ \|\boldsymbol{e}_{k}\|^{2} \end{bmatrix}
=
\frac{1}{2\left[\boldsymbol{e}_{i}~\boldsymbol{e}_{j}~\boldsymbol{e}_{k}\right]}
\begin{bmatrix}
\boldsymbol{e}_{j}\times\boldsymbol{e}_{k} & \boldsymbol{e}_{k}\times\boldsymbol{e}_{i} & \boldsymbol{e}_{i}\times\boldsymbol{e}_{j}
\end{bmatrix}
\begin{bmatrix} \|\boldsymbol{e}_{i}\|^{2} \\ \|\boldsymbol{e}_{j}\|^{2} \\ \|\boldsymbol{e}_{k}\|^{2} \end{bmatrix}
\\
=
\frac{
\|\boldsymbol{e}_{i}\|^{2}\boldsymbol{e}_{j}\times\boldsymbol{e}_{k}
+\|\boldsymbol{e}_{j}\|^{2}\boldsymbol{e}_{k}\times\boldsymbol{e}_{i}
+\|\boldsymbol{e}_{k}\|^{2}\boldsymbol{e}_{i}\times\boldsymbol{e}_{j}
}
{2\left[\boldsymbol{e}_{i}~\boldsymbol{e}_{j}~\boldsymbol{e}_{k}\right]}
}

ただし，$\left[\boldsymbol{e}_{i}~\boldsymbol{e}_{j}~\boldsymbol{e}_{k}\right]=\boldsymbol{e}_{i}^{\mathrm{T}}(\boldsymbol{e}_{j}\times\boldsymbol{e}_{k})=\boldsymbol{e}_{j}^{\mathrm{T}}(\boldsymbol{e}_{k}\times\boldsymbol{e}_{i})=\boldsymbol{e}_{k}^{\mathrm{T}}(\boldsymbol{e}_{i}\times\boldsymbol{e}_{j})$はスカラー三重積と呼ばれます．中心座標は$\boldsymbol{p}_{\mathrm{C}}=\boldsymbol{p}_{l}+\boldsymbol{e}_{\mathrm{C}}$となります．
半径は，前のケースと同様に$r=\|\boldsymbol{p}_{i}-\boldsymbol{p}_{\mathrm{C}}\|$と求まります．

点の数N≧5のときの最小包含球

$N\geq 5$のときは，次のような再帰的関数を定義します．ただし，$\mathcal{P}_{n}=\left\{\boldsymbol{p}_{i}|i=1\sim n\right\}$とおいています（$\mathcal{P}_{n}=\mathcal{P}_{n-1}\cup\left\{\boldsymbol{p}_{n}\right\}$，$\mathcal{P}=\mathcal{P}_{N}$です）．

image.png
image.png

ただし，$\mathcal{B}(\boldsymbol{p},r)$は中心座標$\boldsymbol{p}$および半径$r$の球，$\left|\mathcal{R}\right|$は$\mathcal{R}$の要素数をそれぞれ意味します．

それぞれの関数は

  • BoundingBall($\mathcal{P}_{n}$)…点群$\mathcal{P}_{n}$の最小包含球を求める関数
  • BoundingBallSub($\mathcal{P}_{n},\mathcal{R}$)…二つの点群$\mathcal{P}_{n}$と$\mathcal{R}$全て（つまり$\mathcal{P}_{n}\cup\mathcal{R}$）の最小包含球を求める関数
  • BoundingBallDirect($\mathcal{R}$)…点群$\mathcal{R}$（高々4個）の最小包含球を直接的に求める関数（前節に書いた計算を行う）

という役割を持ちます．

BoundingBallSub($\mathcal{P}_{n},\mathcal{R}$)が鍵で，点群$\mathcal{R}$は，$\mathcal{P}_{n}$と$\mathcal{R}$の最小包含球を同時に考えたときに必ず球の表面に含まれる点群になっています．

BoundingBall($\mathcal{P}_{n}$)の中身から説明します．再帰関数ですので，$\mathcal{P}_{n}$から$\boldsymbol{p}_{n}$だけを取り出した$\mathcal{P}_{n-1}$の最小包含球$\mathcal{BB}(\mathcal{P}_{n-1})=$BoundingBall($\mathcal{P}_{n-1}$)が既に求まっているとしましょう．$\boldsymbol{p}_{n}\in\mathcal{BB}(\mathcal{P}_{n-1})$ならば明らかに$\mathcal{BB}(\mathcal{P}_{n})=\mathcal{BB}(\mathcal{P}_{n-1})$です．$\boldsymbol{p}_{n}\notin\mathcal{BB}(\mathcal{P}_{n-1})$のとき，$\boldsymbol{p}_{n}$は必ず$\mathcal{BB}(\mathcal{P}_{n})$の表面上の点となるはずです．別の言い方をすれば，このとき $\boldsymbol{p}_{n}$は$\mathcal{BB}(\mathcal{P}_{n})$の中心からの最遠点ということになります．そこで，$\mathcal{R}=\{\boldsymbol{p}_{n}\}$とおいてBoundingBallSub($\mathcal{P}_{n-1},\mathcal{R}$)を呼びます．

BoundingBallSub($\mathcal{P}_{n},\mathcal{R}$)も再帰関数ですので，まず$\mathcal{BB}(\mathcal{P}_{n-1}\cup\mathcal{R})=$BoundingBallSub($\mathcal{P}_{n-1},\mathcal{R}$)が既に求まっているとしましょう．$\boldsymbol{p}_{n}\in\mathcal{BB}(\mathcal{P}_{n-1}\cup\mathcal{R})$ならば明らかに$\mathcal{BB}(\mathcal{P}_{n}\cup\mathcal{R})=\mathcal{BB}(\mathcal{P}_{n-1}\cup\mathcal{R})$です．$\boldsymbol{p}_{n}\notin\mathcal{BB}(\mathcal{P}_{n-1}\cup\mathcal{R})$のとき，$\boldsymbol{p}_{n}$は必ず$\mathcal{BB}(\mathcal{P}_{n}\cup\mathcal{R})$の表面上の点となるはずです．$\mathcal{BB}(\mathcal{P}_{n})$の中心からの最遠点は既に$\mathcal{R}$に含まれているので，このときの$\boldsymbol{p}_{n}$は $\mathcal{R}$に含まれる点を除く$\mathcal{BB}(\mathcal{P}_{n})$の中心からの最遠点ということになります．そこでこれを$\mathcal{R}$に追加し，改めて$\mathcal{BB}(\mathcal{P}_{n-1}\cup\mathcal{R})$を求めます．

命題$\star$より， 中心から遠い順に4つ点が見つかれば，それら4点の最小包含球は$\mathcal{BB}(\mathcal{P}_{n}\cup\mathcal{R})$に一致すると言えます．もちろん，4つ見つかる前に点が枯渇してしまった（つまり$\mathcal{P}_{n}=\emptyset$）ならば$\mathcal{BB}(\mathcal{P}_{n}\cup\mathcal{R})=\mathcal{BB}(\mathcal{R})$です．これが上記のアルゴリズムの根拠になっています．

アルゴリズムの計算量オーダー

前節で紹介したアルゴリズムの計算量オーダーを調べてみましょう．BoundingBall$(\mathcal{P}_{n})$の計算量を$c(n)$，BoundingBallSub$(\mathcal{P}_{n-1},\mathcal{R})$の計算量を$c_{\mathrm{s}}(n)$とそれぞれおくと，

\displaylines{
c(n)=O(1)+c(n-1)+\pi(\boldsymbol{p}_{n}\notin\mathcal{BB}(\mathcal{P}_{n-1}))c_{\mathrm{s}}(n-1) \\
c_{\mathrm{s}}(n)=O(1)+c_{\mathrm{s}}(n-1)+\pi(\boldsymbol{p}_{n}\notin\mathcal{BB}(\mathcal{P}_{n-1}\cup\mathcal{R}))c_{\mathrm{s}}(n-1)
}

ただし，$\pi(\boldsymbol{p}\notin\mathcal{BB}(\mathcal{Q}))$は，点$\boldsymbol{p}$が$\mathcal{BB}(\mathcal{Q})$に含まれない確率です．$\pi(\boldsymbol{p}_{n}\notin\mathcal{BB}(\mathcal{P}_{n-1}))$が少なくとも$4/n$以下であることは，直感的に分かります．$\pi(\boldsymbol{p}_{n}\notin\mathcal{BB}(\mathcal{P}_{n-1}\cup\mathcal{R}))$および$c_{\mathrm{s}}(n-1)$は，$\mathcal{R}$の要素数によって異なります．すなわち命題$\star$により$\mathcal{R}$の要素数は高々4ですが，

  1. $\mathcal{R}$の要素数が4のとき，$\pi(\boldsymbol{p}_{n}\notin\mathcal{BB}(\mathcal{P}_{n-1}\cup\mathcal{R}))=0$かつ$c_{\mathrm{s}}(n-1)=0$
  2. $\mathcal{R}$の要素数が3のとき，$\pi(\boldsymbol{p}_{n}\notin\mathcal{BB}(\mathcal{P}_{n-1}\cup\mathcal{R}))=1/n$かつ$c_{\mathrm{s}}(n-1)=O(1)$
  3. $\mathcal{R}$の要素数が2のとき，$\pi(\boldsymbol{p}_{n}\notin\mathcal{BB}(\mathcal{P}_{n-1}\cup\mathcal{R}))=2/n$かつ$c_{\mathrm{s}}(n-1)=O(1)$
  4. $\mathcal{R}$の要素数が1のとき，$\pi(\boldsymbol{p}_{n}\notin\mathcal{BB}(\mathcal{P}_{n-1}\cup\mathcal{R}))=3/n$かつ$c_{\mathrm{s}}(n-1)=O(1)$

なので結局，$c_{\mathrm{s}}(n)$は$\mathcal{R}$の要素数によらず$O(1)$であると分かります．したがって，上記漸化式より$c(n)=O(n)$となります．

Zeoでの実装

Zeoでは，zeo_bv3d_bball.cに最小包含球を求めるアルゴリズムを実装しています．関数は

int zBBall3DPL(zSphere3D *bb, zVec3DList *pl, zVec3D **vp);
int zBBall3D(zSphere3D *bb, zVec3DArray *pa, zVec3D **vp);

で，どちらも第1引数にアドレスが与えられたzSphere3Dインスタンスを，求まった最小包含球とします．第2引数で点群を与えます．前者は点群をzVec3DListで，後者はzVec3DArrayでそれぞれ与えるようになっています．第3引数はzVec3Dポインタへのポインタとなっており，これがNULLでない場合，最小包含球の表面に乗る点（高々4個）へのポインタ配列が格納されます．返り値はいずれもその表面に乗っている点の数です．

1
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
1
0