%?http://www.personal.kent.edu/~rmuhamma/Compgeometry/MyCG/CG-Applets/Center/centercli.htm
%?算法思路：
%?1.?在點集中任取3點A,B,C。?
%?2.?作壹個包含A,B,C三點的最小圓,圓周可能通過這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;?
endend;
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;
%總結：根據算法我寫的這個程序有個隱藏的問題，由於要看比賽了，沒時間再糾正這個問題了。?
-------------------------------------------------------------------------------
function?R=minCirclePoints3(A,B,C)
X=[A(1)?B(1)?C(1)];
Y=[A(2)?B(2)?C(2)];
%計算三邊的長度AB?BC?CA
len=[sqrt((X(1)-X(2))^2+(Y(1)-Y(2))^2)?sqrt((X(2)-X(3))^2+(Y(2)-Y(3))^2)?sqrt((X(3)-X(1))^2+(Y(3)-Y(1))^2)];
%在非特殊情況下計算三角形三角的余弦值?cosA,cosB,cosC
if(sum(len>0)==3)
abc=[cosABC(len(2),len(1),len(3))?cosABC(len(3),len(1),len(2))?cosABC(len(1),len(2),len(3))];
end
%兩點重合、三點重合、三點***線
if(len(1)==len(2)+len(3))
r=len(1)/2;
a=(X(1)+X(2))/2;
b=(Y(1)+Y(2))/2;
R=[a?b?r];
elseif(len(2)==len(1)+len(3))
r=len(2)/2;
a=(X(2)+X(3))/2;
b=(Y(2)+Y(3))/2;
R=[a?b?r];
elseif(len(3)==len(1)+len(2))
r=len(3)/2;
a=(X(1)+X(3))/2;
b=(Y(1)+Y(3))/2;
R=[a?b?r];
%--------------------------------------------------------------------------
else
tmp=(abc<=0);
if(tmp(1))
r=len(2)/2;
a=(X(2)+X(3))/2;
b=(Y(2)+Y(3))/2;
R=[a?b?r];
elseif(tmp(2))
r=len(3)/2;
a=(X(1)+X(3))/2;
b=(Y(1)+Y(3))/2;
R=[a?b?r];
elseif(tmp(3))
r=len(1)/2;
a=(X(1)+X(2))/2;
b=(Y(1)+Y(2))/2;
R=[a?b?r];
elseif(sum(tmp)==0)
a=(((X(1)^2-X(2)^2+Y(1)^2-Y(2)^2)*(Y(2)-Y(3)))-((X(2)^2-X(3)^2+Y(2)^2-Y(3)^2)*(Y(1)-Y(2))))/(2*(X(1)-X(2))*(Y(2)-Y(3))-2*(X(2)-X(3))*(Y(1)-Y(2)));
b=(((X(1)^2-X(2)^2+Y(1)^2-Y(2)^2)*(X(2)-X(3)))-((X(2)^2-X(3)^2+Y(2)^2-Y(3)^2)*(X(1)-X(2))))/(2*(Y(1)-Y(2))*(X(2)-X(3))-2*(Y(2)-Y(3))*(X(1)-X(2)))?;
r=sqrt((X(1)-a)^2+(Y(1)-b)^2);
R=[a?b?r];
end
end
%d=linspace(0,2*pi,100);
%plot(a+r*sin(d),b+r*cos(d),'-',X(1),Y(1),'ro',X(2),Y(2),'bo',X(3),Y(3),'ko',a,b,'.')
%axis([0?10?0?10])
function?c=cosABC(x,y,z)
c=(z^2+y^2-x^2)/(2*z*y);
end
end