最小覆蓋圓演算法

標籤：style   blog   http   color   os   io   for   art   

最小圓覆蓋，很經典的問題。題目大概是，平面上n個點，求一個半徑最小的圓，能夠覆蓋所有的點。

 

演算法有點難懂，於是講講我的理解。

如果要求一個最小覆蓋圓，這個圓至少要由三個點確定。有一種演算法就是任意取三個點作圓，然後判斷距離圓心最遠的點是否在圓內，若在，則完成；若不在則用最遠點更新這個圓。

這裡介紹的演算法是，先任意選取兩個點，以這兩個點的連線為直徑作圓。再以此判斷剩餘的點，看它們是否都在圓內（或圓上），如果都在，說明這個圓已經找到。如果沒有都在：假設我們用的最開始的兩個點為p[1],p[2]，並且找到的第一個不在圓內（或圓上）的點為p[i]，於是我們用這個點p[i]去尋找覆蓋p[1]到p[i-1]的最小覆蓋圓。

那麼，過確定點p[i]的從p[1]到p[i-1]的最小覆蓋圓應該如何求呢？

我們先用p[1]和p[i]做圓，再從2到i-1判斷是否有點不在這個圓上，如果都在，則說明已經找到覆蓋1到i-1的圓。如果沒有都在：假設我們找到第一個不在這個圓上的點為p[j]，於是我們用兩個已知點p[j]與p[i]去找覆蓋1到j-1的最小覆蓋圓。

而對於兩個已知點p[j]與p[i]求最小覆蓋圓，只要從1到j-1中，第k個點求過p[k],p[j],p[i]三個點的圓，再判斷k+1到j-1是否都在圓上，若都在，說明找到圓；若有不在的，則再用新的點p[k]更新圓即可。

於是，這個問題就被轉化為若干個子問題來求解了。

由於三個點確定一個圓，我們的過程大致上做的是從沒有確定點，到有一個確定點，再到有兩個確定點，再到有三個確定點來求圓的工作。

關於正確性的證明以及複雜度的計算這裡就不介紹了，可以去看完整的演算法介紹：http://wenku.baidu.com/view/162699d63186bceb19e8bbe6.html

恩。關於細節方面。

a.通過三個點如何求圓？

   先求叉積。

   若叉積為0，即三個點在同一直線，那麼找到距離最遠的一對點，以它們的連線為直徑做圓即可；

   若叉積不為0，即三個點不共線，那麼就是第二個問題，如何求三角形的外接圓？

b.如何求三角形外接圓？

   假設三個點(x1,y1),(x2,y2),(x3,y3)；

   設過(x1,y1),(x2,y2)的直線l1方程為Ax+By=C，它的中點為(midx,midy)=((x1+x2)/2,(y1+y2)/2)，l1中垂線方程為A1x+B1y=C1；則它的中垂線方程中A1=-B=x2-x1，B1=A=y2-y1,C1=-B*midx+A*midy=((x2^2-x1^2)+(y2^2-y1^2))/2；

   同理可以知道過(x1,y1),(x3,y3)的直線的中垂線的方程。

   於是這兩條中垂線的交點就是圓心。

c.如何求兩條直線交點？

   設兩條直線為A1x+B1y=C1和A2x+B2y=C2。

   設一個變數det=A1*B2-A2*B1；

   如果det=0，說明兩直線平行；若不等於0，則求交點：x=(B2*C1 -B1*C2)/det，y=(A1*C2-A2*C1)/det；

d.於是木有了。。

 1 #include<stdio.h> 2 #include<math.h> 3 struct   TPoint   4 {   5          double x,y;   6 }; 7 TPoint a[1005],d; 8 double r; 9 10 double   distance(TPoint   p1,   TPoint   p2)   //兩點間距離11 {  12     return (sqrt((p1.x-p2.x)*(p1.x -p2.x)+(p1.y-p2.y)*(p1.y-p2.y)));      13 }14 double multiply(TPoint   p1,   TPoint   p2,   TPoint   p0)  15 {  16     return   ((p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y));          17 }  18 void MiniDiscWith2Point(TPoint   p,TPoint   q,int   n)19 {20     d.x=(p.x+q.x)/2.0;21     d.y=(p.y+q.y)/2.0;22     r=distance(p,q)/2;23     int k;24     double c1,c2,t1,t2,t3;25     for(k=1;k<=n;k++)26     {27         if(distance(d,a[k])<=r)continue;28         if(multiply(p,q,a[k])!=0.0)29         {30             c1=(p.x*p.x+p.y*p.y-q.x*q.x-q.y*q.y)/2.0;31             c2=(p.x*p.x+p.y*p.y-a[k].x*a[k].x-a[k].y*a[k].y)/2.0;32 33             d.x=(c1*(p.y-a[k].y)-c2*(p.y-q.y))/((p.x-q.x)*(p.y-a[k].y)-(p.x-a[k].x)*(p.y-q.y));34             d.y=(c1*(p.x-a[k].x)-c2*(p.x-q.x))/((p.y-q.y)*(p.x-a[k].x)-(p.y-a[k].y)*(p.x-q.x));35             r=distance(d,a[k]);36         }37         else38         {39             t1=distance(p,q);40             t2=distance(q,a[k]);41             t3=distance(p,a[k]);42             if(t1>=t2&&t1>=t3)43             {d.x=(p.x+q.x)/2.0;d.y=(p.y+q.y)/2.0;r=distance(p,q)/2.0;}44             else if(t2>=t1&&t2>=t3)45             {d.x=(a[k].x+q.x)/2.0;d.y=(a[k].y+q.y)/2.0;r=distance(a[k],q)/2.0;}46             else47             {d.x=(a[k].x+p.x)/2.0;d.y=(a[k].y+p.y)/2.0;r=distance(a[k],p)/2.0;}48         }49     }50 }51 52 void MiniDiscWithPoint(TPoint   pi,int   n)53 {54     d.x=(pi.x+a[1].x)/2.0;55     d.y=(pi.y+a[1].y)/2.0;56     r=distance(pi,a[1])/2.0;57     int j;58     for(j=2;j<=n;j++)59     {60     if(distance(d,a[j])<=r)continue;61         else62         {63             MiniDiscWith2Point(pi,a[j],j-1);64         }65     }66 }67 int main()68 {69     int i,n;70     while(scanf("%d",&n)&&n)71     {72         for(i=1;i<=n;i++)73         {74             scanf("%lf %lf",&a[i].x,&a[i].y);75         }76         if(n==1)77         { printf("%.2lf %.2lf 0.00\n",a[1].x,a[1].y);continue;}78         r=distance(a[1],a[2])/2.0;79         d.x=(a[1].x+a[2].x)/2.0;80         d.y=(a[1].y+a[2].y)/2.0;81         for(i=3;i<=n;i++)82         {83             if(distance(d,a[i])<=r)continue;84             else85             MiniDiscWithPoint(a[i],i-1);86         }87         printf("%.2lf %.2lf %.2lf\n",d.x,d.y,r);88     }89     return 0;90 }

View Code