Minimal Circle

http://icpc.upc.edu.cn/problem.php?cid=1440&pid=3

第 1 步. 在点集中任取 3 点 A , B , C .

第 2 步. 作一个包 含 A , B , C 三点的最小 圆. 圆周可 能通过这 3 点( 如图 1 所示) , 也 可能只通过
其中两点, 但包含第 3 点. 后一种情况圆周上的两点一定是 位于圆的一条直径的两端.

第 3 步. 在点集中找 出距离第 2 步所建圆圆心最远的点 D . 若 D 点已在圆内或圆周上, 则该
圆即为所求的圆, 算法结束. 否则, 执行第 4 步.

第 4 步. 在 A , B , C , D 中 选 3 个点, 使 由它们生成的一个 包含这 4 点的圆为最 小. 这 3 点成
为新的 A , B 和 C , 返回执 行第 2 步 .

若在第 4 步生成的圆的圆周只通过 A , B , C , D 中的两点, 则圆周上的两点取成新的 A 和 B ,
从另两点中任取一点作为新的 C .

但是我读的时候感觉其中最关键的第四步说的最为模糊了。

下面是我自己的一些关于求四边形即四个点的最小圆覆盖的一些想法。

1. 若此四边形中存在大于等于180度的内角时，则可看作求三角形的最小覆盖圆（钝角三角形为以长边做直径的圆，否则为该三角形外接圆）。

2. 若存在两个锐角，两个非锐角，则必为两个锐角顶点连线为直径的圆。

3. 若有三个锐角，一个非锐角，则必为三个锐角顶点组成的三角形的外接圆。

4. 若有四个直角，则任选其中三点组成三角形求其外接圆即可。

解决了第四步，剩下的就是按照上述算法求解了。
 

 

题意：

给你 N 个点，求最小的覆盖所有点的圆。

输出圆的坐标和半径。

 

 

算法：计算几何 

                         向量叉积，向量点积，两点间的距离，三角形外心

                        

                        一点定圆，两点定圆问题。

 

 

思路： 

 直接贴我在网上看的资料吧，感觉自己不能比他解释的更好： 资料出处

这里介绍的算法是，先任意选取两个点，以这两个点的连线为直径作圆。再以此判断剩余的点，看它们是否都在圆内（或圆上），如果都在，说明这个圆已经找到。如果没有都在：假设我们用的最开始的两个点为p1,p[2]，并且找到的第一个不在圆内（或圆上）的点为p[i]，于是我们用这个点p[i]去寻找覆盖p1到p[i-1]的最小覆盖圆。

那么，过确定点p[i]的从p1到p[i-1]的最小覆盖圆应该如何求呢？

我们先用p1和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]更新圆即可。

于是，这个问题就被转化为若干个子问题来求解了。

由于三个点确定一个圆，我们的过程大致上做的是从没有确定点，到有一个确定点，再到有两个确定点，再到有三个确定点来求圆的工作。

关于正确性的证明以及复杂度的计算这里就不介绍了，可以去看完整的算法介绍：ACM 计算几何 最小圆覆盖算法

 

关于细节方面：PS:细节部分推荐路过的 acmer 们搞自己的模板比较好。

 

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=-Bmidx+Amidy=((x2^2-x1^2)+(y2^2-y1^2))/2；
同理可以知道过(x1,y1),(x3,y3)的直线的中垂线的方程。
于是这两条中垂线的交点就是圆心。

c.如何求两条直线交点？
设两条直线为A1x+B1y=C1和A2x+B2y=C2。
设一个变量det=A1B2-A2B1；
如果det=0，说明两直线平行；若不等于0，则求交点：x=(B2C1 -B1C2)/det，y=(A1C2-A2C1)/det；

d.于是木有了。。

 

关于细节的处理更改：

 求三角形的外接圆：推荐资料

三角形垂心， 重心，外心坐标公式

 

P  = （a^2 + b^2 + c^2) / 2;

Q = 1 / (1 / (P - a^2) + 1 / (P - b^2) + 1 / （P - c^2));

λ1=Q/(P-a^2),λ2=Q/(P-b^2),λ3=Q/(P-c^2)；(注:λ1+λ2+λ3=1) 

垂心:H(x,y)=λ1*A(x,y)+λ2*B(x,y)+λ3*C(x,y), 

重心:G(x,y)=1/3*A(x,y)+1/3*B(x,y)+1/3*C(x,y), 

外心:O(x,y)=(1-λ1)/2*A(x,y)+(1-λ2)/2*B(x,y)+(1-λ3)/2*C(x,y);

 

 

模板化：外接圆的圆心就是各边中垂线的交点也就是所谓的外心。

//三角形的外接圆圆心,先前已经保证了 A,B,C 三点不共线
void TriangleCircumCircle(Point A, Point B, Point C)
{
    Vector a = C-B, b = A-C, c = B-A;
    double da = pow(dist(a), 2), db = pow(dist(b), 2), dc = pow(dist(c), 2);
    double px = (da + db + dc) / 2.0;
    double Q = 1 / (1/(px-da) + 1/(px-db) + 1/(px-dc));
    //t1+t2+t3 = 1
    double t1 = Q/(px-da), t2 = Q/(px-db), t3 = Q/(px-dc);
    Point O = A*((1-t1)/2.0) + B*((1-t2)/2.0)  + C*((1-t3)/2.0)  ; //外心

    //Point H = A*t1 + B*t2 + C*t3; //垂心
    //Point G = A*(1/3) + B*(1/3) + C*(1/3); //重心

    center = O;
}

