最小半径圆

题目描述：

平面上有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);
    }
}