用MATLAB实现最小圆覆盖问题!

用MATLAB实现最小圆覆盖问题!
我用 MALTAB实现了这样一个图形,怎样用MALTAB语句构建一个最小圆,能覆盖上面所有的九个点.
这个图的语句：
x=[22 8 4 51 38 17 81 18 62]
y=[38 13 81 32 11 12 63 45 12]
plot(x,y,'*')%
grid on%

%  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;            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;%总结：根据算法我写的这个程序有个隐藏的问题,由于要看比赛了,没时间再纠正这个问题了. -------------------------------------------------------------------------------function R=minCirclePoints3(A,B,C)X=[A(1) B(1) C(1)];Y=[A(2) B(2) C(2)];%计算三边的长度AB BC CAlen=[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,cosCif(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];    endend%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);endend
再问： 真的真的非常谢谢您！不得不承认，牛人还是多啊 不过，我在实现时，出现了这样的提错语句： ??? Error: File: zuixiaoyuan.m Line: 69 Column: 1 Function definitions are not permitted at the prompt or in scripts. 是什么意思？ 我自行调节了下，可还是没出结果 能不能麻烦您帮忙解释下 随便问下，、可以把圆心的坐标和半径都显示出来吗？ 不胜感激！