source: prextra/modeclust.m @ 23

%MODECLUST Clustering by mode-seeking in feature space
% 
%       [LAB,J] = MODECLUST(A,K)
% 
% INPUT
%   A       Dataset
%   K       Number of neighbours to search for local mode (default: 10)
%
% OUTPUT
%   LAB     Cluster assignments, 1..K
%   J       Indices of modal samples
%
% DESCRIPTION
% A K-NN modeseeking method is used to assign objects to their nearest mode.
% The nearest neighbor implementation uses the ANN Matlab Wrapper package. 
% It should be in the path. If needed download it from
% http://webscripts.softpedia.com/scriptDownload/ANN-MATLAB-Wrapper-Download-33976.html
%
% LITERATURE
% Cheng, Y. "Mean shift, mode Seeking, and clustering", IEEE Transactions
% on Pattern Analysis and Machine Intelligence, vol. 17, no. 8, pp. 790-799,
% 1995.
% 
% SEE ALSO
% MAPPINGS, DATASETS, KMEANS, HCLUST, KCENTRES, PROXM

% Copyright: R.P.W. Duin, r.p.w.duin@prtools.org
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands


function [assign,J] = modeclust(a,k)

        prtrace(mfilename);
        
        % ANNQUERY needs to be in the path
        annquerycheck;

        % prepare data
        if (nargin < 2), k = 10; end
        m = size(a,1);
        b = (+a)';
        
        % Run the k-NN search, indices in J, distances in d
        [J,d]= annquery(b,b,k);

        % density estimate from distance to most remote neighbor
        f = 1./(max(d,[],1)+realmin)';
        
        % Find local indices of local density maxima in neighbourhood.
        [dummy,I] = max(reshape(f(J),size(J)));

        % Translate back to indices in all the data. N will contain the
        % dataset index of their initial mode estimate in the K-neighbourhood.
        N = J(I+[0:k:k*(m-1)]);

        % Climb the mode
        % Re-assign samples to the sample temporily mode is assigned to.
        % Iterate until assignments don't change anymore. Samples that then point 
        % to themselves are modes; all other samples point to the closest mode.

        M = N(N);
        while (any(M~=N))
                N = M; M = N(N);
        end

        % Use renumlab to obtain assignments 1, 2, ... and the list of unique
        % assignments (the modes).
        [assign,J] = renumlab(M');

return