hdu 3007 Buried memory （最小覆盖圆）

http://acm.hdu.edu.cn/showproblem.php?pid=3007

　　前几天认真看了一下最小圆覆盖的论文，对于随机点集，只需要按顺序插入。每次插入分两种情况，一种是点在当前的圆中，另一种是点在园外。在圆中可以直接不理，在圆外的就要枚举所有构成最小覆盖圆的可能了。当发现点在圆外的时候，这个点必定是当前要找到的最小圆经过的点之一，于是就可以枚举出两个点和三个点的情况。看上去是一个O(n^3)的算法，但是对于随机的插入的点，剪枝可以使得复杂度降到期望的O(n)的规模。

代码如下：

#include <cstdio>
#include <cmath>
#include <iostream>
#include <algorithm>
#include <cstring>

using namespace std;

const double EPS = 1e-10;
inline int sgn(double x) { return (x > EPS) - (x < -EPS);}
struct Point {
    double x, y;
    Point() {}
    Point(double x, double y) : x(x),y(y) {}
    bool operator < (Point a) const { return sgn(x - a.x) < 0 || sgn(x - a.x) == 0 && sgn(y - a.y) < 0;}
    bool operator == (Point a) const { return sgn(x - a.x) == 0 && sgn(y - a.y) == 0;}
    Point operator + (Point a) const { return Point(x + a.x, y + a.y);}
    Point operator - (Point a) const { return Point(x - a.x, y - a.y);}
    Point operator * (double p) const { return Point(x * p, y * p);}
    Point operator / (double p) const { return Point(x / p, y / p);}
} ;
typedef Point Vec;
inline double crossDet(Vec a, Vec b) { return a.x * b.y - a.y * b.x;}
inline double crossDet(Point o, Point a, Point b) { return crossDet(a - o, b - o);}
inline double dotDet(Vec a, Vec b) { return a.x * b.x + a.y * b.y;}
inline double vecLen(Vec x) { return sqrt(dotDet(x, x));}
inline Point normal(Vec x) { return Point(-x.y, x.x) / vecLen(x);}

Point lineIntersect(Point P, Vec v, Point Q, Vec w) {
    Vec u = P - Q;
    double t = crossDet(w, u) / crossDet(v, w);
    return P + v * t;
}
inline Point getMid(Point a, Point b) { return (a + b) / 2.0;}

struct Circle {
    Point c;
    double r;
    Circle() {}
    Circle(Point c, double r) : c(c), r(r) {}
} ;

Circle getCircle(Point a, Point b, Point c) {
    Vec v1 = b - a, v2 = c - a;
    if (sgn(dotDet(b - a, c - a)) <= 0) return Circle(getMid(b, c), vecLen(b - c) / 2.0);
    if (sgn(dotDet(a - b, c - b)) <= 0) return Circle(getMid(a, c), vecLen(a - c) / 2.0);
    if (sgn(dotDet(a - c, b - c)) <= 0) return Circle(getMid(a, b), vecLen(a - b) / 2.0);
    Point ip = lineIntersect(getMid(a, b), normal(v1), getMid(a, c), normal(v2));
    return Circle(ip, vecLen(ip - a));
}

int andrew(Point *pt, int n, Point *ch) {
    sort(pt, pt + n);
    int m = 0;
    for (int i = 0; i < n; i++) {
        while (m > 1 && sgn(crossDet(ch[m - 2], ch[m - 1], pt[i])) <= 0) m--;
        ch[m++] = pt[i];
    }
    int k = m;
    for (int i = n - 2; i >= 0; i--) {
        while (m > k && sgn(crossDet(ch[m - 2], ch[m - 1], pt[i])) <= 0) m--;
        ch[m++] = pt[i];
    }
    if (n > 1) m--;
    return m;
}

const int N = 555;
Point pt[N], ch[N];
int rnd[N];

void randPoint(Point *pt, int n) {
    for (int i = 0; i < n; i++) rnd[i] = (rand() % n + n) % n;
    for (int i = 0; i < n; i++) swap(pt[i], pt[rnd[i]]);
}

inline bool inCircle(Point p, Circle C) { return sgn(vecLen(C.c - p) - C.r) <= 0;}

int main() {
//    freopen("in", "r", stdin);
    int n;
    while (cin >> n && n) {
        for (int i = 0; i < n; i++) scanf("%lf%lf", &pt[i].x, &pt[i].y);
        n = andrew(pt, n, ch);
        randPoint(ch, n);
        Circle ans = Circle(ch[0], 0.0), tmp;
        for (int i = 0; i < n; i++) {
            if (inCircle(ch[i], ans)) continue;
            ans = Circle(ch[i], 0.0);
            for (int j = 0; j < i; j++) {
                if (inCircle(ch[j], ans)) continue;
                ans = Circle(getMid(ch[i], ch[j]), vecLen(ch[i] - ch[j]) / 2.0);
                for (int k = 0; k < j; k++) {
                    if (inCircle(ch[k], ans)) continue;
                    ans = getCircle(ch[i], ch[j], ch[k]);
                }
            }
        }
        printf("%.2f %.2f %.2f\n", ans.c.x, ans.c.y, ans.r);
//        printf("%.2f\n", ans.r + 0.5);
    }
    return 0;
}

——written by Lyon