参考资料：

楼上的思路出处：点击打开链接 
这也是我第一次去 德问 感觉这个社区超级棒 ！

楼上的思路证明：点击打开链接 个人感觉已经很清楚了，不用看这个证明。

三角形的五心问题：百度百科之三角形的五心

关于三角形五心的公式：ACM三角形五心公式 自己还没有证明，等待自己的模板ing...

 

PS: 咋一看，一篇东拼西凑的题解

 

开始找思路的时候看到的一篇题解，当时没看懂，现在看来思路和这个一样，代码比较简洁，这里推荐一下：点击打开链接

比赛时错误的分析：

 

凸包的直径就是最小覆盖圆的直径，明显错了，一个等边三角形的情况就直接否决。

但是提供了一个思路，枚举的时候如果点非常多，是否可以从直径开始找最小覆盖圆，效率是否会更优？

考虑到这里的点非常少，确定凸包再用旋转卡壳确定凸包直径显得复杂了点。

那么这里有一个问题：确定凸包的直径的端点是否一定在最小覆盖圆的圆圈上。。。。

 

#include <iostream>
#include <cstring>
#include <cstdlib>
#include <cstdio>
#include <queue>
#include <cmath>
#include <algorithm>
 
#define LL long long
#define PI (acos(-1.0))
#define EPS (1e-10)
 
using namespace std;
 
struct P
{
    double x,y,a;
}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;
}
 
//规范相交
bool 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)) );
}
bool 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);
    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;
}
 
bool 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);
}
bool 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);
}
bool Is_Acute_Triangle(P p1,P p2,P p3)
{
    //三点共线
    if(fabs(X_Mul(p1,p2,p1,p3)) < EPS)
    {
        return false;
    }
    double a = Cal_Angle_Three_Point(p1,p2,p3);
    if(a > PI || fabs(a-PI) < EPS)
    {
        return false;
    }
    a = Cal_Angle_Three_Point(p2,p1,p3);
    if(a > PI || fabs(a-PI) < EPS)
    {
        return false;
    }
    a = Cal_Angle_Three_Point(p3,p1,p2);
    if(a > PI || fabs(a-PI) < EPS)
    {
        return false;
    }
    return true;
}
bool 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]),tc1;
    double dis,temp,rad = Cal_Point_Dis(p[0],center),tra1;
    int i,site,a = 0,b = 1,c = 2,d,s1,s2,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;
        }
    }
 
    while( Is_In_Circle(center,rad,p[site]) == false )
    {
        d = site;
       // printf("a = %d b = %d c = %d d = %d\n",a,b,c,d);
 
        //printf("x = %.2lf  y = %.2lf\n",center.x,center.y);
 
       // printf("rad = %.2lf\n\n",rad);
       // getchar();
 
        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]);
 
        if((l1 > l2 || fabs(l1-l2) < EPS) && (l1 > l3 || fabs(l1-l3) < EPS))
        {
            s1 = a,s2 = d,tp = 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]) == false || Is_In_Circle(center,rad,p[b]) == false || Is_In_Circle(center,rad,p[c]) == false )
        {
            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]) == false)
            {
                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]) == false)
            {
                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,n;
    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);
    }
}

第 1 步. 在点集中任取 3 点 A , B , C .

第 2 步. 作一个包 含 A , B , C 三点的最小 圆. 圆周可 能通过这 3 点( 如图 1 所示) , 也 可能只通过
其中两点, 但包含第 3 点. 后一种情况圆周上的两点一定是 位于圆的一条直径的两端.

第 3 步. 在点集中找 出距离第 2 步所建圆圆心最远的点 D . 若 D 点已在圆内或圆周上, 则该
圆即为所求的圆, 算法结束. 否则, 执行第 4 步.

第 4 步. 在 A , B , C , D 中 选 3 个点, 使 由它们生成的一个 包含这 4 点的圆为最 小. 这 3 点成
为新的 A , B 和 C , 返回执 行第 2 步 .

若在第 4 步生成的圆的圆周只通过 A , B , C , D 中的两点, 则圆周上的两点取成新的 A 和 B ,
从另两点中任取一点作为新的 C .

但是我读的时候感觉其中最关键的第四步说的最为模糊了。

下面是我自己的一些关于求四边形即四个点的最小圆覆盖的一些想法。

1. 若此四边形中存在大于等于180度的内角时，则可看作求三角形的最小覆盖圆（钝角三角形为以长边做直径的圆，否则为该三角形外接圆）。

2. 若存在两个锐角，两个非锐角，则必为两个锐角顶点连线为直径的圆。

3. 若有三个锐角，一个非锐角，则必为三个锐角顶点组成的三角形的外接圆。

4. 若有四个直角，则任选其中三点组成三角形求其外接圆即可。

解决了第四步，剩下的就是按照上述算法求解了。