「算法笔记」最小圆覆盖

最小圆覆盖利用随机增量法,是一种看似暴力但实际上很高效的算法。


概述

最小圆覆盖算法用来求:给定 $n$ 个点求覆盖所有点的最小的圆。


实现

首先我们对这 $n$ 个点随机打乱,可以利用 $\text{C++}$ 中的 $\text{random_shuffle}$ 函数。

首先令当前最小圆的圆心为点 $1$,半径为 $0$;接下来分成 $3​$ 重循环。

  1. 枚举点 $i\in [1,n]$,如果点 $i$ 不在圆中则把圆心定为点 $i$ 并进入下层循环。
  2. 枚举点 $j\in[1,i)$,如果点 $j$ 不在圆中则用点 $i,j$ 构造圆并进入下层循环。
  3. 枚举点 $k\in[1,j)$,如果点 $k$ 不在圆中则用点 $i,j,k$ 构造圆(显然三点定圆),圆心使用原方程随便推一下就可以计算了。

时间复杂度:$\mathcal O(n)$,详见下方证明。


复杂度

上述过程看起来复杂度为 $\mathcal O(n^3)$,接下来进行严格分析。

由于这 $n$ 个点的最小覆盖圆只需要 $3$ 个点就可以确定了,因此每个点参与确定最小圆的概率为 $\frac{3}{n}$。所以第 $1,2$ 层循环中在第 $i$ 个点进入下一层循环的概率约为 $\frac{3}{i}$ 和 $\frac{2}{i}$。

设该算法的三个循环从外到内的复杂度分别为 $T_1(n)​$;则有:

$$ \begin{aligned} T_3(n)&=\mathcal O(n) \\ T_2(n)&=\mathcal O(n)+\sum_{i=1}^n\frac{i-2}{i}\cdot \mathcal O(1)+\frac{2}{i}\cdot T_3(i)=\mathcal O(n) \\ T_1(n)&=\mathcal O(n)+\sum_{i=1}^n\frac{i-3}{i}\cdot \mathcal O(1)+\frac{3}{i}\cdot T_2(i)=\mathcal O(n) \end{aligned} $$

故整个算法的复杂度为 $\mathcal O(n)$!


代码

#include <cstdio>
#include <cmath>
#include <cstdlib>
#include <ctime>
#include <algorithm>

const int N = 1e6 + 5;
const double eps = 1e-12;

int n;
double r;

struct Point {
    double x, y;
    Point(double _x = 0, double _y = 0) {
        x = _x, y = _y;
    }
} p[N], o;

double sqr(double x) {
    return x * x;
}
double dis(Point a, Point b) {
    return sqrt(sqr(a.x - b.x) + sqr(a.y - b.y));
}
bool cmp(double a, double b) {
    return fabs(a - b) < eps;
}
bool inCir(Point a, Point o, double r) {
    double d = dis(a, o);
    return d < r || cmp(d, r);
}
void calc(Point a, Point b, Point c, Point &o, double &r) {
    double a1 = 2 * (a.x - b.x), b1 = 2 * (a.y - b.y), c1 = sqr(a.x) + sqr(a.y) - sqr(b.x) - sqr(b.y);
    double a2 = 2 * (a.x - c.x), b2 = 2 * (a.y - c.y), c2 = sqr(a.x) + sqr(a.y) - sqr(c.x) - sqr(c.y);
    o.x = (c1 * b2 - c2 * b1) / (a1 * b2 - a2 * b1);
    o.y = (c1 * a2 - c2 * a1) / (b1 * a2 - b2 * a1);
    r = dis(a, o);
}
void getCir(Point *p, int n, Point &o, double &r) {
    srand(time(0) + *(new char));
    std::random_shuffle(p + 1, p + n + 1);
    o = p[1], r = 0;
    for (int i = 1; i <= n; i++) {
        if (inCir(p[i], o, r)) continue;
        o = p[i], r = 0;
        for (int j = 1; j < i; j++) {
            if (inCir(p[j], o, r)) continue;
            o.x = (p[i].x + p[j].x) / 2;
            o.y = (p[i].y + p[j].y) / 2;
            r = dis(o, p[i]);
            for (int k = 1; k < j; k++) {
                if (inCir(p[k], o, r)) continue;
                calc(p[i], p[j], p[k], o, r);
            }
        }
    }
}
int main() {
    scanf("%d", &n);
    for (int i = 1; i <= n; i++) {
        scanf("%lf%lf", &p[i].x, &p[i].y);
    }
    getCir(p, n, o, r);
    printf("%.10lf %.10lf %.10lf\n", o.x, o.y, r);
    return 0;
}

习题

发表评论