图的最小顶点覆盖，用matlab实现

求一个图的最小顶点覆盖集Red，图用关联矩阵B表示，若删除掉集合Red中一个元素仍然能覆盖整个图，则可以通过删除该元素得到新的最简Red集。若不能，则保留原来的Red集。
例：图为B,原有的属性集为A=[1,2,3,4,5]，但通过别的算法可以得到约简为Red=[1,2,3,4];现要对已求得的Red集再进行约简得到最小的约简集。
B=[0,1,0,0,0; （对应边e1）
0,0,1,0,0; (对应边e2)
0,1,1,0,0; （对应边e3）
1,0,0,1,0; (对应边e4)
1,0,0,0,1]; (对应边e5)
通过（最小约简）步骤后，应该得到更小的属性约简集为Red=[1,2,3].
设计代码
clear
clc
B=[0,1,0,0,0;0,0,1,0,0;0,1,1,0,0;1,0,0,1,0;1,0,0,0,1];
Red=[1,2,3,4];
[m,n]=size(B);
for k=1:length(Red(:))
a=Red;
a(k)=[];
clear A;
for i=1:m
[~,j]=find(B(i,:)>0); %找到第i行不为0的元素坐标
A=intersect(a,j);
if isempty(A) %若交集为空时,跳出循环
break
end
disp(a) %若对于B中从1-m行每行与Red交不为空，则输出新的Red集。（此处为难点）
end
end
以上代码实现不了，难点处设计应该有误，不知如何改正，能实现结果。

1. 理论依据
   定理：如果点p不在集合S的最小覆盖球内，则p一定在S∪{p}的最小覆盖球上。
   根据这个定理，我们可以分三次确定前i个点的最小覆盖圆。
2. 令前i−1个点的最小覆盖球为C
3. 如果第i个点在C内，则前i个点的最小覆盖球也是C
4. 如果不在，那么第i个点一定在前i个点的最小覆盖球上，接着确定前i−1个点中还有哪两个在最小覆盖球上。因此，设当前球心为Pi，半径为0，做固定了第i个点的前i个点的最小球覆盖。
5. 固定了一个点：不停地在范围内找到第一个不在当前最小球上的点Pj，设当前球心为(Pi+Pj)/2，半径为∣PiPj∣/2，做固定了两个点的，前j个点外加第i个点的最小球覆盖。
6. 固定了两个点：不停地在范围内找到第一个不在当前最小球上的点Pk，设当前球为Pi,Pj,Pk的外接球。
7. 伪代码
   圆 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循环就可以搞定，还有一个问题是如何求外接球。
8. 外接球算法分析
   对于已知的三个点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|。
9. 复杂度分析
   由于一堆点最多只有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=(2A-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+4R^2<=(norm(z-x))^2
ooo=(z+x)/2;
R=norm(z-ooo);
elseif norm(z-x)^2+4R^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=Rq+o(1);
W=Rw+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. 在点集中任取3点A,B,C。

% 2. 作一个包含A,B,C三点的最小圆,圆周可能通过32313133353236313431303231363533e58685e5aeb931333264646533这3点，也可能只通过其中两点,但包含第3点.后一种情况圆周上的两点一定是位于圆的一条直径的两端。

% 3. 在点集中找出距离第2步所建圆圆心最远的D点，若D点已在圆内或圆周上，则该圆即为所求的圆，算法结束.否则执行第4步。

% 4. 在A,B,C,D中选3个点,使由它们生成的一个包含这4个点的圆为最小，这3 点成为新的A,B,C，返回执行第2步。

% 若在第4步生成的圆的圆周只通过A,B,C,D 中的两点，则圆周上的两点取成新的A和B,从另两点中任取一点作为新的C。

clear all;close all;clc;

x=[22 8 4 51 38 17 81 18 62]

y=[38 13 81 32 11 12 63 45 12]

plot(x,y,'*');hold on;

grid on%

set_3P=nchoosek(1:length(x),3);

AI=set_3P(1,1);

BI=set_3P(1,2);

CI=set_3P(1,3);

A=[x(AI) y(AI)];

B=[x(BI) y(BI)];

C=[x(CI) y(CI)];

while 1

R=minCirclePoints3(A,B,C);

cr=[R(1),R(2)];

r=zeros(1,length(x));

for i=1:length(x)

r(i)=sqrt((x(i)-cr(1))^2+(y(i)-cr(2))^2);

end;

maxValue=max(r);    %或者N=max(r(:))

[mc]=find(maxValue==r);

if r(mc)<=R(3)%没有点在圆外，结束回家吃饭去

alpha=0:pi/20:2*pi;%角度[0,2*pi]

plot(cr(1)+R(3)*cos(alpha),cr(2)+R(3)*sin(alpha),'--r');%中心点在(R(1),R(2))半径为R(3)的圆

axis equal;

break;%所有点都被圆覆盖

else

%距离圆心最远的点在圆外

end;

D=[x(mc),y(mc)];

P=[A;B;C;D];%保存这四个点的坐标

DI=mc;

set_3P=nchoosek([AI,BI,CI,DI],3);

rSet=[];

for i=1:length(set_3P)

A=[x(set_3P(i,1)) y(set_3P(i,1))];

B=[x(set_3P(i,2)) y(set_3P(i,2))];

C=[x(set_3P(i,3)) y(set_3P(i,3))];

R=minCirclePoints3(A,B,C);

rSet=[rSet;[R,i]];%每行：圆心坐标,半径,第几组(每组包括随机的三个点)

end;

rSet=sortrows(rSet,3);%按照半径排序

%   在四个圆中找一个最小半径圆包含这四个点

for i=1:size(rSet,1)

for j=1:4

if sqrt((rSet(i,1)-(P(j,1) ))^2+ ( rSet(i,2)-(P(j,2)))^2) >rSet(i,3)%这个圆不行

break;

end

end;

if j>4%第i组三个点产生的圆可行--必然可以找到一个

break;

end;

end;

mc=rSet(i,4);

A=[x(set_3P(mc,1)) y(set_3P(mc,1))];

B=[x(set_3P(mc,2)) y(set_3P(mc,2))];

C=[x(set_3P(mc,3)) y(set_3P(mc,3))];

end;