はじめに

点群$\mathcal{P}=\{\mathrm{P}_{i}|i=1\sim N\}$に含まれる$N$個の点を全て内部および表面上に含む球のうち，半径が最小のものを最小包含球(Bounding Ball, Smallest Enclosing Ball)と呼びます．最も単純な包含立体の一つであり，領域判定や干渉判定などのいろいろな場面で有用なものです．

以下，中心点が$\mathrm{P}_{\mathrm{C}}(\boldsymbol{p}_{\mathrm{C}})$かつ半径が$r$である球を$\mathcal{B}(\boldsymbol{p}_{\mathrm{C}},r)$，点群$\mathcal{P}$の最小包含球を$\mathcal{BB}(\mathcal{P})=\mathcal{B}(\boldsymbol{p}_{\mathrm{BB}},r_{\mathrm{BB}})$と表すことにします．$\boldsymbol{p}_{\mathrm{BB}}$と$r_{\mathrm{BB}}$を求める計算は，次のように最適化問題として表現できます．

\left\{\boldsymbol{p}_{\mathrm{BB}},r_{\mathrm{BB}}\right\}=\mathop{\mathrm{arg~min}}_{\boldsymbol{p}_{\mathrm{C}},r}\left\{r
\left|\|\boldsymbol{p}-\boldsymbol{p}_{\mathrm{C}}\|\leq r, \forall\boldsymbol{p}\in\mathcal{P}\right.
\right\}

これを解く効率の良いアルゴリズムは，幾つか知られています．本記事ではその中でもエレガントな，次の論文で提案されているものを解説します．

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.

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

ただし後に書くように，実用上は再帰呼び出しの深さが問題になるので，再帰関数を用いずに実装する方法が必要です．本記事ではこれも示します．¹

点の数が4以下のときの最小包含球

N≦4ならば最小包含球は直接求まる

点の数$N=1$のときは，$\mathcal{BB}(\mathcal{P})$は$\mathcal{P}$に含まれるただ一つの点を中心とする半径$0$の球となります．すなわち$\boldsymbol{p}_{\mathrm{BB}}=\boldsymbol{p}_{1}$，$r_{\mathrm{BB}}=0$です．

最小包含球が実用上の意味を持つのは$N\geq 2$の場合で，このとき$\mathcal{BB}(\mathcal{P})$は，次の3つの条件のうちいずれかを必ず満たします．

$\mathcal{P}$に含まれるある2点を直径の両端に持つ．
$\mathcal{P}$に含まれるある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$通りに絞られるので，その中から最小半径のものを選べば良いわけです．

具体的にそれぞれの球を求めてみましょう．

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点を頂点とする三角形の外接円を最大断面とする球

$\mathrm{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={\displaystyle \sum_{m=i,j,k}\frac{\|\boldsymbol{p}_{m}-\boldsymbol{p}_{\mathrm{C}}\|}{3}}$とします．

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

$\mathrm{P}_{\mathrm{C}}$は四面体の外心となります．これは，四面体の全ての辺の垂直二等分面の交点です．四面体の辺は6本ありますが，実際は全ての辺について考える必要は無く，ある一点を共有する3辺のみ考えれば十分です（四面体の外心の性質より，残り3辺の垂直二等分面も同一点で交わるため）．ここでは点$\mathrm{P}_{l}$を共有する3辺$\mathrm{P}_{i}\mathrm{P}_{l}$，$\mathrm{P}_{j}\mathrm{P}_{l}$，$\mathrm{P}_{k}\mathrm{P}_{l}$を考えましょう．上図はちょっと分かりづらいですが，3つの垂直二等分面による四面体の断面（三角形）の交点が，求めたい$\mathrm{P}_{\mathrm{C}}$であることを表しています．

$\mathrm{P}_{\mathrm{C}}$が辺$\mathrm{P}_{i}\mathrm{P}_{l}$の垂直二等分面上の点であるということは，$\mathrm{P}_{\mathrm{C}}$を辺$\mathrm{P}_{i}\mathrm{P}_{l}$上に射影したときに辺$\mathrm{P}_{i}\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}
}

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

\displaylines{
\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}_{i}$，$\boldsymbol{e}_{j}$，$\boldsymbol{e}_{k}$が互いに独立ならば，上式は$\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})$はスカラー三重積またはGrassmann積と呼ばれます．

