function [triout, xout, yout, zout] = triage(tri, x, y, z)

% triage -- Triangle repair.
%  [triout, xout, yout, zout] = triage(tri, x, y, z) performs
%   triage on the triangles "tri", whose vertices lie at (x, y, z),
%   so that the Matlab "tsearch" function can be applied, despite
%   concavities.  The algorithm is:
%      1. Cast existing "tri" to counter-clockwise sense.
%      2. Fill all concavities and holes with clockwise
%            triangles that have at least one NaN vertex.
%   If a "z" vector is given, it will be made compatible
%   with the (xout, yout) values.  Added boundary points
%   will be interpolated from the original z, and all
%   supplemental interior points will be set to NaN.
%   The output consists of the original triangles and
%   data at the top, followed by the supplemental items
%   at the bottom.
%  triage(N) demonstrates itself with N vertices (default = 32).
 
% 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 09-Jul-1999 04:43:15.
% Updated    23-Jul-1999 16:22:12.

if nargin < 1, tri = 'demo'; end
if isequal(tri, 'demo'), tri = 32; end
if ischar(tri), tri = eval(tri); end

if length(tri) == 1
	nvertices = tri;
	scale = nvertices^(1/4);
	rad = 1 + rand(nvertices, 1) / scale;
	ang = sort(rand(size(rad)) * 2 * pi);
	xy = rad .* exp(ang*sqrt(-1));
	x = real(xy); y = imag(xy);
	
	[tt, xx, yy, ii] = tripoly(x, y);
	
	new = ones(size(x));
	new(ii) = 0;
	new = find(new);
	xnew = xx(new); ynew = yy(new);
	
	s = trisense(tt, xx, yy);
	ccw = find(s > 0);
	cw = find(s < 0);
	if any(cw), tt(cw, :) = []; end
	[to, xo, yo] = feval(mfilename, tt, xx, yy);
	delete(get(gca, 'Children'))
	hold off
	h = tripatch(to, xo, yo);
	theBDF = [mfilename ' ' int2str(nvertices)];
	set(h, 'ButtonDownFcn', theBDF)
	title(theBDF)
	hold on, plot(x, y, 'o', xnew, ynew, '+'), hold off
	set(gca, 'XLim', [-2 2], 'YLim', [-2 2])
	figure(gcf)
	zoomsafe
	return
end

if nargin < 3, help(mfilename), return, end
if nargin < 4, z = zeros(size(x)); end

x = x(:); y = y(:); z = z(:);

trinew = [];

s = trisense(tri, x, y);
cw = (s < 0); ccw = (s > 0);
tri(cw, :) = fliplr(tri(cw, :));

%	disp(' ## trihull...')
[h, overlapping] = trihull(tri);

if ~overlapping
%	disp(' ## tripoly...')
	h = [NaN; h(:); NaN];
	f = find(isnan(h));
	
	kstart = 2;
	kstop = length(f);
	
	for k = kstart:kstop
		
		remaining = length(f)-k+1;
		if 0 & k == 2 | rem(remaining, 10) == 0
			disp([' ## Remaining: ' int2str(remaining) ' segments.'])
		end
		
		hi = h(f(k-1)+1:f(k)-1);
		xhi = x(hi); yhi = y(hi);
		[ti, xi, yi, ind] = tripoly(xhi(:), yhi(:));

% Renumber and rearrange with old on top and new below.

		inserted = length(xi) - length(xhi);
		if any(inserted)
			disp([' ## "tripoly" has inserted ' int2str(inserted) ' points.'])
		
			old = ind;
		
			new = ones(size(xi));
			new(old) = 0;
			new = find(new);
		
			htemp = [old(:); new(:)];
			xi = xi(htemp);
			yi = yi(htemp);
			
			[ignore, hinv] = sort(htemp);
			ti = hinv(ti);
			
			old = 1:length(old);
			new = length(old) + (1:length(new));
			old = old(:);
			new = new(:);
		
			hi_old = hi;
			hi_new = length(x) + [1:length(new)];
			hi = [hi_old(:); hi_new(:)];
		
			xi_new = xi(new);
			yi_new = yi(new);
			x = [x; xi_new(:)];
			y = [y; yi_new(:)];
		end

		ti = hi(ti);   % Map to original indices.
		
		[m, n] = size(ti);
		if n == 1, ti = ti.'; end
		s = trisense(ti, x, y);
		cw = (s < 0); ccw = (s > 0);
		a(k-1) = areaxy(xi, yi);
		if a < 0   % cw
			ti = ti(ccw, :);
		else   % ccw.
			ti = ti(cw, :);
		end
		trinew = [trinew; ti];
	end
end

% Interpolate z for any new (x, y) points.

if nargin > 3 & length(z) < length(x)
	inserted = (length(z)+1):length(x);
	z(length(x)) = 0;
	z(inserted) = triterp(trinew, x, y, z, x(inserted), y(inserted));
end

xnew = []; ynew = []; znew = [];

if ~isempty(trinew)
	[m, n] = size(trinew);
	if n == 1, trinew = trinew.'; end
	xnew = x(trinew);
	ynew = y(trinew);
	[m, n] = size(xnew);
	if n == 1, xnew = xnew.'; ynew = ynew.'; end
	xnew = xnew * [1;1;1] / 3;
	ynew = ynew * [1;1;1] / 3;
	znew = zeros(size(xnew)) + NaN;
end

xout = [x; xnew];
yout = [y; ynew];
if nargin > 3, zout = [z; znew]; end

if ~isempty(trinew)
	k = (length(x) + (1:length(xnew))).';
	trinew = [[trinew(:, [1 2]); trinew(:, [2 3]); trinew(:, [3 1])] [k; k; k]];
	s = trisense(trinew, xout, yout);
	ccw = (s > 0);
	if any(ccw)
		trinew(ccw, :) = fliplr(trinew(ccw, :));
	end
end

if nargout > 0
	triout = [tri; trinew];
end