Dyd's Blog

He who has a strong enough why can bear almost any how.

最小圆覆盖

Posted on In
Symbols count in article: 2.9k Reading time ≈ 3 mins.

玄学 $O(n)$

最小圆覆盖

问题

给出 $n$ 个点，让你画一个最小的包含所有点的圆

性质

  1. 最小覆盖圆是唯一的

    简证：若存在两个不同的最小覆盖圆，则点一定都在两圆交集内，以相交弦为直径的圆也覆盖所有点且直径更小，矛盾

  2. 若点 $p$ 不在点集 $S$ 的最小覆盖圆内（包括边上），则 $p$ 一定在点集 $\{p\} \vee S$ 的最小覆盖圆的边上

    简证：先画一个很大的圆 $P_1$ 包括 $\{p\} \vee S$ ，然后在保证包括 $S$ 的前提下把圆缩小至不能再缩小，此时圆 $P_2$ 就是 $S$ 的最小覆盖圆，因为 $p \not\in P_2, p \in P_1$ ，所以在缩小过程中必然有一个圆 $P_3$ 使得 $p$ 在其边上，当圆的半径小于 $P_3$ 时， $p$ 就到圆外了，所以 $P_3$ 就是 $\{p\} \vee S$ 的最小覆盖圆

  3. 最后得到的最小覆盖圆的边上一定要么至少有三个点，要么只有两个点且两点连线为圆的直径

    简证：如果圆上只有一个点，则圆无法确定，肯定可以通过平移加缩小找到一个更小的圆包住所有点，与最小不符，而三个点或一条直径确定了一个圆

过程

现在有了如上性质，我们考虑如何求出最小覆盖圆

由性质 $3$ ，我们考虑如何求出圆上的三个点，具体地：

  1. 先将点随机化，先让第一个点作为一个“圆”
  2. 枚举 $i$ 从 $2$ 到 $n$ ，如果点 $i$ 在当前圆内（包括边）就继续；如果不在，进入步骤 $3$
  3. 由性质 $2$ ， $i$ 一定在点 $1 \sim i$ 的最小覆盖圆边上，所以我们先将圆改成点 $i$ ，然后进入步骤 $4$
  4. 枚举 $j : 1 \to i - 1$ ，如果点 $j$ 在当圆内（包括边）就继续；如果不在，将圆改成以点 $i, j$ 为直径的圆，进入步骤 $5$
  5. 枚举 $k : 1 \to j - 1$ ，如果点 $k$ 在当圆内（包括边）就继续；如果不在，将圆改成以点 $i, j, k$ 确定的圆

时间

看似有三层循环，时间为 $O(n^3)$ ，但我们进入下一层循环是有条件的，即该点不在圆上，而由性质 $1$ ，最小覆盖圆唯一的，又因为三个点确定一个圆，所以在这个圆被确定前，每个点有 $\frac{3}{n}$ 的几率是在最小覆盖圆的边上的，这也是进入下一层的条件，所以二、三层的时间复杂度为 $O(n) * \frac{3}{n} = O(3)$ ，故总时间为 $(n) * O(3) * O(3) = O(n)$

代码

板子

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
#include <bits/stdc++.h>
#define STC static
#define LL long long
#define DB double
#define DRG default_random_engine
using namespace std;
const int N = 1e5 + 5;
const DB PI = acos(-1), eps = 1e-12;
DRG e{20051221};
struct Point
{
    double x, y;
    Point operator + (const Point &t) const
    {
        return {x + t.x, y + t.y};
    }
    Point operator - (const Point &t) const
    {
        return {x - t.x, y - t.y};
    }
    Point operator * (const DB &t) const
    {
        return {x * t, y * t};
    }
    DB operator * (const Point &t) const
    {
        return x * t.y - y * t.x;
    }
    Point operator / (const DB &t) const
    {
        return {x / t, y / t};
    }
} ;
struct Circle
{
    Point p;
    DB r;   
} ;
int cmp(DB x, DB y)
{
    if (fabs(x - y) < eps)
        return 0;
    return x > y ? 1 : -1;
}
Point rotate(Point x, DB th)
{
    return {x.x * cos(th) + x.y * sin(th), -x.x * sin(th) + x.y * cos(th)};
}
DB dist(Point x, Point y)
{
    Point t = x - y;
    return sqrt(t.x * t.x + t.y * t.y);
}
Point jiao(Point p, Point u, Point q, Point v)
{
    DB t = (v * (p - q)) / (u * v);
    return p + u * t;
}
pair<Point, Point> get_line(Point x, Point y) //找中垂线
{
    return {(x + y) / 2, rotate(y - x, PI / 2)};
}
Circle get_c(Point x, Point y, Point z)
{
    auto u = get_line(x, y), v = get_line(x, z);
    auto p = jiao(u.first, u.second, v.first, v.second);
    return {p, dist(p, x)};
}
int main()
{
    int n;
    scanf("%d", &n);
    STC Point q[N];
    for (int i = 0; i < n; ++i)
        scanf("%lf %lf", &q[i].x, &q[i].y);
    shuffle(q, q + n, e);
    Circle c({q[0], 0});
    for (int i = 1; i < n; ++i)
        if (cmp(c.r, dist(c.p, q[i])) < 0)
        {
            c = {q[i], 0};
            for (int j = 0; j < i; ++j)
                if (cmp(c.r, dist(c.p, q[j])) < 0)
                {
                    c = {(q[i] + q[j]) / 2, dist(q[i], q[j]) / 2};
                    for (int k = 0; k < j; ++k)
                        if (cmp(c.r, dist(c.p, q[k])) < 0)
                            c = get_c(q[i], q[j], q[k]);
                }
        }
    printf("%.10lf\n", c.r);
    printf("%.10lf %.10lf\n", c.p.x, c.p.y);
    return 0;
}