MATLAB求解最小球覆盖问题
一、题目描述
最小圆覆盖是寻找能够覆盖平面给定的n个离散点的最小圆。该问题存在线性时间算法，即复杂度是O(n)。提供参考文献。编程实现三维最小球覆盖的一个算法。
二、算法分析
1. 理论依据
定理：如果点p不在集合S的最小覆盖球内，则p一定在S∪{p}的最小覆盖球上。
根据这个定理，我们可以分三次确定前i个点的最小覆盖圆。
- 1.令前i−1个点的最小覆盖球为C
- 2.如果第i个点在C内，则前i个点的最小覆盖球也是C
- 3.如果不在，那么第i个点一定在前i个点的最小覆盖球上，接着确定前i−1个点中还有哪两个在最小覆盖球上。因此，设当前球心为Pi，半径为0，做固定了第i个点的前i个点的最小球覆盖。
- 4.固定了一个点：不停地在范围内找到第一个不在当前最小球上的点Pj，设当前球心为(Pi+Pj)/2，半径为∣PiPj∣/2，做固定了两个点的，前j个点外加第i个点的最小球覆盖。
- 5.固定了两个点：不停地在范围内找到第一个不在当前最小球上的点Pk，设当前球为Pi,Pj,Pk的外接球。
2. 伪代码
圆 C;
for(i=1 to n)
{
if(P[i] 不在 C 内)
{
C = {P[i], 0};
for(j=1 to i-1)
{
if(P[j] 不在 C 内)
{
C = {0.5*(P[i]+P[j]), 0.5*dist(P[i], P[j])};
for(k=1 to j-1)
{
if(P[k] 不在 C 内)
{
C = 外接球(P[i], P[j], P[k]);
}
}
}
}
}
}
对于这个算法只需要三个模式完全相同的for循环就可以搞定，还有一个问题是如何求外接球。
3. 外接球算法分析
对于已知的三个点A,B,C，需求其最小外接球。如果三角形ABC为钝角三角形或直角三角形，球心M即为最长边的中点，半径R为最长边的一半。
如果三角形ABC为锐角三角形，其球心M为任意两边中垂线的交点。
求锐角三角形的球心M，需要得到中垂线的直线方程。
设A,B,C均为1*3的行向量，则AB中点为P=(A+B)/2，BC中点为Q=(B+C)/2。
设平面ABC的一个法向量为L，且L=(A-B)×(B-C)。(两向量的叉乘是两向量的法向量)
设PF=L×(A-B)，则PF与AB垂直且PF平行于平面ABC；设QF=L×（B-C），则QF与BC垂直且QF平行于平面ABC。
所以PF与QF分别为AB和BC在平面ABC内的中垂线的平行向量。
因此可以得到中垂线方程：
(Mx-Px)/PFx=(My-Py)/PFy=(Mz-Pz)/PFz
(Mx-Qx)/QFx=(My-Qy)/QFy=(Mz-Qz)/QFz
可以得到球心M(Mx，My，Mz)，半径R=|MA|。
4. 复杂度分析
由于一堆点最多只有3个点确定了最小覆盖求，因此n个点中每个点参与确定最小覆盖圆的概率不大于3/n
所以，每一层循环在第i个点处调用下一层的概率不大于3/i
那么设算法的三个循环的复杂度分别为T1(n),T2(n),T3(n)，则有：
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-wgZBp0Jo-1588649026859)(C:\Users\crjia\AppData\Roaming\Typora\typora-user-images\image-20200504224806385.png)]
不难解得，T1(n)=T2(n)=T3(n)=O(n)
三、代码展示
-
min_ball.m(框架部分)
function min_ball
syms a b c;
ooo=[a b c];
o=[0 0 0]; %球心
R=0; %球的半径
n=input(\'请输入你想要生成的离散点个数:\');
t=input(\'请输入你想要生成的点的范围(上限）:\');
A=rand(n,3);
B=(2*A-1)*t;
%随机产生n个离散点，储存在pot里面
pot=fix(B)+1;
%生成最小球
for i=1:n
if norm(pot(i,:)-o)>R
o = pot(i,:);
R=0;
for j=1:i-1
if norm(pot(j,:)-o)>R
o = (pot(i,:)+pot(j,:))/2;
R=norm(pot(j,:)-o);
for k=1:j-1
if norm(pot(k,:)-o)>R
%最小外接球(补充ball函数）
%ball函数输入三个点坐标，返回这三个点的最小球的球心坐标
x=pot(i,:);
y=pot(j,:);
z=pot(k,:);
if norm(z-y)^2+4*R^2<=(norm(z-x))^2
ooo=(z+x)/2;
R=norm(z-ooo);
elseif norm(z-x)^2+4*R^2<=(norm(z-y))^2
ooo=(z+y)/2;
R=norm(z-ooo);
else %三个点构成锐角三角形的情况
ooo=ballcenter(x,y,z);
R=norm(z-ooo);
o = ooo;
R=norm(pot(k,:)-o);
end
end
end
end
end
end
end
o=double(o); %保证小数防止mesh出bug
R=double(R);
[q, w, e]=sphere(30);
Q=R*q+o(1);
W=R*w+o(2);
E=R*e+o(3);
subplot(1,3,1)
mesh(Q,W,E)
hold on
%画n个点
x1=pot(:,1);
y1=pot(:,2);
z1=pot(:,3);
subplot(1,3,2)
scatter3(x1,y1,z1,\'r\');
hold on
subplot(1,3,3)
mesh(Q,W,E);
alpha(0.8) %设置透明度
shading flat %去掉那道些线
hold on
scatter3(x1,y1,z1,\'r\');
-
ballcenter.m(求最小球球心)
function p = ballcenter(x, y, z)
syms a b c;
% 圆的法向量
pf= cross(x-y, x-z);
p12 = (x + y)/2;
p23 = (y + z)/2;
% 求两条在三角形面内的中垂线的向量
p12f = cross(pf, x-y);
p23f = cross(pf, y-z);
eq1=(a-p12(1))*p12f(2)-(b-p12(2))*p12f(1);
eq2=(a-p12(1))*p12f(3)-(c-p12(3))*p12f(1);
eq3=(a-p23(1))*p23f(2)-(b-p23(2))*p23f(1);
eq4=(a-p23(1))*p23f(3)-(c-p23(3))*p23f(1);
[a,b,c]=solve(eq1,eq2,eq3,eq4,a,b,c);
p=[a b c];
end
四、结果展示
20个点，范围[-10,10]
30个点，范围[-20,20]
30个点，范围[-50,50]
40个点，范围[-100,100]
五、参考资料
[1] https://www.luogu.com.cn/problemnew/solution/P1742
[2] LINEAR-TIME ALGORITHMS FOR LINEAR PROGRAMMING IN R3 AND RELATED PROBLEMS*
NIMROD MEGIDDOt