%
% fuzzy c-means clustering
%
%   [ xhat, u ] = fcmeans( data, mpow, nc [, c] )
%
%  xhat = cluster centres
%  u    = membership functions for each point
%  data = data points
%  mpow = degree of fuzziness (0-1)
%  nc   = number of clusters
%  c    = initial cluster centres (optional)
%
% bob newell, june 1997
%
%-----------------------------------------------------
function [ xhat, u ] = fcmeans( data, mpow, nc, c )
%
[ n, m ] = size( data );
d = zeros( n, nc );
u = zeros( n, nc );
%
% given or random centres
if nargin > 3,
  xhat = c;
else
  notuniq = 1;
  while notuniq,
    sel = round( rand( 1, nc ) * ( n - 1 ) ) + 1;
    notuniq = 0;
    for i = 2:nc,
      for j = 1:i-1,
        if sel(i) == sel(j), notuniq = 1; end;
      end;
    end;
  end;
  xhat = data( sel, : ) * 1.1;
end;
%
% iterate
crit = 1;
while crit > 0.0001,
%
% calculate the distances
  for i = 1:nc,
    for k = 1:n,
      d(k,i) = norm( data(k,:) - xhat(i,:) );
    end;
  end;
%
% calculate membership matrix
  for i = 1:nc,
    for k = 1:n,
      s = 0.0;
      for j = 1:nc,
        s = s + ( d(k,i) / d(k,j) )^( ( 2 / mpow ) - 1 );
      end;
      u(k,i) = 1 / s;
    end;
  end;
%
% calculate new centres
  crit = 0;
  for i = 1:nc,
    sn = zeros( 1, m ); sd = 0;
    for k = 1:n,
      sn = sn + data(k,:) .* u(k,i) ^ mpow;
      sd = sd + u(k,i) ^ mpow;
    end;
    vi = sn ./ sd;
    crit = crit + norm( xhat(i,:) - vi );
    xhat(i,:) = vi;
  end; 
%
end;
%
% calculate the distances
  for i = 1:nc,
    for k = 1:n,
      d(k,i) = norm( data(k,:) - xhat(i,:) );
    end;
  end;
%
% calculate membership matrix
  for i = 1:nc,
    for k = 1:n,
      s = 0.0;
      for j = 1:nc,
        s = s + ( d(k,i) / d(k,j) )^( ( 2 / mpow ) - 1 );
      end;
      u(k,i) = 1 / s;
    end;
  end;
% 
