function [tout, xo, yo, io] = triclean(tri, x, y, ind, verbose)

% triclean -- Triangle clean-up.
%  triclean(tri, x, y, ind) cleans-up the output of "tripoly"
%   by combining triangles and flipping quadralaterals that
%   involve "fake" points, until the interior triangles no
%   longer have any such points.  The triangle vertices that
%   are indexed in "tri" correspond to the polygon perimeter
%   (x, y); the original polygon points (as output by "tripoly")
%   are indexed in "ind".  The same information is returned,
%   now cleaned-up.  Only the counter-clockwise triangles remain.
 
% 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 26-Jul-1999 15:24:23.
% Updated    29-Jul-1999 14:07:29.

if nargin < 1, help(mfilename), return, end
if nargin < 5, isVerbose = 0; end

result = [];

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

if any(cw), tri(cw, :) = []; end

to = tri; xo = x; yo = y; io = ind;

tic

done = 0;
while ~done
	done = 1;

% In the first pass, we squeeze adjacent triangles
%  together, if they have a simple geometry.  The
%  slowest lines are "f = find(to...)" and the next.
%  Also the calls to "trisort" and "uniq".

to = [to(:, [1 2 3]); to(:, [1 3 2]); to(:, [2 1 3]); ...
		to(:, [2 3 1]); to(:, [3 1 2]); to(:, [3 2 1])];

	squeezed = 1;
	while any(squeezed)
		squeezed = 0;
		fake = ones(size(xo));
		fake(io) = 0;
		fake = find(fake);
		if isempty(fake), break, end
		if rem(length(fake), 10) == 0
			disp([' ## Remaining: ' int2str(length(fake))])
		end
		for k = 1:length(fake)
			node = fake(k);
			before = rem(node - 2 + length(xo), length(xo)) + 1;
			after = rem(node, length(xo)) + 1;
			f = find(to(:, 1) == before & to(:, 2) == node);   % Slow 23-31%.
			g = find(to(:, 1) == node & to(:, 2) == after);   % Slow 22-30%.
			for j = 1:length(g)
				for i = 1:length(f)
					if to(f(i), 3) == to(g(j), 3)
						opposite = to(f(i), 3);
						temp = sort([before opposite after]);
						indices = [1 2 3; 1 3 2; 2 1 3; ...
									2 3 1; 3 1 2; 3 2 1];
						tnew = temp(indices);
						if isVerbose
							disp([' ## Squeezed: ' int2str(node) ...
								' ' mat2str(to(f(i), :)) ...
								' ' mat2str(to(g(j), :)) ...
								' ' mat2str(tnew(1, :))])
						end
						h = find(any(to.' == node));
						to(h(1:6), :) = tnew;
						to(h(7:end), :) = [];
						xo(node) = [];
						yo(node) = [];
						io(io == node) = [];
						io(io > node) = io(io > node) - 1;
						to(to > node) = to(to > node) - 1;
						squeezed = 1;
						done = 0;
						break
					end
				end
				if squeezed, break, end
			end
			if squeezed, break, end
		end
	end
	
	to = trisort(to);   % Slow 7-13%.
	to = uniq(to);   % Slow 7-12%.
	
	s = trisense(to, xo, yo);
	cw = (s < 0);
	if any(cw), to(cw, :) = fliplr(to(cw, :)); end

% In the second pass, we flip one convex quadralateral
%  that involves a fake point.  Then, we go back to
%  pass #1.

	[m, n] = size(to);
	to = [to(:, [1 2 3]); to(:, [1 3 2]); to(:, [2 1 3]); ...
			to(:, [2 3 1]); to(:, [3 1 2]); to(:, [3 2 1])];
					
	indices = (1:m).';
	indices = indices * ones(1, 6);
	to = [to indices(:)];
	
	flipped = 1;
	
	while any(flipped)
		flipped = 0;
		flag = 0;
		fake = ones(size(xo));
		fake(io) = 0;
		fake = find(fake);
		if isempty(fake), break, end
		for k = 1:length(fake)
			node = fake(k);
			f = find(to(:, 1) == node);
			for j = 1:length(f)-1
				for i = j+1:length(f)
					if to(f(i), 2) == to(f(j), 2)
						a = node;
						b = to(f(i), 3);
						c = to(f(i), 2);
						d = to(f(j), 3);
						abcd = [a b c d];
						xtemp = xo(abcd);
						ytemp = yo(abcd);
						if isconvex(xtemp, ytemp)
							if isVerbose
								disp([' ## Flipped: ' mat2str([a b c d])])
							end
							triangle = to(f(i), 4);
							h = find(to(:, 4) == triangle);
							to(h(1), :) = [a b d triangle];
							triangle = to(f(j), 4);
							h = find(to(:, 4) == triangle);
							to(h(1), :) = [b c d triangle];
							flipped = 1;
							done = 0;
							break
						end
					end
				end
				if flipped, break, end
			end
			if flipped, break, end
		end
		flipped = 0;   % Do just one.
	end
	
	to = to(1:m, 1:3);
	to = trisort(to);
	to = uniq(to);
	
	s = trisense(to, xo, yo);
	cw = (s < 0);
	if any(cw), to(cw, :) = fliplr(to(cw, :)); end
end

toc

if nargout > 0
	tout = to;
else
	assignin('caller', 'ans', to)
end

fake = ones(length(x), 1);
fake(ind) = 0;
fake = find(fake);

subplot(1, 2, 1)
delete(get(gca, 'Children'))
hold off
axis normal
tripatch(tri, x, y)
hold on
% plot(x(ind), y(ind), 'o-')
plot(x(fake), y(fake), '+')
hold off
zoomsafe out

fake = ones(length(xo), 1);
fake(io) = 0;
fake = find(fake);

subplot(1, 2, 2)
delete(get(gca, 'Children'))
hold off
axis normal
tripatch(to, xo, yo)
hold on
% plot(xo(io), yo(io), 'o-')
plot(xo(fake), yo(fake), '+')
hold off
zoomsafe out

zoomsafe all