最小圆覆盖（经典算法【三点定圆）

刚刚学了一些基础的三维计算几何

接触到了增量法——一种看似暴力，实际睿智的算法

下面就是增量法在另一类问题上的展现

算法原文

问题描述

给定n个点，用一个最小的圆把这些点全部覆盖，求这个圆的圆心半径

算法

① 将所有点随机排布（这样可以保证算法的复杂度）

② 初始随意找到两点，设为P1,P2P1,P2，以P1P2P1P2为直径得到初始圆，设为C2C2（CiCi表示包含前i个点的最小圆）

③ 按顺序依次加点，设当前点为PiPi：若PiPi在当前圆Ci−1Ci−1内，Ci=Ci−1Ci=Ci−1；否则进入④

④ 一旦进入④，就说明我们需要在构造一个新的圆

显然插入点PiPi一定在新圆的边界上

简单的，我们直接以P1PiP1Pi为直径暂且得到一个CiCi

⑤ 新得到的CiCi不一定能包含1~i所有的点

我们找到不在CiCi中的一点Pj(j<i)Pj(j<i)，那么Pi,PjPi,Pj一定在更新的圆的边界上

现在为止，我们能确定有两个点（Pi,PjPi,Pj）在更新的圆的边界上

因此，简单的，我们直接以P1PjP1Pj为直径暂且得到一个CjCj

⑥ 同样的，新得到的CjCj不一定能包含1~j所有的点

我们找到不在CjCj中的一点Pk(k<j<i)Pk(k<j<i)，那么Pi,Pj,PkPi,Pj,Pk一定在更新的圆的边界上

现在我们能确定有三个点（Pi,Pj，PkPi,Pj，Pk）在更新的圆的边界上

因为三点确定一个圆，Pi,Pj，PkPi,Pj，Pk构成了新的圆，一定能覆盖前ii个点



上图中展示的就是一个简单维护过程

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

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

时间复杂度：O(N)

空间复杂度：O(N)

小细节

Q1.

过三点如何求圆？

A1.

先求叉积

若叉积为0，即三个点在同一直线，那么找到距离最远的一对点，以它们的连线为直径做圆即可；

若叉积不为0，即三个点不共线，那么就求三角形的外接圆

Q2.

如何求三角形外接圆？

A1.

设三个点(x1,y1),(x2,y2),(x3,y3)(x1,y1),(x2,y2),(x3,y3)

设过(x1,y1),(x2,y2)(x1,y1),(x2,y2)的直线l1l1方程为Ax+By=CAx+By=C

它的中点为(xmid,ymid)=(x1+x22,y1+y22)(xmid,ymid)=(x1+x22,y1+y22)

方法一：

l1l1中垂线方程为A1x+B1y=C1A1x+B1y=C1

则它的中垂线方程中：

A1=−B=x2−x1A1=−B=x2−x1

B1=A=y2−y1B1=A=y2−y1

C1=−B∗xmid+A∗ymid=(x22−x21)+(y22−y21)2C1=−B∗xmid+A∗ymid=(x22−x12)+(y22−y12)2

方法二：

如果我们把直线用“点+向量”的形式记录，求中垂线就简单多了

旋转原直线方向向量九十度即可

向量旋转

x’=xcosθ-ysinθ

y’=xsinθ+ycosθ

同理可以知道过(x1,y1),(x3,y3)(x1,y1),(x3,y3)的直线的中垂线的方程

于是这两条中垂线的交点就是圆心

Q3.

如何求两条直线交点？

A3.

方法一：

设两条直线为A1x+B1y=C1A1x+B1y=C1和A2x+B2y=C2A2x+B2y=C2

求delta=A1∗B2−A2∗B1delta=A1∗B2−A2∗B1（类似叉积）

如果delta=0delta=0，说明两直线平行；

若不等于0，则求交点：x=B2∗C1−B1∗C2delta，y=A1∗C2−A2∗C1deltax=B2∗C1−B1∗C2delta，y=A1∗C2−A2∗C1delta

方法二：

如果我们把直线用“点+向量”的形式记录，这个问题也很简单

node jiao(node p,node v,node q,node w)
//p+tv
//q+tw
{
node u=p-q;
double t=Cross(w,u)/Cross(v,w);
return p+v*t;
}

Code

思路基本上就是这样，代码实现可能不是这么显然

瞧好了您内

const double Pi=acos(-1.0);

int dcmp(double x)
{
if (fabs(x)<eps) return 0;
else if (x<0) return -1;
else return 1;
}

double lenth(node a) {return sqrt(Dot(a,a));}

node rotate(node a,double t)    //向量旋转
{
return node(a.x*cos(t)-a.y*sin(t),a.x*sin(t)+a.y*cos(t));
}

node jiao(node p,node v,node q,node w)
//p+tv
//q+tw
{
node u=p-q;
double t=Cross(w,u)/Cross(v,w);
return p+v*t;
}

node get_c(node a,node b,node c)
{
node p=(a+b)/2;    //ad中点
node q=(a+c)/2;    //ac中点
node v=rotate(b-a,Pi/2.0),w=rotate(c-a,Pi/2.0);   //中垂线的方向向量
if (dcmp(lenth(Cross(v,w)))==0)    //平行
{
if (dcmp(lenth(a-b)+lenth(b-c)-lenth(a-c))==0)
return (a+c)/2;
if (dcmp(lenth(b-a)+lenth(a-c)-lenth(b-c))==0)
return (b+c)/2;
if (dcmp(lenth(a-c)+lenth(c-b)-lenth(a-b))==0)
return (a+b)/2;
}
return jiao(p,v,q,w);
}

void min_circular
{
random_shuffle(P+1,P+n+1);    //随机化
c=P[1],r=0;
//c 圆心
//r 半径
for (int i=2;i<=n;i++)
if (dcmp(lenth(c-P[i])-r)>0)    //不在圆内
{
c=P[i],r=0;
for (int j=1;j<i;j++)
if (dcmp(Lenth(c-P[j])-r)>0)
{
c=(P[i]+P[j])/2.0;
r=lenth(c-P[i]);
for (int k=1;k<j;k++)
if (dcmp(lenth(c-P[k])-r)>0)
{
c=get_c(P[i],P[j],P[k]);
r=lenth(c-P[i]);
}
}
}
}