多点寻找最小外接圆

#include <ansi_c.h>
 
  
#define LL long long  
#define PI (acos(-1.0))  
#define EPS (1e-10)  
  

  
typedef struct 
{  
    double x,y,a;  
}P;
P p[1100],cp[1100]; 
  
struct L  
{  
    P p1,p2;  
} line[50];  
  
double X_Mul(P a1,P a2,P b1,P b2)  
{  
    P v1 = {a2.x-a1.x,a2.y-a1.y},v2 = {b2.x-b1.x,b2.y-b1.y};  
    return v1.x*v2.y - v1.y*v2.x;  
}  

//¼ÆËãÁ½µã¼ä¾à
double Cal_Point_Dis(P p1,P p2)  
{  
    return sqrt((p1.x-p2.x)*(p1.x-p2.x) (p1.y-p2.y)*(p1.y-p2.y));  
}  

//¼ÆËãÁ½µã¼ä¾àµÄƽ·½
double Cal_Point_Dis_Pow_2(P p1,P p2)  
{  
    return (p1.x-p2.x)*(p1.x-p2.x) (p1.y-p2.y)*(p1.y-p2.y);  
}  
  
P Cal_Segment_Cross_Point(P a1,P a2,P b1,P b2)  
{  
    double t = X_Mul(b1,b2,b1,a1) / X_Mul(a1,a2,b1,b2);  
    P cp = {a1.x (a2.x-a1.x)*t, a1.y (a2.y-a1.y)*t};  
    return cp;  
}  
  
//¹æ·¶Ïཻ  
int Is_Seg_Cross_Without_Point(P a1,P a2,P b1,P b2)  
{  
    double xm1 = X_Mul(a1,a2,a1,b1);  
    double xm2 = X_Mul(a1,a2,a1,b2);  
    double xm3 = X_Mul(b1,b2,b1,a1);  
    double xm4 = X_Mul(b1,b2,b1,a2);  
  
    return (xm1*xm2 < (-EPS) && xm3*xm4 < (-EPS));  
}  
  
//ÏòÁ¿abÓëXÖáÕý·½ÏòµÄ¼Ð½Ç  
double Cal_Angle(P t,P p)  
{  
    return ((t.x-p.x)*(t.x-p.x) 1 - (t.x-p.x-1)*(t.x-p.x-1))/(2*sqrt((t.x-p.x)*(t.x-p.x) (t.y-p.y)*(t.y-p.y)));  
}  
  
//¼ÆËã ¡Ïb2.a.b1 µÄ´óС  
double Cal_Angle_Three_Point(P a,P b1,P b2)  
{  
    double l1 = Cal_Point_Dis_Pow_2(b1,a);  
    double l2 = Cal_Point_Dis_Pow_2(b2,a);  
    double l3 = Cal_Point_Dis_Pow_2(b1,b2);  
  
    return acos( (l1 l2-l3)/(2*sqrt(l1*l2)) );  
}  
  
int cmp_angle(P p1,P p2)  
{  
    return (p1.a < p2.a || (fabs(p1.a-p2.a) < EPS && p1.y < p2.y) ||(fabs(p1.a-p2.a) < EPS && fabs(p1.y-p2.y) < EPS && p1.x < p2.x)   );  
}  
  
//pΪµã¼¯  Ó¦¸Ã±£Ö¤ÖÁÉÙÓÐÈýµã²»¹²Ïß  
//n Ϊµã¼¯´óС  
//·µ»ØֵΪ͹°üÉϵ㼯ÉϵĴóС  
//͹°üÉϵĵã´æ´¢ÔÚcpÖÐ  
int Creat_Convex_Hull(P *p,int n)  
{  
    int i,top;  
    P re; //re Ϊ½¨Á¢¼«½ÇÅÅÐòʱµÄ²Î¿¼µã  ÏÂ×óµã  
  
    for(re = p[0],i = 1; i < n; i)  
    {  
        if(re.y > p[i].y || (fabs(re.y-p[i].y) < EPS && re.x > p[i].x))  
            re = p[i];  
    }  
  
    for(i = 0; i < n; i)  
    {  
        if(p[i].x == re.x && p[i].y == re.y)  
        {  
            p[i].a = -2;  
        }  
        else p[i].a = Cal_Angle(re,p[i]);  
    }  
  
    //sort(p,p n,cmp_angle);  
//Sort(input_array, total_num, ANALYSIS_SORT_ASCENDING, output_array)
  
    top =0 ;  
    cp[top ] = p[0];  
    cp[top ] = p[1];  
  
    for(i = 2 ; i < n;)  
    {  
        if(top < 2 || X_Mul(cp[top-2],cp[top-1],cp[top-2],p[i]) > EPS)  
        {  
            cp[top ] = p[i ];  
        }  
                    else  
        {  
            --top;  
        }  
    }  
  
    return top;  
}  
  
