function tri = delaunay(x,y,flag)
%DELAUNAY Delaunay triangulation.
%   TRI = DELAUNAY(X,Y) returns a set of triangles such that no data
%   points are contained in any triangle's circumcircle. Each row of
%   the M-by-3 matrix TRI defines one such triangle and contains
%   indices into the vectors X and Y.  
%
%   To avoid the degeneracy of collinear data, DELAUNAY adds some random
%   fuzz to the data before calling DELAUNAYC. The default fuzz standard
%   deviation is chosen to maintain about 7 digits of accuracy in the data. 
%
%   TRI = DELAUNAY(X,Y,FUZZ) uses the specified value for the fuzz
%   standard deviation instead of the default value of 4*sqrt(eps).  It is
%   possible that there is no value of FUZZ that will produce a correct
%   triangulation.  In this unlikely situation, you will need to
%   preprocess your data to avoid collinear or nearly collinear data.
%
%   TRI  = DELAUNAY(X,Y,'sorted') assumes that the points X and Y are
%   sorted first by  Y and then by X, duplicate points have been
%   eliminated, and there are no collinear or nearly collinear data.  No
%   fuzz is added to the data.
%
%   The Delaunay triangulation is used with: GRIDDATA (to interpolate
%   scattered data), CONVHULL, VORONOI (to compute the VORONOI diagram),
%   and is useful by itself to create a triangular grid for scattered data
%   points.
%
%   The functions DSEARCH and TSEARCH search the triangulation to
%   find nearest neighbor points or enclosing triangles, respectively.
%   
%   Based on sweepline Voronoi code by Steven Fortune.
%   
%   See also DELAUNAYC, VORONOI, TRIMESH, TRISURF, GRIDDATA, CONVHULL, 
%            DSEARCH, TSEARCH, DELAUNAYN, QHULL.

%   Clay M. Thompson 7-8-95, 6-5-98
%   Copyright 1984-2000 The MathWorks, Inc.
%   $Revision: 1.18 $  $Date: 2000/03/17 16:14:10 $

tol = 4*sqrt(eps);
if nargin==3 & ~ischar(flag),
  tol = flag;
end

if nargin==2 | (nargin==3 & ~ischar(flag))
  S = rand('state'); % Remember state
  rand('state',0) % Always use the same random sequence

  x = x(:); y = y(:);
  xmin = min(x); xmax = max(x);
  ymin = min(y); ymax = max(y);
  
  % Remove duplicate points
  [xf,xe] = log2(x);
  [yf,ye] = log2(y);  
  [xy,ndx] = unique([xf xe yf ye],'rows');
  xf = xy(:,1); xe = xy(:,2);
  yf = xy(:,3); ye = xy(:,4);
  xy = [];

  % Add fuzz to mantissa
  %
  % This approach can fail if there are points so close together that
  % the fuzz changes the relative position of the points.  In these cases,
  % you may have to specify TOL to be smaller or do something else
  % special to make sure that there are no collinear points in the dataset.
  xf = xf + tol*(rand(size(xf))-.5)*2;
  yf = yf + tol*(rand(size(xf))-.5)*2; 
  xx = pow2(xf,xe);
  yy = pow2(yf,ye);

  rand('state',S); % reset state  

  % Offset data so it is centered around zero
  xx = xx - (xmin+xmax)/2;
  yy = yy - (ymin+ymax)/2;

  % Now do the triangulation
  tri = delaunayc(xx,yy);
  tri = reshape(ndx(tri),size(tri)); % Map to original indices

  % Remove any zero area triangles from list
  Area = abs((x(tri(:,2))-x(tri(:,1))) .* (y(tri(:,3))-y(tri(:,1))) - ...
             (x(tri(:,3))-x(tri(:,1))) .* (y(tri(:,2))-y(tri(:,1))))/2;
  tri((Area == 0),:) = [];
  
elseif nargin==3 & ischar(flag)
  tri = delaunayc(x,y,flag);
else
  error('Wrong number of input arguments.');
end
