最小覆盖模型matlab_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]

五、参考资料

[2] LINEAR-TIME ALGORITHMS FOR LINEAR PROGRAMMING IN R3 AND RELATED PROBLEMS*

NIMROD MEGIDDOt