int Is_Seg_Parallel(P a1,P a2,P b1,P b2)  
{  
    double xm1 = X_Mul(a1,a2,a1,b1);  
    double xm2 = X_Mul(a1,a2,a1,b2);  
  
    return (fabs(xm1) < EPS && fabs(xm2) < EPS);  
}  
  
int Point_In_Seg(P m,P a1,P a2)  
{  
    double l1 = Cal_Point_Dis(m,a1);  
    double l2 = Cal_Point_Dis(m,a2);  
    double l3 = Cal_Point_Dis(a1,a2);  
  
    return (fabs(l1 l2-l3) < EPS);  
}  
  
//¼ÆËãÈý½ÇÐÎÍâ½ÓÔ²Ô²ÐÄ  
P Cal_Triangle_Circumcircle_Center(P a, P b, P c)  
{ 
    P mp1 = {(b.x a.x)/2,(b.y a.y)/2}, mp2 = {(b.x c.x)/2,(b.y c.y)/2};  
    P v1 = {a.y-b.y, b.x-a.x}, v2 = {c.y-b.y, b.x-c.x};  
  
    P p1 = {mp1.x v1.x, mp1.y v1.y}, p2 = {mp2.x v2.x,mp2.y v2.y};  
  
    return Cal_Segment_Cross_Point(mp1,p1,mp2,p2);  
}  
  
int Is_Acute_Triangle(P p1,P p2,P p3)  
{  
    //Èýµã¹²Ïß  
    if(fabs(X_Mul(p1,p2,p1,p3)) < EPS)  
    {  
        return 0;  
    }  
  
    double a = Cal_Angle_Three_Point(p1,p2,p3);  
    if(a > PI || fabs(a-PI) < EPS)  
    {  
        return 0;  
    }  
    a = Cal_Angle_Three_Point(p2,p1,p3);  
    if(a > PI || fabs(a-PI) < EPS)  
    {  
        return 0;  
    }  
    a = Cal_Angle_Three_Point(p3,p1,p2);  
    if(a > PI || fabs(a-PI) < EPS)  
    {  
        return 0;  
    }  
    return 1;  
}  

//ÅжϵãÊÇ·ñÔÚÔ²ÄÚ
int Is_In_Circle(P center,double rad,P p)  
{  
    double dis = Cal_Point_Dis(center,p);  
    return (dis < rad || fabs(dis-rad) < EPS);  
}  

//¼ÆËãÖеã
P Cal_Seg_Mid_Point(P p2, P p1)  
{  
    P mp = {(p2.x p1.x) / 2, (p1.y p2.y) / 2};  
    return mp;  
}  

