Welzlの最小球アルゴリズム（理屈編）

　沢山の点（点群）を囲う最も小さい円（3Dなら球）を求めたいという要求は色々な局面で見られます。例えば距離センサーが返す大量の点群を覆う境界球を知りたいとか、ゲーム内にあるメッシュの衝突判定の第1カリング（詳細な衝突を行う前の粗いカリング）のために境界球を求めたい場合など。

　この割と古典的な問題に対して様々なアルゴリズムが提案されてきました。その中でWelzlが1991年に発表した論文「Smallest enclosing disks (balls and ellipsoids)」で紹介されているアルゴリズムは、シンプルかつ高い効率で最小境界円（球）を求められます。論文は検索すると多分出てきますので、見たい方は調べてみて下さい。

　今回はこのアルゴリズムの理屈を論文を元にしながら見て行こうと思います。

　ただ先に注意しておくことがありまして、Welzlのアルゴリズムは再帰関数を使用している為、論文に掲載されている筋道のまま実装するとすぐにスタックオーバーフローを起こしてしまいます。つまり実用化するには再帰関数を展開して非再帰化する必要があります。それについては次回考えてみます。

1点を取り除き、残りの最小円を求める

　まず最初にminiSphereという関数の仕様を定めます。miniSphere関数は引数に渡された点群Pの最小円（球）を求めます。つまり今回の目的をそのまま果たしてくれる関数です（注：論文内ではこの関数名をminidiskとしていますが、ここではより一般的な意味でminiSphereとします）。関数に渡される点群を点群Pと呼ぶ事にします。

　関数の最初では、点群Pから点を一つ取り除きます。その除去点以外の点群P'を囲う最小円を算出します。この1点を取り除くというのが論文の肝であり上手い考え方なんです：

　上図のオレンジの点が点群P'で、その最小円（暫定最小円）が求められています。この時、青い点に対して左と右の2つのパターンが考えられます。左の場合、青い点は暫定最小円に含まれています。つまりこの1点は暫定最小円に何も関与しない、あっても無くても変わらないので、暫定最小円が点群Pの最小円と確定します。

　一方右側のように除いた1点が暫定最小円の外にある場合、暫定最小円は点群Pの最小円ではないとわかります。よって青い点を含んだ最小円を改めて求め直す必要があります。しかしこの時「除いた1点は点群Pの最小円の円周上に絶対ある！」という事が確定します。これが超大事です：

　感覚的には除いた1点が対象点群の大外に位置する事になるわけです。外側に出っ張っている点というイメージでしょうか。ここから右のパターンだった場合「この青い点を最小円の円周（境界）に含む」という条件を満たす最小円を求めるフェーズに移行します。そういう条件付きの最小円を求める事ができる関数をb_miniSphre関数とする事にしましょう。

　先の除去点のように円周上にある事が確定している点群を点群Rとします。除去点が暫定最小円に含まれない場合、その除去点は点群Rに追加される事になるわけです。そして点群P'と点群Rの両方を用いて改めてb_miniSphere関数に渡して最小円（球）を求めてもらいます。

　ここまでを疑似的なコードで書くとこういう感じになります：

Sphere miniSphere( List P ) {

	// 点群Pから1点を取り除く
	Point removeP = P.removePoint();
	
	// 取り除いた点以外の点群P'の最小球を算出
	Sphere smallestSph = miniSphere( P );
	
	// 除いた点が暫定最小円の外にある場合点群Rに追加して
	// Rが円周上にある条件で最小球を求める
	if ( !smallestSph.contain( removeP ) ) {
		List R;
		R.add( removeP );
		smallestSph = b_miniSphere( P, R )
	}
	
	return smallestSph;
}

　上の疑似コードで実装がまだ不明なのは点群Rを使うb_miniSphere関数のみです。これさえあれば、最小円（球）は原理的に求まる事になりますね。次からはb_miniSphere関数について考えてみましょう。

やっぱり1点を除いて最小円を求める

　b_miniSphere関数は点群Pと最小円の円周（境界）上にある事が確定している点群Rの2つの点群を引数に取ります。

　b_miniSphere関数の実装はminiSphere関数と凄く似ています。点群Pから1点を取り除いて（下図の赤い点）その最小円を求めます。ただしminiSphere関数と違いその最小円は点群R（青い点）を円周上に含まないといけないという条件が追加されています。それを求めてくれるのはb_miniSphere関数自体ですよね。なのでここに再帰が発生するんです：

赤いのが取り除いた点。青い点は点群Rに所属する点で、最小円の周（境界）上になければならない（という条件）

