凸包和最小覆盖圆问题

给出平面上的一些点，求覆盖这些点的最小圆。

具体问题可以见hdu 2215。

具体解法是，先求凸包，然后枚举凸包上任意3个点，若枚举的三个点构成钝角三角形，则最大半径为最长边的一半

否则，半径r=a*b*c/(4*s)其中s是面积，具体面积可以用叉乘求得，s=(向量a叉乘向量b)的绝对值的一半。

View Code

  1 #include<iostream>
  2 #include<string>
  3 #include<cmath>
  4 #include<algorithm>
  5 using namespace std;
  6 
  7 struct node
  8 {
  9     double x;
 10     double y;
 11 };
 12 
 13 node point[110];
 14 int n;
 15 
 16 double dist(node a,node b)
 17 {
 18     return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
 19 }
 20 
 21 double mul(node p,node a,node b)
 22 {
 23     return (a.x-p.x)*(b.y-p.y)-(b.x-p.x)*(a.y-p.y);
 24 }
 25 
 26 int cmp(node a,node b)
 27 {
 28     double d=mul(point[0],a,b);
 29     if(d>0)
 30         return 1;
 31     if(d==0 && dist(point[0],a)<=dist(point[0],b))
 32         return 1;
 33     return 0;
 34 }
 35 
 36 int stack[110],top;
 37 
 38 void tubao()
 39 {
 40     stack[0]=0;
 41     stack[1]=1;
 42     top=2;
 43     for(int i=2;i<n;i++)
 44     {
 45         while(top>1 && mul(point[stack[top-2]],point[stack[top-1]],point[i])<=0)
 46             top--;
 47         stack[top++]=i;
 48     }
 49 }
 50 
 51 double max(double a,double b)
 52 {
 53     return a>b?a:b;
 54 }
 55 
 56 int main()
 57 {
 58     int i,j,k;
 59     double a,b,c,d;
 60     freopen("in.txt","r",stdin);
 61     while(scanf("%d",&n)==1 && n)
 62     {
 63         for(i=0;i<n;i++)
 64             scanf("%lf%lf",&point[i].x,&point[i].y);
 65         if(n==1)
 66         {
 67             printf("0.50\n");
 68             continue;
 69         }
 70         /**********************************///凸包
 71         k=0;
 72         for(i=1;i<n;i++)
 73         {
 74             if(point[i].x<point[k].x)
 75                 k=i;
 76             else if(point[i].x==point[k].x && point[i].y<point[k].y)
 77                 k=i;
 78         }
 79         swap(point[0],point[k]);
 80         sort(point+1,point+n,cmp);
 81         tubao();
 82         /**********************************/
 83         double ans=0,temp,ares;
 84         if(top<3)
 85         {
 86             a=dist(point[stack[0]],point[stack[1]]);
 87             printf("%0.2lf\n",a/2+0.5);
 88             continue;
 89         }
 90         for(i=0;i<top;i++)
 91         {
 92             for(j=i+1;j<top;j++)
 93             {
 94                 for(k=j+1;k<top;k++)
 95                 {
 96                     a=dist(point[stack[i]],point[stack[j]]);
 97                     b=dist(point[stack[i]],point[stack[k]]);
 98                     c=dist(point[stack[j]],point[stack[k]]);
 99                     d=max(max(a,b),c);
100                     ares=fabs(mul(point[stack[i]],point[stack[j]],point[stack[k]])); //面积的2倍
101                     if(a==d)
102                     {
103                         a=c;
104                         c=d;
105                     }
106                     else if(b==d)
107                     {
108                         b=c;
109                         c=d;
110                     }
111                     if(a*a+b*b<c*c) //钝角三角形
112                     {
113                         temp=d/2;
114                     }
115                     else
116                     {
117                         temp=a*b*c/(2*ares);
118                     }
119                     ans=max(ans,temp);
120                 }
121             }
122         }
123         printf("%.2lf\n",ans+0.5);
124     }
125     return 0;
126 }