hdu 3007 Buried memory （最小覆盖圆）

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

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

代码如下：

  1 #include <cstdio>
  2 #include <cmath>
  3 #include <iostream>
  4 #include <algorithm>
  5 #include <cstring>
  6 
  7 using namespace std;
  8 
  9 const double EPS = 1e-10;
 10 inline int sgn(double x) { return (x > EPS) - (x < -EPS);}
 11 struct Point {
 12     double x, y;
 13     Point() {}
 14     Point(double x, double y) : x(x),y(y) {}
 15     bool operator < (Point a) const { return sgn(x - a.x) < 0 || sgn(x - a.x) == 0 && sgn(y - a.y) < 0;}
 16     bool operator == (Point a) const { return sgn(x - a.x) == 0 && sgn(y - a.y) == 0;}
 17     Point operator + (Point a) const { return Point(x + a.x, y + a.y);}
 18     Point operator - (Point a) const { return Point(x - a.x, y - a.y);}
 19     Point operator * (double p) const { return Point(x * p, y * p);}
 20     Point operator / (double p) const { return Point(x / p, y / p);}
 21 } ;
 22 typedef Point Vec;
 23 inline double crossDet(Vec a, Vec b) { return a.x * b.y - a.y * b.x;}
 24 inline double crossDet(Point o, Point a, Point b) { return crossDet(a - o, b - o);}
 25 inline double dotDet(Vec a, Vec b) { return a.x * b.x + a.y * b.y;}
 26 inline double vecLen(Vec x) { return sqrt(dotDet(x, x));}
 27 inline Point normal(Vec x) { return Point(-x.y, x.x) / vecLen(x);}
 28 
 29 Point lineIntersect(Point P, Vec v, Point Q, Vec w) {
 30     Vec u = P - Q;
 31     double t = crossDet(w, u) / crossDet(v, w);
 32     return P + v * t;
 33 }
 34 inline Point getMid(Point a, Point b) { return (a + b) / 2.0;}
 35 
 36 struct Circle {
 37     Point c;
 38     double r;
 39     Circle() {}
 40     Circle(Point c, double r) : c(c), r(r) {}
 41 } ;
 42 
 43 Circle getCircle(Point a, Point b, Point c) {
 44     Vec v1 = b - a, v2 = c - a;
 45     if (sgn(dotDet(b - a, c - a)) <= 0) return Circle(getMid(b, c), vecLen(b - c) / 2.0);
 46     if (sgn(dotDet(a - b, c - b)) <= 0) return Circle(getMid(a, c), vecLen(a - c) / 2.0);
 47     if (sgn(dotDet(a - c, b - c)) <= 0) return Circle(getMid(a, b), vecLen(a - b) / 2.0);
 48     Point ip = lineIntersect(getMid(a, b), normal(v1), getMid(a, c), normal(v2));
 49     return Circle(ip, vecLen(ip - a));
 50 }
 51 
 52 int andrew(Point *pt, int n, Point *ch) {
 53     sort(pt, pt + n);
 54     int m = 0;
 55     for (int i = 0; i < n; i++) {
 56         while (m > 1 && sgn(crossDet(ch[m - 2], ch[m - 1], pt[i])) <= 0) m--;
 57         ch[m++] = pt[i];
 58     }
 59     int k = m;
 60     for (int i = n - 2; i >= 0; i--) {
 61         while (m > k && sgn(crossDet(ch[m - 2], ch[m - 1], pt[i])) <= 0) m--;
 62         ch[m++] = pt[i];
 63     }
 64     if (n > 1) m--;
 65     return m;
 66 }
 67 
 68 const int N = 555;
 69 Point pt[N], ch[N];
 70 int rnd[N];
 71 
 72 void randPoint(Point *pt, int n) {
 73     for (int i = 0; i < n; i++) rnd[i] = (rand() % n + n) % n;
 74     for (int i = 0; i < n; i++) swap(pt[i], pt[rnd[i]]);
 75 }
 76 
 77 inline bool inCircle(Point p, Circle C) { return sgn(vecLen(C.c - p) - C.r) <= 0;}
 78 
 79 int main() {
 80 //    freopen("in", "r", stdin);
 81     int n;
 82     while (cin >> n && n) {
 83         for (int i = 0; i < n; i++) scanf("%lf%lf", &pt[i].x, &pt[i].y);
 84         n = andrew(pt, n, ch);
 85         randPoint(ch, n);
 86         Circle ans = Circle(ch[0], 0.0), tmp;
 87         for (int i = 0; i < n; i++) {
 88             if (inCircle(ch[i], ans)) continue;
 89             ans = Circle(ch[i], 0.0);
 90             for (int j = 0; j < i; j++) {
 91                 if (inCircle(ch[j], ans)) continue;
 92                 ans = Circle(getMid(ch[i], ch[j]), vecLen(ch[i] - ch[j]) / 2.0);
 93                 for (int k = 0; k < j; k++) {
 94                     if (inCircle(ch[k], ans)) continue;
 95                     ans = getCircle(ch[i], ch[j], ch[k]);
 96                 }
 97             }
 98         }
 99         printf("%.2f %.2f %.2f\n", ans.c.x, ans.c.y, ans.r);
100 //        printf("%.2f\n", ans.r + 0.5);
101     }
102     return 0;
103 }

View Code

——written by Lyon