/*
	tri.c
	Written by Joseph O'Rourke
	orourke@cs.smith.edu
	Triangulation code for "Computational Geometry in C."
        Assumes polygon vertices are entered in ccw order. 
*/
#include	<stdio.h>
#include	<math.h>
#define	ERROR	1
#define	X	0
#define	Y	1
typedef	enum { FALSE, TRUE }	bool;

#define DIM     2               /* Dimension of points */
typedef int     tPointi[DIM];   /* type integer point */
typedef double  tPointd[DIM];   /* type double point */

#define PMAX    1000            /* Max # of pts in polygon */
typedef tPointi tPolygoni[PMAX];/* type integer polygon */

int     Area2( tPointi a, tPointi b, tPointi c );
int     AreaPoly2( int n, tPolygoni P );
bool 	Xor( bool x, bool y );
bool    IntersectProp( tPointi a, tPointi b, tPointi c, tPointi d );
bool    Left( tPointi a, tPointi b, tPointi c );
bool    LeftOn( tPointi a, tPointi b, tPointi c );
bool    Collinear( tPointi a, tPointi b, tPointi c );
void	SubVec( tPointi a, tPointi b, tPointi c );
int     Dot( tPointi a, tPointi b );
int     Length2( tPointi p );
bool    Between( tPointi a, tPointi b, tPointi c );
bool    Intersect( tPointi a, tPointi b, tPointi c, tPointi d );
bool	Diagonalie( int i, int j, int n, tPolygoni P );
bool	Diagonal( int i, int j, int n, tPolygoni P );
void	PointAssign( tPointi a, tPointi b );
void	Triangulate1( int n, tPolygoni P );
void	Triangulate( int n, tPolygoni P );
void	TriRecurse( int n, tPolygoni P, int labels[] );
bool	InCone( int i, int j, int n, tPolygoni P );
void	ClipEar1( int i, int n, tPolygoni P, int labels[] );
void	ClipEar( int i, int n, tPolygoni P );
int	ReadPoints( tPolygoni P );
void	PrintPoly( int n, tPolygoni P, int labels[] );
void	PrintPoint( tPointi p );

main()
{
	int	n;
	tPolygoni	P;

	n = ReadPoints( P );
	Triangulate1( n, P );
	Triangulate( n, P );
}

/* 
        Returns twice the signed area of the triangle determined by a, b, c,
        positive if a, b, c are oriented ccw, and negative if cw.
*/
int     Area2( tPointi a, tPointi b, tPointi c )
{
/*	The text has this:
        return
                a[0] * b[1] - a[1] * b[0] +
                a[1] * c[0] - a[0] * c[1] +
                b[0] * c[1] - c[0] * b[1];

	The following computation is algebraically equivalent but
	uses four fewer multiplications.  It is obtained by shifting the
	coordinate system so that point a is the origin.
*/

	return
		(b[0] - a[0]) * (c[1] - a[1]) -
		(c[0] - a[0]) * (b[1] - a[1]);
}

/*      
        Returns twice the area of polygon P.
*/
int     AreaPoly2( int n, tPolygoni P )
{
        int     i;
        int     sum = 0;        /* Partial area sum */

        for (i = 1; i < n-1; i++)
                sum += Area2( P[0], P[i], P[i+1] );
        return  sum;
}

/*
	Exclusive or: true iff exactly one argument is true.
	The arguments are negated to ensure that they are 0/1
	values.  Then the bitwise Xor operator may apply.
	(This idea is due to Michael Baldwin.)
*/
bool 	Xor( bool x, bool y )
{
	return	!x ^ !y;
}
/*
	Returns true iff ab properly intersects cd: they share
	a point interior to both segments.  The properness of the
	intersection is ensured by using strict leftness.
*/
bool    IntersectProp( tPointi a, tPointi b, tPointi c, tPointi d )
{
	/* Eliminate improper cases. */
	if (
		Collinear(a,b,c) ||
		Collinear(a,b,d) ||
		Collinear(c,d,a) ||
		Collinear(c,d,b)
	   )
		return FALSE;

        return
                Xor( Left(a,b,c), Left(a,b,d) )
        &&      Xor( Left(c,d,a), Left(c,d,b) );
}