　結果はやはり2パターンに分かれます。暫定最小円に除いた点（赤い点）が含まれていれば、最小円確定です。含まれていない場合、miniSphere関数で話した事と全く同様で、青い点と共にその赤点もまた最小円の円周上に含まれる事になります。ですから右のパターンの場合は赤い点も点群Rに追加し、再びb_miniSphere関数で最小円を求める事になります。

　この段階でのb_miniSphere関数を疑似コードで書いてみます：

Sphere b_miniSphere( List P, List R ) {

	// 点群Pから1点を取り除く
	Point removeP = P.removePoint();
	
	// 取り除いた点以外の点群の最小球を算出
	Sphere smallestSph = b_miniSphere( P, R );
	
	// 最小球に除いた点が含まれなければ点群Rに追加して
	// Rが境界にある条件で最小球を求める
	if ( !smallestS.contain( removeP ) ) {
		R.add( removeP );
		smallestSph = b_miniSphere( P, R )
	}
	
	return smallestSph;
}

　miniSphere関数とほとんど一緒です。違うのはb_miniSphere関数自体を再帰的に利用している点です。

　ただし、この段階のコードではb_miniSphere関数は完全には動きません。そりゃそうで、肝心の最小円を求める具体的な部分がどこにもないからです。b_miniSphereはいつ、どのように最小円を求めるのか？それを次に見て行きましょう。

b_miniSphere関数での最小円の算出

　b_miniSphere関数の冒頭に注目です。1点取り除いた点群P'と点群Rがセットになってまたb_miniSphere関数に放り込まれています。つまりb_miniSphere関数が呼ばれるたびに点群Pから点が一つずつ失われていき、底の底まで行くと最終的に点群P内の点がすべて無くなってしまう事になります。

　b_miniSphere関数の目的は引数の点群PとRから最小円（球）を求める事でした。その引数に渡された点群Pが空っぽの場合、残りの点群Rが円周上にある円を求めるしかありません。そう、ここが最小円算出のタイミングです。

　点群Pが空の状態で点群Rに含まれる点の数は、円の場合は0～3点、球の場合は0～4点です。円の場合は4パターン、3D球だと5パターンを検討すれば良い事になります。

　円のケースで具体的にどうするか列挙してみます。
　点群R内の点が0個の場合、まぁどうしようもないので原点中心で半径0の仮最小円を返します。アルゴリズムの初期段階でこの状況は必ず起きるのですが、元の点群Pに点が1つ以上あればこの仮の最小円が採択される事はありません。
　点が1つの場合はその点が中心で半径0が最小円です。
　点が2つの場合はその2点を結ぶ線分が直径になります（3Dの球でも同様）。なので中心点はその中点で、半径は中心点と2点のどちらかまでの距離となります。
　点が3つあった場合はその3点からなる三角形の外心が中心点になります。外心を求めるのはやや大変なので、ここだけちょっと掘り下げて説明します。

三角形の外心を3点の座標から求める

　三角形の外心とは3点からの距離が等しい点です。それすなわちその3点を含む円（外接円）の中心点に他ならないですよね。幾何学的には「各辺の垂直二等分線の交点」が外心になります：

　上図の青い点q01、q02がそれぞれ中点で、ベクトルv01、v02は各辺の垂線（法線）です。そのベクトルを持つ2直線の交点cが外心となります。各中点は、

と両端の座標の平均で求められます。法線ベクトルv01はp1-p0で求められるベクトルu01の成分を入れ替えると出てきます：

v02も同様です。

　さてここから2本の直線の交点cを求めるため、次の2つの式を立てます：

　上式はq01を通る垂直線上の点をtで表現する陽関数で、下は辺p0p2と平行なベクトルu02と括弧内のベクトルとが垂直関係になるという性質を表した陰関数です。あるtにおいてxは交点cになるはずです。そこで上式を下に代入して整理すると、

このように双方を満たすtを算出できます。ここから外心cは、

と3点の座標から直接計算できます。

　上の最終式の分母がゼロになる事がありえます。それは三角形が縮退している、すなわち点群Rに含まれる3点が1直線状に並んでいる場合です。これは例外として扱う必要があります。具体的には3点の内2点を選び、その2点を通る最小円を算出します。その円の中に残りの1点が含まれている物を採択します。

　という事で、b_miniSphere関数は大きく2つのフェーズに分かれる事になります。一つは点群Pが空っぽの時に最小円（球）を求めるフェーズ。もう一つは1点を取り除いた暫定最小円（球）が最小円になれるかどうかを判定するフェーズです。疑似コードを書いてみましょう：

