bzoj 1336 最小圆覆盖

最小圆覆盖

问题：给定平面上的一个点集，求半径最小的一个圆，使得点集中的点都在其内部或上面。

随机增量算法：

　　定义：点集A的最小圆覆盖是Circle(A)

　　定理：如果Circle(A)=C₁，且a不被C₁覆盖，那么a在Circle(AU{a})的边界上。

　　证明：换一种找最小圆覆盖的思路，我们初始化一些圆，圆心为A中的点，半径为0，并且让半径慢慢变大，必定存在一个时刻，所有圆的交集由空变为非空，那个最开始的非空交集是一个点，并且就是我们最小圆覆盖的圆心位置。当A中的所有点代表的圆有交集时，点a代表的圆还没有到达那个点（否则点a就被C₁覆盖掉了），我们让半径继续增大，必然会有一个时刻点a代表的圆与A代表的圆的公共区域相交，这个点就是AU{a}的最小圆覆盖的圆心，它到点a的距离就是半径。

　　算法：

 1 c = ( p[1] )
 2 for i = 2 to n
 3     if p[i] in c then continue
 4     c = ( p[i] )
 5     for j = 1 to i-1
 6         if p[j] in c then continue
 7         c = ( p[i], p[j] )
 8         for k = 1 to j-1
 9             if p[k] in c then continue
10             c = ( p[i], p[j], p[k] )

　　（(p[i],p[j])代表包含这两个点的最小的圆）

　　第一层循环的循环不变量是：c是p[1],p[2],...,p[i-1]的最小圆覆盖。

　　第二层循环的循环不变量是：c是p[1],p[2],...,p[j-1]和p[i]的最小圆覆盖。

　　第一层循环的循环不变量是：c是p[1],p[2],...,p[k-1],p[i]和p[j]的最小圆覆盖。

　　转移用上面的定理证明。

　　有个性质：上面伪代码的第10行中的三个点不可能共线（只需分别证明三个点中的一个不会在另外两个代表的线段上就行了）。

 1 /**************************************************************
 2     Problem: 1336
 3     User: idy002
 4     Language: C++
 5     Result: Accepted
 6     Time:372 ms
 7     Memory:2372 kb
 8 ****************************************************************/
 9  
10 #include <cstdio>
11 #include <cmath>
12 #include <algorithm>
13 #define line(a,b) ((b)-(a))
14 #define N 100010
15 #define eps 1e-10
16 using namespace std;
17  
18 int sg( double x ) { return (x>-eps)-(x<eps); }
19 struct Vector {
20     double x, y;
21     Vector(){}
22     Vector( double x, double y ):x(x),y(y){}
23     Vector operator+( const Vector &b ) const { return Vector(x+b.x,y+b.y); }
24     Vector operator-( const Vector &b ) const { return Vector(x-b.x,y-b.y); }
25     Vector operator*( double b ) const { return Vector(x*b,y*b); }
26     Vector operator/( double b ) const { return Vector(x/b,y/b); }
27     double operator^( const Vector &b ) const { return x*b.y-y*b.x; }
28     double len() { return sqrt(x*x+y*y); }
29     Vector normal() { return Vector(-y,x); }
30 };
31 typedef Vector Point;
32 Point inter( Point P, Vector u, Point Q, Vector v ) {
33     return P+u*((line(P,Q)^v)/(u^v));
34 }
35 struct Circle {
36     Point o;
37     double r;
38     Circle(){}
39     Circle( Point &a ) {
40         o = a;
41         r = 0;
42     }
43     Circle( Point &a, Point &b ) {
44         o = (a+b)/2;
45         r = (a-b).len()/2;
46     }
47     Circle( Point &a, Point &b, Point &c ) {
48         Point P=(a+b)/2, Q=(b+c)/2;
49         Vector u=(a-b).normal(), v=(b-c).normal();
50         o = inter(P,u,Q,v);
51         r = (a-o).len();
52     }
53     bool contain( Point &a ) {
54         return sg( line(a,o).len()-r ) <= 0;
55     }
56 };
57  
58 int n;
59 Point pts[N];
60 Circle cir;
61  
62 int main() {
63     scanf( "%d", &n );
64     for( int i=1; i<=n; i++ ) {
65         double x, y;
66         scanf( "%lf%lf", &x, &y );
67         pts[i] = Point(x,y);
68     }
69     random_shuffle( pts+1, pts+1+n );
70     cir = Circle(pts[1]);
71     for( int i=2; i<=n; i++ ) {
72         if( cir.contain(pts[i]) ) continue;
73         cir = Circle(pts[i]);
74         for( int j=1; j<i; j++ ) {
75             if( cir.contain(pts[j]) ) continue;
76             cir = Circle(pts[i],pts[j]);
77             for( int k=1; k<j; k++ ) {
78                 if( cir.contain(pts[k]) ) continue;
79                 cir = Circle(pts[i],pts[j],pts[k]);
80             }
81         }
82     }
83     printf( "%.10lf
", cir.r );
84     printf( "%.10lf %.10lf
", cir.o.x, cir.o.y );
85 }
86

View Code