NX二次开发Welzl算法求点云的最小外接圆（最小覆盖圆）

一、概述

        最小圆问题（也称为最小覆盖圆问题、边界圆问题、最小边界圆问题、最小封闭圆问题）是可以在线性时间复杂度内求出覆盖n个点的最小圆。在一次学习的过程中遇到了这个问题，网上有好多种现成的方法，直接抄写代码即可，但是个人认为有必要搞清楚一下原理，所以抽取半天时间进行整理，个人认为如果一味的抄袭只会让自己的思路越来越窄，我们可以抄，但是一定要清楚其必要的原理，我们可以不用，但不能不会。

二、思路

2.1Welzl’s算法：

假设我们已经知道第P(n-1)个点的最小外接圆D内或者D上，那么第Pn个点有两种情况

（1）Pn在D里面或者上面，则最小外接圆没有变化，还是D；

（2）Pn不在D里面或者上面，需要更新一个新的外接圆，因为最终要使所有点都在最小圆内，因此我们要计算出最小圆为D’。

上述就是主要的Welzl’s算法思想，通过以下三个定理就可以计算出点云的最小外接圆：

设P是n个点的集合，P不是空的，p是P中的一个点。R也是一个点的集合，最小外接圆的边界点。

（1）如果存在一包含P且在其边界上有R的位置外接圆，则D(P,R)是唯一的；

（2）如果p不在D(P-{p},R)中，则p在D(P,R)上，只要存在，即D(P,R)=D(P-{p},Ru{p});

（3）如果D(P,R)存在，则P中存在一个由最多max{0,3-|R|}个点组成的集合S(除了边界点集合围成最小覆盖原所用到的点)，使得D(P,R)=D(S,R)。这意味着D(P)由P中最多3个点决定边界。（至于为什么三个点，个人人为最多三个点就可以得到一个外接圆，四个点在圆上也只要求三个就可以）。

2.2算法步骤

①首先现将所有点随机排列；

②按顺序把点一个一个的加入（一步一步的求前i个点的最小覆盖圆），每加入一个点就进入③；

③如果发现当前i号点在当前的最小圆的外面，那么说明点i一定在前i个点的最小覆盖圆边界上，我们转到④来进一步确定这个圆，否则前i个点的最小覆盖圆与前i-1个点的最小覆盖圆是一样的，则不需要更新，返回②；

④此时已经确认点i一定在前i个点的最小覆盖圆的边界上了，那么我们可以把当前圆的圆心设为第i个点，半径为0，然后重新把前i-1个点加入这个圆中（类似上面的步骤，只不过这次我们提前确定了点i在圆上，目的是一步一步求出包含点i的前j个点的最小覆盖圆），每加入一个点就进入⑤；

⑤如果发现当前j号点在当前的最小圆的外面，那么说明点j也一定在前j个点（包括i）的最小覆盖圆边界上，我们转到⑥来再进一步确定这个圆，否则前j个点（包括i）的最小覆盖圆与前i-1个点（包括i）的最小覆盖圆是一样的，则不需要更新，返回④；

⑥此时已经确认点i，j一定在前j个点（包括i）的最小覆盖圆的边界上了，那么我们可以把当前圆的圆心设为第i个点与第j的点连线的中点，半径为到这两个点的距离（就是找一个覆盖这两个点的最小圆），然后重新把前j-1个点加入这个圆中（还是类似上面的步骤，只不过这次我们提前确定了两个点在圆上，目的是求出包含点i，j的前k个点的最小覆盖圆），每加入一个点就进入⑦；

⑦如果发现当前k号点在当前的最小圆的外面，那么说明点k也一定在前k个点（包括i，j）的最小覆盖圆边界上，我们不需要再进一步确定这个圆了（因为三个点能确定一个圆！），直接求出这三点共圆，否则前k个点（包括i，j）的最小覆盖圆与前k-1个点（包括i，j）的最小覆盖圆是一样的，则不需要更新。

三、用到的数学知识

3.1两点求最小覆盖圆的中心和半径

已知有两点分别为p1(x,y，p2(x,y)

        最小覆盖圆P的圆心：P.X=(p1.x+p2.x)/2；        P.Y=(p1.y+p2.y)/2;       

        半径：R=

3.2三点求最小覆盖圆的中心和半径

三点共线：找到距离最远的一对点，以它们的连线为直径做圆即可（这里我们进行点云求最小覆盖圆数据已经处理了所以不会出现）

三个点不共线，那么就求三角形的外接圆：

设圆心坐标O为（x0,y0）,半径为r；三个点的坐标为分别A(x1,y1)，B(x2,y2)，C(x3,y3)。

三个点到圆心的距离相等：

化简得到：

使用克拉默法则对行列式求解

令：

四、代码模板

代码网上有好多写法，但是思路都是一样的，下面分享一下大佬写过的模板。

const int N = 100010;

struct P_Nod
{
	double x, y;
}b[N];
double sqr(double x)
{
	return x * x;
}

double dis(P_Nod x, P_Nod y)
{
	return sqrt(sqr(x.x - y.x) + sqr(x.y - y.y));
}
bool incircle(P_Nod x, P_Nod y, double R)
{
	if (dis(y, x) <= R)
	{
		return true;
	}
	return false;
}
P_Nod solve(double a, double b, double c, double d, double e, double f)
{
	double y = (f * a - c * d) / (b * d - e * a);
	double x = (f * b - c * e) / (a * e - b * d);
	P_Nod node1 = { x, y };
	return node1;
}
void GetSmallestCircle(vector<Point3d> vecPoint3d, double circleData[3], double& R)
{
	int n = vecPoint3d.size();
	if (!n)
	{
		return;
	}
	int i, j, k;
	for (i = 0; i < n; i++)
	{
		b[i].x = vecPoint3d[i].X;
		b[i].y = vecPoint3d[i].Y;
	}
	R = 0;
	P_Nod O;
	O.x = b[0].x;
	O.y = b[0].x;
	for (i = 0; i < n; i++)
	{
		if (!incircle(b[i], O, R))
		{
			O.x = b[i].x; O.y = b[i].y; R = 0;
			for (j = 0; j < i; j++)
			{
				if (!incircle(b[j], O, R))
				{
					O.x = (b[i].x + b[j].x) / 2;
					O.y = (b[i].y + b[j].y) / 2;
					R = dis(O, b[i]);
					for (k = 0; k < j; k++)
					{
						if (!incircle(b[k], O, R))
						{
							O = solve(b[i].x - b[j].x, b[i].y - b[j].y, (sqr(b[j].x) + sqr(b[j].y) - sqr(b[i].x) - sqr(b[i].y)) / 2,
								b[i].x - b[k].x, b[i].y - b[k].y, (sqr(b[k].x) + sqr(b[k].y) - sqr(b[i].x) - sqr(b[i].y)) / 2);
							R = dis(b[i], O);
						}
					}
				}
			}
		}
	}
	circleData[0] = O.x;	//圆心 X 坐标
	circleData[1] = O.y;	//圆心 Y 坐标
	circleData[2] = R;		//圆半径 R
}