中心座標は$\boldsymbol{p}_{\mathrm{C}}=\boldsymbol{p}_{l}+\boldsymbol{e}_{\mathrm{C}}$となります．半径は，$r=\|\boldsymbol{e}_{\mathrm{C}}\|$と求まりますが，前のケースと同様に$r={\displaystyle \sum_{m=i,j,k,l}\frac{\|\boldsymbol{p}_{m}-\boldsymbol{p}_{\mathrm{C}}\|}{4}}$としても良いでしょう．

10通りの可能性を総当りで調べる必要は無い

前述の通り，$N=4$ならば最小包含球は高々${}_{4}C_{2}+{}_{4}C_{3}+{}_{4}C_{4}=10$個の可能性のうちの一つとなるわけですが，最初から10通りの計算をして総当りで比較する必要はありません．

1つ目の可能性，「$\mathcal{P}$に含まれるある2点（$\mathrm{P}_{i}$，$\mathrm{P}_{j}$とします）を直径の両端に持つ」球が他の2点（$\mathrm{P}_{k}$，$\mathrm{P}_{l}$とします）を内部に含むとき，「$\mathrm{P}_{i}$，$\mathrm{P}_{k}$を直径の両端に持つ球」「$\mathrm{P}_{i}$，$\mathrm{P}_{l}$を直径の両端に持つ球」「$\mathrm{P}_{k}$，$\mathrm{P}_{l}$を直径の両端に持つ球」はいずれも$\mathrm{P}_{j}$を包含しません．
また，この球の半径は，「$\mathrm{P}_{i}$，$\mathrm{P}_{j}$，$\mathrm{P}_{k}$を通る円を最大断面とする球」「$\mathrm{P}_{i}$，$\mathrm{P}_{j}$，$\mathrm{P}_{l}$を通る円を最大断面とする球」「$\mathrm{P}_{j}$，$\mathrm{P}_{k}$，$\mathrm{P}_{l}$を通る円を最大断面とする球」のいずれの半径よりも小さいことが保証されます（証明は省略します）．

同様のことが2つ目の可能性についても言えます．すなわち，「$\mathcal{P}$に含まれるある3点（$\mathrm{P}_{i}$，$\mathrm{P}_{j}$，$\mathrm{P}_{k}$とします）を通る円を最大断面とする」球が残りの1点（$\mathrm{P}_{l}$とします）を内部に含むとき，「$\mathrm{P}_{i}$，$\mathrm{P}_{j}$，$\mathrm{P}_{l}$を通る円を最大断面とする球」は$\mathrm{P}_{k}$を包含しません．
また，この球の半径は，「$\mathrm{P}_{i}$，$\mathrm{P}_{j}$，$\mathrm{P}_{k}$，$\mathrm{P}_{l}$を通る球」の半径よりも小さいことが保証されます（こちらも証明は省略します）．

以上より，

まず4点から2点を順番に選んでいき，それを直径の両端とする球が残り2点を包含するものが見つかったら，その時点でそれが最小包含球と分かる
上記のような球が存在しない場合，4点から3点を順番に選んでいき，それらを通る円を最大断面とする球が残り1点を包含するものが見つかったら，その時点でそれが最小包含球と分かる
上記のような球が存在しない場合，4点を通る球が最小包含球と分かる

ということが言えます．

この部分のアルゴリズムは，下記のAlgorithm 1のようになります．

ただし，$|\mathcal{P}|$は$\mathcal{P}$の要素数を意味します．また，BoundingBall2$(\boldsymbol{p}_{i},\boldsymbol{p}_{j})$，BoundingBall3$(\boldsymbol{p}_{i},\boldsymbol{p}_{j},\boldsymbol{p}_{k})$，BoundingBall4$(\boldsymbol{p}_{i},\boldsymbol{p}_{j},\boldsymbol{p}_{k},\boldsymbol{p}_{l})$はそれぞれ，上で説明した「$\boldsymbol{p}_{i}$，$\boldsymbol{p}_{j}$を直径の両端とする球」「$\boldsymbol{p}_{i}$，$\boldsymbol{p}_{j}$，$\boldsymbol{p}_{k}$を通る円を最大断面とする球」「$\boldsymbol{p}_{i}$，$\boldsymbol{p}_{j}$，$\boldsymbol{p}_{k}$，$\boldsymbol{p}_{l}$を通る球」を求める計算に対応します．

点の数が5以上のときの最小包含球

