isodata.m

模式识别较全的讲解

function ISODATA
%%%%%%%%%%%%%ISODATA算法%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 输出为每次迭代之后得到的聚类中心以及最终得到的聚类中心
% 并用图表示最后得到的聚类结果
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc;
clear;

x=[0 3 2 1 5 4 6 5 6 7;
   0 8 2 1 3 8 3 4 4 5];
Nc=1; %初始聚类中心数目
K= 2; %预期的聚类中心数目；
theta_N = 1; %每一个聚类域中最少的样本数目，即如少于此数就不作为一个独立的聚类；
theta_S = 1; %一个聚类域中样本距离分布的标准差；
theta_c = 4; %两聚类中心之间的最小距离，如小于此数，两个聚类进行合并；
L = 0; %在一次迭代运算中可以合并的聚类中心的最多对数；
I = 4; %迭代运算的次数序号

z=[];
x_num=length(x);
label=zeros(x_num,1);
for i=1:Nc
    z=[z,x(:,i)];
end
count=0;
while(1)
    count=count+1;
    %%%%%将N个模式样本分给最近的聚类%%%%%%%
    while(1)
        for i=1:x_num;
            label(i)=dist_min(z,x(:,i));
        end
        z_before=z;
        Nc_before=Nc;
        for i=1:Nc_before
            if length(find(label==i))<theta_N
                z=[z(:,1:i-1-(Nc_before-Nc)),z(:,i+1-(Nc_before-Nc):Nc)];
                Nc=Nc-1;
            end
        end
        if(z_before==z)
            break;
        end
    end

    %%%%%%%%%%%修正各聚类中心值%%%%%%%%%%%%%
    z=zeros(2,Nc);
    for i=1:x_num
        z(:,label(i))=z(:,label(i))+x(:,i);
    end

    for i=1:Nc
        z(:,i)=z(:,i)/length(find(label==i));
    end

    fprintf('The %dth iteration:\n',count);
    z
    %%%%%%%%%%计算各聚类中心间的平均距离%%%%%%%%%
    D_total=zeros(Nc,1);
    D=zeros(Nc,1);
    for i=1:x_num
        D_total(label(i))=D_total(label(i))+norm(x(:,i)-z(:,label(i)));
    end
    for i=1:Nc
        D(i)=D_total(i)/length(find(label==i));
    end

    %%%%%%%%%%%%%计算全部模式样本对其相应聚类中心的总平均距离%%%%%%%%%%%%
    D_avg=sum(D_total)/x_num;
    %%%%%%%%%%%%%判别分裂、合并及迭代运算等步骤%%%%%%%%%%%%%%%%%%
    flag=0;                                 %判断分裂是否成功的标志位
    if count==I                             %最后一次迭代，算法结束
        theta_c=0;
        break;
    elseif (mod(count,2)~=0 & Nc<2*K)| Nc<=K/2       %进行分裂处理
        sigma=zeros(2,Nc);
        for i=1:x_num
            sigma(1,label(i))=sigma(1,label(i))+(x(1,i)-z(1,label(i)))^2;
            sigma(2,label(i))=sigma(2,label(i))+(x(2,i)-z(2,label(i)))^2;
        end
        for i=1:Nc
            sigma(:,i)=sqrt(sigma(:,i)/length(find(label==i)));
        end

        sigma_max=zeros(2,Nc);
        for i=1:Nc
            [sigma_max(1,i),sigma_max(2,i)]=max(sigma(:,i));
        end

        Nc_temp = Nc;
        for j=1:Nc_temp                          %将所有满足条件的都分裂
            if sigma_max(1,j)>theta_S
                if (D(j)>D_avg & length(find(label==j))>2*(theta_N+1)) | Nc<=K/2
                    z_temp=zeros(2,1);
                    z_temp(sigma_max(2,j))=0.5*sigma_max(1,j);
                    z=[z(:,1:j-1+Nc-Nc_temp),z(:,j+Nc-Nc_temp)+z_temp,z(:,j+Nc-Nc_temp)-z_temp,z(:,j+1+Nc-Nc_temp:Nc)];
                    Nc=Nc+1;
                    flag=1;
                else
                    continue;
                end
            end
        end
    end
    if flag==1                                   %分裂成功，跳到第二步
        continue;
    end
    if mod(count,2)==0 | Nc>=2*K | flag==0        %进行合并处理
        if L==0
            continue;
        end
        D_ij=[];
        lookup_table=[];
        ij=[];
        for i=1:Nc-1
            for j=i+1:Nc
                lookup_table=[lookup_table,[i;j]];
                D_ij=[D_ij,norm(z(:,i)-z(:,j))];
            end
        end
        [D_ij,idx]=sort(D_ij);
        z_temp=[];
        Nc_temp=Nc;
        for i=1:L
            if D_ij(i)>=theta_c
                break;
            else
                m = lookup_table(1,idx(i));
                n = lookup_table(2,idx(i));
                if length(find(ij==m))==0 & length(find(ij==n))==0
                    ij=[ij,m,n];
                    z_temp=[z_temp,(z(:,m)*length(find(label==m))+z(:,n)*length(find(label==n)))/(length(find(label==m))+length(find(label==n)))];
                    Nc_temp=Nc_temp-1;
                end
            end
        end
        z_unite=[];
        for i=1:Nc
            if length(find(ij==i))==0
                z_unite=[z_unite,z(:,i)];
            end
        end
        z_unite=[z_unite,z_temp];
        z=z_unite;
        Nc=Nc_temp;

    end

end
fprintf('The iteration is finished.\n');
fprintf('The result is:\n');
z

%%%%%%%%%%%绘出最终聚类结果%%%%%%%%%%%%%%%
figure(1);
hold on;
for i=1:x_num
    if label(i)==1;
        plot(x(1,i),x(2,i),'b+');
    elseif label(i)==2;
        plot(x(1,i),x(2,i),'rx');
    else
        plot(x(1,i),x(2,i),'g*');
    end
end
hold off;


function y=dist_min(z,x);
for i=1:size(z,2)
    if i==1
        dist=norm(z(:,i)-x);
        y=i;
    else
        if norm(z(:,i)-x)<dist
            dist=norm(z(:,i)-x);
            y=i;
        end
    end
end