THESE/codes/kmeans/kmeans.m

   1 function [mu,mask]=kmeans(ima,k)\r
   2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\r
   3 %\r
   4 %   kmeans image segmentation\r
   5 %\r
   6 %   Input:\r
   7 %          ima: grey color image\r
   8 %          k: Number of classes\r
   9 %   Output:\r
  10 %          mu: vector of class means \r
  11 %          mask: clasification image mask\r
  12 %\r
  13 %   Author: Jose Vicente Manjon Herrera\r
  14 %    Email: jmanjon@fis.upv.es\r
  15 %     Date: 27-08-2005\r
  16 %\r
  17 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\r
  18 \r
  19 % check image\r
  20 ima=double(ima);\r
  21 copy=ima;         % make a copy\r
  22 ima=ima(:);       % vectorize ima\r
  23 mi=min(ima);      % deal with negative \r
  24 ima=ima-mi+1;     % and zero values\r
  25 \r
  26 s=length(ima);\r
  27 \r
  28 % create image histogram\r
  29 \r
  30 m=max(ima)+1;\r
  31 h=zeros(1,m);\r
  32 hc=zeros(1,m);\r
  33 \r
  34 for i=1:s\r
  35   if(ima(i)>0) h(ima(i))=h(ima(i))+1;end;\r
  36 end\r
  37 ind=find(h);\r
  38 hl=length(ind);\r
  39 \r
  40 % initiate centroids\r
  41 \r
  42 mu=(1:k)*m/(k+1);\r
  43 \r
  44 % start process\r
  45 \r
  46 while(true)\r
  47   \r
  48   oldmu=mu;\r
  49   % current classification  \r
  50  \r
  51   for i=1:hl\r
  52       c=abs(ind(i)-mu);\r
  53       cc=find(c==min(c));\r
  54       hc(ind(i))=cc(1);\r
  55   end\r
  56   \r
  57   %recalculation of means  \r
  58   \r
  59   for i=1:k, \r
  60       a=find(hc==i);\r
  61       mu(i)=sum(a.*h(a))/sum(h(a));\r
  62   end\r
  63   \r
  64   if(mu==oldmu) break;end;\r
  65   \r
  66 end\r
  67 \r
  68 % calculate mask\r
  69 s=size(copy);\r
  70 mask=zeros(s);\r
  71 for i=1:s(1),\r
  72 for j=1:s(2),\r
  73   c=abs(copy(i,j)-mu);\r
  74   a=find(c==min(c));  \r
  75   mask(i,j)=a(1);\r
  76 end\r
  77 end\r
  78 \r
  79 mu=mu+mi-1;   % recover real range\r
  80 \r