% 聚類方法：C－均值算法
clear;clc; % 清內存，工作空間及命令行空間

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