最小半径圆

题目描述：

平面上有n个点（0<=n<=1000），每个点用一对坐标（x，y）表示。其中x，y分别为点的x轴和y轴坐标。同时约定0<=x<=10000，0<=y<=10000，且x，y为整数。

当n个点的坐标给出后，试找出一个半径最小的圆，将n个点全部包围，点可以在圆周上。

输入格式：

第一行包含一个整数m（m<=10）表示有m组测试数据。每组测试数据的第一行为一个整数n，表示当前这组测试数据包含n个点，接下来的n行中每行有两个点，表示点的坐标。

输出格式：

对于每组数据，在一行内输出三个数，分别表示圆心坐标和半径，精确到小数点后两位。

输入输出样例：

Sample input	Output for the input
3 3 0 0 1 0 0 1 4 0 0 0 1 1 1 1 0 4 0 0 2 0 4 0 2 2	0.50 0.50 0.71 0.50 0.50 0.71 2.00 0.00 2.00

 

解：

一个比较直接的想法就是枚举所有点的情况，这样就可以找到最小半径的外接圆。但是这样做的时间复杂度太高，我们需要寻求其它的方法：

给定点集A，记其最小外接圆为mincircle（A），显然对A来说mincircle（A）是存在且唯一的。还要注意一些特殊情况：当A为空集时，mincircle（A）为空；当A中只有一个点时，mincircle（A）的圆心就为该点，并且半径为0；当A中只有两个点时，mincircle（A）的圆心为两点连线的中点，半径为两点距离的一半。

显然，mincircle（A）可由A上最多三点确定，也就是说存在点集合B，|B|<=3，且B包含于A，有mincircle（A）＝mincircle（B）。所以如果点a不属于集合B，那么mincircle（A－{a}）＝mincircle（B）；如果mincircle（A－{a}）<> mincircle（B），那么点a一定属于集合B。

所以我们从一个空集R开始，不断向里面添加点，并维护R的外接圆最小。这样就可以得到解决该题的算法。

 

代码：

 

#include < stdio.h >
#include < math.h >
#include < string .h >
#include < stdlib.h >


#define     sqr(a)    ((a)*(a))
#define  MAX 1000
#define  EPS 1E-6


struct  TPoint{     // 点
     double  x,y;
}point[MAX];

struct  TCircle{     // 圆
    TPoint centre;
     double  r;
}circle;

struct  TTriangle{     // 三角形
    TPoint p0,p1,p2;
}circleedge;

int  casenum,pointnum;

inline  double  
distance(TPoint p1,TPoint p2){     // 计算两点间的距离
     return  sqrt(sqr(p1.x - p2.x) + sqr(p1.y - p2.y));
}

inline  double
triangleArea(TTriangle t){     // 计算三角形面积
     return  fabs(t.p0.x * t.p1.y + t.p1.x * t.p2.y + t.p2.x * t.p0.y  -  t.p1.x * t.p0.y - t.p2.x * t.p1.y - t.p0.x * t.p2.y) / 2 ;
}

inline TCircle
circumcircle(TTriangle t){     // 计算三角形的外接圆
    TCircle temp;
     double  a,b,c,c1,c2;
     double  xa,ya,xb,yb,xc,yc;
    a = distance(t.p0,t.p1);
    b = distance(t.p1,t.p2);
    c = distance(t.p2,t.p0);
    temp.r = a * b * c / triangleArea(t) / 4 ;
    xa = t.p0.x;    ya = t.p0.y;
    xb = t.p1.x;    yb = t.p1.y;
    xc = t.p2.x;    yc = t.p2.y;
    c1 = (sqr(xa) + sqr(ya) - sqr(xb) - sqr(yb)) / 2 ;
    c2 = (sqr(xa) + sqr(ya) - sqr(xc) - sqr(yc)) / 2 ;
    temp.centre.x = (c1 * (ya - yc) - c2 * (ya - yb)) / ((xa - xb) * (ya - yc) - (xa - xc) * (ya - yb));
    temp.centre.y = (c1 * (xa - xc) - c2 * (xa - xb)) / ((ya - yb) * (xa - xc) - (ya - yc) * (xa - xb));
     return  temp;
}

inline  bool
incircle(TPoint p,TCircle c){     // 判断点p是否在圆c内,如果在圆内放回true否则放回false
     return  distance(p,c.centre) < c.r + EPS ?   true : false ;
}


inline TCircle
minCircle( int  pcount,TTriangle ce){     // 计算最小外接圆
    TCircle temp;
    memset( & temp, 0 , sizeof (temp));
     switch (pcount){
     case   0 :         // 没有任何点时半径等于负数
        temp.r =- 2 ;    
         break ;
     case   1 :         // 只有一个点时这个点就是圆心，半径等于0
        temp.centre = ce.p0;    
         break ;
     case   2 :         // 如果有两个点，圆心为两点连线的中点，半径为两点距离的一半
        temp.r = distance(ce.p0,ce.p1) / 2 ;    
        temp.centre.x = (ce.p0.x + ce.p1.x) / 2 ;
        temp.centre.y = (ce.p0.y + ce.p1.y) / 2 ;
         break ;
     default :     // 三个点的情况
        temp = circumcircle(ce);
         break ;
    }
     return  temp;
}


void  
caculate( int  t, int  ecount,TTriangle ce){
     int     i;
    circle = minCircle(ecount,ce);
     if (ecount == 3 )  return ;
     for (i = 0 ;i < t;i ++ ){
         // 当新的点到圆心的距离大于半径时，重新计算重新计算外接圆
         if (distance(point[i],circle.centre) > circle.r + EPS){    
             switch (ecount){
             case   0 : ce.p0 = point[i];     break ;
             case   1 : ce.p1 = point[i];     break ;
             case   2 :    ce.p2 = point[i];     break ;
            }
            caculate(t - 1 ,ecount + 1 ,ce);
        }
    }
}


void
main(){

     int  i;
    TTriangle ce;
    scanf( " %d " , & casenum);
     while (casenum -- ){
        scanf( " %d " , & pointnum);
         for (i = 0 ;i < pointnum;i ++ )
            scanf( " %lf%lf " , & point[i].x, & point[i].y);
        caculate(pointnum, 0 ,ce);
        printf( " %.2lf %.2lf %.2lf " ,circle.centre.x,circle.centre.y,circle.r);
    }
}