MATLAB求解最小球覆盖问题（复杂度O(N))

数学软件大作业一

一、题目描述

最小圆覆盖是寻找能够覆盖平面给定的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)

三、代码展示

1. 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');

2. 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