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

刚刚学了一些基础的三维计算几何
接触到了增量法——一种看似暴力，实际睿智的算法
下面就是增量法在另一类问题上的展现

算法原文

问题描述

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

算法

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

② 初始随意找到两点，设为$P_1,P_2$，以$P_1{P}_2$为直径得到初始圆，设为$C_2$（$C_i$表示包含前i个点的最小圆）

③ 按顺序依次加点，设当前点为$P_i$：若$P_i$在当前圆$C_{i-1}$内，$C_i=C_{i-1}$；否则进入④

④ 一旦进入④，就说明我们需要在构造一个新的圆
显然插入点$P_i$一定在新圆的边界上
简单的，我们直接以$P_1{P}_i$为直径暂且得到一个$C_i$

⑤ 新得到的$C_i$不一定能包含1~i所有的点
我们找到不在$C_i$中的一点$P_j(j<i)$，那么$P_i,P_j$一定在更新的圆的边界上
现在为止，我们能确定有两个点（$P_i,P_j$）在更新的圆的边界上
因此，简单的，我们直接以$P_i{P}_j$为直径暂且得到一个$C_j$

⑥ 同样的，新得到的$C_j$不一定能包含1~j所有的点
我们找到不在$C_j$中的一点$P_k(k<j<i)$，那么$P_i,P_j,P_k$一定在更新的圆的边界上

现在我们能确定有三个点（$P_i,P_j，P_k$）在更新的圆的边界上
因为三点确定一个圆，$P_i,P_j，P_k$构成了新的圆，一定能覆盖前$i$个点

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

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

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

时间复杂度：O(N)

空间复杂度：O(N)

小细节

Q1.

过三点如何求圆？

A1.

先求叉积

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

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

Q2.

如何求三角形外接圆？

A1.

设三个点$(x_1,y_1),(x_2,y_2),(x_3,y_3)$
设过$(x_1,y_1),(x_2,y_2)$的直线$l_1$方程为$Ax+By=C$
它的中点为$(x_{mid},y_{mid})=(\frac{x_1+x_2}{2},\frac{y_1+y_2}{2})$

方法一：
$l_1$中垂线方程为$A_{1}x+B_{1}y=C_1$
则它的中垂线方程中：
$A_1=-B=x_2-x_1$
$B_1=A=y_2-y_1$
$C_1=-B\ast x_{mid}+A\ast y_{mid}=\frac{(x_2^2-x_1^2)+(y_2^2-y_1^2)}{2}$

方法二：
如果我们把直线用“点+向量”的形式记录，求中垂线就简单多了
旋转原直线方向向量九十度即可

向量旋转
x’=xcosθ-ysinθ
y’=xsinθ+ycosθ

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

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

Q3.

如何求两条直线交点？

A3.

方法一：
设两条直线为$A_{1}x+B_{1}y=C_1$和$A_{2}x+B_{2}y=C_2$
求$delta=A_1\ast B_2-A_2\ast B_1$（类似叉积）
如果$delta=0$，说明两直线平行；
若不等于0，则求交点：$x=\frac{B2\ast C1-B1\ast C2}{delta}，y=\frac{A1\ast C2-A2\ast C1}{delta}$

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

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]);
						}
				}
		} 
}