N≧5ならば最小包含球は再帰的に求まる

$N\geq 5$のときは，次のように考えます．

まず，$N$個のうち$1\sim n$番目（$n\leq N$）の点のみを集めた部分点群を$\mathcal{P}_{n}=\{\mathrm{P}_{i}|i=1\sim n\}$と表すことにします．今，$\mathcal{P}_{n-1}$の最小包含球$\mathcal{BB}(\mathcal{P}_{n-1})$が求まっているとして，

$\mathrm{P}_{n}$が$\mathcal{BB}(\mathcal{P}_{n-1})$内部（境界含む）の点ならば，明らかに$\mathcal{BB}(\mathcal{P}_{n})=\mathcal{BB}(\mathcal{P}_{n-1})$
$\mathrm{P}_{n}$が$\mathcal{BB}(\mathcal{P}_{n-1})$外部の点ならば，$\mathrm{P}_{n}$は明らかに$\mathcal{BB}(\mathcal{P}_{n})$の境界上の点

が言えます．

ここで，境界上の点の候補のみを選り分けた点群$\mathcal{R}$を作り，$\mathcal{P}_{n}$と$\mathcal{R}$を併せた点群の最小包含球$\mathcal{BB}(\mathcal{P}_{n}\cup\mathcal{R})$を求めることを考えましょう．前節で説明したことより，$\left|\mathcal{R}\right|$は高々4ですので，もし，$\left|\mathcal{R}\right|=4$である，あるいは$\mathcal{P}_{n}=\emptyset$であるならば，$\mathcal{BB}(\mathcal{P}_{n}\cup\mathcal{R})=\mathcal{BB}(\mathcal{R})$のはずです．$\mathcal{BB}(\mathcal{R})$は前節で示した関数BoundingBallDirectで求めることが出来ます．

今，$\mathcal{P}_{n-1}\cup\mathcal{R}$の最小包含球$\mathcal{BB}(\mathcal{P}_{n-1}\cup\mathcal{R})$が求まっているとして，

$\mathrm{P}_{n}$が$\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})$
$\mathrm{P}_{n}$が$\mathcal{BB}(\mathcal{P}_{n-1}\cup\mathcal{R})$外部の点ならば，$\mathrm{P}_{n}$は明らかに$\mathcal{BB}(\mathcal{P}_{n}\cup\mathcal{R})$境界上の点，すなわち$\mathcal{BB}(\mathcal{P}_{n}\cup\mathcal{R})=\mathcal{BB}(\mathcal{P}_{n-1}\cup(\mathcal{R}\cup\left\{\mathrm{P}_{n}\right\}))$

が言えます．以上のことから，次のAlgorithm 2，Algorithm 3が導けます．

これら二つの再帰的手続きはよく似ていますが，Algorithm 2は$|\mathcal{P}_{n}|$が4以下にまで減ったら最小包含球の候補を更新するのに対し，Algorithm 3は$|\mathcal{R}|$が4まで増えたら最小包含球の候補を更新する，という点で性質が異なります．

ここで次のような思考実験をしてみましょう．仮に$\mathcal{P}_{n}=\mathcal{P}$，$\mathcal{R}=\emptyset$としてBoundingBallInternalを呼ぶと，まず$\mathcal{P}$から一つずつ点を減らして空集合になったところでBoundingBallDirectが呼ばれます．このときは$\mathcal{R}$も空集合ですので，最初の最小包含球の候補$\mathcal{B}_{0}$は数学的に意味を持たないものとなります．しかしこの$\mathcal{B}_{0}$に対して点の内外判定処理が定義されている（どのような点もこの球の外部の点と判定される）ならば，$\boldsymbol{p}_{1}\notin\mathcal{B}_{0}$ですので$\mathcal{R}\leftarrow\left\{\mathrm{P}_{1}\right\}$と更新され，この後は意味のある計算がなされることになります．

すなわち，BoundingBall$(\mathcal{P}_{n})$は実質的にBoundingBallInternal$(\mathcal{P}_{n},\emptyset)$と等価な計算であることが分かります．したがって，Algorithm 2はAlgorithm 3に統合できる（最小包含球はBoundingBallInternal$(\mathcal{P},\emptyset)$より求まる）わけです．

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

BoundingBallInternal$(\mathcal{P}_{n})$の計算量オーダーを$c(n)$とおくと，

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