/*
	Returns true iff c is strictly to the left of the directed
	line through a to b.
*/
bool    Left( tPointi a, tPointi b, tPointi c )
{
        return  Area2( a, b, c ) > 0;
}

bool    LeftOn( tPointi a, tPointi b, tPointi c )
{
        return  Area2( a, b, c ) >= 0;
}

bool    Collinear( tPointi a, tPointi b, tPointi c )
{
        return  Area2( a, b, c ) == 0;
}

/*
        Returns T iff (a,b,c) are collinear and point c lies 
	on the closed segement ab.
*/
bool    Between( tPointi a, tPointi b, tPointi c )
{
        tPointi ba, ca;
        if ( ! Collinear( a, b, c ) )
                return  FALSE;

	/* If ab not vertical, check betweenness on x; else on y. */
	if ( a[X] != b[X] ) 
		return	((a[X] <= c[X]) && (c[X] <= b[X])) ||
			((a[X] >= c[X]) && (c[X] >= b[X]));
	else
		return	((a[Y] <= c[Y]) && (c[Y] <= b[Y])) ||
			((a[Y] >= c[Y]) && (c[Y] >= b[Y]));
}

/*
        Returns true iff segments ab and cd intersect, properly or improperly.
*/
bool    Intersect( tPointi a, tPointi b, tPointi c, tPointi d )
{
        if ( IntersectProp( a, b, c, d ) )
                return  TRUE;
        else if (       Between( a, b, c )
                     || Between( a, b, d )
                     || Between( c, d, a )
                     || Between( c, d, b )
                )
                return  TRUE;
        else    return  FALSE;
}

/*
	Returns T iff (v_i, v_j) is a proper internal *or* external
	diagonal of P, *ignoring edges incident to v_i and v_j*.
*/
bool	Diagonalie( int i, int j, int n, tPolygoni P )
{
	int	k;
	int	k1;

	/* For each edge (k,k+1) of P */
	for ( k = 0; k < n; k++ ) {
		k1 = (k+1) % n;
		/* Skip edges incident to i or j */
		if ( ! (
				( k == i ) || ( k1 == i )
			||	( k == j ) || ( k1 == j )
		       )
		   )
			if ( Intersect( P[i], P[j], P[k], P[k1] ) )
				return	FALSE;
	}
	return	TRUE;
}

/*
	Implements a = b, assignment of points/vectors.
	Assignment between arrays is not possible in C.
*/
void	PointAssign( tPointi a, tPointi b )
{
	int	i;

	for ( i = 0; i < DIM; i++ )
		a[i] = b[i];
}

/*
	Prints out n-3 diagonals (as pairs of integer indices)
	which form a triangulation of P.
*/
void	Triangulate1( int n, tPolygoni P )
{
	tPolygoni	Pt;
	int	labels[PMAX];
	int	i;

	for ( i = 0; i < n; i++ ){
		PointAssign( Pt[i], P[i] );
		labels[i] = i;
	}

	TriRecurse( n, Pt, labels );
}

void	TriRecurse( int n, tPolygoni P, int labels[] )
{
	int	i, i1, i2;

	printf("TriRecurse: n = %d\n", n);
	/*PrintPoly( n, P, labels );*/
	if ( n > 3 )
		for ( i = 0; i < n; i++ ) {
			i1 = (i+1) % n;
			i2 = (i+2) % n;
			if ( Diagonal( i, i2, n, P ) ) {
				printf("%d %d\n", labels[i], labels[i2]);
				ClipEar1( i1, n, P, labels );
				TriRecurse( n-1, P, labels );
				break;
			}
		}
}

void	PrintPoint( tPointi p )
{
	int	i;

	putchar('(');
	for ( i = 0; i < DIM; i++ ) {
		printf("%d", p[i]);
		if ( i != DIM-1 ) putchar(',');
	}
	putchar(')');
}