//¼ÆËã×îСÍâ½ÓÔ²
void Cal_Min_Circle(P *p, int n)  
{  
    if(n == 0)  
        return ; 
//Ö»ÓÐÒ»¸öµã
    if(n == 1)  
    { 
        printf("%.2lf %.2lf %.2lf\n", p[0].x, p[0].y, 0.00);  
        return ;  
    }
//Ö»ÓÐÁ½¸öµã
    if(n == 2)  
    {  
        printf("%.2lf %.2lf %.2lf\n",(p[1].x p[0].x) / 2,(p[0].y p[1].y) / 2, Cal_Point_Dis(p[0], p[1]) / 2);  
        return ;  
    }  
  
//³¬¹ýÁ¬¸öµã,ÏȽ«Á½µãÖÐÐÄÉèÖÃΪԲÐÄ
    P center = Cal_Seg_Mid_Point(p[0], p[1]);
P tc1;  
  
//¼ÆËã°ë¾¶
    double rad = Cal_Point_Dis(p[0],center);
double dis = -1, temp = 0.0, tra1;  
    int a = 0, b = 1, c = 2;
int i = 0, site = 0, d = 0, s1 = 0, s2 = 0, tp = 0;  
  
//¼ÆËã×î´ó°ë¾¶¼°±êʶ
    for(dis = -1, site = 0, i = 0; i < n; i)  
    {  
        temp = Cal_Point_Dis(center, p[i]);  
        if(temp > dis)  
        {  
            dis = temp;  
            site = i;  
        }  
    }  
  
//Èôijµã²»ÔÚÔ²ÉÏ
    while( Is_In_Circle(center, rad, p[site]) == 0 )  
    {  
        d = site;  
       
//¼ÆËãµã¾à
        double l1 = Cal_Point_Dis(p[a], p[d]);  
        double l2 = Cal_Point_Dis(p[b], p[d]);  
        double l3 = Cal_Point_Dis(p[c], p[d]);  
  
//aµã¾àÀë×îÔ¶
        if((l1 > l2 || fabs(l1 - l2) < EPS) && (l1 > l3 || fabs(l1 - l3) < EPS))  
        {  
            s1 = a, s2 = d, tp = b;  
        } 
//bµã¾àÀë×îÔ¶
        else if((l2 > l1 || fabs(l1 - l2) < EPS) && (l2 > l3 || fabs(l2 - l3) < EPS))  
        {  
            s1 = b, s2 = d, tp = c;  
        }  
        else  
        {  
            s1 = c, s2 = d, tp = a;  
        }  
  
//¼ÆËãÖеã
        center = Cal_Seg_Mid_Point(p[s1], p[s2]);  
//¼ÆËã¼ä¾à
        rad = Cal_Point_Dis(center, p[s1]);  
  
//ÆäÓàÈý¸öµãÈÎÒ»¸öµã²»ÔÚÔ²ÉÏ»òÕßÔ²ÄÚ
        if( Is_In_Circle(center,rad,p[a]) == 0 || Is_In_Circle(center,rad,p[b]) == 0 || Is_In_Circle(center,rad,p[c]) == 0 )  
        {  
//¼ÆËãÈý½ÇÐÎÍâ½çÔ²ÐÎ
            center = Cal_Triangle_Circumcircle_Center(p[a],p[c],p[d]);  
            rad = Cal_Point_Dis(p[d],center);  
            s1 = a,s2 = c,tp = d;  
  
//ÅжϵãÊÇ·ñÔÚÔ²ÄÚ  ²»ÔÚÔ²ÄÚ
            if(Is_In_Circle(center,rad,p[b]) == 0)  
            {  
                center = Cal_Triangle_Circumcircle_Center(p[a],p[b],p[d]);  
                rad = Cal_Point_Dis(p[d],center);  
                s1 = a,s2 = b,tp = d;  
            }  
            else  
            {  
                tc1 = Cal_Triangle_Circumcircle_Center(p[a],p[b],p[d]);  
                tra1 = Cal_Point_Dis(p[d],tc1);  
  
                if(tra1 < rad && Is_In_Circle(tc1,tra1,p[c]))  
                {  
                    rad = tra1,center = tc1;  
                    s1 = a,s2 = b,tp = d;  
                }  
            }  
  
            if(Is_In_Circle(center,rad,p[c]) == 0)  
            {  
                center = Cal_Triangle_Circumcircle_Center(p[c],p[b],p[d]);  
                rad = Cal_Point_Dis(center,p[d]);  
                s1 = c,s2 = b,tp = d;  
            }  
            else  
            {  
                tc1 = Cal_Triangle_Circumcircle_Center(p[c],p[b],p[d]);  
                tra1 = Cal_Point_Dis(p[d],tc1);  
                if(tra1 < rad && Is_In_Circle(tc1,tra1,p[a]))  
                {  
                    rad = tra1,center = tc1;  
                    s1 = b,s2 = c,tp = d;  
                }  
            }  
        }  
  
        a = s1, b = s2, c = tp;  
  
//¼ÆËã×î´ó°ë¾¶¼°±êʶ 
        for(dis = -1, site = 0, i = 0; i < n; i)  
        {  
            temp = Cal_Point_Dis(center, p[i]);  
            if(temp > dis)  
            {  
                dis = temp;  
                site = i;  
            }  
        }  
    } 

    printf("%.2f %.2f %.2f\n",center.x, center.y, rad);  
}  

//Ö÷º¯Êý
int main()  
{  
    int i = 0, n = 0;  
    while(scanf("%d", &n) && n)  
    {  
//ÔØÈëÊý¾Ý
        for(i = 0; i < n; i)  
        {  
            scanf("%lf %lf",&p[i].x,&p[i].y);  
        }

//¼ÆËã×îСÍâ½ÓÔ²
        Cal_Min_Circle(p,n);  
    }  
}