%
% multidimensional ellipsoid clustering algorithm
%
%  [ xhat, shape, alloc ] = 
%            ellclust( data, dmin, dmax, pflag )
%
%   xhat  = cluster centres (nclust x ndim)
%   shape = [ centre radii angle ]
%               (nclust x 3*ndim)
%   alloc = cluster allocation of data (npts x 1)
%   data  = data points (npts x ndim)
%           (best normalise all dimensions to 0-1)
%   dmin  = minimum agglomeration distance
%   dmax  = maximum agglomeration distance
%   pflag = plot points and clusters if nonzero
%
% bob newell, june 1997
%
%--------------------------------------------------------
function [ xhat, shape, alloc ] = ...
                    ellclust( data, dmin, dmax, pflag )
%
[ n, m ] = size( data );
alloc = zeros( n, 1 );
dagg = dmin;
minpts = 2;
dmult = 1.2;
if pflag & m == 2,
  plot( data(:,1), data(:,2), '+' );
end;
%
% find clusters encompassing all points
k = 0;
a = find( alloc == 0 );
while ~isempty(a),
  k = k + 1;
%
% select an unallocated point
  xhat(k,:) = data( a(1),: );
%
% agglomerate
  more = 1;
  while more & ~isempty(a),
    more = 0;
    for i = 1:length(a),
      d2 = (data(a(i),:) - xhat(k,:)) ...
                          * (data(a(i),:) - xhat(k,:))';
      if d2 <= dagg,
        alloc(a(i)) = k; more = 1;
      end;
    end;
%
% calculate ellipsoid
    a = find( alloc == k );
    if length(a) == 1,
      xhat(k,:) = data(a,:);
    else
      xhat(k,:) = mean( data(a,:) );
    end;
    a = find( alloc == 0 );
  end;
end;
%
% deallocate small clusters
[ c m ] = size( xhat );
for k = c:-1:1,
  a = find( alloc == k );
  if length(a) < minpts,
    alloc(a) = zeros(size(a));
    for i = 1:n,
      if alloc(i) > k,
        alloc(i) = alloc(i) - 1;
      end;
    end;
    xhat(k,:) = [];
  end;
end;
%
% try secondary allocation of stray points
[ c m ] = size( xhat );
a = find( alloc == 0 );
while length(a) > 0 & dagg < dmax,
  for i = 1:length(a),
    for k = 1:c,
      d2 = (data(a(i),:) - xhat(k,:)) ...
                          * (data(a(i),:) - xhat(k,:))';
      if d2 <= dagg,
        alloc(a(i)) = k;
      end;
    end;
  end;
  a = find( alloc == 0 );
  dagg = dagg * dmult;
end;
%
% final calculation of position and shape
shape = [];
xhat = [];
for k = 1:c,
  a = find( alloc == k );
  if length(a) == 1,
    xhat = [ xhat ; data(a,:) ];
  else
    xhat = [ xhat ; mean( data(a,:) ) ];
  end;
  [ vec, val ] = eig( cov( data(a,:) ) );
  shape = [ shape ; xhat(k,:) ...
                  1.5*sqrt(diag(val)') vec(:,1)' ];
  if pflag & m == 2,
    xc = shape(k,1); yc = shape(k,2);
    rx = shape(k,3); ry = shape(k,4);
    xr = shape(k,5); yr = shape(k,6);
    hold on;
    plot( xc, yc, 'o' );
    b = atan2( yr, xr );
    if b < 0, b = b + 2 * pi; end;
    x1 = rx*cos(b)+xc; y1 = rx*sin(b)+yc;
    for a = 0.1:0.1:2*pi,
      x2 = rx * cos(a); y2 = ry * sin(a);
      rr = sqrt(x2*x2+y2*y2); ar = atan2(y2,x2);
      if ar < 0, ar = ar + 2*pi; end;
      x2 = rr*cos(ar+b)+xc; y2 = rr*sin(ar+b)+yc;
      line( [x1 x2], [y1 y2]);
      x1 = x2; y1 = y2;
    end;
    x2 = rx*cos(b)+xc; y2 = rx*sin(b)+yc;
    line( [x1 x2], [y1 y2]);
    hold off; 
  end;
end;
%
