function t = delaunayn(x)
%DELAUNAYN  N-D Delaunay tessellation.
%   T = DELAUNAYN(X) computes a set of simplices such that no data 
%   points of X are contained in any circumspheres of the simplices. 
%   The set of simplices forms the Delaunay tessellation. X is 
%   an m-by-n array representing m points in n-D space. T is a 
%   numt-by-(n+1) array where each row is the indices into X of the
%   vertices of the corresponding simplex.
%
%   DELAUNAYN is based on Qhull.  
%
%   See also QHULL, VORONOIN, CONVHULLN, DELAUNAY, DELAUNAY3.

%   Copyright 1984-2000 The MathWorks, Inc.
%   $Revision: 1.11 $ $Date: 2000/08/28 15:53:27 $

[m,n] = size(x);
if m < n+1,
  error('Not enough points to do tessellation.');
end
if any(isinf(x(:)) | isnan(x(:)))
  error('Data containing Inf or NaN cannot be tessellated.');
end
if m == n+1
  t = 1:n+1;
  return;
end

% qhull needs the sum of squares at the end of the points
x = [x sum(x.*x,2)];
t = qhullmx(x', 'd ', 'QJ ', 'Pp');
x(:,end) = [];

% try to get rid of zero volumn simplices. They are generated
% because of the fuzzy jiggling.

[m, n] = size(t);
v = logical(ones(m,1));

seps = sqrt(eps)*max(abs(x(:)));
for i=1:m
  if abs(det(x(t(i,1:n-1),:)-x(t(i,n)*ones(n-1,1),:))) < seps
    v(i) = 0;
  end
end

t = t(v,:);
