離散點的最小包圍圓

% 算法思路：

% 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。

 

matlab代碼：

close all;
clc;

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

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

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

grid on%

set_3P=nchoosek(1:length(x),3);%從1到9，每次取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)];%初始化為9個點中選擇3個點的全排列的第一種

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)]; %三個點的x坐標

Y=[A(2) B(2) C(2)]; %三個點的y坐標

%計算三邊的長度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

 

 

 

 

來源：https://zhidao.baidu.com/question/308990622.html