function [theHull, nonUnique] = trihull(theTri)

% trihull -- Indices of boundary of triangulation.
%  trihull(theTri) returns the indices of the bounding
%   polygon associated with theTri triangulation.  This
%   corresponds to the convex-hull if the triangles form
%   a filled convex set.  For triangulations that contain
%   holes, the individual segments are separated by NaNs.
%  [theHull, nonUnique] = trihull(theTri) also returns
%   nonUnique = TRUE (non-zero) if any point appears in
%   more than one segment.
 
% Copyright (C) 1999 Dr. Charles R. Denham, ZYDECO.
%  All Rights Reserved.
%   Disclosure without explicit written consent from the
%    copyright owner does not constitute publication.
 
% Version of 02-Jul-1999 17:02:50.
% Updated    09-Jul-1999 15:31:56.

if nargin < 1
	help(mfilename)
	theTri = 'demo';
end

if isequal(theTri, 'demo'), theTri = 10; end
if ischar(theTri), theTri = eval(theTri); end

if length(theTri) == 1
	nPoints = theTri;
	x = rand(nPoints, 1); y = rand(size(x));
	tri = delaunay(x, y);
	h = feval(mfilename, tri);
	trimesh(tri, x, y, x*0)
	view(2)
	px = x(h); py = y(h);
	px(end+1) = px(1); py(end+1) = py(1);
	hold on
	plot(px, py, '-bo')
	hold off
	f = findobj(gcf, 'Type', 'patch');
	theBDF = [mfilename ' ' int2str(nPoints)];
	title(theBDF)
	set(f, 'buttondownfcn', theBDF)
	set(gcf, 'Name', [mfilename ' Demo'])
	zoomsafe
	figure(gcf)
	if nargout > 0, theHull = h; end
	return
end

t = [theTri(:, [1 2]); theTri(:, [2 3]); theTri(:, [3 1])];

m = max(max(t)); n = m;

s = sparse(t(:, 1), t(:, 2), ones(size(t(:, 1))), m, n);

t = sort(t.').';
s1 = sparse(t(:, 1), t(:, 2), ones(size(t(:, 1))), m, n);

s1 = s1 + s1.';

[i, j] = find(s == 1 & s1 == 1);   % Too slow.

[i, indices] = sort(i);
j = j(indices);

connect(i) = j;

h = zeros(size(i));
f = find(connect);
f = f(1);
h(1) = connect(f);
connect(f) = 0;
starts(1) = 1;

for k = 2:length(h)
%	remaining = length(h)-k+1;
%	if rem(remaining, 1000) == 0 & k > 2
%		disp([' ## Remaining: ' int2str(remaining) ' edges.'])
%	end
	if connect(h(k-1))
		h(k) = connect(h(k-1));
		connect(h(k-1)) = 0;
	else
		f = find(connect);
		if any(f)
			f = f(1);
			h(k) = connect(f);
			connect(f) = 0;
			starts(end+1) = k;
		else
			break
		end
	end
end

if 0 & length(starts) > 1
	disp([' ## Segments: ' int2str(length(starts))])
end

% Check perimeters for points represented more than once.

htemp = h(h ~= 0);
s = sparse(htemp, htemp, ones(size(htemp)));
[i, j] = find(s > 1);
nonUnique = any(i);
if nonUnique & nargout < 2
	disp([' ## ' mfilename ': Some perimeters overlap.'])
end

if length(starts) > 1
	htemp = [];
	starts(end+1) = length(h) + 1;
	for k = 1:length(starts)-1
		i = starts(k):starts(k+1)-1;
		if k > 1, htemp = [htemp; NaN]; end
		htemp = [htemp; h(i)];
	end
	h = htemp;
end

if nargout > 0
	theHull = h;
else
	assignin('caller', 'ans', h)
end