function [mu,mask]=kmeans1(ima,k) 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% 
%   kmeans image segmentation 
% 
%   Input: 
%          ima: grey color image 
%          k: Number of classes 
%   Output: 
%          mu: vector of class means  
%          mask: clasification image mask 
% 
%   Author: Jose Vicente Manjon Herrera 
%    Email: jmanjon@fis.upv.es 
%     Date: 27-08-2005 
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
% check image 
mask=ima; 
ima=double(ima); 
copy=ima;         % make a copy 
ima=ima(:);       % vectorize ima 圖像轉成一個向量 
mi=min(ima);      % deal with negative 求圖像中的最小灰度值 
ima=ima-mi+1;     % and zero values 圖像所有值向下浮動至最小值為1 
 
s=length(ima);    % 計算向量總元素數 
 
% create image histogram 
 
m=max(ima)+1;    %取最大像素值 
h=zeros(1,m);    %聚類前圖像向量中各灰度出現的次數統計 
hc=zeros(1,m);   %聚類后圖像向量所屬類編號 
 
for i=1:s 
  if(ima(i)>0)  
      h(ima(i))=h(ima(i))+1; 
  end 
end 
ind=find(h);    %找到圖像向量中所有出現次數非0的值 
hl=length(ind); %圖像向量中所有非0值的個數 
 
% initiate centroids 
mu=(1:k)*m/(k+1);%確定初始聚類中心（k個） 
 
% start process 
 
while(true) 
   
  oldmu=mu; %記錄上一次的聚類中心 
   
  % current classification   
  
  for i=1:hl 
      c=abs(ind(i)-mu);     %計算圖像向量中每一個分量與所有聚類中心的距離 
      cc=find(c==min(c));   %找到圖像向量中每一個分量與聚類中心距離中最短的那一個聚類中心 
      hc(ind(i))=cc(1);     %記錄圖像向量中每一個分量所屬的類號 
  end 
   
  %recalculation of means   
  %更新聚類中心 
  for i=1:k  
      a=find(hc==i);        %找到圖像向量中所有屬于第i類的分量 
      if sum(h(a))~=0 
          mu(i)=sum(a.*h(a))/sum(h(a)); %更新第i類的聚類中心 
      end 
  end 
   
  %聚類中心不再變化，說明聚類結束，退出 
  if(mu==oldmu)  
      break; 
  end 
   
end 
 
% calculate mask 
s=size(copy); 
mask=zeros(s); 
for i=1:s(1) 
    for j=1:s(2) 
        c=abs(copy(i,j)-mu); 
        a=find(c==min(c));   
        mask(i,j)=a(1); 
    end 
end