function [tout, xout, yout, ind] = tripoly(x, y, doPlot)

% tripoly -- Triangulation of a polygon.
%  [tout, xout, yout, ind] = tripoly(x, y, doPlot) modifies polygon
%   points (x, y), such that its triangulation tout consists solely
%   of strictly interior and exterior triangles, identifiable by their
%   counter- clockwise or clockwise sense, respectively.  Use the
%   "trisense" function to distinguish between them.  The returned
%   ind vector contains indices for reconstructing the original
%   polygon, using x = xout(ind) and y = yout(ind).  Points not
%   listed in ind are "fake", i.e. added to achieve the goal of
%   this routine.
%  tripoly(nVertices, doPlot) demonstrates itself with nVertices
%   (random).  Clicking on any triangle invokes the demonstration
%   again with different data.  The "fake" points are marked by
%   a solid symbol.  The plot is zoomable via "zoomsafe".
 
% 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 07-Jul-1999 09:15:01.
% Updated    29-Jul-1999 08:51:57.

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

if nargin < 3, doPlot = 0; end

if length(x) == 1
	if nargin > 1
		doPlot = y;
	elseif nargout < 1
		doPlot = 1;
	else
		doPlot = 0;
	end
	nVertices = x;
	scale = nVertices^(1/3);
	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);
end

xi = x(:); yi = y(:);
[zi, indices] = uniq([xi yi]);
xi = zi(:, 1); yi = zi(:, 2);

if length(xi) ~= length(x)
	warning(' ## Some data were culled because of non-uniqueness.')
end

% Make the polygon counter-clockwise.
%  Perhaps this step should not be done.

if (0)
	if areaxy(xi, yi) < 0
		xi = flipud(xi);
		yi = flipud(yi);
	end
end

% Algorithm: after triangulating the (x, y) points,
%  search for polygon edges that do not appear in
%  the triangulation itself.  Bisect those orphaned
%  edges, re-triangulate, then test again, iterating
%  until no such edges remain.  This scheme will
%  converge so long as the poly-line does not intersect
%  itself.

% Eventually, we will stamp out the supplemental points,
%  ideally producing a triangularization from the original
%  points alone.  This will be done by combining adjacent
%  triangles into a larger one where possible.  Where not,
%  we will flip diagonals in existing quadralaterals, thereby
%  eliminating all edges that refer to one or more supplemental
%  points.  Those points can then be deleted.

% disp(' ')

% tic

added = zeros(length(xi), 1);

orphans = 1;

iteration = 0;
while any(orphans)
	iteration = iteration + 1;
	tri = delaunay(xi, yi);
	tri = sort(tri.').';
	edges = [tri(:, [1 2]); tri(:, [1 3]); tri(:, [2 3]); [length(xi) length(xi)]];
	s = sparse(edges(:, 1), edges(:, 2), ones(size(edges(:,1))));
	orphans = zeros(1, length(xi));
	k = [1:length(xi) 1];
	len = length(k);
	for i = len:-1:2
		if i == len
			found = s(k(i), k(i-1));
		else
			found = s(k(i-1), k(i));
		end
		orphans(k(i-1)) = ~found;
	end
	
	f = find(orphans);
	if any(f)
		added(end+1) = 0;
		added = [added ones(size(added))].';
		k = 1:length(xi)+1;
		ki = zeros(2, length(xi)+1);
		ki(1, :) = 1:length(xi)+1;
		ki(2, f) = f+0.5;
		added = added(logical(ki));
		added = added(:);
		added(end) = [];
		ki = ki(logical(ki));
		xtemp = xi(:); xtemp(end+1) = xi(1);
		ytemp = yi(:); ytemp(end+1) = yi(1);
		xi = interp1(k, xtemp, ki, '*linear');
		yi = interp1(k, ytemp, ki, '*linear');
		xi(end) = [];
		yi(end) = [];
		xi = xi(:);
		yi = yi(:);
	end
	
	if any(orphans) & rem(iteration, 10) == 0
		disp([' ## ' mfilename ': orphans: ' int2str(sum(orphans))...
				'; length: ' int2str(length(xi))]);
	end
end

tri = trisort(tri);

% To polish the solution, traverse around the original polygon,
%  seeking each "added" point in turn.  If its geometry is
%  "simple", then the two adjacent interior triangles can be
%  squeezed together.  Iterate until no such situations remain.
%  If we are lucky, the interior will thereby become devoid
%  of added triangles.

% disp([' ## Elapsed time: ' num2str(toc) ' seconds.'])

if nargout > 0
	tout = tri; xout = xi; yout = yi; ind = find(~added);
end

if ~doPlot, return, end

disp([' ## Added ' int2str(sum(~~added)) ' points.'])

red = [1 0 0];
green = [0 1 0];
blue = [0 0 1];
white = red + green + blue;
black = 0*white;

delete(get(gca, 'Children'))
hold off
h = tripatch(tri, xi, yi);
theBDF = [mfilename ' ' int2str(length(x))];
set(h(~isnan(h)), 'ButtonDownFcn', theBDF, 'EdgeColor', white)
hold on
f = find(added);
plot(xi(f), yi(f), 'k+', 'MarkerFaceColor', 'k');
hold off
view(2)
title(theBDF)
set(gca, 'XLim', [-2 2], 'YLim', [-2 2])
figure(gcf)
zoomsafe

if (0)

theSense = trisense(tri, xi, yi);
[m, n] = size(tri);

delete(get(gca, 'Children'))
hold off

red = [1 0 0];
green = [0 1 0];
blue = [0 0 1];
white = red + green + blue;
black = 0*white;

for i = 1:m
	switch theSense(i)
	case -1
		theColor = red;
	otherwise
		theColor = green;
	end
	h(i) = patch(xi(tri(i, :)), yi(tri(i, :)), theColor);
	hold on
end

theBDF = [mfilename ' ' int2str(length(x))];
set(h(~isnan(h)), 'ButtonDownFcn', theBDF, 'EdgeColor', white)
f = find(added);
h = plot(xi(f), yi(f), 'ko', 'MarkerFaceColor', 'k');
set(h, 'LineWidth', 1.5)
f = find(~added);
f(end+1) = 1;
h = plot(xi(f), yi(f), '-ko');
set(h, 'LineWidth', 1.5)
hold off
view(2)
title(theBDF)
set(gca, 'XLim', [-2 2], 'YLim', [-2 2])
figure(gcf)
zoomsafe

end