c－cluster.m

一种基于matlab语言编程实现的数据聚类算法

disp(strvcat('基于C－均值算法的特征点聚类',strcat(datestr(now)))); 
 
% S=[-1,0,0,1,3,5,5,4;-1,1,0,0,5,4,5,5] 
% S=input('S=') 
S=[0,0,4,4,5,5,1;0,1,4,5,4,5,0] 
sx=S(1,:); 
sy=S(2,:); 
len = length(sx);  
group = zeros(len,1); % group存储特征点与类相对应的情况 
c_num = 0; % c_num存储聚类次数 
 
%yesno='n'; 
%while yesno=='n' | yesno=='N' 
disp('请输入类别数：'); C=input('C='); % 聚类中心数 
 
tic; % 计时开始 
 
% i是特征点控制变量 
  % j,k是聚类中心控制变量 
  mx = zeros(C,1); my = zeros(C,1); % 初始化聚类中心坐标 
  mx(1,1) = sx(1,1); 
  my(1,1) = sy(1,1); % 第一个聚类中心的初始坐标为源图中第一个特征点的坐标 
  
 for k = 2:1:C % 计算第k个聚类中心的初始坐标 
maxd = zeros(len,1); 
for j = 1:1:k-1 % 求每个特征点与前k-1个聚类中心的距离之和 
i = 1:1:len;  
d = (sx(i)-mx(j)).^2+(sy(i)-my(j)).^2; 
maxd = maxd+d'; 
end 
d=find(maxd==max(maxd)); % 第k个聚类中心是与前k-1个聚类中心相距最远的特征点 
mx(k,1)=sx(1,d); 
my(k,1)=sy(1,d); 
  end 
  disp('所有开始聚类中心为:') 
 disp([mx,my]); 
 
  mxtemp = zeros(C,1); mytemp = zeros(C,1); N = zeros(C,1); 
  for i=1:1:len % 基于初始聚类中心，进行初始划分 
k=1:1:C; 
d=(sx(i)-mx(k)).^2+(sy(i)-my(k)).^2; 
g=find(d==min(d));  
group(i)=g(1);  %根据最小距离分类 
N(g(1),1) = N(g(1),1)+1; %计算各个类别元素个数 
mxtemp(g(1),1) = mxtemp(g(1),1)+sx(i); %求同一类别中元素x值之和 
mytemp(g(1),1) = mytemp(g(1),1)+sy(i); %求同一类别中元素y值之和 
end 
 
mx(:)=0; my(:)=0; N(find(N==0))=1; 
  mx=mxtemp./N; my=mytemp./N; % 基于初始划分，重新计算聚类中心 
  Je=0; 
  for i=1:1:len % 基于初始划分，计算误差平方和 
Je = Je+(sx(i)-mx(group(i))).^2+(sy(i)-my(group(i))).^2; 
  end 
   
  no_change =1; 
  while no_change==1 % 连续迭代thre次，Je不变，则停止聚类 
ifbreak = 0; 
for i=1:1:len 
if N(group(i))~=0 
k=1:1:C; % 计算p 
d=(N(k).*((sx(i)-mx(k)).^2+(sy(i)-my(k)).^2))./(N(k)+1); 
d(group(i))=(N(group(i)).*((sx(i)-mx(group(i))).^2+(sy(i)-my(group(i))).^2))./(N(group(i))-1);  
 
g=find(d==min(d)); 
if g(1)~=group(i) % 如果pk<=pj,把特征点从i类移到k类 
mxtemp(group(i))=mxtemp(group(i))-sx(i); % 把特征点从i类中移出 
mytemp(group(i))=mytemp(group(i))-sy(i); 
N(group(i))=N(group(i))-1; 
 
group(i)=g(1); % 修改特征点的类别归属 
mxtemp(g(1))=mxtemp(g(1))+sx(i); % 把特征点移入k类 
mytemp(g(1))=mytemp(g(1))+sy(i); 
N(g(1))=N(g(1))+1; 
       
mx(:)=0; my(:)=0; N(find(N==0))=1; 
mx=mxtemp./N; my=mytemp./N; % 重新计算聚类中心 
          Je=0; 
          for i=1:1:len % 计算误差平方和 
Je = Je+(sx(i)-mx(group(i))).^2+(sy(i)-my(group(i))).^2; 
          end 
ifbreak=1; 
end 
end 
    end 
no_change=ifbreak==1; 
 disp('所有最终聚类中心为:') 
 disp([mx,my]); 
end 
 
  time=toc; 
  figure; hold on; %  
 color='.r.g.b.m.c.y';  col2='rgbmcy';%'gbrymc'; 
  for k=1:1:C 
xtemp=sx(find(group==k)); ytemp=sy(find(group==k)); 
plot(xtemp,ytemp,[color((k-1)*2+1),color((k-1)*2+2)],'MarkerSize',30); 
plot(mx(k),my(k),'p','MarkerSize',10) 
%rectangle('Position',[round(mx(k)),round(my(k)),3,3],'Curvature',[1,1],'FaceColor',col2(k)); 
disp(k) 
disp('类元素有:') 
disp([xtemp',ytemp']) 
  end 
  title(['类别数为',num2str(C),', 聚类时间',num2str(time),'s']); 
  hold off; 
 
c_num=c_num+1; 
judge(c_num,1)=C; judge(c_num,2)=Je; % judge存储类别数和误差平方和 
%figure;  
%plot(judge(:,1),judge(:,2),'-s','MarkerFaceColor','g'); 
%disp('是否结束聚类?'); yesno=input('Y/N? ','s'); 
end