算法 | 最小椭圆覆盖 · RedKold的小站

众所周知，南京大学 计算机学院非常喜欢考察计算几何（~~和一些奇怪东西~~），2025 年开放日 jyy 出了一道最小椭圆覆盖，哀鸿遍野。在 2024 年，则考了一道图形识别题（判断三角形/椭圆/矩形）。

这类题往往八仙过海，突破了传统算法八股，对码力的要求仍然很高，必须做好准备。

~~这里祝福笔者在后天(9.5)上午的机考中获得成功~~

最小圆覆盖

问

问题阐述
输入 $N$ 个点的坐标 $(x,y)$ ，你需要输出一个最小的覆盖这些点的圆（输出圆心坐标和半径大小即可）。

解决这道题我们使用随机增量法。

注

随机增量法 (OI Wiki)
随机增量算法是计算几何的一个重要算法，它对理论知识要求不高，算法时间复杂度低，应用范围广大．

增量法 (Incremental Algorithm) 的思想与第一数学归纳法类似，它的本质是将一个问题化为规模刚好小一层的子问题．解决子问题后加入当前的对象．写成递归式是：

$T(n)=T(n-1)+g(n)$
增量法形式简洁，可以应用于许多的几何题目中．

增量法往往结合随机化，可以避免最坏情况的出现．

对于此题的做法：

假设 $C_{i-1}$ 是前 $i-1$ 个点所确定下来的最小圆（覆盖），则考虑新加入的点 $P_{i}$

若 $P_{i}$ 在 $C_{i-1}$ 中，则什么都不做
否则，加入这个点，具体来说是以 $P_{i}$ $P i $ 为基础，重复以上流程加入第 $j$ $j$ 个点 $(j<i)$ $(j < i)$ 。（这里使用两个点 $p[i]$ $p [i]$ 和 $p[j]$ $p [j]$ 构造圆）
- 构造完毕后，由于 $p[i]$ 永远是直径的一个端点，所以 $p[i]$ 一定在圆上。
- 但是可能会有 $p[j]$ 在这个过程偏移开来。所以需要再做一次。
- 对所有的 $p[i],p[j]$ ，我们遍历第 $k$ 个点 $p[k]$ 满足（ $k<j$ ）。由于三个点能唯一确定一个圆，所以三层遍历实际遍历了所有的三点组合找到了 所有最小圆覆盖 中最大的一个（或者说最小上界？），可以覆盖所有点。（又是一个 $\max\min$ 类型问题）

Code

# 万能头不是好习惯，但是先这样吧..
#include <bits/stdc++.h>
using namespace std;

int N;

#define eps 1e-9

struct Point {
    double x;
    double y;
    Point(double x, double y) : x(x), y(y) {}
};

double distance(Point u, Point v) {
    auto dx = abs(u.x - v.x);
    auto dy = abs(u.y - v.y);
    return sqrt(dx * dx + dy * dy);
}

struct Circle {
    Point o = Point(0, 0);

    double r;

    Circle(){
        r = 0;
    }

    Circle(Point o, double r) : o(o), r(r) {}

    Circle circle_from(Point p1, Point p2) {
        // two get one;
        auto _o = Point((p1.x + p2.x) / 2, (p1.y + p2.y) / 2);
        auto _r = distance(p1, p2) / 2;

        return Circle(_o, _r);
    }

    static Circle circle_from(Point u, Point v, Point w) {
        auto a1 = 2 * (v.x - u.x);
        auto b1 = 2 * (v.y - u.y);
        auto c1 = v.x * v.x + v.y * v.y - u.x * u.x - u.y * u.y;
        auto a2 = 2 * (w.x - v.x);
        auto b2 = 2 * (w.y - v.y);
        auto c2 = w.x * w.x + w.y * w.y - v.x * v.x - v.y * v.y;

        // 克拉默法则推导出来
        auto x = ((c1 * b2) - c2 * b1) / ((a1 * b2) - (a2 * b1));
        auto y = ((a1 * c2) - (a2 * c1)) / ((a1 * b2) - (a2 * b1));
        auto o = Point(x, y);
        // 圆心到任一顶点的距离就是外接圆半径，不需要再除以 2
        auto r = distance(o, u);

        return Circle(o, r);
    }
};

bool inside(Point p, Point o, double r) {
    return distance(p, o) <= r + eps;
}

Circle smallest_circle_cover(vector<Point>& p) {
    Point o = p[0];
    double r = 0;

    for (int i = 0; i < N; i++) {
        if (inside(p[i], o, r)) continue;

        // p[i] 在圆外，新圆必经过 p[i]；
        // 先取过 p[i]、p[0] 的两点圆作为候选
        o.x = (p[i].x + p[0].x) / 2;
        o.y = (p[i].y + p[0].y) / 2;
        r = distance(p[i], p[0]) / 2;

        for (int j = 1; j < i; j++) {
            if (inside(p[j], o, r)) continue;

            // p[j] 也在圆外，改取过 p[i]、p[j] 的两点圆
            o.x = (p[i].x + p[j].x) / 2;
            o.y = (p[i].y + p[j].y) / 2;
            r = distance(p[i], p[j]) / 2;

            for (int k = 0; k < j; k++) {
                if (inside(p[k], o, r)) continue;

                // p[k] 也在圆外，取 p[i]、p[j]、p[k] 的外接圆，
                // 并把它同步为当前圆
                Circle c = Circle::circle_from(p[i], p[j], p[k]);
                o = c.o;
                r = c.r;
            }
        }
    }

    return Circle(o, r);
}

int main() {
    cin >> N;
    vector<Point> ps;
    ps.reserve(N);
    double x, y;
    for (int i = 0; i < N; i++) {
        cin >> x >> y;
        ps.emplace_back(x, y);
    }

    // 打乱（Fisher-Yates），下标必须落在 [0, N)
    srand(static_cast<unsigned>(time(nullptr)));
    for (int i = N - 1; i > 0; i--) {
        swap(ps[i], ps[rand() % (i + 1)]);
    }

    auto c = smallest_circle_cover(ps);

    printf("%.9f\n%.9f %.9f\n", c.r, c.o.x, c.o.y);
}

后记

结构体是一个伟大的发明，多用结构体和函数来写更好读的代码。

你又不是 OIer 那样的高手，写稳点呗

最小圆覆盖

Code

后记

评 论

评论