が言えます．ただし，$\pi(\boldsymbol{p}\notin\mathcal{B})$は点$\boldsymbol{p}$が球$\mathcal{B}$に含まれない確率です．

$\pi(\boldsymbol{p}_{n}\notin\mathcal{BB}(\mathcal{P}_{n-1}\cup\mathcal{R}))$は，$\mathcal{R}$の要素数によって異なります．すなわち

$\left|\mathcal{R}\right|=4$のとき，$\pi(\boldsymbol{p}_{n}\notin\mathcal{BB}(\mathcal{P}_{n-1}\cup\mathcal{R}))=0$
$\left|\mathcal{R}\right|=3$のとき，$\pi(\boldsymbol{p}_{n}\notin\mathcal{BB}(\mathcal{P}_{n-1}\cup\mathcal{R}))=1/n$
$\left|\mathcal{R}\right|=2$のとき，$\pi(\boldsymbol{p}_{n}\notin\mathcal{BB}(\mathcal{P}_{n-1}\cup\mathcal{R}))=2/n$
$\left|\mathcal{R}\right|=1$のとき，$\pi(\boldsymbol{p}_{n}\notin\mathcal{BB}(\mathcal{P}_{n-1}\cup\mathcal{R}))=3/n$

となります．

ここで，

$kO(n)$=$O(kn)=O(n)$（$k$は$n$によらない定数）
$O(n)/n$=$O(n)-1$
$O(1)-1$=$O(1)$

ですから，上記漸化式は結局$\left|\mathcal{R}\right|$によらず

c(n)=O(1)+c(n-1)

となります．したがって，$c(N)=NO(1)=O(N)$であると分かります．

実用上の問題と実装

再帰的呼び出しを用いない

上記のアルゴリズムはとてもエレガントで効率が良いのですが，点を一つ取り出すごとに1回再帰的呼び出しをするので，$N$が数千〜数万になると大抵のシステムではスタック上限に引っかかってしまいます．

そこで，再帰的呼び出しを使う代わりに，自前のスタックをループの中でpush/popすることで，等価な処理を行わせることを考えます．

BoundingBallInternalは非常にシンプルな構造をしており，次のように機能分解出来ます．

$\mathcal{P}_{n}=\emptyset$または$|\mathcal{R}|=4$のときに最小包含球の候補を直接求める部分
$\mathcal{P}_{n}$から$\mathrm{P}_{n}$を取り出し，再帰的に$\mathcal{BB}(\mathcal{P}_{n-1}\cup\mathcal{R})$を求める部分
$\mathrm{P}_{n}$が最新の最小包含球候補の内部にあるか否かを判定し，
1. 内部にあれば，最新候補をそのまま返す部分
2. 外部にあれば，再帰的に$\mathcal{BB}(\mathcal{P}_{n-1}\cup(\mathcal{R}\cup\left\{\mathrm{P}_{n}\right\}))$を求める部分

上記のうち「再帰的に」と付いている部分は，点群を更新する（$\mathcal{P}_{n}$から$\mathrm{P}_{n}$を除く，$\mathcal{R}$に$\mathrm{P}_{n}$を加える）処理をしてから頭に戻る，と読み替えられます．また，この処理の後は点群を復元する（$\mathcal{P}_{n-1}$に$\mathrm{P}_{n}$を戻す，$\mathcal{R}$から$\mathrm{P}_{n}$を除く）処理が求められます．したがって，これらの処理（アクション）のうち頭に戻ってから次にすべきは何（どれ）かを正しく記憶しておけば，ループでも等価な処理が行えるわけです．

アクションは次の3つにまとめられます．

SAVE $\mathrm{P}_{n}$を一時的に$\mathcal{P}$から除く処理（最初の再帰呼び出し直前）
CHECK $\mathrm{P}_{n}$の内外判定を行う処理（最初の再帰呼出し直後）
REVERT $\mathrm{P}_{n}$を$\mathcal{P}$に戻し，同時に$\mathcal{R}$から除く処理（最後の再帰呼び出し直後）

これらを用いれば，BoundingBallInternal$(\mathcal{P},\mathcal{R})$と等価なアルゴリズムが次のAlgorithm 4のように書けます．

ただし，

