说在前面
这个算法还是挺简单的…
然而me并不能理解其正确性
网上的博客也都没有证明
于是手动模拟算法用心感知正确性hhhhh
题目
题面
给出N个点,让你画一个最小的包含所有点的圆。
点数不超过100000
输入输出格式
输入格式:
第一行一个整数N,表示点的个数
接下来每行两个实数(x,y)表示一个点的坐标
输出格式:
输出共两行,第一行为半径,第二行为圆心坐标,有SPJ
解法
就是一个板子题,关于算法流程,最小圆覆盖算法 by commonc
me感性理解(口胡)的正确性
对一个点集作最小圆覆盖,等于对这个点集的凸包作最小圆覆盖
假设已经获得了包含前i个点的最小圆,现在要让它包含第i+1个点。如果第i+1个点本来就在圆内那么不用管。不然的话,这个点一定在原来的圆可以覆盖的凸包外,那么新的圆一定要过这个点。(这一个for是求包含前i+1个点的最小圆覆盖)
现在求包含第i+1个点和前j个点的最小圆覆盖,对于第二个点,按照for j = 1 to i的顺序依次加入,如果第j个点在圆外,同理,新圆一定要过第j个点。
现在求包含第i+1个点和第j个点和前k个点的最小圆覆盖,对于第三个点,按照for k = 1 to j-1的顺序依次加入,原理同上。并且一定不会出现有两个点u,v(i>j>u,v)不能被同时包含的情况(如果出现,那么圆i,u,v应该包含了点j,且先被枚举到,即在第二层循环时j已在圆内,根本就不会再进第三层循环)。
那么一定可以求出「包含第i+1个点和第j个点和前k个点的最小圆覆盖」,当k=j-1时求出「包含第i+1个点和前j个点的最小圆覆盖」,当j=i时求出「包含前i+1个点的最小圆覆盖」。所以最后能够求出所有点的最小圆覆盖
时间复杂度是期望 Θ(n) Θ ( n ) 的,上面那篇博客有证明
下面是自带大长度的代码
#include <cmath>
#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std ;
int N ;
double R , PI = 3.1415926535897932384 , eps = 1e-10 ;
struct Vector{
double x , y ;
double len(){
return sqrt( x * x + y * y ) ;
}
Vector rotate( double rad ){
return ( Vector ){ x*cos(rad) - y*sin(rad) , x*sin(rad) + y*cos(rad) } ;
}
}p[100005] , O ;
typedef Vector Point ;
Vector operator + ( const Vector &A , const Vector &B ){ return ( Vector ){ A.x + B.x , A.y + B.y } ; }
Vector operator - ( const Vector &A , const Vector &B ){ return ( Vector ){ A.x - B.x , A.y - B.y } ; }
Vector operator * ( const Vector &A , const double &p ){ return ( Vector ){ A.x * p , A.y * p } ;}
Vector operator / ( const Vector &A , const double &p ){ return ( Vector ){ A.x / p , A.y / p } ;}
double cross( const Vector &A , const Vector &B ){ return A.x * B.y - A.y * B.x ; }
double dot( const Vector &A , const Vector &B ){ return A.x * B.x + A.y * B.y ; }
double dis ( const Point &A , const Point &B ){
double t1 = A.x - B.x , t2 = A.y - B.y ;
return sqrt( t1 * t1 + t2 * t2 ) ;
}
bool inside( const Point &A ){ return dis( A , O ) <= R ; }
double dcmp( double x ){
if( x > -eps && x < eps ) return 0 ;
return x > eps ? 1 : -1 ;
}
Point glt( const Point &P , const Vector &u , const Point &Q , const Vector &v ){
double t = cross( Q - P , v ) / cross( u , v ) ;
return P + u * t ;
}
Point cal( Point &A , Point &B , Point &C ){
Vector a = A - B , b = B - C ;
if( dcmp( cross( a , b ) ) == 0 ){
double t1 = dis( A , B ) , t2 = dis( B , C ) , t3 = dis( A , C ) ;
if( t1 >= t2 && t1 >= t3 ) return ( A + B ) / 2 ;
if( t2 >= t1 && t2 >= t3 ) return ( B + C ) / 2 ;
if( t3 >= t1 && t3 >= t2 ) return ( A + C ) / 2 ;
} return glt( ( A + B ) / 2 , a.rotate( PI / 2 ) , ( B + C ) / 2 , b.rotate( PI / 2 ) ) ;
}
void solve(){
for( int i = 1 ; i <= N ; i ++ )
if( !inside( p[i] ) ){
O = p[i] ; R = 0 ;
for( int j = 1 ; j < i ; j ++ )
if( !inside( p[j] ) ){
O = ( p[i] + p[j] ) / 2 ;
R = dis( p[i] , O ) ;
for( int k = 1 ; k < j ; k ++ )
if( !inside( p[k] ) ){
O = cal( p[i] , p[j] , p[k] ) ;
R = dis( p[i] , O ) ;
}
}
}
printf( "%f\n" , R ) ;
printf( "%f %f\n" , O.x , O.y ) ;
}
int main(){
srand( 20031109 ) ;
scanf( "%d" , &N ) ;
for( int i = 1 ; i <= N ; i ++ )
scanf( "%lf%lf" , &p[i].x , &p[i].y ) ;
random_shuffle( p + 1 , p + N + 1 ) ;
solve() ;
}