一、定义
什么是最小圆覆盖?其实和最小矩形覆盖定义是类似的,给出一个点集,求能覆盖住所有点的最小圆。
二、两种算法
求最小圆覆盖有两种算法,分别是增量法和模拟退火,个人推荐增量法,它的精度更高一些,且时间复杂度是稳定的线性级(点的顺序打乱后),所以下面也主要介绍增量法的原理。
增量法
前置知识
1. 圆上三点确定唯一的一个圆。
这个道理很简单,考虑三角形外接圆就行。
2. 若已有某个点集的最小圆覆盖,向该点集中再加入一个圆外的点,这时最小圆会被更新,这个点一定出现在新的最小圆的边界上。
对于任意覆盖点集的圆,如果新加入的点不在边界上,那么总可以通过平移或缩放得到一个过该点还更小的圆。如下图:
算法步骤
1. 先对点集乱序处理,不然会被特殊数据卡成O(n^3),打乱顺序可以用这个函数random_shuffle,用法和sort一样。
2. 第一层循环i取1~n,分别求前i个点的最小圆覆盖。考虑若已知前i-1个点的最小圆覆盖,第i个点位于这个圆内,那前i-1个点的最小圆覆盖就是前i个点的最小圆覆盖,往后递推即可。若第i个点在圆外,那这个圆需要被更新,由前置知识2可知新的圆一定过第i个点,此时为了求该圆需要执行第二层循环。
3. 第二层循环j取1~i-1,分别求在过第i个点的前提下前j个点的最小圆覆盖。和上面的步骤类似,如果已知前j-1个点的最小圆覆盖,第j个点位于圆内,那前j-1个点的最小圆覆盖就是前j个点的最小圆覆盖。如果第j个点在圆外,那同样需要更新圆,且这个圆一定过第i个点和第j个点,此时为了求该圆需要执行第三层循环。
4. 第三层循环k取1~j-1,分别求在过第i个点和第j个点的前提下前k个点的最小圆覆盖。和上面的步骤类似,如果已知前k-1个点的最小圆覆盖,第k个点位于圆内,那前k-1个点的最小圆覆盖就是前k个点的最小圆覆盖。若第k个点在圆外,该圆需要被更新,新圆一定过第i个点、第j个点和第k个点,此时由前置知识1可知该圆已经被确定了,对三点求外接圆就是待求圆。
示意图
总结
1. 其实整个算法有点类似递归,一般情况下都是到达最后一层循环时才能够真正得到一个圆,之后向上回溯使用这个圆。
2. 步骤中加粗的字对于理解该算法非常重要。
3. 第一层循环开始前可以初始化圆心为第1个点,半径为0,之后相当于向圆内依次加入n个点。第二层循环开始前可以初始化圆心为第i个点(前提条件),半径为0,之后相当于向圆内依次加入前i-1个点。第三层循环开始前可以初始化圆心为第i个点和第j个点的中点(前提条件),半径为圆心到第i个点或第j个点的距离,也就是这两点的最小圆覆盖,之后相当于向圆内依次加入前j-1个点。
4. 最小圆覆盖问题 算法步骤与证明+代码模板,最后推荐一篇博客,写的非常好。
模拟退火
在确定规模下增量法的时间复杂度和精度都是有保证的,但模拟退火的做法则是不断向精确解逼近,其精度和时间复杂度是正相关的,一般情况下要想获得更高的精度就需要运行更久的时间。
三、模板
//最小圆覆盖(增量法)
//别忘了对点随机排序
//时间复杂度O(n)
circle minCircleCover(polygon A)
{
circle c;
random_shuffle(A.p, A.p+A.n);
c.p = A.p[0];
c.r = 0.0;
for(int i = 0; i < A.n; i++)
if(!c.relation(A.p[i]))
{
c.p = A.p[i];
c.r = 0.0;
for(int j = 0; j <= i-1; j++)
if(!c.relation(A.p[j]))
{
c.p = (A.p[i]+A.p[j])/2;
c.r = c.p.distance(A.p[i]);
for(int k = 0; k <= j-1; k++)
if(!c.relation(A.p[k]))
c = circle(A.p[i], A.p[j], A.p[k]);
}
}
return c;
}
//最小圆覆盖(模拟退火)
//精度不如增量法高,慎用!
circle minCircleCover(polygon A, int t)//t只是为了函数重载
{
circle c;
c.p.x = c.p.y = c.r = 0.0;
for(int i = 0; i < A.n; i++)
c.p.x += A.p[i].x, c.p.y += A.p[i].y;
c.p.x /= A.n, c.p.y /= A.n;
for(int i = 0; i < A.n; i++)
c.r = max(c.r, c.p.distance(A.p[i]));
double delta = 100.0;//delta越大结果越精确
while(delta > eps)
{
int id = 0;
double maxDis = c.p.distance(A.p[id]);
for(int i = 1; i < A.n; i++)
{
double tempDis = c.p.distance(A.p[i]);
if(tempDis > maxDis)
{
maxDis = tempDis;
id = i;
}
}
c.r = min(c.r, maxDis);
c.p.x += (A.p[id].x-c.p.x)/maxDis*delta;
c.p.y += (A.p[id].y-c.p.y)/maxDis*delta;
delta *= 0.98;
}
return c;
}