function [f,v] = fisofast(x,f0,m,e,s)
% FISODATA Fuzzy isodata clustering algorithm (faster version, with
%          a fixed distance measure and center calculation).
%
% [F V] = FISODATA(X,F0,m,e,s)
%     Input:
%       X  ... M by N data matrix, M is the number of data points
%              and N the dimension of X
%       F0 ... either an initial fuzzy partition matrix, or the
%              number of clusters. In the latter case, a default
%              partition matrix is generated.
%       m  ... (optional) parameter m > 1, determines the "fuzziness"
%              of clustering, for m close to 1 clusters become crisp
%              default value is 2.
%       e  ... (optional) termination tolerance, the algorithm stops
%              when |F(k-1) - F(k)| <= e, default 1e-3.
%       s  ... (optional) parameter for plotting intermediate results
%              (only for 1D and 2D data), default 0 - i.e. no plot, set
%              to 1 to view the clustering process on-line.
%     Output:
%       F  ... fuzzy partition matrix
%       V  ... centroid prototype matrix

% (c) Robert Babuska, June 1993

[mx,nx] = size(x);

if (nargin < 3), m = 2; elseif isempty(m), m = 2; end;
if (nargin < 4), e = 1e-3; elseif isempty(e), e = 1e-3; end;
if (nargin < 5), s = 0; elseif isempty(e), s = 0; end;

[mf0,c] = size(f0);
if (c==mf0 & c==1)
   c = f0;
% Initialize fuzzy partition matrix
   f0  = zeros(c,mx);
   for j=1:mx
     cc = 1 + fix((j-1)/round(mx/c));
     while (cc > c), cc = cc - 1; end
     f0(cc,j) = 0.6+0.4*rand;
     for i=1:c
       sumf0 = sum(f0(:,j));
       if (i ~= cc)
         f0(i,j) = 1-sumf0;
         if (i ~= c), f0(i,j) = f0(i,j)*rand; end;
       end
       if (c ~= 1 & cc==c & i==c), f0(i-1,j) = f0(i-1,j) + 1 - sumf0; end
     end
   end
   f0 = (f0.*(f0 > 0))';
end
f = zeros(mx,c);

if (s ~= 0)  &  (nx < 3),
    v = ones(c,1)*mean(x);
    if (nx == 1) x1 = [x ones(mx,1)]; v1 = [v ones(c,1)];
    else x1 = x; v1 = v; end;
    H = plot(x1(:,1),x1(:,2),'b.',v1(:,1),v1(:,2),'r+','EraseMode','xor','markersize',15,'linewidth',3);
    minx = min(x1(:,1)); maxx = max(x1(:,1));
    miny = min(x1(:,2)); maxy = max(x1(:,2));
    ma = 0.2*max(max(abs(x1)));
    axis([minx-ma maxx+ma miny-ma maxy+ma]);
    set(H(2),'LineWidth',3);
    set(gcf,'UserData',H);
    axis([-4 4 -2 2]);
    set(gca,'xticklabel',[],'yticklabel',[])
    set(gca,'xtick',[],'ytick',[])
    xlabel('x','fontname','times','fontsize',14); 
    ylabel('y','fontname','times','fontsize',14);
end;

% Iterate
iter = 0;
while  max(max( (f0-f) > e ))
   iter = iter + 1;
  f = f0;
% Calculate centers
  fm = f.^m;
  if (nx == 1)
    v = (sum(x*ones(1,c).*fm) ./ sum(fm))';
  else
    v = zeros(c,nx);
    for i = 1 : nx
        v(:,i) = (sum(x(:,i)*ones(1,c).*fm) ./ sum(fm))';
    end;
  end;
% Calculate distances
  if (nx == 1)
    d = (x*ones(1,c) - ones(mx,1)*v').^2;
  else
    d = zeros(mx,c);
    for i = 1 : c
        d(:,i) = sum((x - ones(mx,1)*v(i,:))'.^2)';
    end;
  end;
% Update f0
  d = d.^(1/(m-1));
  sumd = sum((d.^(-1))')'*ones(1,c);
  f0 = (d .* sumd).^(-1);
% Plot intermediate results
  if (s ~= 0), if (nx < 3)
    if (nx == 1) x1 = [x ones(mx,1)]; v1 = [v ones(c,1)];
    else x1 = x; v1 = v; end;
    H = get(gcf,'UserData');
    set(H(2),'xdata',v1(:,1),'ydata',v1(:,2)); drawnow;
  end; end;
if iter == 1, 
  disp('Hit any key to start ...');
  pause;
else
  pause(1)
end;
end
f = f0;

