function k = convhull(x,y,tri)
%CONVHULL Convex hull.
%   K = CONVHULL(X,Y) returns indices into the X and Y vectors of the
%   points on the convex hull.  
%
%   K = CONVHULL(X,Y,TRI) uses the triangulation (as obtained from
%   DELAUNAY) instead of computing it.
%
%   See also DELAUNAY, VORONOI, POLYAREA.

%   Clay M. Thompson 8-15-95
%   Copyright (c) 1984-98 by The MathWorks, Inc.
%   $Revision: 1.10 $  $Date: 1998/06/09 20:37:06 $

tol = 8*sqrt(eps);

isrow = ndims(x)==2 & size(x,1)==1 & size(x,2)~=1;
[msg,x,y] = xychk(x(:),y(:)); % Treat these a columns.
error(msg);

if isempty(x) | isempty(y), kk = []; return, end

if nargin==2, 
  tri = delaunay(x,y);
end

n = prod(size(x));
ntri = size(tri,1);
t = (1:ntri)';

% Triangle edge if T(i,j) is non-zero
T = sparse(tri(:,[1 1 2 2 3 3]),tri(:,[2 3 1 3 1 2]),1,n,n);
[i,j] = find(triu(T==1));
k = [i;j];

% Sort the points
xc = x(k); yc = y(k);
xm = mean(xc); ym = mean(yc);
phi = angle((xc-xm) + (yc-ym)*sqrt(-1));
[dum,ind] = unique(phi);
k = k(ind);

% Remove any collinear points on the convex hull
if ~isempty(k)
  xy = [x(k(end)) y(k(end));x(k) y(k);x(k(1)) y(k(1))];
  a = xy(2:end-1,:) - xy(1:end-2,:);
  b = xy(3:end  ,:) - xy(2:end-1,:);
  mag = sqrt(sum(a .* a,2) .* sum(b .* b,2));

  % a dot b - norm(a)*norm(b) < tol when a and b are aligned
  ind = abs(sum(a .* b,2) - mag) < tol & sum(a.*a,2) > tol;
  k(ind) = []; % Remove collinear points
end

if ~isempty(k), k = [k; k(1)]; end
if isrow, k = k'; end
