최소 외접원 (Minimum Enclosing Circle, MEC) 은 평면 위 N 개의 점을 모두 포함하는 가장 작은 반지름의 원을 찾는 문제. 대표 알고리즘은 Emo Welzl (1991) 의 randomized incremental algorithm, 기댓값 O(N).
문제 상황과 동기
모든 점을 하나의 원으로 커버하면서 반지름을 최소화.
Naive 접근: 점 2개 또는 3개로 결정되는 모든 원을 시도하며 O(N^4).
핵심 통찰: MEC 는 항상 2개 또는 3개의 support point 에 의해 결정. 점을 무작위 순서로 추가하며 원 밖의 점만 경계로 편입.
MEC 의 성질:1. MEC 는 점 2개가 직경을 결정 (지름) 하거나, 3개 점이 원주 위에 있음.2. MEC 내부에 있는 점은 제거해도 원이 변하지 않음.3. 새로운 점 p 가 현재 원 밖에 있으면, p 는 새 원의 경계 위에 있어야 함.Welzl 의 재귀: welzl(P, R): // P: 남은 점, R: 경계 점 if P == {} or |R| == 3: return trivial_circle(R) p = random(P) C = welzl(P - {p}, R) if p in C: return C return welzl(P - {p}, R ∪ {p})
알고리즘
welzl(P, R): if P is empty or |R| == 3: return circle_from(R) // R 크기에 따라 0/1/2/3 점 처리 p = random element from P C = welzl(P \ {p}, R) if p inside C: return C R.push(p) C = welzl(P \ {p}, R) R.pop() return Ccircle_from(R): if |R| == 0: return zero-radius at origin if |R| == 1: return zero-radius at R[0] if |R| == 2: return circle with R[0],R[1] as diameter if |R| == 3: return circumcircle of R[0],R[1],R[2]
구현
// Welzl's algorithm, expected O(N)#include <bits/stdc++.h>using namespace std;struct Pt { double x, y; };struct Cir { Pt c; double r; };double dist(const Pt& a, const Pt& b) { return hypot(a.x-b.x, a.y-b.y);}bool inside(const Pt& p, const Cir& cir) { return dist(p, cir.c) <= cir.r + 1e-9;}Cir c2(const Pt& a, const Pt& b) { return {{(a.x+b.x)/2, (a.y+b.y)/2}, dist(a,b)/2};}Cir c3(const Pt& a, const Pt& b, const Pt& c) { double d = 2*(a.x*(b.y-c.y)+b.x*(c.y-a.y)+c.x*(a.y-b.y)); if (abs(d) < 1e-12) return c2(a, b); double ux = ((a.x*a.x+a.y*a.y)*(b.y-c.y) + (b.x*b.x+b.y*b.y)*(c.y-a.y) + (c.x*c.x+c.y*c.y)*(a.y-b.y)) / d; double uy = ((a.x*a.x+a.y*a.y)*(c.x-b.x) + (b.x*b.x+b.y*b.y)*(a.x-c.x) + (c.x*c.x+c.y*c.y)*(b.x-a.x)) / d; return {{ux, uy}, dist({ux,uy}, a)};}Cir welzl(vector<Pt>& P, vector<Pt>& R, int n) { if (n == 0 || R.size() == 3) { if (R.empty()) return {{0,0}, 0}; if (R.size() == 1) return {R[0], 0}; if (R.size() == 2) return c2(R[0], R[1]); return c3(R[0], R[1], R[2]); } int idx = rand() % n; swap(P[idx], P[n-1]); Pt p = P[n-1]; Cir C = welzl(P, R, n-1); if (inside(p, C)) return C; R.push_back(p); C = welzl(P, R, n-1); R.pop_back(); return C;}int main() { int N; cin >> N; vector<Pt> P(N); for (auto& p : P) cin >> p.x >> p.y; Cir ans = welzl(P, *new vector<Pt>(), N); cout << fixed << setprecision(10); cout << ans.c.x << " " << ans.c.y << " " << ans.r << "\n";}
import math, random, sysdef dist(a, b): return math.hypot(a[0]-b[0], a[1]-b[1])def inside(p, cx, cy, r): return dist(p, (cx, cy)) <= r + 1e-9def c2(a, b): return ((a[0]+b[0])/2, (a[1]+b[1])/2, dist(a,b)/2)def c3(a, b, c): d = 2*(a[0]*(b[1]-c[1]) + b[0]*(c[1]-a[1]) + c[0]*(a[1]-b[1])) if abs(d) < 1e-12: return c2(a, b) ux = ((a[0]**2+a[1]**2)*(b[1]-c[1]) + (b[0]**2+b[1]**2)*(c[1]-a[1]) + (c[0]**2+c[1]**2)*(a[1]-b[1])) / d uy = ((a[0]**2+a[1]**2)*(c[0]-b[0]) + (b[0]**2+b[1]**2)*(a[0]-c[0]) + (c[0]**2+c[1]**2)*(b[0]-a[0])) / d return (ux, uy, dist((ux,uy), a))def welzl(P, R): if not P or len(R) == 3: if len(R) == 0: return (0.0, 0.0, 0.0) if len(R) == 1: return (R[0][0], R[0][1], 0.0) if len(R) == 2: return c2(R[0], R[1]) return c3(R[0], R[1], R[2]) idx = random.randrange(len(P)) P[idx], P[-1] = P[-1], P[idx] p = P.pop() cx, cy, r = welzl(P, R) if inside(p, cx, cy, r): P.append(p) return (cx, cy, r) R.append(p) cx, cy, r = welzl(P, R) P.append(p) R.pop() return (cx, cy, r)N = int(sys.stdin.readline())pts = [tuple(map(float, sys.stdin.readline().split())) for _ in range(N)]cx, cy, r = welzl(pts, [])print(f"{cx:.10f} {cy:.10f} {r:.10f}")
import java.util.*;import java.io.*;public class Main { static class Pt { double x, y; Pt(double x, double y) { this.x=x; this.y=y; } } static class Cir { Pt c; double r; Cir(Pt c, double r) { this.c=c; this.r=r; } } static double d(Pt a, Pt b) { return Math.hypot(a.x-b.x, a.y-b.y); } static boolean in(Pt p, Cir cir) { return d(p, cir.c) <= cir.r+1e-9; } static Cir c2(Pt a, Pt b) { return new Cir(new Pt((a.x+b.x)/2,(a.y+b.y)/2), d(a,b)/2); } static Cir c3(Pt a, Pt b, Pt c) { double dd = 2*(a.x*(b.y-c.y)+b.x*(c.y-a.y)+c.x*(a.y-b.y)); if (Math.abs(dd)<1e-12) return c2(a,b); double ux = ((a.x*a.x+a.y*a.y)*(b.y-c.y)+(b.x*b.x+b.y*b.y)*(c.y-a.y)+(c.x*c.x+c.y*c.y)*(a.y-b.y))/dd; double uy = ((a.x*a.x+a.y*a.y)*(c.x-b.x)+(b.x*b.x+b.y*b.y)*(a.x-c.x)+(c.x*c.x+c.y*c.y)*(b.x-a.x))/dd; return new Cir(new Pt(ux,uy), d(new Pt(ux,uy), a)); } static Cir welzl(Pt[] P, List<Pt> R, int n, Random rnd) { if (n==0||R.size()==3) { if (R.isEmpty()) return new Cir(new Pt(0,0),0); if (R.size()==1) return new Cir(R.get(0),0); if (R.size()==2) return c2(R.get(0),R.get(1)); return c3(R.get(0),R.get(1),R.get(2)); } int idx = rnd.nextInt(n); Pt tmp = P[idx]; P[idx]=P[n-1]; P[n-1]=tmp; Pt p = P[n-1]; Cir C = welzl(P, R, n-1, rnd); if (in(p, C)) return C; R.add(p); C = welzl(P, R, n-1, rnd); R.remove(R.size()-1); return C; } public static void main(String[] args) throws IOException { BufferedReader br = new BufferedReader(new InputStreamReader(System.in)); int N = Integer.parseInt(br.readLine()); Pt[] P = new Pt[N]; for (int i=0; i<N; i++) { StringTokenizer st = new StringTokenizer(br.readLine()); P[i] = new Pt(Double.parseDouble(st.nextToken()), Double.parseDouble(st.nextToken())); } Cir ans = welzl(P, new ArrayList<>(), N, new Random()); System.out.printf("%.10f %.10f %.10f%n", ans.c.x, ans.c.y, ans.r); }}
stdin
30 01 00 1
결과
0.5000000000 0.5000000000 0.7071067812
stdin
40 02 02 20 2
결과
1.0000000000 1.0000000000 1.4142135624
복잡도
항목
값
시간 (기댓값)
O(N)
시간 (최악)
O(N^2) (매우 드묾)
공간
O(N) (재귀 스택)
변형 / 활용
Minimum enclosing ball: 고차원 공간으로 일반화. Welzl 은 임의 차원에 동작.
Convex hull 최적화: MEC 는 항상 볼록 껍질 위의 점에 의해 결정. 껍질만 추려낸 후 Welzl 적용 가능.
응용: 클러스터링 반경 최소화, GPS 오차 보정, 충돌 감지.
함정
1. 3 점 공선 (collinear)
세 점이 일직선이면 외접원이 정의되지 않음. circle_from_3 에서 분모 d = 0 이면 circle_from_2 로 fallback.
2. 실수 오차
eps = 1e-9 로 inside 판정. 특히 공선 판정에서 오차 보정이 중요.
3. 재귀 깊이 / 복사
Welzl 의 naive 복사 구현은 O(N^2). 인덱스 기반 swap 으로 제자리 처리하여 O(N) 유지.
💬 댓글