NextAction$(\mathcal{S})$は，$\mathcal{S}$の先頭アクションを参照する（$\mathcal{S}=\emptyset$ならばSAVEを返す）処理
ChangeNextAction$(\mathcal{S},a)$は，$\mathcal{S}$の先頭アクションを$a$に上書きする処理
Save$(\mathcal{P},\bar{\mathcal{P}},\mathcal{S},\bar{\mathcal{S}},a)$は，$\mathcal{P}$からpopした点$\mathrm{P}$を$\bar{\mathcal{P}}$にpushし，
同時に$\mathcal{S}$からpopしたアクションを$a\in\{$SAVE,CHECK,REVERT$\}$に上書きして$\bar{\mathcal{S}}$にpushする処理
Revert$(\mathcal{P},\bar{\mathcal{P}},\mathcal{S},\bar{\mathcal{S}})$は，$\bar{\mathcal{P}}$からpopした点$\mathrm{P}$を$\mathcal{P}$にpush（$\bar{\mathcal{P}}=\emptyset$ならばnilを返す）し，
同時に$\bar{\mathcal{S}}$からpopしたアクションを$\mathcal{S}$にpushする処理
PushExtreme$(\mathrm{P},\mathcal{R})$は，$\mathrm{P}$を$\mathcal{R}$にpushする処理
PopExtreme$(\mathcal{R})$は，$\mathcal{R}$の先頭要素をpopする処理

です．$\mathcal{P}$，$\bar{\mathcal{P}}$，$\mathcal{R}$，$\mathcal{S}$，$\bar{\mathcal{S}}$は全てスタックであることを前提としています．

Algorithm 4は基本的に，$\mathcal{P}=\emptyset$または$|\mathcal{R}|=4$のときには最小包含球の候補を更新する，そうでないときには次のアクションに基づいてスタックを更新する，という構造ですが，後者のうちREVERTだけは，最小包含球の更新よりも前にPopExtremeを正しく発動する必要があるため，先に行っています．

元のアルゴリズムが低確率で失敗する問題について

筆者も理由を理解出来ていないのですが，元のアルゴリズムは点をSAVE/REVERTする順番によっては失敗することがあるようです（経験的に，発生確率は1%未満です）．これについては，Welzl氏自身が共著となっている関連論文があります．

Jiri Matousek, Micha Sharir, and Emo Welzl, A subexponential bound for linear programming, Algorithmica, Vol. 16, No. 4-5, pp. 498--516, 1996.

これも著者自身の公開しているものが読めます．

実用的には，上記アルゴリズムで最小包含球を求めた後に，全点がそれの内部または境界上に含まれているか改めてチェックし，含まれていないものがあればそれを点群先頭に移して再度計算し直す，という方法で回避できます．このチェックも$O(N)$ですので，全体の計算量オーダーは変わらず$O(N)$のままです．

Zeoでの実装

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

int zVec3DDataBoundingBallRecursive(zVec3DData *data, zSphere3D *bb, zVec3D **vp);
int zVec3DDataBoundingBall(zVec3DData *pointdata, zSphere3D *bb, zVec3D **vp);

で，前者が再帰的呼び出しを使ってオリジナルアルゴリズムをそのまま実装したもの，後者が再帰的呼び出しを使わず実装したものです．どちらも第1引数に与えられた点群dataの最小包含球を求め，第2引数にアドレスが与えられたzSphere3Dインスタンスに結果を代入します．第3引数はzVec3Dポインタへのポインタとなっており，これがNULLでない場合，境界（最小包含球の表面）上に乗る点（高々4個）へのポインタ配列が格納されます（したがって，元の配列サイズは4以上である必要があります）．返り値はいずれもその表面に乗っている点の数です．

2次元版最小包含球（最小包含円）も，同様にzeo_bv2d_boundingball.cに

int zVec2DDataBoundingBallRecursive(zVec2DData *data, zDisk2D *bb, zVec2D **vp);
int zVec2DDataBoundingBall(zVec2DData *pointdata, zDisk2D *bb, zVec2D **vp);

として実装しています．引数の意味は3次元版と同じです．
3次元の場合は球構造体zSphere3D，2次元の場合は円構造体zDisk2Dをそれぞれ使っていることに注意して下さい．

これらの実体は，zeo_bvxd_boundingball.hに定義したマクロによって共通化されています（Zeoでは同様の書き方をしているものが幾つかあります．我ながらエキセントリックな書き方だと思います）．

同じ問題意識を持って書かれたこちらの記事は大変貴重で，参考にさせて頂きました． ↩