void	Triangulate( int n, tPolygoni P )
{
	int	i, i1, i2;

	if ( n > 3 )
		for ( i = 0; i < n; i++ ) {
			i1 = (i+1) % n;
			i2 = (i+2) % n;
			if ( Diagonal( i, i2, n, P ) ) {
				PrintPoint( P[i]  ); putchar('\t');
				PrintPoint( P[i2] ); putchar('\n');
				ClipEar( i1, n, P );
				Triangulate( n-1, P );
				break;
			}
		}
}

/*
	Returns true iff the diagonal (i,j) is striclty internal to the 
	polygon P in the neighborhood of the i endpoint.
*/
bool	InCone( int i, int j, int n, tPolygoni P )
{
	int	i1;	/* i + 1 */
	int	in1;	/* i - 1 */

	i1  = (i + 1)     % n;
	in1 = (i + n - 1) % n;

	/* If P[i] is a convex vertex [ i+1 left or on (i-1,i) ]. */
	if ( LeftOn( P[in1], P[i], P[i1] ) )
		return	Left( P[i], P[j], P[in1] )
		     &&	Left( P[j], P[i], P[i1]  );

	/* Assume (i-1,i,i+1) not collinear. */
	/* else P[i] is reflex. */
	else
		return	!(	LeftOn( P[i], P[j], P[i1]  )
		            &&	LeftOn( P[j], P[i], P[in1] ) );
}

/*
	Returns T iff (v_i, v_j) is a proper internal
	diagonal of P.
*/
bool	Diagonal( int i, int j, int n, tPolygoni P )
{
	return	InCone( i, j, n, P ) && Diagonalie( i, j, n, P );
}


/*
	Removes P[i] by copying P[i+1]...P[n-1] left one index.
*/
void	ClipEar1( int i, int n, tPolygoni P, int labels[] )
{
	int	k;
	
	for ( k = i; k < n-1; k++ ) {
		PointAssign( P[k], P[k+1] );
		labels[k] = labels[k+1];
	}
}

void	ClipEar( int i, int n, tPolygoni P )
{
	int	k;
	
	for ( k = i; k < n-1; k++ )
		PointAssign( P[k], P[k+1] );
}

int	ReadPoints( tPolygoni P )
/*
        Reads in the coordinates of the vertices of a polygon from stdin,
        puts them into P, and returns n, the number of vertices.
        The input is assumed to be pairs of whitespace-separated coordinates,
        one pair per line.  The number of points is not part of the input.
*/
{
	int	n = 0;

	printf("Polygon:\n");
	printf("  i   x   y\n");
	while ( (n < PMAX) && 
		(scanf("%d %d",&P[n][0],&P[n][1]) != EOF) ) {
		printf("%3d%4d%4d\n", n, P[n][0], P[n][1]);
		++n;
	}
	if (n < PMAX)
		printf("n = %3d vertices read\n",n);
	else	printf("Error in ReadPoints:\
			too many points; max is %d\n", PMAX);
	putchar('\n');

	return	n;
}

void	PrintPoly( int n, tPolygoni P, int labels[] )
{
	int	i;

	printf("Polygon:\n");
	printf("  i   l   x   y\n");
	for( i = 0; i < n; i++ )
		printf("%3d%4d%4d%4d\n", 
			i, labels[i], P[i][0], P[i][1]);
}
/*
        Puts a - b into c.
*/
void	SubVec( tPointi a, tPointi b, tPointi c )
{
        int     i;
        
        for( i = 0; i < DIM; i++ )
                c[i] = a[i] - b[i];
}

/*
        Returns the dot product of the two input vectors.
*/
int     Dot( tPointi a, tPointi b )
{
        int     i;
        int     sum = 0;

        for( i = 0; i < DIM; i++ )
                sum += a[i] * b[i];
        
        return  sum;
}

/*
        Returns the square of the length of the vector.
*/
int     Length2( tPointi p )
{
        return  Dot( p, p );
}