xx=[148723.8400,160336.6900;
146444.3400,161044.5800;
149877.2200,161091.0400;
156650.5300,160402.7300;
154830.0900,160255.0500;
155780.6100,161155.0600;
155973.6900,161520.8100;
155757.6200,161698.6300;
153146.9700,160216.4700;
155477.3700,162651.1300;
164429.8700,163489.9700;
158431.2500,162322.4000;
158487.3400,164413.3800;
169251.0600,165333.1600;
165295.9600,166307.4800;
157354.7000,166956.6600;
163428.0200,166293.3500;
161442.7300,165833.1700;
150413.2200,159364.6600;
146987.6800,159418.2100;
148375.7200,159294.2600;
149100.7800,159241.6700;
147221.0600,157904.4100;
149007.8800,158255.2300;
147900.6400,158283.6700;
150211.9400,158547.6400;
148886.1800,156757.9800;
149951.9000,157738.2800;
150870.2500,156890.8400;
149088.1600,155678.1900;
152100,165900;
145693.6500,160063.0700;
145614.2600,159480.6800;
145488.5900,158833.2500;
139881.7800,156200.2800;
137348.1400,155938.4600;
139202.4500,156601.0100;
141521.1700,156506.1000;
143425.7500,157339.4000;
144718.6800,157602.5600;
137889.6800,157781.6600;
136656.2300,157553.5500;
139713.8100,158198.0500;
142527.7600,158312.9700;
144512.6700,158820.0300;
142990.5500,159397.0900;
144348.8600,159828.8200;
143385.2200,160463.4800;
141629.7100,160163.6500;
140522.9900,160396.7500;
137404.8000,160392.0300;
139356.8500,160405.1700;
145961.7700,157814.6700;
143191.8600,155025.2100;
144403.2900,155930.0100;
146013.2900,155683.1200;
146935.2600,156555.2700;
147771.4400,156163.0400;
142144.7900,154597.2600;
147066.4000,155149.2900;
144690.9800,154316.9900;
134100,160300;
138760.2800,153074.0100;
138393.6800,154240.8900;
139465.5100,154309.7300;
137872.6200,155016.2100;
139064.1200,155217.4300;
138015.1800,155514.9600;
140500.8400,155007.6800;
138793.5200,155638.1900;
139842.8200,155778.7000;
138989.3200,156131.4500;
141339.4500,153959.7200;
140309.1500,153473.7800;
140416.5700,152146.0800;
141612.0200,152424.7800;
145313.2100,152711.9500;
142147.4400,152788.0900;
143839.8000,153020.8300;
142733.3500,153099.5900;
143678.0600,153444.6400;
145035.9900,153315.7000;
144302.0500,152239.5000;
140808.2000,152432.4700;
143106.6700,152356.1300;
143800.1000,151838.3900;
142199.9400,151897.7700;
145146.7700,151489.5600;
143976.3000,151159.1600;
142646.6000,151242.9600;
142863.2700,150030.9400;
143903.2700,150127.8400;
139800,149000;
149026.5600,154716.6800;
146543.6800,153913.7800;
151449.4400,153764.2500;
149366.8900,153648.0500;
147188.7900,152823.1500;
151500.7200,152128.3200;
149214.8100,151472.2300;
147951.6800,151173.0700;
146813.3800,151435.8200;
149198.4600,149918.9200;
145062.6800,150092.5500;
146277.0000,149377.0500;
148037.4200,149443.9000;
146343.9800,150237.0900;
155100,150000];
z=xx(:,1)
y=xx(:,2)
plot(z,y,'*');hold on;
grid on%
if length(z)==2
cc=sqrt(((z(1)-z(2))^2+(y(1)-y(2))^2));
r=cc/2;
a=(z(1)+z(2))/2;
b=(y(1)+y(2))/2;
R=[a b r]
alpha=0:pi/20:2*pi;%角度[0,2*pi]
plot(a+r*cos(alpha),b+r*sin(alpha),'--r');%中心点在（a,b）半径为r的圆
axis equal;
end
if length(z)>=3
[s,d]=size(xx);
for i=1:1:s
for j=i:1:s
summ=0;
for k=1:1:d
summ=summ+(xx(i,k)-xx(j,k))^2; %计算数据集中的点间距，
end
DistanceBO(i,j)=sqrt(summ); %构造出点间距矩阵DistanceBO
DistanceBO(j,i)=sqrt(summ);
end
end
for i=1:1:s %找到各个点间距离
for j=i:1:s
DD(i,j)=DistanceBO(i,j);
end;
end;
aaa=DD(:);
TempSortDistanceBO=sort(aaa);
bbb=fliplr(TempSortDistanceBO');
%下边找出3个距离最远的点
maxfrist=[];
maxsecond=[];
maxtwo=[];
allmaxfrist=[];
finalpoints=[];
for i=1:1:s %找到各个点间距离
for j=i:1:s
if DD(i,j)==bbb(1)
maxfrist=[i j];
end;
end;
end;
for i=1:1:s %找到各个点间距离
for j=i:1:s
if i==maxfrist(1)&&j==maxfrist(2)
d=0;
elseif bbb(2)~=bbb(1)
if DD(i,j)==bbb(2)
maxsecond=[i j];
if maxsecond(1)==maxfrist(1)||maxsecond(1)==maxfrist(2)
maxtwo=maxsecond(2);
elseif maxsecond(2)==maxfrist(1)||maxsecond(2)==maxfrist(2)
maxtwo=maxsecond(1);
else maxtwo=maxsecond;
end;
finalpoints=[maxfrist maxtwo];
end;
elseif bbb(2)==bbb(1)
if DD(i,j)==bbb(2)
bb=[i j];
allmaxfrist=[allmaxfrist bb];
%maxtwo=[allmaxfrist(2)];
% if allmaxfrist(1)-allmaxfrist(2)==0||allmaxfrist(1)-allmaxfrist(3)==0||allmaxfrist(1)-allmaxfrist(4)==0
% finalpoints=[allmaxfrist(2) allmaxfrist(3) allmaxfrist(4)] ;
% elseif allmaxfrist(2)-allmaxfrist(3)==0||allmaxfrist(2)-allmaxfrist(4)==0
% finalpoints=[allmaxfrist(1) allmaxfrist(3) allmaxfrist(4)] ;
% elseif allmaxfrist(4)-allmaxfrist(3)==0
% finalpoints=[allmaxfrist(1) allmaxfrist(2) allmaxfrist(4)] ;
% else finalpoints=[allmaxfrist(1) allmaxfrist(2) allmaxfrist(3) allmaxfrist(4)];
if allmaxfrist(1)==maxfrist(1)||allmaxfrist(1)==maxfrist(2)
finalpoints=[maxfrist allmaxfrist(2)];
elseif allmaxfrist(2)==maxfrist(1)||allmaxfrist(2)==maxfrist(2)
finalpoints=[maxfrist allmaxfrist(1)];
else finalpoints=[maxfrist allmaxfrist(1)];%此时最终三个点取前三个点
end;
end;
end;
end;
end;
%finalpoints=[maxfrist maxtwo];
%set_3P=nchoosek(1:length(finalpoints),3);
%AI=set_3P(1,1);
%BI=set_3P(1,2);
%CI=set_3P(1,3);
%while 1
A=[z(finalpoints(1)) y(finalpoints(1))];
B=[z(finalpoints(2)) y(finalpoints(2))];
C=[z(finalpoints(3)) y(finalpoints(3))];
R=minCirclePoints3(A,B,C);
cr=[R(1),R(2)];
r=zeros(1,length(z));
for i=1:length(z)
r(i)=sqrt((z(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)<=0.000001%没有点在圆外，结束回家吃饭去
alpha=0:pi/20:2*pi;%角度[0,2*pi]
plot(R(1)+R(3)*cos(alpha),R(2)+R(3)*sin(alpha),'--r');%中心点在（R(1),R(2)）半径为R(3)的圆
axis equal;
break; %所有点都被圆覆盖
else
%距离圆心最远的点在圆外
end;
if r(mc)-R(3)>0.000001
D=[z(mc),y(mc)];
P=[A;B;C;D];%保存这四个点的坐标
DI=mc;
set_3P=nchoosek([finalpoints(1),finalpoints(2),finalpoints(3),DI],3);
rSet=[];
for i=1:length(set_3P)
A=[z(set_3P(i,1)) y(set_3P(i,1))];
B=[z(set_3P(i,2)) y(set_3P(i,2))];
C=[z(set_3P(i,3)) y(set_3P(i,3))];
R=minCirclePoints3(A,B,C);
alpha=0:pi/20:2*pi;%角度[0,2*pi]
% plot(R(1)+R(3)*cos(alpha),R(2)+R(3)*sin(alpha),'--r');%中心点在（R(1),R(2)）半径为R(3)的圆
axis equal;
rSet=[rSet;[R,i]];%每行：圆心坐标,半径,第几组(每组包括随机的三个点)
end;
rSet=sortrows(rSet,3);%按照半径排序
% 在四个圆中找一个最小半径圆包含这四个点
for i=1:1:size(rSet,1)
k=1;
for j=1:1:4
if sqrt((rSet(i,1)-(P(j,1) ))^2+ ( rSet(i,2)-(P(j,2)))^2) -rSet(i,3)>0.00001%这个圆不行
break;
else
k=k+1;
end
end;
if k>4%第i组三个点产生的圆可行--必然可以找到一个
break;
评论0