Day 154-Q4 | Welzl's Algorithm（最小包含円）

問題

平面上の $N$ 個の点が与えられる。すべての点を内部または境界に含む最小の円（最小包含円）の中心座標と半径を求めよ。

$S \subseteq D$ を「その円の境界上に乗っている点の集合」とすると、最小包含円は $|S| \le 3$ の点だけで一意に決まることが知られている。この性質を使い、点を1つずつ処理しながら「これまで処理した点をすべて含む最小円」を維持し、新しい点がその円の外にあれば境界点集合を更新して再構成する（Welzlのアルゴリズム）。

入力形式

N
x_1 y_1
...
x_N y_N

制約

$1 \le N \le 500$

$-10^4 \le x_i, y_i \le 10^4$

絶対誤差または相対誤差 $10^{-6}$ 以内で正解

入出力例

入力例1

出力例1

2.000000 1.500000 2.500000

概念図：直角三角形の斜辺=直径

ヒント（段階的開示）

ヒント1（方向性）

すべての3点の組み合わせで外接円を計算し他の全点を含むか確認する $O(N^4)$ の全探索でも答えは出るが、$N$ が大きいと間に合わない。「境界を決める点は高々3個」という事実を、点を後ろから確定させていく再帰に落とし込む。

ヒント2（アプローチ）

再帰関数 mec(i, boundary) を「点 $0,\dots,i-1$ と境界候補集合 boundary を必ず含む最小円」を返すものとして定義する。点 $i$ を追加したとき、$i-1$ 番目までの円に $i$ が入っていればそのまま、入っていなければ「$i$ は必ず境界上に乗る」ことが証明できるので boundary に $i$ を加えて再帰する。boundary が3点集まった時点でその3点の外接円が確定版として返る。

ヒント3（誘導）

def mec(i, boundary):
    if i == len(pts) or len(boundary) == 3:
        return circle_from(boundary)
    p = pts[i]
    c = mec(i + 1, boundary)
    if in_circle(c, p):
        return c
    return mec(i + 1, boundary + [p])

模範解答 (Python)

import sys
input = sys.stdin.readline
EPS = 1e-7


def circle_from_2(a, b):
    cx = (a[0] + b[0]) / 2
    cy = (a[1] + b[1]) / 2
    r = ((a[0] - b[0]) ** 2 + (a[1] - b[1]) ** 2) ** 0.5 / 2
    return (cx, cy, r)


def circle_from_3(a, b, c):
    ax, ay = a; bx, by = b; cx, cy = c
    d = 2 * (ax * (by - cy) + bx * (cy - ay) + cx * (ay - by))
    if abs(d) < 1e-12:
        cands = [circle_from_2(a, b), circle_from_2(a, c), circle_from_2(b, c)]
        best = None
        for cc in cands:
            if in_circle(cc, a) and in_circle(cc, b) and in_circle(cc, c):
                if best is None or cc[2] < best[2]:
                    best = cc
        return best
    ux = ((ax**2 + ay**2) * (by - cy) + (bx**2 + by**2) * (cy - ay) + (cx**2 + cy**2) * (ay - by)) / d
    uy = ((ax**2 + ay**2) * (cx - bx) + (bx**2 + by**2) * (ax - cx) + (cx**2 + cy**2) * (bx - ax)) / d
    r = ((ux - ax) ** 2 + (uy - ay) ** 2) ** 0.5
    return (ux, uy, r)


def in_circle(circle, p):
    cx, cy, r = circle
    return (p[0] - cx) ** 2 + (p[1] - cy) ** 2 <= r * r + EPS


def circle_from(boundary):
    if len(boundary) == 0:
        return (0.0, 0.0, 0.0)
    if len(boundary) == 1:
        return (boundary[0][0], boundary[0][1], 0.0)
    if len(boundary) == 2:
        return circle_from_2(boundary[0], boundary[1])
    return circle_from_3(*boundary)


def welzl(pts):
    sys.setrecursionlimit(5000)

    def mec(i, boundary):
        if i == len(pts) or len(boundary) == 3:
            return circle_from(boundary)
        p = pts[i]
        c = mec(i + 1, boundary)
        if c and in_circle(c, p):
            return c
        return mec(i + 1, boundary + [p])

    return mec(0, [])


def main():
    N = int(input())
    pts = [tuple(map(float, input().split())) for _ in range(N)]
    cx, cy, r = welzl(pts)
    print(f"{cx:.6f} {cy:.6f} {r:.6f}")


main()

Step-by-Step 解説

1幾何的な下見
入力例の3点 $(0,0),(4,0),(0,3)$ は直角三角形（直角は原点）をなし、斜辺 $(4,0)$–$(0,3)$ の長さは $\sqrt{16+9}=5$。直角三角形の外接円は斜辺を直径とする円になる（タレスの定理の逆）ので中心は斜辺の中点 $(2,1.5)$、半径は $2.5$。

24点すべてがこの円に収まるか検証
$(0,0),(4,0),(0,3)$ はいずれも中心から距離 $2.5$（境界上）。$(2,1)$ は距離 $0.5$（内部）。4点すべてが収まり、3点が既に境界上にあるためこれ以上小さくはできない。

3アルゴリズムの動きとの対応
最初の数点だけでは円が小さすぎて後の点が外側にはみ出すたびにboundaryに点が追加される。最終的に $(0,0),(4,0),(0,3)$ が境界集合として確定し、circle_from_3 がちょうど中心 $(2.0,1.5)$・半径 $2.5$ を返す。$(2,1)$ はどの段階でも既存の円の内部に収まるため境界集合には加わらない。

4出力誤差の許容
浮動小数点演算（除算・平方根）を含むため、厳密一致でなくても絶対誤差または相対誤差 $10^{-6}$ 以内であれば正解とする。

よくあるミス

ミス	原因	正しい書き方
3点が共線のとき外接円公式でゼロ除算する	分母d=0のケースを見落とす	dがほぼ0なら共線とみなし、最も離れた2点を直径とする円で残り1点を含むかチェックするフォールバックを用意する
ptsを毎回シャッフルしないと正しさ自体が壊れると誤解する	「乱択アルゴリズム=ランダム性が正しさの前提」という思い込み	シャッフルは期待計算量を保証するためのもので、正しさは処理順序によらず成り立つ
in_circleの判定に誤差マージンを入れず厳密比較する	浮動小数点誤差で境界上の点が「外」と誤判定される	r*r + EPSのように許容誤差を持たせて判定する
再帰の終了条件をlen(boundary)==3だけにする	全点処理し終わったときの終了条件を見落とす	i==len(pts) または len(boundary)==3 のどちらかで確定円を返す
半径0の円（1点のみ）のケースを特別扱いし忘れる	2点・3点前提のロジックのみ実装	circle_fromに0点・1点のケース（半径0）も明示的に用意する

次のステップ

発展: 3次元空間での最小包含球（Miniball問題）に拡張する。境界点集合は最大4点になる。
次回予告: Simulated Annealing（焼きなまし法）

自己評価

完全理解おおむね理解要復習

Day 154-Q4 — Welzl's Algorithm（最小包含円・乱択増分法）

問題