function idx = grabk(xx, K)
%GRABK points using scheme for farthest-point clustering
%
% idx = grabk(xx, K);
%        xx DxN input points
%         K 1x1 number of clusters
%
%       idx 1xK indices of points
%
% Update: FPC, a farthest-point clustering function, can be used instead of
% GRABK. The new function can also label points (optional second output).
% This function may be deprecated and eventually removed from the imurray-matlab
% toolbox.

% Iain Murray, August 2008

[D, N] = size(xx);

% Naive O(NK) algorithm. Good enough for many purposes.
centers = zeros(D, K);
idx = ones(1, K);

sq_length = sum(xx.*xx, 1)'; % Nx1

idx(1) = ceil(rand*N);
centers(:, 1) = xx(:, idx(1));
sq_dist_to_any = bsxfun(@minus, sq_length, xx'*(2*xx(:, idx(1)))) + sq_length(idx(1));

for kk = 2:K
    [dummy, idx(kk)] = max(sq_dist_to_any);
    centers(:, kk) = xx(:, idx(kk));
    sq_dist_to_new = bsxfun(@minus, sq_length, xx'*(2*xx(:, idx(kk)))) + sq_length(idx(kk));
    sq_dist_to_any = min(sq_dist_to_new, sq_dist_to_any);
end