Sphere b_miniSphere( List P, List R ) {

	// 点群Pが空の場合は最小円を算出する
	if ( P.isEmpty() ) {
		if ( R.size() == 0 ) return Sphere();
		if ( R.size() == 1 ) return Sphere( R[0], 0.0 );
		if ( R.size() == 2 ) return calcSmallest2( R );
		if ( R.size() == 3 ) return calcSmallest3( R );
		return Sphere();
	}
	
	// 点群Pから1点を取り除く
	Point removeP = P.removePoint();
	
	// 取り除いた点以外の点群の最小球を算出
	Sphere smallestSph = b_miniSphere( P, R );
	
	// 最小球に除いた点が含まれなければ点群Rに追加して
	// Rが境界にある条件で最小球を求める
	if ( !smallestSph.contain( removeP ) ) {
		R.add( removeP );
		smallestSph = b_miniSphere( P, R )
	}
	
	return smallestSph;
}

　これで動作するb_miniSphere関数が完成です。そして、頭で定義したminiSphere関数は、いまや次のようにb_miniSphere関数に置き換えられる事がわかります：

Sphere miniSphere( List P ) {
    // 点群Rを空にしてb_miniSphereを呼べばよい
	return b_miniSphere( P, List() );
}

これでWelzlの最小円のアルゴリズム完成です！

3次元の場合

　3次元の場合もアルゴリズムは基本変わりません。変わるのはm_miniSphere関数で最小球を求める所です：

Sphere b_miniSphere( List P, List R ) {

	// 点群Pが空の場合は最小円を算出する
	if ( P.isEmpty() ) {
		if ( R.size() == 0 ) return Sphere();
		if ( R.size() == 1 ) return Sphere( R[0], 0.0 );
		if ( R.size() == 2 ) return calcSmallest2( R );
		if ( R.size() == 3 ) return calcSmallest3( R );
		if ( R.size() == 3 ) return calcSmallest4( R );  // ここ
		return Sphere();
	}
	
	// 点群Pから1点を取り除く
	Point removeP = P.removePoint();
	
	// 取り除いた点以外の点群の最小球を算出
	Sphere smallestSph = b_miniSphere( P, R );
	
	// 最小球に除いた点が含まれなければ点群Rに追加して
	// Rが境界にある条件で最小球を求める
	if ( !smallestSph.contain( removeP ) ) {
		R.add( removeP );
		smallestSph = b_miniSphere( P, R )
	}
	
	return smallestSph;
}

　点群Rに点が4つある場合、4点を通る唯一の球を求める事になります。これは4点から3つ点を選び、三角形を2つ作ります。次にその三角形の外心、つまり三角形の外接球の中心点を求めます。球上の3点を結んでできる三角形の外心から垂線を伸ばすと、その垂線上には球の中心点が存在します。三角形が2つあれば、その垂線の交点が中心点になるわけです。

　ただしこれにもやはり例外が発生します。4点が同一平面上にあると球は定まらなくなってしまいます。その場合は3点を通る球を4パターン算出し、残りの1点が含まれる球を採用します。

　4点を通る球を求める方法を別記事にまとめましたのでよろしければご参照下さい：

そのまま実装するとスタックオーバーフロー

　さてWelzlの最小円（球）アルゴリズムはこれで動くわけですが、冒頭にも述べたようにこの再帰的なアルゴリズムをそのまま実装するとあっと言う間にスタックオーバーフローを起こしてしまいます。

　「スタックオーバーフローとは？」という方のために少しお話しします。関数を呼び出すとスタックメモリという一時メモリに引数とローカル変数を扱う領域が確保されます。これは関数を抜けると解放されるのですが、再帰関数のように入れ子関数がどんどんコールされ続けると、関数から抜ける事無くスタックにどんどん値が積まれ続けてしまいます。スタックメモリのサイズは上限がある為（通常数M～数十MB）、極端に深くまで関数が呼ばれるとスタックメモリが溢れてしまい、アプリケーションがクラッシュしてしまいます。これをスタックオーバーフローと言います。

　今時のPCならスタックメモリを多く設定しても問題ありませんが、点群数は不明であるわけで、潜在的なスタックオーバーフローはそれでは解決しません。

　そのため、Welzlのアルゴリズムは非再帰化しないと実用するにちょっと厳しいんです。そこで次回はWelzlのアルゴリズムの非再帰関数版を考えてみようと思います。