計算幾何經典問題：最小覆蓋圓算法 (Welzl 算法) 完整解析與 C++ 實現

問題定義

給定平面上 n 個點，求包含所有點的最小圓 (Minimum Enclosing Circle, MEC)，即圓心 (cx, cy) 與半徑 r，使得所有點到圓心距離 ≤ r，且 r 最小。此圓唯一存在。

Welzl 算法：期望線性時間 O(n)

隨機增量構建，遞歸處理邊界點集 R (|R| ≤ 3)。

核心遞歸關係

MEC(P, R):
    if P 為空 或 |R| == 3:
        return trivial(R) // 直接由 R 決定圓
    隨機選取 p ∈ P
    D = MEC(P - {p}, R)
    if p 在 D 內:
        return D
    else:
        return MEC(P - {p}, R ∪ {p})

基礎情況 trivial(R)

C++ 完整實現

#include 
using namespace std;

struct Point { double x, y; };
struct Circle { Point c; double r; };

const double EPS = 1e-9;

double dist(const Point& a, const Point& b) {
    return hypot(a.x - b.x, a.y - b.y);
}

bool contains(const Circle& C, const Point& p) {
    return dist(C.c, p) <= C.r + EPS;
}

// 兩點定圓
Circle circleFrom2(const Point& a, const Point& b) {
    return {{(a.x+b.x)/2, (a.y+b.y)/2}, dist(a,b)/2};
}

// 三點定圓 (外接圓)
Circle circleFrom3(const Point& a, const Point& b, const Point& c) {
    double A = b.x - a.x, B = b.y - a.y;
    double C = c.x - a.x, D = c.y - a.y;
    double E = A*(a.x+b.x) + B*(a.y+b.y);
    double F = C*(a.x+c.x) + D*(a.y+c.y);
    double G = 2*(A*(c.y-b.y) - B*(c.x-b.x));
    if (fabs(G) < EPS) // 共線
        return { {0,0}, -1 }; // 標記無效
    Point o = { (D*E - B*F)/G, (A*F - C*E)/G };
    return { o, dist(o, a) };
}

// trivial(R)
Circle trivial(vector& R) {
    if (R.empty()) return {{0,0}, 0};
    if (R.size() == 1) return {R[0], 0};
    if (R.size() == 2) return circleFrom2(R[0], R[1]);
    // 3 點：檢查任意兩點定圓是否包含第三點
    Circle c = circleFrom3(R[0], R[1], R[2]);
    if (c.r >= 0) return c;
    // 共線情況：取最遠兩點
    double d01 = dist(R[0], R[1]);
    double d02 = dist(R[0], R[2]);
    double d12 = dist(R[1], R[2]);
    if (d01 >= d02 && d01 >= d12) return circleFrom2(R[0], R[1]);
    if (d02 >= d01 && d02 >= d12) return circleFrom2(R[0], R[2]);
    return circleFrom2(R[1], R[2]);
}

// Welzl 遞歸
Circle welzl(vector& P, vector R, int n) {
    if (n == 0 || R.size() == 3) return trivial(R);
    int idx = rand() % n;
    swap(P[idx], P[n-1]);
    Point p = P[n-1];
    Circle D = welzl(P, R, n-1);
    if (contains(D, p)) return D;
    R.push_back(p);
    return welzl(P, R, n-1);
}

Circle minEnclosingCircle(vector points) {
    shuffle(points.begin(), points.end(), mt19937(random_device{}()));
    return welzl(points, {}, points.size());
}

int main() {
    ios::sync_with_stdio(false);
    cin.tie(nullptr);
    int n; cin >> n;
    vector pts(n);
    for (auto& p : pts) cin >> p.x >> p.y;
    Circle ans = minEnclosingCircle(pts);
    cout << fixed << setprecision(3);
    cout << ans.c.x << " " << ans.c.y << " " << ans.r << "\n";
    return 0;
}

關鍵點解析

1. 隨機化：shuffle 確保期望線性時間，避免最惡 O(n⁴)。
2. 邊界集 R：最多 3 個點 (平面上圓由 3 點決定)，遞歸深度 ≤ 3 + n。
3. 浮點數精度：EPS = 1e-9 判斷包含關係，避免精度問題導致錯誤。
4. 共線處理：三點共線時外接圓不存在，退化為最遠兩點定圓。

應用場景

• 無線感測器網絡覆蓋優化
• 數據聚類初始化 (k-means++ 類似思想)
• 碰撞檢測包圍球
• 地理圍欄生成

擴展：最小覆蓋球 (高維)

Welzl 算法自然推廣到 d 維，邊界集大小 ≤ d+1，期望時間 O(d n)。但數值穩定性隨維度急劇下降，高維常用近似算法 (Core-set, Ritter's algorithm)。