/***********************************************************************/
/* Open Visualization Data Explorer                                    */
/* (C) Copyright IBM Corp. 1989,1999                                   */
/* ALL RIGHTS RESERVED                                                 */
/* This code licensed under the                                        */
/*    "IBM PUBLIC LICENSE - Open Visualization Data Explorer"          */
/***********************************************************************/

#include <dxconfig.h>


#define	TRIANGULATE_ALWAYS		1
#define	TRIANGULATE_DEBUG		1
#define	TRIANGULATE_MAIN		0
#define	FIX_TRIANGLE_ORIENTATION	1


#include <stdio.h>
#include <math.h>
#include <dx/dx.h>


#ifndef TRUE
#define	TRUE	1
#define	FALSE	0
#endif

#define	LOCAL_LOOPS	256

typedef struct
{
    float		u;
    float		v;
} Point2;

typedef struct
{
    Point2		ll;
    Point2		ur;
} BBox;

typedef struct BVertex
{
    struct BVertex	*next;		/* boundary pointers         */
    struct BVertex	*prev;
    struct Loop		*lp;		/* initial parent loop       */
    int			id;		/* point id                  */
    Point2		pos;		/* projection of vertex      */
    int			tried;		
    int			pass;
} BVertex;

typedef struct
{
    int			*faces;		/* the input faces		*/
    int			*loops;		/* the input loops		*/
    int			*edges;		/* the input edges		*/
    float		*points;	/* the input points		*/
    int			*map;		/* the faces to triangles map   */
    int			nfaces;		/* # of faces			*/
    int			nloops;		/* # of loops			*/
    int			nedges;		/* # of edges			*/
    int			npoints;	/* # of points			*/
    int			nDim;		/* dimensionality of points     */
} FLEP;

typedef struct Loop
{
    struct Loop		*next;
    struct Loop		*prev;
    int			nvert;		/* number of vertices        */
    int			*vids;		/* pointer to vertices       */
    BVertex		*base;		/* base of linked list       */
    int			merged;		/* this loop has been merged */
    int			sign;		/* direction of rotation     */
    BBox		box;		/* loop bounding box         */
} Loop;


/*
 * For the tree representing the loop nesting...
 */
typedef struct treenode *TreeNode;

struct treenode
{
    Loop 		*loop;
    TreeNode 		sibling;
    TreeNode 		children;
    int 		nChildren;
    int 		level;		
};

/*
 * To make this more efficient declare a register float _det_ in any
 * routine that uses this macro.
 */

static double	_det_;
#define DetSign(p0,p1,p2)\
(\
    _det_ = (double)(p0.u * ((double)p1.v - p2.v)) +\
	    (double)(p1.u * ((double)p2.v - p0.v)) +\
	    (double)(p2.u * ((double)p0.v - p1.v)),\
    _det_ < (float) 0.0 ? -1 : _det_ > (float) 0.0 ? 1 : 0\
)

#define Determinant(p0,p1,p2)\
    ( p0.u * (p1.v - p2.v) + p1.u * (p2.v - p0.v) + p2.u * (p0.v - p1.v))

int	_dxf_TriangleCount (int *f, int nf, int *l, int nl, int *e, int ne);
int	_dxf__TriangleCount (FLEP *in);
Error	_dxf_TriangulateFLEP (int *f, int nf, int *l, int nl, int *e, int ne,
		float *p, int np, Vector *normal, Triangle *tris, int *ntri,
		int *map, int nDim);
Error	_dxf__TriangulateFLEP (FLEP *in, Vector *normal, Triangle *tris, int *ntri,
			  FLEP *valid);



static void	FlipLoop	 (Loop *loops);
static int	IntersectsEdge	 (Loop *loops, BVertex *p, BVertex *q);
static void	LoopSign	 (Loop *lp);
static void	LoopBBox	 (Loop *lp);
static int	PointInTri	 (Point2 p, Point2 plast, Point2 pnext,
				  Point2 v0, Point2 v1, Point2 v2, int tsign);
static void	ProjectPoints	 (Point *normal, float *vertices,
			          int nloops, Loop *loops, int nDim);
static int	AddLinkNested	 (Loop *loops, BVertex *p, BVertex *spares,
				  int *nvert, BVertex **first, BVertex **end);
static int	AnyInsideNested (TreeNode holes, BVertex *p0, BVertex *p1,
				  BVertex *p2, int psign);
static Error	InitNestedLoops	 (int *nloops, Loop **loops, int np);
static Error	TriangulateNestedLoops (int nloops, Loop *loops,
				  float *vertices,
				  Point *normal, int np,
				  Triangle *tptr, int *rntri, int nDim);
static void	LoopSignNested	 (Loop *lp);
static int      AddLink          (Loop *loops, BVertex *p, BVertex **spares,
				  int *nvert, BVertex **first, BVertex **end);
static int      AnyInside        (Loop *loops, BVertex *p0, BVertex *p1,
                                  BVertex *p2, int psign);
static Error    InitLoops        (int *nloops, Loop **loops, BVertex *spares,
				  int np);
static Error	TriangulateLoops (int nloops, Loop *loops, float *vertices,
				  Point *normal, int np,
				  Triangle *tptr, int *rntri, int nDim);

extern Array    _dxf_FLE_Normals (int *f, int nf, int *l, int nl,
				  int *e, int ne, float *v, int nd);

#if 0
void ShowLoop(BVertex *);
#endif

static int srprint = 0;
#if TRIANGULATE_DEBUG
#define srprintf	if (srprint) printf
#else
#define	srprintf
#endif


int
_dxf_TriangleCount (int *f, int nf, int *l, int nl, int *e, int ne)
{
    FLEP	in;

    in.faces	= f;
    in.loops	= l;
    in.edges	= e;
    in.nfaces	= nf;
    in.nloops	= nl;
    in.nedges	= ne;

    return (_dxf__TriangleCount (&in));
}


int
_dxf__TriangleCount (FLEP *in)
{
    FLEP	lin;
    int		f;
    int		l, ls, le;
    int		e, es, ee;
    int		ntri	= 0;
    int		nvert;

    lin = *in;

    for (f = 0; f < lin.nfaces; f++)
    {
	nvert = 0;

	ls = lin.faces[f];
	le = f < lin.nfaces - 1 ? lin.faces[f + 1] : lin.nloops;

	for (l = ls; l < le; l++)
	{
	    es = lin.loops[l];
	    ee = l < lin.nloops - 1 ? lin.loops[l + 1] : lin.nedges;
	    nvert += ee - es;
	}

	nvert += 2 * (le - ls);
	ntri  += nvert - 2;
    }

    return (ntri);
}



extern Error
_dxfTriangulateField(Field field)
{
    Error	ret	= ERROR;
    Array	af, al, ae, ap, an;
    int		*f, *l, *e;
    int		nf, nl, ne, np;
    float 	*p;
    Vector	*n;
    Triangle	*tri	= NULL;
    int		ntri;
    int		newtri	= 0;
    Array	arr;
    Array	col;
    int		size;
    int		i, *map = NULL;
    Triangle	*src;
    Triangle	t;
    Array       sA=NULL, dA=NULL;
    char	*name;
    int		nDim;

    af = (Array) DXGetComponentValue (field, "faces");
    al = (Array) DXGetComponentValue (field, "loops");
    ae = (Array) DXGetComponentValue (field, "edges");
    ap = (Array) DXGetComponentValue (field, "positions");

    if (af == NULL || al == NULL || ae == NULL || ap == NULL)
    {
	DXSetError (ERROR_BAD_PARAMETER,
	      "missing at least one of faces/loops/edges/positions");
	goto cleanup;
    }

    f = (int    *) DXGetArrayData (af);
    l = (int    *) DXGetArrayData (al);
    e = (int    *) DXGetArrayData (ae);
    p = (float  *) DXGetArrayData (ap);

    DXGetArrayInfo (af, &nf,  NULL, NULL, NULL, NULL);
    DXGetArrayInfo (al, &nl,  NULL, NULL, NULL, NULL);
    DXGetArrayInfo (ae, &ne,  NULL, NULL, NULL, NULL);
    DXGetArrayInfo (ap, &np,  NULL, NULL, NULL, &nDim);

    if (nDim == 2 || nDim == 3)
    {
	an = (Array) DXGetComponentValue(field, "normals");
	if (! an)
	{
	    an = _dxf_FLE_Normals(f, nf, l, nl, e, ne, p, nDim);
	    if (! an)
		goto cleanup;
	    
	    if (! DXSetComponentValue(field, "normals", (Object)an))
	    {
		DXDelete((Object)an);
		goto cleanup;
	    }

	    if (! DXSetStringAttribute((Object)an, "dep", "faces"))
	       goto cleanup;
	}

	n = (Vector *) DXGetArrayData(an);
    }
    else
    {
	DXSetError(ERROR_INVALID_DATA, 
		"Faces/Loops/Edges vertices must be 2 or 3D");
	goto cleanup;
    }

    ntri = _dxf_TriangleCount (f, nf, l, nl, e, ne);
    size = ntri * sizeof (Triangle);
    tri = (Triangle *) DXAllocateLocal (size);
    if (tri == NULL)
	goto cleanup;

    map = (int *)DXAllocate(nf*sizeof(int));
    if (! map)
	goto cleanup;

    if (! _dxf_TriangulateFLEP (f, nf, l, nl, e, ne, 
				p, np, n, tri, &newtri, map, nDim))
	goto cleanup;
    
    if (newtri == 0)
    {
	ret = OK;
	goto cleanup;
    }

#define	STRING(_x)	((Object) DXNewString (_x))
    arr = DXNewArray (TYPE_INT, CATEGORY_REAL, 1, 3);
    if (arr == NULL)
	goto cleanup;
    if (! DXAddArrayData (arr, 0, newtri, (Pointer) tri))
	goto cleanup;
    if (! DXSetComponentValue (field, "connections", (Object) arr))
	goto cleanup;
    arr = NULL;
    if (! DXSetComponentAttribute (field, "connections", "element type",
				 STRING ("triangles")))
	goto cleanup;
    
    i = 0;
    while (NULL != (sA = (Array)DXGetEnumeratedComponentValue(field, i++, &name)))
    {
	if (! strcmp(name, "invalid faces"))
	{
	    int i, j, k;
	    InvalidComponentHandle ifaces, itris;

	    ifaces = DXCreateInvalidComponentHandle((Object)field,
							NULL, "faces");
	    itris  = DXCreateInvalidComponentHandle((Object)field,
							NULL, "connections");
	    
	    if (! ifaces || ! itris)
		goto cleanup;
	    
	    DXSetAllValid(itris);

	    for (i = k = 0; i < nf; i++)
	    {
		if (DXIsElementInvalid(ifaces, i))
		{
		    for (j = 0; j < map[i]; j++)
			DXSetElementInvalid(itris, k++);
		}
		else
		    k += map[i];
	    }

	    if (! DXSaveInvalidComponent(field, itris))
	    {
		DXFreeInvalidComponentHandle(itris);
		DXFreeInvalidComponentHandle(ifaces);
		goto cleanup;
	    }

	    DXFreeInvalidComponentHandle(itris);
	    DXFreeInvalidComponentHandle(ifaces);
	}
	else
	{
	    int itemSize;
	    Type type;
	    Category cat;
	    int rank, shape[30];
	    Object attr = DXGetComponentAttribute(field, name, "dep");

	    dA = NULL;

	    if (! attr)
		continue;
	    
	    if (strcmp(DXGetString((String)attr), "faces"))
		continue;
	    
	    DXGetArrayInfo(sA, NULL, &type, &cat, &rank, shape);
	    itemSize = DXGetItemSize(sA);

	    if (DXQueryConstantArray(sA, NULL, NULL))
	    {
		dA = (Array)DXNewConstantArrayV(newtri, 
			    DXGetConstantArrayData(sA), type, cat, rank, shape);
		if (! dA)
		    goto cleanup;
	    }
	    else
	    {
		int i, j;
		char *srcPtr;
		char *dstPtr;

		dA = DXNewArrayV(type, cat, rank, shape);
		if (! dA)
		    goto cleanup;
		
		if (! DXAddArrayData(dA, 0, newtri, NULL))
		    goto cleanup;
		
		dstPtr = (char *)DXGetArrayData(dA);
		srcPtr = (char *)DXGetArrayData(sA);

		for (i = 0; i < nf; i++)
		{
		    for (j = 0; j < map[i]; j++)
		    {
			memcpy(dstPtr, srcPtr, itemSize);
			dstPtr += itemSize;
		    }
		    srcPtr += itemSize;
		}
	    }

	    if (! DXSetComponentValue(field, name, (Object)dA))
		goto cleanup;
	    dA = NULL;

	    if (! DXSetComponentAttribute(field, name, "dep", 
					    (Object)DXNewString("connections")))
		goto cleanup;

	}
    }

    /*
     * Get rid of the old stuff
     */

    DXDeleteComponent (field, "invalid faces");
    DXDeleteComponent (field, "faces");
    DXDeleteComponent (field, "loops");
    DXDeleteComponent (field, "edges");
    if (! DXEndField (field))
	goto cleanup;

    ret = OK;

cleanup:
    if (newtri == 0 && ret == OK)
    {
	char *name;

	while (DXGetEnumeratedComponentValue (field, 0, &name))
	    DXDeleteComponent (field, name);
    }

    DXDelete((Object) dA);
    DXFree ((Pointer) tri);
    DXFree ((Pointer) map);
    return (ret);
}

Error
_dxf_TriangulateFLEP (int *f, int nf, int *l, int nl, int *e, int ne,
		float *p, int np, Vector *normal, Triangle *tris, int *ntri,
		int *map, int nDim)
{
    FLEP	in;

    in.faces	= f;
    in.loops	= l;
    in.edges	= e;
    in.points	= p;
    in.map	= map;
    in.nfaces	= nf;
    in.nloops	= nl;
    in.nedges	= ne;
    in.npoints	= np;
    in.nDim     = nDim;

    return (_dxf__TriangulateFLEP (&in, normal, tris, ntri, NULL));
}


Error
_dxf__TriangulateFLEP (FLEP *in, Vector *normal, Triangle *tris, int *ntri,
		  FLEP *valid)
{
    Vector      tri_direction;
    int         temp, k;
    float       dir;
    FLEP	lin;
    FLEP	lvalid;
    Error	ret	= ERROR;
    Loop	*loops	= NULL;
    int		size;
    Triangle	*ltris, *tmp_tris;
    int		lntri;
    int		ftri;
    int		f;
    int		l, ls, le;
    int		e, es, ee;
    int		skips;

    lin = *in;
    if (valid)
	lvalid = *valid;

    size  = lin.nloops * sizeof (Loop);
    loops = (Loop *) DXAllocateLocal (size);
    if (loops == NULL)
	goto cleanup;
    
    ltris = tris;
    lntri = 0;

    for (f = 0; f < lin.nfaces; f++)
    {
	if (valid && ! lvalid.faces[f])
	    continue;

	skips = 0;
	ls = lin.faces[f];
	le = f < lin.nfaces - 1 ? lin.faces[f + 1] : lin.nloops;

	memset (loops, NULL, (le - ls) * sizeof (Loop));

	for (l = ls; l < le; l++)
	{
	    if (valid && ! lvalid.loops[l])
	    {
		skips++;
		continue;
	    }

	    es = lin.loops[l];
	    ee = l < lin.nloops - 1 ? lin.loops[l + 1] : lin.nedges;

	    loops[l - ls - skips].nvert = ee - es;
	    loops[l - ls - skips].vids  = lin.edges + es;
	}

	if (getenv("DX_NESTED_LOOPS"))
	{
	    if (! TriangulateNestedLoops (le - ls - skips, loops, lin.points,
		    (normal) ? normal + f : NULL, lin.npoints, ltris, &ftri, lin.nDim))
	    {
#if TRIANGULATE_ALWAYS
		if (DXGetError () == ERROR_INTERNAL)
		    DXResetError ();
		else
#endif
		    goto cleanup;
	    }
	}
	else
	{
	    if (! TriangulateLoops (le - ls - skips, loops, lin.points,
		    (normal) ? normal + f : NULL, lin.npoints, ltris, &ftri, lin.nDim))
	    {
#if TRIANGULATE_ALWAYS
		if (DXGetError () == ERROR_INTERNAL)
		    DXResetError ();
		else
#endif
		    goto cleanup;
	    }
	}




#if FIX_TRIANGLE_ORIENTATION
        for (k=0, tmp_tris=ltris; k<ftri; k++) {
	  if (lin.nDim == 2)
	  {
	      float *p = lin.points + 2*tmp_tris->p;
	      float *q = lin.points + 2*tmp_tris->q;
	      float *r = lin.points + 2*tmp_tris->r;
	      float x0 = q[0] - p[0];
	      float y0 = q[1] - p[1];
	      float x1 = r[0] - q[0];
	      float y1 = r[1] - q[1];

	      tri_direction.x = 0.0;
	      tri_direction.y = 0.0;
	      tri_direction.z = x0*y1 - x1*y0;
	  }
	  else
	  {
	      Point *pts = (Point *)lin.points;

	      tri_direction = DXCross(DXSub(pts[tmp_tris->q], pts[tmp_tris->p]),
                                  DXSub(pts[tmp_tris->r], pts[tmp_tris->q]));

	  }

	  dir = DXDot(tri_direction, normal[f]);

          if (dir < 0) {
            temp = tmp_tris->r;
            tmp_tris->r = tmp_tris->q;
            tmp_tris->q = temp;
          }
          tmp_tris++;
        }

#endif

	
	lntri += ftri;
	ltris += ftri;
	in->map[f] = ftri;
    }

    ret = OK;

cleanup:
    if (ret == OK && ntri)
	*ntri = lntri;

    DXFree ((Pointer) loops);
    return (ret);
}


#define Shift \
{ \
    p0 = p1; \
    p1 = p2; \
    p2 = p2->next; \
    continue; \
}

#define EmitTri012 \
{ \
    i0 = p0->id; \
    i1 = p1->id; \
    i2 = p2->id; \
    if (i0 != i1 && i0 != i2 && i1 != i2) \
    { \
	tptr->p = p0->id; \
	tptr->q = p1->id; \
	tptr->r = p2->id; \
	tptr++; \
	ntri++; \
    } \
    p0->next = p2; \
    p2->prev = p0; \
}


/*
 * $$$$$ NOTE:  This code, as did the previous version from Jarek et al,
 * $$$$$        assumes that the first loop is the outer loop and that
 * $$$$$        all subsequent loops are holes _dxf_inside of it.
 */

#if 0
static int dumpLoop = 0;
static int loopLimit = 99999;

static void
DumpLoop(BVertex *p0)
{
    FILE *fd;
    int n = 0;
    BVertex *b0, *b1;
    
    fd = fopen("dump.out", "w");

    b0 = b1 = p0;
    do
    {
	n++;
	b0 = b0->next;
    } while (b0 != b1);
    fprintf(fd, "%d\n", n);
    b0 = b1 = p0;
    do
    {
	fprintf(fd, "    %d\n", b0->id);
	b0 = b0->next;
    } while (b0 != b1);

    fclose(fd);
}
#endif


static Error
TraverseTree(TreeNode node, int *rntri, Triangle *tbase)
{
    TreeNode child;

    if (! node)
	return OK;

    if ((node->level & 0x1) != 0)
    {
	/* 
	 * Then its an outer contour
	 */
	if (! TriangulateLoopWithHoles(node, rntri, tbase))
	    goto error;
    }

    for (child = node->children; child; child = child->sibling)
	if (! TraverseTree(child, rntri, tbase))
	    goto error;

    return OK;

error:
    return ERROR;
}

static Error
TriangulateLoopWithHoles(TreeNode node, int *rntri, Triangle *tbase)
{
    Loop *loop = node->loop;
    int nvert = loop->nvert, not_done;
    BVertex *p0 = loop->base, *p1, *p2;
    BVertex *end = NULL;
    BVertex *first = NULL;
    TreeNode child;
    int	i0, i1, i2;
    int nholes;
    TreeNode holes;
    BVertex *spares;
    int n;
    int tsign, psign = loop->sign;
    Triangle *tptr = tbase + *rntri;
    int ntri = *rntri;
    int pass = 0;
    Loop *l;

    if (! psign)
	return OK;

    /*
     * Get spares required for linking.  While you're at it, make sure
     * the holes are oriented opposite to the outer contour.  Also - create
     * linked list of loops
     */
    node->loop->next = node->loop->prev = node->loop;
    for (holes = node->children, nholes = 0; holes; holes = holes->sibling)
    {
	nholes ++;
	holes->loop->next = node->loop->next;
	node->loop->next = holes->loop;
	if (holes->loop->sign == psign)
	    FlipLoop(holes->loop);
    }
    
    /*
     * and doubly link the list 
     */
    l = node->loop;
    do {
	l->next->prev = l;
    } while ((l = l->next) != node->loop);

    spares = (BVertex *)DXAllocateZero(nholes * 2 * sizeof(BVertex));
    if (! spares)
	goto error;

    n = 0;
    while (1)
    {
	do
	{
	    not_done = 0;
	    pass ++;

	    if (first != NULL)
		p0 = first;
	    
	    if (first == end)
		end = NULL;

	    while (p0->pass < pass && p0 != end)
	    {
		p0->pass = pass;

		p1 = p0->next;
		p2 = p1->next;

		if (!p0->tried)
		{
		    srprintf ("DXTri: %2d %2d %2d (%3d)",
			   p0->id, p1->id, p2->id, nvert);

		    if (p0->id == p1->id || p1->id == p2->id ||
			(tsign = DetSign(p0->pos, p1->pos, p2->pos)) == 0)
		    {
			/*
			 * degenerate triangle.  eliminate p1 from the loop.
			 */
			p0->next = p2;
			p2->prev = p0;
			not_done++;
			nvert--;
		    }
		    else if (tsign != psign ||
			AnyInsideNested (node->children, p0, p1, p2, psign))
		    {
			p0->tried = 1;
		    }
		    else
		    {
			not_done++;
			p0->prev->tried = 0;
			if (p0 == first)
			    first = first->prev;
			if (p2 == end) 
			{
			    if (end == first)
				end = NULL;
			    else
				end = end->next;
			}

			if (p1 == first)
			    first = p0;

			EmitTri012;
			nvert--;

			while (nvert > 2)
			{
			    p1 = p0->next;
			    p2 = p1->next;

			    if (DetSign(p0->pos, p1->pos, p2->pos) == 0)
			    {
				p0->next = p2;
				p2->prev = p0;
				if (p1 == first)
				    first = p0;
				nvert --;
			    }
			    else
				break;
			}
		    }
		}

		p0 = p0->next;
	    }

	} while (not_done && nvert >= 3);

	if (! AddLinkNested(node->loop, p0, spares+2*n, &nvert, &first, &end))
	    break;
	
	n ++;
    }

    if (nvert == 3)
	EmitTri012;

    *rntri = ntri;

    DXFree((Pointer)spares);
    return OK;

error:
    DXFree((Pointer)spares);
    return ERROR;
}

static TreeNode CreateNestingTree(int nLoops, Loop *loops);

static Error TriangulateNestedLoops (int nloops, Loop *loops, float *vertices,
			       Point *normal, int np,
			       Triangle *tptr, int *rntri, int nDim)
{
    Error	ret	= ERROR;	/* return status                    */
    BVertex	*vbase	= NULL;		/* base of list blocks              */
    int		nvert;			/* total vertex count               */
    int		size;			/* size temporary                   */
    int		i, j;
    BVertex	*bptr;
    TreeNode    tree = NULL;

    if (rntri)
	*rntri = 0;

    for (nvert = 0, i = 0; i < nloops; i++)
	nvert += loops[i].nvert;

    size = nvert * sizeof (BVertex);
    vbase = (BVertex *) DXAllocateLocal (size);
    if (vbase == NULL)
	goto error;

    /*
     * For each loop:
     *   Initialize its boundary vertex list as a circularly linked list.
     *   For now we'll just put in the vertex id.  We'll fill in the 
     *   vertex positions later.  We separate this for two reasons.  First,
     *   it involves random accessing and we don't want to guarantee that
     *   we destroy any savings we get from the PVS's line cache, and 
     *   second we want to use tight loops for the projection from 3-space
     *   to 2-space.
     */
    bptr = vbase;
    for (j = 0; j < nloops; j++)
    {
	loops[j].base = bptr;

	for (i = 0, bptr = loops[j].base; i < loops[j].nvert; i++, bptr++)
	{
	    bptr->next   = bptr + 1;
	    bptr->prev   = bptr - 1;
	    bptr->lp     = loops + j;
	    bptr->id     = *loops[j].vids++;
	}

	/*
	 * Circularly link the list
	 */

	loops[j].base[loops[j].nvert - 1].next = loops[j].base;
	loops[j].base->prev = loops[j].base + (loops[j].nvert - 1);
    }

    ProjectPoints (normal, vertices, nloops, loops, nDim);

    /*
     * Get rid of colinear loops and loops with fewer than 3 points and
     * link all of the others together.
     */

    if (! InitLoopsNested (&nloops, &loops, np))
	goto error;

    tree = CreateNestingTree(nloops, loops);
    if (! tree)
	goto error;
    
    AssignTreeLevel(tree, 0);

    if (! TraverseTree(tree, rntri, tptr))
	goto error;

    ret = OK;

error:
cleanup:
    DXFree((Pointer)tree);
    DXFree ((Pointer) vbase);
    return (ret);
}


#if 1
static int nShow = 0;

void 
ShowLoop(BVertex *first)
{
    int i, n;
    BVertex *p0;
    char buf[32];
    FILE *foo;

    sprintf(buf, "bad%d.dx", nShow++);

    foo = fopen(buf, "w");
    if ( foo )
    {
	n = 0;
	p0 = first;
	do
	{
	    n++;
	} while ((p0 = p0->next) != first);

	fprintf(foo,
    "object \"faces\" class array type int rank 0 items 1 data follows\n0\nattribute \"ref\" string \"loops\"\n");

	fprintf(foo,
    "object \"loops\" class array type int rank 0 items 1 data follows\n0\nattribute \"ref\" string \"edges\"\n");

	fprintf(foo,
    "object \"edges\" class array type int rank 0 items %d data follows\n", n);
	for (i = 0; i < n; i++)
	    fprintf(foo, "%d\n", i);
	fprintf(foo, "attribute \"ref\" string \"positions\"\n");

	fprintf(foo,
    "object \"positions\" class array type float rank 1 shape 3 items %d data follows\n", n);
	for (i = 0, p0 = first; i < n; i++, p0 = p0->next)
	    fprintf(foo, "%f %f 0\n", p0->pos.u, p0->pos.v);
	
	fprintf(foo,
    "object \"data\" class array type int rank 0 items %d data follows\n", n);
	for (i = 0, p0 = first; i < n; i++, p0 = p0->next)
	    fprintf(foo, "%d\n", p0->id);
	fprintf(foo, "attribute \"dep\" string \"positions\"\n");
	
	fprintf(foo, "object \"fle\" class field\n");
	fprintf(foo, "    component \"faces\" \"faces\"\n");
	fprintf(foo, "    component \"loops\" \"loops\"\n");
	fprintf(foo, "    component \"edges\" \"edges\"\n");
	fprintf(foo, "    component \"positions\" \"positions\"\n");
	fprintf(foo, "    component \"data\" \"data\"\n");
	fprintf(foo, "end\n");

	fclose(foo);
    }
    
}
#endif

/*
 * Figure out what the smallest dimension of the polygon is, e.g. the
 * dimension that contains the portion of the normal vector with the
 * largest magnitude.  Now copy in the vertices, projecting them from
 * 3D space to the appropriate 2D space.
 */

static void ProjectPoints (Point *normal, float *vertices,
			   int nloops, Loop *loops, int nDim)
{
    Point	lnormal;
    int		elim;
    int		i, j;
    BVertex	*bptr;

    if (nDim == 2)
    {
	float	*pptr;
	Loop	*loop;

	for (j = 0, loop = loops; j < nloops; j++, loop++)
	    for (i = 0, bptr = loop->base; i < loop->nvert; i++, bptr++)
	    {
		pptr = vertices + 2*bptr->id;
		bptr->pos.u = pptr[0];
		bptr->pos.v = pptr[1];
	    }
    }
    else
    {
	Point   *vptr = (Point *)vertices;
	Point	*pptr;
	Loop	loop;

	lnormal = *normal;
	if (lnormal.x < (float) 0.0) lnormal.x = -lnormal.x;
	if (lnormal.y < (float) 0.0) lnormal.y = -lnormal.y;
	if (lnormal.z < (float) 0.0) lnormal.z = -lnormal.z;

	elim = 0;
	if (lnormal.y > lnormal.x)
	    elim = 1;
	if (lnormal.z > (elim == 0 ? lnormal.x : lnormal.y))
	    elim = 2;

	for (j = 0; j < nloops; j++)
	{
	    loop = loops[j];

	    switch (elim)
	    {
		case 0:	/* get rid of X */
		    for (i = 0, bptr = loop.base; i < loop.nvert; i++, bptr++)
		    {
			pptr = vptr + bptr->id;
			bptr->pos.u = pptr->y;
			bptr->pos.v = pptr->z;
		    }
		    break;

		case 1:	/* get rid of Y */
		    for (i = 0, bptr = loop.base; i < loop.nvert; i++, bptr++)
		    {
			pptr = vptr + bptr->id;
			bptr->pos.u = pptr->x;
			bptr->pos.v = pptr->z;
		    }
		    break;

		case 2:	/* get rid of Z */
		    for (i = 0, bptr = loop.base; i < loop.nvert; i++, bptr++)
		    {
			pptr = vptr + bptr->id;
			bptr->pos.u = pptr->x;
			bptr->pos.v = pptr->y;
		    }
		    break;
	    }
	}
    }
}


static PseudoKey
HashPoint(Key k)
{
    PseudoKey pk;

    Point2 *p = &((*(BVertex **)k)->pos);
    pk = (PseudoKey)(*(int *)((&p->u))*17 + *((int *)(&p->v))*53);

    return pk;
}

static int
ComparePoint(Key k0, Key k1)
{
    Point2 *p0 = &((*(BVertex **)k0)->pos);
    Point2 *p1 = &((*(BVertex **)k1)->pos);

    return ! (p0->u == p1->u && p0->v == p1->v);
}


#define A_INSIDE_B  	 1
#define B_INSIDE_A  	 2
#define A_OUTSIDE_B 	 3
#define AB_CROSS    	 4
#define AB_DISJOINT 	 5
#define AB_INDETERMINATE 6


static Error
_CheckNesting(Loop *A, Loop *B, int *result)
{
  BVertex *i, *o;
  int j, knt, inside = 0, outside = 0;
  float ui, uo0, uo1;
  float vi, vo0, vo1;
  float u, t;
  float mU, MU, mV, MV;

  /*
   * Determine the overlap region
   */
  if (A->box.ll.u > B->box.ll.u) mU = A->box.ll.u;
  else				 mU = B->box.ll.u;
  if (A->box.ll.v > B->box.ll.v) mV = A->box.ll.v;
  else				 mV = B->box.ll.v;
  if (A->box.ur.v < B->box.ur.v) MU = A->box.ur.u;
  else				 MU = B->box.ur.u;
  if (A->box.ur.v < B->box.ur.v) MV = A->box.ur.v;
  else				 MV = B->box.ur.v;


  for (i = A->base, j = 0; j < A->nvert; j++, i = i->next)
  {
    ui = i->pos.u;
    vi = i->pos.v;

    if (ui < mU || ui > MU || vi < mV || vi > MV)
    {
      outside ++;
      continue;
    }

    o = B->base;
    knt = 0;
    do
    {
      uo1 = o->pos.u;
      uo0 = o->next->pos.u;
      vo1 = o->pos.v;
      vo0 = o->next->pos.v;

      /*
       * Orient the edge endpoints so that vo0 <= vo1
       */
      if (vo0 > vo1)
      {
	t = vo0; vo0 = vo1; vo1 = t;
	t = uo0; uo0 = uo1; uo1 = t;
      }

      /*
       * if the edge is horizontal, forget it.
       */
      if (vo0 == vo1)
	continue;
	  
      /* 
       * if both endpoints of the edge lie to the left of
       * the point, forget it.
       */
      if (ui > uo0 && ui > uo1)
	continue;

      /*
       * If the edge lies entirely above or below the point,
       * forget it.
       */
      if (vi > vo1 || vi < vo0)
	continue;
      
      /*
       * If the v exactly coincides with the upper endpoint, count
       * it a cross.  If it exactly coincides with the lower, skip it.
       */
      if (vi == vo1 && ui < uo1)
      {
	knt ++;
	continue;
      }
      else if (vi == vo0)
	  continue;
      
      /*
       * Get the edge value of u at v = vi. If thats to the right
       * of the test point, count it a cross.
       */
      t = (vi - vo0) / (vo1 - vo0);
      u = uo0 + t*(uo1 - uo0);
      if (u > ui)
      {
	knt ++;
	continue;
      }
    }
    while ((o = o->next) != B->base);

    if (knt & 0x1)
      inside ++;
    else
      outside ++;
  }

  if (! inside && ! outside)
  {
      DXSetError(ERROR_INTERNAL, "empty loop?");
      return ERROR;
  }

  if (inside && !outside)
    *result = A_INSIDE_B;
  else if (!inside && outside)
    *result = A_OUTSIDE_B;
  else if (!inside && !outside)
    *result = AB_CROSS;

  return OK;
}

static Error
_CheckNestingBoxes(Loop *A, Loop *B, int *result)
{
    if (A->box.ll.u > B->box.ur.u) *result = AB_DISJOINT;
    else if (A->box.ll.v > B->box.ur.v) *result = AB_DISJOINT;
    else if (A->box.ur.u < B->box.ll.u) *result = AB_DISJOINT;
    else if (A->box.ur.v < B->box.ll.v) *result = AB_DISJOINT;
    else *result = AB_INDETERMINATE;

    return OK;
}


static Error
CheckNesting(Loop *A, Loop *B, int *result)
{
  int ab;
  
   
  /*
   * If boxes don't overlap, loops are disjoint
   */
  if (! _CheckNestingBoxes(A, B, result))
      return ERROR;
  else if (*result == AB_DISJOINT)
      return OK;

  /*
   * Check points of A against B loop.  If all are inside,
   * then A is inside B and we're done.  Also - if points
   * of A are both inside and outside B, then we're also done.
   */
  if (!_CheckNesting(A, B, result))
    return ERROR;

  if (*result == A_INSIDE_B || *result == AB_CROSS)
    return OK;
 
  /*
   * All points of A are outside B.  Check to see if A is nested
   * inside B or whether they are disjoint.
   */
  if (! _CheckNesting(B, A, result))
    return ERROR;

  if (*result == A_INSIDE_B)
  {
    *result = B_INSIDE_A;
    return OK;
  }
  else if (*result == AB_CROSS)
  {
    DXSetError(ERROR_INTERNAL, "B x A but A !x B");
    return ERROR;
  }
  else
  {
    *result = AB_DISJOINT;
    return OK;
  }
}

static Error
InsertTreeNode(TreeNode current, TreeNode new)
{
    TreeNode last = NULL, next;

    last = NULL;
    next = current->children;
    while (next)
    {
	int ab;
	
	if (! CheckNesting(new->loop, next->loop, &ab))
	    return ERROR;

	if (ab == A_INSIDE_B)
	{
	    if (new->children)
	    { 
		DXSetError(ERROR_INTERNAL, "topology error");
		return ERROR;
	    }
	    else
		return InsertTreeNode(next, new);
	}
	else if (ab == B_INSIDE_A)
	{
	    TreeNode tmp = next;

	    if (last)
		last->sibling = next->sibling;
	    else
		current->children = next->sibling;

	    next = next->sibling;

	    tmp->sibling = new->children;
	    new->children = tmp;
	}
	else if (ab == AB_CROSS)
	{ 
	      DXSetError(ERROR_INTERNAL, "topology error");
	      return ERROR;
	}
	else
	{
	    last = next;
	    next = next->sibling;
	}
    }

    new->sibling = current->children;
    current->children = new;

    return OK;
}

static Error
AssignTreeLevel(TreeNode node, int level)
{
    TreeNode next;

    if (! node)
	return OK;
      
    node->level = level;
    
    for (next = node->children; next; next = next->sibling)
    {
	if (! AssignTreeLevel(next, level+1))
	    return ERROR;
    }
    
    return OK;
}

static TreeNode
CreateNestingTree(int nLoops, Loop *loops)
{
    int i;
    struct treenode *root;
    struct treenode *nodes;

    nodes = (struct treenode *)DXAllocate((nLoops+1) * sizeof(struct treenode));
    if (! nodes)
	goto error;
    
    root = nodes + 0;

    root->sibling = NULL;
    root->children = NULL;
    root->nChildren = 0;
    root->loop = NULL;

    for (i = 0; i < nLoops; i++)
    {
	/* 
	 * remember that forst node is the root
	 */
	TreeNode node = nodes + (i+1);

	node->loop = loops + i;
	node->sibling = NULL;
	node->children = NULL;
	node->nChildren = 0;

	if (! InsertTreeNode(root, node))
	    goto error;
    }

    return root;
  
error:
    DXFree((Pointer)nodes);
    return NULL;
}

static Error
InitLoopsNested (int *nloopsPtr, Loop **loopsPtr, int np)
{
    int		nloops = *nloopsPtr;
    Loop 	*loops = *loopsPtr;
    int		i, j, k;
    BVertex	*p0, *p1, *p2;
    int		nvert;
    int		maxLoops = *nloopsPtr;
    HashTable   ptHash =
	DXCreateHash(sizeof(BVertex *), HashPoint, ComparePoint);
    if (! ptHash)
	goto error;

    for (j = 0; j < nloops; j++)
    {
	p0 = loops[j].base;
	do
	{
	    BVertex **q;
	    q = (BVertex **)DXQueryHashElement(ptHash, (Key)&p0);
	    if (q)
		p0->id = (*q)->id;
	    else
		if (! DXInsertHashElement(ptHash, (Element)&p0))
		    goto error;
	    p0 = p0->next;
	}
	while (p0 != loops[j].base); 
    }

    DXDestroyHash(ptHash);
    ptHash = NULL;

    for (j = 0; j < nloops; j++)
    {
	p0 = loops[j].base;
	do 
	{
	    for (p1 = loops[j].base; p1 != p0; p1 = p1->next)
		if (p1->id == p0->id)
		{
		    if (p0->next == p1)
		    {
			p0->next = p1->next;
			p0->next->prev = p0;
			loops[j].nvert -= 1;
		    }
		    else
		    {
			int p1loop, p0loop;
			
			for (p1loop = 0, p2 = p1; p2 != p0; p2 = p2->next, p1loop++);
			p0loop = loops[j].nvert - p1loop;

			p2 = p0->prev;

			p0->prev = p1->prev;
			p0->prev->next = p0;

			p1->prev = p2;
			p1->prev->next = p1;

			if (p0loop > 2 && p1loop <= 2)
			{
			    loops[j].nvert = p0loop;
			    loops[j].base  = p0;
			}
			else if (p0loop <= 2 && p1loop > 2)
			{
			    loops[j].nvert = p1loop;
			    loops[j].base  = p1;
			}
			else if (p0loop > 2 && p1loop > 2)
			{
			    loops[j].nvert = p0loop;
			    loops[j].base  = p0;

			    if (nloops == maxLoops)
			    {
				Loop *newloops = (Loop *)
					DXAllocate(2*maxLoops*sizeof(Loop));
				if (! newloops)
				    goto error;
				
				memcpy(newloops, loops, nloops*sizeof(Loop));

				if (loops != *loopsPtr)
				    DXFree((Pointer)loops);
				
				loops = newloops;
			    
				maxLoops *= 2;
			    }
		    
			    loops[nloops] = loops[j];

			    loops[nloops].base = p1;
			    loops[nloops].nvert = p1loop;
			    nloops ++;

			    /*
			     * If its the outer loop, then one may lie inside
			     * the other. If so, its a hole.  Otherwise, 
			     * it should be handled as a single outer loop.
			     * So check.
			     */
			}
			else
			    loops[j].nvert = 0;
		    }

		    p1 = p0 = loops[j].base;
		    break;
		}
	    
	    p0 = p0->next;

	} while (p0 != loops[j].base);
    }

    for (j = 0; j < nloops; j++)
    {
	p0 = loops[j].base;
	if (p0)
	{
	    int nv = loops[j].nvert;

	    if (nv)
	    {
		do 
		{
		    if (nv > 3 && p0->id == p0->next->next->id)
		    {
			nv -= 2;
			loops[j].base = p0;
			p0->next = p0->next->next->next;
			p0->next->prev = p0;
		    }
		    p0 = p0->next;
		} while (p0 != loops[j].base);

		loops[j].nvert = nv;
	    }
	}
    }

    for (i = j = 0; j < nloops; j++)
	if (loops[j].nvert != 0)
	{
	    if (j != 0)
		loops[i] = loops[j];
	    i++;
	}

    nloops = i;

    for (i = 0; i < nloops; i++)
    {
	p0 = loops[i].base;
	do
	{
	    p0->tried = 0;
	    p0->pass = 0;
	    p0 = p0->next;
	} while (p0 != loops[i].base);
    }

    *nloopsPtr = nloops;
    *loopsPtr = loops;

    for (j = 0; j < nloops; j++)
    {
	LoopBBox (loops + j);
	LoopSignNested (loops + j);
    }

    for (j = 0; j < nloops; j++)
	if (loops[j].sign == 0 || loops[j].nvert < 3)
	    loops[j].merged = TRUE;

    return OK;

error:
    if (ptHash)
	DXDestroyHash(ptHash);
    return ERROR;
}


static void FlipLoop (Loop *loops)
{
    Loop	loop;
    int		i;
    BVertex	*curr;
    BVertex	*next;

    loop = *loops;

    for (i = 0, curr = loop.base; i < loop.nvert; i++, curr = next)
    {
	next = curr->next;
	curr->next = curr->prev;
	curr->prev = next;
    }

    loops->sign *= -1;
}


static void
LoopBBox (Loop *lp)
{
    BVertex	*p;
    Point2	pos;
    Point2	ll;		
    Point2	ur;

    if (! lp->base)
    {
        lp->box.ll.u = lp->box.ll.v = 0.0;
        lp->box.ur.u = lp->box.ur.v = 0.0;
	return;
    }

    ur = ll = lp->base->pos;
    for (p = lp->base->next; p != lp->base; p = p->next)
    {
	pos = p->pos;

	if (pos.u < ll.u)      ll.u = pos.u;
	if (pos.u > ur.u) ur.u = pos.u;
	
	if (pos.v < ll.v)      ll.v = pos.v;
	if (pos.v > ur.v) ur.v = pos.v;
    }

    lp->box.ll = ll;
    lp->box.ur = ur;

    return;
}

static void LoopSignNested (Loop *lp)
{
    BVertex	*p;
    BVertex	*p0, *p1, *p2;
    BVertex	*tmp;
    Point2	pos;
    Point2	ll;			/* the lower left */
    int		psign;

    /*
     * Find the "leftmost bottom-most" vertex for use as our starting midpoint.
     * While were at it, find the rest of the bounding box too.
     */

    p = lp->base;
    p1 = lp->base;
    ll = p1->pos;
    for (p = lp->base->next; p != lp->base; p = p->next)
    {
	pos = p->pos;

	if (pos.u < ll.u)
	{
	    ll.u = pos.u;
	    p1   = p;
	}
	else if (pos.u == ll.u && pos.v < ll.v)
	{
	    ll.v = pos.v;
	    p1   = p;
	}
    }

    /*
     * Walk around the polygon until we either get a non-zero sign for the
     * determinant or until we've come all the way around.  If we do get
     * a sign then that is the characteristic direction for the polygon
     * and we can return it immediately.  Otherwise all of the points in
     * the polygon are colinear.
     */

    tmp = p1;
    p0  = p1->prev;
    p2  = p1->next;

    do
    {
	psign = DetSign (p0->pos, p1->pos, p2->pos);
	if (psign)
	{
	    lp->sign = psign;
	    return;
	}

	p0 = p1;
	p1 = p2;
	p2 = p2->next;

    } while (p2 != tmp);

    lp->sign = 0;
}


#define CheckPoint()\
{\
    ip = p->id;\
    if (ip == i0 || ip == i1 || ip == i2)\
	continue;\
    ppos = p->pos;\
    if (ppos.u < tribox.ll.u || ppos.v < tribox.ll.v ||\
        ppos.u > tribox.ur.u || ppos.v > tribox.ur.v)\
	continue;\
    if (PointInTri (ppos, p->prev->pos, p->next->pos, v0, v1, v2, psign))\
	return (TRUE);\
}

static int AnyInsideNested (TreeNode holes, BVertex *p0,
			BVertex *p1, BVertex *p2, int psign)
{
    BVertex	*p;
    int		i0, i1, i2;
    Point2	v0, v1, v2;
    Point2	ppos;
    BBox	tribox, loopbox;
    int		ip;
    Loop	*lp;
    int		i;
    int		nvert;
    TreeNode    hole;

    i0 = p0->id;  i1 = p1->id;  i2 = p2->id;
    v0 = p0->pos; v1 = p1->pos; v2 = p2->pos;

    tribox.ll.u = v0.u<v1.u ? v0.u<v2.u ? v0.u : v2.u : v1.u<v2.u ? v1.u : v2.u;
    tribox.ll.v = v0.v<v1.v ? v0.v<v2.v ? v0.v : v2.v : v1.v<v2.v ? v1.v : v2.v;
    tribox.ur.u = v0.u>v1.u ? v0.u>v2.u ? v0.u : v2.u : v1.u>v2.u ? v1.u : v2.u;
    tribox.ur.v = v0.v>v1.v ? v0.v>v2.v ? v0.v : v2.v : v1.v>v2.v ? v1.v : v2.v;

    /*
     * Check the boundary
     */

    for (p = p2->next; p != p0; p = p->next)
	CheckPoint ();

    /*
     * Check the other loops.
     */

    for (hole = holes; hole; hole = hole->sibling)
    {
	lp = hole->loop;
	loopbox = lp->box;
	if (loopbox.ur.u < tribox.ll.u || loopbox.ll.u > tribox.ur.u ||
	    loopbox.ur.v < tribox.ll.v || loopbox.ll.v > tribox.ur.v)
	    continue;

	for (i = 0, nvert = lp->nvert, p = lp->base;
	     i < nvert;
	     i++, p = p->next)
	    CheckPoint ();
    }

    return (FALSE);
}


#if 0
    foreach point P on the boundary
      foreach unconnected loop L
	foreach point Q on an unconnected loop
	    if P or Q has been split then continue
	    make the edge PQ
		foreach vertex V on the boundary and all loops
		    W = V->next;
		    if (PQ intersects VW)
			continue on to next Q or P
		no intersections so connect PQ
    could not form a link return FALSE

    NOTE: the above has been changed slightly.  It now finds
    the first inside loop visible from some point on the outside
    loop.  Then of all the possible pairs between these two,
    it chooses the closest.

#endif

static int
AddLinkNested(Loop *loops, BVertex *p0, BVertex *spares, int *nvert,
				BVertex **first, BVertex **end)
{
    BVertex	*sp;
    BVertex	*p, *q;
    BVertex	*np, *nq;
    Loop	*lp;
    int		i, maxi;
    int		j, maxj;
    int		k;

    sp = spares;

    for (j = 0, maxj = *nvert, p = p0; j < maxj; j++, p = p->next)
    {
	for (k = 1, lp = loops->next; lp != loops; lp = lp->next, k++)
	{
	    if (lp->merged)
		continue;

	    for (i = 0, maxi = lp->nvert, q = lp->base;
		 i < maxi;
		 i++, q = q->next)
	    {
		if (p->id == q->id)
		    continue;

		if (IntersectsEdge (loops, p, q))
		    continue;

		srprintf ("LINKING at %2d - %2d\n", p->id, q->id);
		np = sp++;
		nq = sp++;
		*np = *p;
		*nq = *q;

		np->pass = 0;
		nq->pass = 0;

		np->next = q;
		nq->next = p;
		np->prev = p->prev;
		nq->prev = q->prev;
		p->prev->next = np;
		q->prev->next = nq;
		p->prev  = nq;
		q->prev  = np;

		lp->merged     = TRUE;
		lp->next->prev = lp->prev;
		lp->prev->next = lp->next;
		lp->next       = NULL;
		lp->prev       = NULL;

		*nvert += lp->nvert + 2;
		loops->nvert = *nvert;
		loops->base  = p0;

		np->prev->tried = 0;
		np->tried = 0;

		*first = np->prev;
		*end   = nq->next;

		return (TRUE);
	    
	    }
	}
    }

    return FALSE;
}

/*
 * If the edge we are looking at contains one of the link points then
 * just move on.  If we find that a vertex has been replicated then we
 * assume an intersection since if we don't we run the risk of introducing
 * a loop that runs counter to the real direction of the boundary.
 * Next we check the rotations of the triangles formed by connecting the
 * endpoints of the segments.  Basically, if we find that both triangles
 * formed using one segment as the base rotate the same way then there
 * isn't an intersection.
 */

#define CheckEdge(_V)\
{\
    v1 = v0->next;\
    srprintf ("   %2d - %2d:  ", v0->id, v1->id);\
    if (v0 == _V || v1 == _V)\
    {\
	srprintf ("point match\n");\
	continue;\
    }\
    v0i  = v0->id;\
    v1i  = v1->id;\
    if (pid == v0i || pid == v1i || qid == v0i || qid == v1i)\
    {\
	srprintf ("ID MATCH\n");\
	return (TRUE);\
    }\
    v0os = v0->pos;\
    v1os = v1->pos;\
    if ((v0os.u < pqbox.ll.u && v1os.u < pqbox.ll.u) ||\
        (v0os.v < pqbox.ll.v && v1os.v < pqbox.ll.v) ||\
        (v0os.u > pqbox.ur.u && v1os.u > pqbox.ur.u) ||\
        (v0os.v > pqbox.ur.v && v1os.v > pqbox.ur.v))\
	{\
	    srprintf ("no bbox overlap\n");\
	    continue;\
	}\
    abc = DetSign (ppos, qpos, v0os);\
    abd = DetSign (ppos, qpos, v1os);\
    cda = DetSign (v0os, v1os, ppos);\
    cdb = DetSign (v0os, v1os, qpos);\
    srprintf ("%2d %2d %2d %2d  ", abc, abd, cda, cdb);\
    if ((abc * abd) == 1 || (cda * cdb) == 1)\
    {\
	srprintf ("no intersection\n");\
	continue;\
    }\
    srprintf ("INTERSECTION\n");\
    return (TRUE);\
}

static int IntersectsEdge (Loop *loops, BVertex *p, BVertex *q)
{
    int		pid, qid;
    Point2	ppos, qpos;
    BVertex	*v0, *v1;
    int		v0i, v1i;
    Point2	v0os, v1os;
    Loop	*lp;
    int		i;
    int		nvert;
    int		abc, abd, cda, cdb;
    BBox	pqbox;
    BBox	loopbox;

    pid  = p->id;
    qid  = q->id;
    ppos = p->pos;
    qpos = q->pos;

    if (ppos.u < qpos.u)
    {
	pqbox.ll.u = ppos.u;
	pqbox.ur.u = qpos.u;
    }
    else
    {
	pqbox.ll.u = qpos.u;
	pqbox.ur.u = ppos.u;
    }

    if (ppos.v < qpos.v)
    {
	pqbox.ll.v = ppos.v;
	pqbox.ur.v = qpos.v;
    }
    else
    {
	pqbox.ll.v = qpos.v;
	pqbox.ur.v = ppos.v;
    }

    /*
     * Check around the boundary
     */

    for (v0 = p->next; v0 != p; v0 = v0->next)
	CheckEdge (p);

    /*
     * Check around each loop
     */

    for (lp = loops->next; lp != loops; lp = lp->next)
    {
	loopbox = lp->box;
	if (loopbox.ur.u < pqbox.ll.u || loopbox.ll.u > pqbox.ur.u ||
	    loopbox.ur.v < pqbox.ll.v || loopbox.ll.v > pqbox.ur.v)
	    continue;

	for (i = 0, nvert = lp->nvert, v0 = lp->base;
	     i < nvert;
	     i++, v0 = v0->next)
	{
	    CheckEdge (q);
	}
    }

    return (FALSE);
}


/*
 * Determines whether a point is _dxf_inside a triangle by examining the
 * rotations in the triangles formed by each pair of vertices in the
 * test triangle and the point.  If all of the rotations are the same
 * then the point is _dxf_inside of the triangle, if one is different then
 * the point is outside, if one doesn't rotate then the point is either
 * on or coincident with one of the triangle's edges.
 */

static int PointInTri (Point2 p, Point2 plast, Point2 pnext,
				Point2 v0, Point2 v1, Point2 v2, int psign)
{
    int		det01p, det12p, det20p;
    int		sum;
    int		np = 0, nn = 0, nz = 0;
    Point2	e0, e1;

    det01p = DetSign (v0, v1, p);
    if (det01p > 0)
        np ++;
    else if (det01p == 0)
    {
	e0 = v0;
	e1 = v1;
        nz ++;
    }
    else 
	nn ++;

    det12p = DetSign (v1, v2, p);
    if (det12p > 0)
        np ++;
    else if (det12p == 0)
    {
	e0 = v1;
	e1 = v2;
        nz ++;
    }
    else 
	nn ++;

    det20p = DetSign (v2, v0, p);
    if (det20p > 0)
        np ++;
    else if (det20p == 0)
    {
	e0 = v2;
	e1 = v0;
        nz ++;
    }
    else 
	nn ++;

    if (nn == 3 || np == 3)
	return TRUE;
    else if ((nn == 0 || np == 0) && nz == 1)
    {
	int p0, p1;

	p0 = DetSign(e0, e1, plast);
	p1 = DetSign(e0, e1, pnext);
	if (p0 == psign || p1 == psign)
	    return TRUE;
	else
	    return FALSE;
    }
    else
	return FALSE;
}


#if TRIANGULATE_MAIN
main ()
{
    Point	pvals [100];
    int		ids[1000];
    Loop	loops[100];
    Triangle	tris[1000];
    int		ntri;
    int		nloops	= 0;
    int		npts	= 0;
    int		i, j, k;
    Point	*pptr;
    Point	p;
    Point	normal;

    p.x = p.y = p.z = 0.0;
    for (i = 0, pptr = pvals; i < 10; i++)
    {
	p.x = (float) i;
	for (j = 0; j < 10; j++)
	{
	    p.y = (float) j;
	    *pptr++ = p;
	}
    }

    normal.x = 0.0;
    normal.y = 0.0;
    normal.z = 1.0;

    for (;;)
    {
	printf ("enter number of loops:  ");
	fflush (stdout);
	scanf ("%d", &nloops);
	if (nloops == 0)
	    exit (0);
	if (nloops < 0)
	{
	    nloops = -nloops;
	    srprint = TRUE;
	}
	for (j = 0; j < nloops; j++)
	{
	    loops[j].vids = ids + npts;
	    printf ("  loop %2d: enter cnt & indices:  ", j);
	    fflush (stdout);
	    scanf ("%d", &loops[j].nvert);
	    for (i = 0; i < loops[j].nvert; i++)
		scanf ("%d", &ids[npts++]);
	}
	printf ("\n");

	TriangulateLoops (nloops, loops, pvals, &normal, 100, tris, &ntri);
	printf ("  %3d triangles in buffer %x\n", ntri, tris);
	for (i = 0; i < ntri; i++)
	    printf ("%3d:  %3d %3d %3d\n", i, tris[i].p, tris[i].q, tris[i].r);
    }
}
#endif


static Error TriangulateLoops (int nloops, Loop *loops, float *vertices,
			       Point *normal, int np,
			       Triangle *tptr, int *rntri, int nDim)
{
    Error	ret	= ERROR;	/* return status                    */
    int		nvert;			/* total vertex count               */
    BVertex	*vbase	= NULL;		/* base of list blocks              */
    BVertex	*spares;		/* spare nodes for splitting        */
    int		size;			/* size temporary                   */
    Loop	loop;			/* local copy of current loop       */
    int		psign;			/* sign of the polygon's  rotation  */
    int		tsign;			/* sign of the triangle's rotation  */
    Point	lnormal;		/* local copy of face normal	    */
    int		elim;			/* dimension to elminate            */
    BVertex	*p0, *p1, *p2;		/* vertices of tri being considered */
    int		i0, i1, i2;
    int		pass;
    int		ntri	= 0;		/* # of triangles generated         */
    int		i, j;
    BVertex	*bptr;
    Loop	*orig_loops = loops;
    BVertex	*first, *end;
    int		not_done;
#if 0
    int		loopknt = 0;
#endif

    if (rntri)
	*rntri = 0;

    /*
     * Since each loop requires a link and each link requires splitting 2
     * vertices we end up with the following for the maximum number of
     * vertices.
     */

    for (nvert = 0, i = 0; i < nloops; i++)
	nvert += loops[i].nvert;
    nvert += 2 * (nloops - 1);

    size = nvert * sizeof (BVertex);
    vbase = (BVertex *) DXAllocateLocal (size);
    if (vbase == NULL)
	goto error;
    spares = vbase;
    memset (vbase, NULL, size);

    /*
     * For each loop:
     *   * Initialize its boundary vertex list as a circularly linked list.
     *   * We first have to chunk off the necessary number of BVertex nodes
     *     for use with this list and set the rest aside as the spares.
     *   * For now we'll just put in the vertex id.  We'll fill in the 
     *     vertex positions later.  We separate this for two reasons.  First,
     *     it involves random accessing and we don't want to guarantee that
     *     we destroy any savings we get from the PVS's line cache, and 
     *     second we want to use tight loops for the projection from 3-space
     *     to 2-space.
     */

    for (j = 0; j < nloops; j++)
    {
	loops[j].base = spares;
	spares += loops[j].nvert;
	loop = loops[j];

	for (i = 0, bptr = loop.base; i < loop.nvert; i++, bptr++)
	{
	    bptr->next   = bptr + 1;
	    bptr->prev   = bptr - 1;
	    bptr->lp     = loops + j;
	    bptr->id     = *loop.vids++;
	}

	/*
	 * Circularly link the list
	 */

	bptr = loop.base + (loop.nvert - 1);
	bptr->next = loop.base;
	loop.base->prev = bptr;
    }

    ProjectPoints (normal, vertices, nloops, loops, nDim);

    /*
     * Get rid of colinear loops and loops with fewer than 3 points and
     * link all of the others together.
     */

    if (! InitLoops (&nloops, &loops, spares, np))
	goto error;

    psign = loops->sign;
    if (psign == 0)
    {
	ret = OK;
	goto error;
    }


    first = loops->base;
    end   = NULL;
    pass  = 0;
    nvert = loops->nvert;
    p0    = loops->base;
    first = NULL;
    while (nloops > 0)
    {
	do
	{
	    if (first != NULL)
		p0 = first;

	    not_done = 0;
	    pass++;
	    
	    if (first == end)
		end = NULL;

	    while (p0->pass < pass && p0 != end)
	    {
		p0->pass = pass;

		p1 = p0->next;
		p2 = p1->next;

		if (!p0->tried)
		{
		    srprintf ("DXTri: %2d %2d %2d (%3d/%3d)",
			   p0->id, p1->id, p2->id, pass, nvert);

		    if (p0->id == p1->id || p1->id == p2->id ||
			(tsign = DetSign(p0->pos, p1->pos, p2->pos)) == 0)
		    {
			/*
			 * degenerate triangle.  eliminate p1 from the loop.
			 */
			p0->next = p2;
			p2->prev = p0;
			not_done++;
			nvert--;
		    }
		    else if (tsign != psign ||
			AnyInside (loops, p0, p1, p2, psign))
		    {
			p0->tried = 1;
		    }
		    else
		    {
			not_done++;
			p0->prev->tried = 0;
			if (p0 == first)
			    first = first->prev;
			if (p2 == end) 
			{
			    if (end == first)
				end = NULL;
			    else
				end = end->next;
			}

#if 0
			if (pass >= which_pass)
			{
			    fprintf(stderr, "output: %d %d %d\n", p0->id, p1->id, p2->id);
			    ShowLoop(p2);
			}
#endif

			EmitTri012;
			nvert--;

			while (nvert > 2)
			{
			    p1 = p0->next;
			    p2 = p1->next;

			    if (DetSign(p0->pos, p1->pos, p2->pos) == 0)
			    {
				p0->next = p2;
				p2->prev = p0;
				nvert --;
			    }
			    else
				break;
			}
		    }
		}

		p0 = p0->next;
	    }

	} while (not_done && nvert >= 3);

	if (nloops > 1)
	    if (! AddLink (loops, p0, &spares, &nvert, &first, &end))
	    {
		DXSetError (ERROR_INTERNAL, "Couldn't add link");
		goto error;
	    }

	nloops--;
    }

    if (nvert > 3)
    {
	DXSetError(ERROR_INVALID_DATA, "topology error");
#if 0
        {
	FILE *foo = fopen("bad.dx", "w");
	if ( foo )
	{
	    if (! first)
		first = p0;

	    for (i = 0, p0 = first; p0->next != first; p0 = p0->next)
		i++;

	    fprintf(foo,
	"object \"faces\" class array type int rank 0 items 1 data follows\n0\nattribute \"ref\" string \"loops\"\n");

	    fprintf(foo,
	"object \"loops\" class array type int rank 0 items 1 data follows\n0\nattribute \"ref\" string \"edges\"\n");

	    fprintf(foo,
	"object \"edges\" class array type int rank 0 items %d data follows\n", nvert);
	    for (i = 0; i < nvert; i++)
		fprintf(foo, "%d\n", i);
	    fprintf(foo,
	"attribute \"ref\" string \"positions\"\n ", nvert);

	    fprintf(foo,
	"object \"positions\" class array type float rank 1 shape 3 items %d data follows\n", nvert);
	    for (i = 0, p0 = first; i < nvert; i++, p0 = p0->next)
		fprintf(foo, "%f %f 0\n", p0->pos.u, p0->pos.v);
	    
	    fprintf(foo, "object \"fle\" class field\n");
	    fprintf(foo, "    component \"faces\" \"faces\"\n");
	    fprintf(foo, "    component \"loops\" \"loops\"\n");
	    fprintf(foo, "    component \"edges\" \"edges\"\n");
	    fprintf(foo, "    component \"positions\" \"positions\"\n");
	    fprintf(foo, "end\n");

	    fclose(foo);
	}
	}
#endif
	goto error;
    }

    if (nvert == 3)
	EmitTri012;

#if 0
    if (dumpLoop || loopknt == loopLimit) 
	DumpLoop(p0);
#endif

    ret = OK;
    goto cleanup;

error:
#if ! TRIANGULATE_ALWAYS
    ntri = 0;
#endif

cleanup:
    if (orig_loops != loops)
       DXFree((Pointer)loops);
    DXFree ((Pointer) vbase);
    if (rntri)
	*rntri = ntri;
    return (ret);
}

static int
AnyInside (Loop *loops, BVertex *p0,
			BVertex *p1, BVertex *p2, int psign)
{
    BVertex	*p;
    int		i0, i1, i2;
    Point2	v0, v1, v2;
    Point2	ppos;
    BBox	tribox, loopbox;
    int		ip;
    Loop	*lp;
    int		i;
    int		nvert;

    i0 = p0->id;  i1 = p1->id;  i2 = p2->id;
    v0 = p0->pos; v1 = p1->pos; v2 = p2->pos;

    tribox.ll.u = v0.u<v1.u?v0.u<v2.u?v0.u : v2.u : v1.u<v2.u?v1.u : v2.u;
    tribox.ll.v = v0.v<v1.v?v0.v<v2.v?v0.v : v2.v : v1.v<v2.v?v1.v : v2.v;
    tribox.ur.u = v0.u>v1.u?v0.u>v2.u?v0.u : v2.u : v1.u>v2.u?v1.u : v2.u;
    tribox.ur.v = v0.v>v1.v?v0.v>v2.v?v0.v : v2.v : v1.v>v2.v?v1.v : v2.v;

    /*
     * Check the boundary
     */

    for (p = p2->next; p != p0; p = p->next)
	CheckPoint ();

    /*
     * Check the other loops.
     */

    for (lp = loops->next; lp != loops; lp = lp->next)
    {
	loopbox = lp->box;
	if (loopbox.ur.u < tribox.ll.u || loopbox.ll.u > tribox.ur.u ||
	    loopbox.ur.v < tribox.ll.v || loopbox.ll.v > tribox.ur.v)
	    continue;

	for (i = 0, nvert = lp->nvert, p = lp->base;
	     i < nvert;
	     i++, p = p->next)
	    CheckPoint ();
    }

    return (FALSE);
}

static Error
InitLoops (int *nloopsPtr, Loop **loopsPtr, BVertex *spares,
		       int np)
{
    int		nloops = *nloopsPtr;
    Loop 	*loops = *loopsPtr;
    int		i, j, k;
    int		flip;
    BVertex	*p0, *p1, *p2;
    int		nvert;
    int		maxLoops = *nloopsPtr;
    HashTable   ptHash =
	DXCreateHash(sizeof(BVertex *), HashPoint, ComparePoint);
    if (! ptHash)
	goto error;

    for (j = 0; j < nloops; j++)
    {
	p0 = loops[j].base;
	do
	{
	    BVertex **q;
	    q = (BVertex **)DXQueryHashElement(ptHash, (Key)&p0);
	    if (q)
		p0->id = (*q)->id;
	    else
		if (! DXInsertHashElement(ptHash, (Element)&p0))
		    goto error;
	    p0 = p0->next;
	}
	while (p0 != loops[j].base); 
    }

    DXDestroyHash(ptHash);
    ptHash = NULL;

    for (j = 0; j < nloops; j++)
    {
	p0 = loops[j].base;
	do 
	{
	    for (p1 = loops[j].base; p1 != p0; p1 = p1->next)
		if (p1->id == p0->id)
		{
		    if (p0->next == p1)
		    {
			p0->next = p1->next;
			p0->next->prev = p0;
			loops[j].nvert -= 1;
		    }
		    else
		    {
			int p1loop, p0loop;
			
			for (p1loop = 0, p2 = p1; p2 != p0; p2 = p2->next, p1loop++);
			p0loop = loops[j].nvert - p1loop;

			p2 = p0->prev;

			p0->prev = p1->prev;
			p0->prev->next = p0;

			p1->prev = p2;
			p1->prev->next = p1;

			if (p0loop > 2 && p1loop <= 2)
			{
			    loops[j].nvert = p0loop;
			    loops[j].base  = p0;
			}
			else if (p0loop <= 2 && p1loop > 2)
			{
			    loops[j].nvert = p1loop;
			    loops[j].base  = p1;
			}
			else if (p0loop > 2 && p1loop > 2)
			{
			    loops[j].nvert = p0loop;
			    loops[j].base  = p0;

			    if (nloops == maxLoops)
			    {
				Loop *newloops = (Loop *)
					DXAllocate(2*maxLoops*sizeof(Loop));
				if (! newloops)
				    goto error;
				
				memcpy(newloops, loops, nloops*sizeof(Loop));

				if (loops != *loopsPtr)
				    DXFree((Pointer)loops);
				
				loops = newloops;
			    
				maxLoops *= 2;
			    }
		    
			    loops[nloops] = loops[j];

			    loops[nloops].base = p1;
			    loops[nloops].nvert = p1loop;
			    nloops ++;
			}
			else
			    loops[j].nvert = 0;
		    }

		    p1 = p0 = loops[j].base;
		    break;
		}
	    
	    p0 = p0->next;

	} while (p0 != loops[j].base);
    }

    for (j = 0; j < nloops; j++)
	loops[j].merged = FALSE;

    for (j = 0; j < nloops; j++)
	LoopSign (loops + j);

    loops->next = loops->prev = loops;
    if (nloops == 1)
	return OK;

    /*
     * If one of the inner loops (a hole) rotates the same way as the
     * outer loop then we need to flip its sense.  We do this by
     * interchanging its next and its prev pointers.  But first, we
     * want to make sure that all of our outer loops always rotate
     * the same way, this way thin faces will have their tops & bottoms
     * (or fronts & backs) triangulated the same way rather than oppositely.
     * This is particularly important for reducing the triangle count when
     * a thin solid collapses to a face during simplification.
     */

    if (loops->sign == 1)
	FlipLoop (loops);
    flip = loops->sign;
    for (j = 1; j < nloops; j++)
	if (loops[j].sign == flip)
	    FlipLoop (loops + j);

    for (j = 0; j < nloops; j++)
	if (loops[j].sign == 0 || loops[j].nvert < 3)
	    loops[j].merged = TRUE;

    for (j = 0; j < nloops; j++)
    {
	p0 = loops[j].base;
	if (p0)
	{
	    do 
	    {
		for (k = j+1; k < nloops; k++)
		    if ((p1 = loops[k].base) != NULL)
			do
			{
			    if (p1->id == p0->id)
			    {
				p2 = p0->next;
				p0->next = p1->next;
				p0->next->prev = p0;
				p1->next = p2;
				p2->prev = p1;
				loops[j].nvert += loops[k].nvert;
				loops[k].nvert = 0;
				loops[k].base = NULL;
				if (loops[j].box.ll.u > loops[k].box.ll.u)
				    loops[j].box.ll.u = loops[k].box.ll.u;
				if (loops[j].box.ll.v > loops[k].box.ll.v)
				    loops[j].box.ll.v = loops[k].box.ll.v;
				if (loops[j].box.ur.u < loops[k].box.ur.u)
				    loops[j].box.ur.u = loops[k].box.ur.u;
				if (loops[j].box.ur.v < loops[k].box.ur.v)
				    loops[j].box.ur.v = loops[k].box.ur.v;
				p1 = NULL;
			    }
			    else
				p1 = p1->next;
			} while (p1 != loops[k].base);
		p0 = p0->next;
	    } while (p0 != loops[j].base);
	}

	p0 = loops[j].base;
	if (p0)
	{
	    int nv = loops[j].nvert;

	    if (nv)
	    {
		do 
		{
		    if (nv > 3 && p0->id == p0->next->next->id)
		    {
			nv -= 2;
			loops[j].base = p0;
			p0->next = p0->next->next->next;
			p0->next->prev = p0;
		    }
		    p0 = p0->next;
		} while (p0 != loops[j].base);

		loops[j].nvert = nv;
	    }
	}
    }

    for (i = j = 0; j < nloops; j++)
	if (loops[j].nvert != 0)
	{
	    if (j != 0)
		loops[i] = loops[j];
	    i++;
	}

    nloops = i;

    for (i = 0; i < nloops; i++)
    {
	p0 = loops[i].base;
	do
	{
	    p0->tried = 0;
	    p0->pass = 0;
	    p0 = p0->next;
	} while (p0 != loops[i].base);
    }

    /*
     * Set up the loop linked lists
     */

    for (j = 0; j < nloops; j++)
    {
	loops[j].next = loops + j + 1;
	loops[j].prev = loops + j - 1;
    }
    loops[nloops - 1].next = loops;
    loops[0].prev = loops + nloops - 1;

    for (j = 0; j < nloops; j++)
    {
	if (loops[j].merged)
	{
	    loops[j].prev->next = loops[j].next;
	    loops[j].next->prev = loops[j].prev;
	}
    }

    *nloopsPtr = nloops;
    *loopsPtr = loops;
    return OK;

error:
    if (ptHash)
	DXDestroyHash(ptHash);
    return ERROR;
}

static int
AddLink(Loop *loops, BVertex *p0, BVertex **spares, int *nvert,
				BVertex **first, BVertex **end)
{
    BVertex	*sp;
    BVertex	*p, *q;
    BVertex	*np, *nq;
    Loop	*lp;
    int		i, maxi;
    int		j, maxj;
    int         k;

    sp = *spares;

    for (j = 0, maxj = *nvert, p = p0; j < maxj; j++, p = p->next)
    {
	for (k = 1, lp = loops->next; lp != loops; lp = lp->next, k++)
	{
	    for (i = 0, maxi = lp->nvert, q = lp->base;
		 i < maxi;
		 i++, q = q->next)
	    {
		if (p->id == q->id)
		    continue;

		if (IntersectsEdge (loops, p, q))
		    continue;

		srprintf ("LINKING at %2d - %2d\n", p->id, q->id);
		np = sp++;
		nq = sp++;
		*np = *p;
		*nq = *q;

		np->pass = 0;
		nq->pass = 0;

		np->next = q;
		nq->next = p;
		np->prev = p->prev;
		nq->prev = q->prev;
		p->prev->next = np;
		q->prev->next = nq;
		p->prev  = nq;
		q->prev  = np;

		lp->merged     = TRUE;
		lp->next->prev = lp->prev;
		lp->prev->next = lp->next;
		lp->next       = NULL;
		lp->prev       = NULL;

		*nvert += lp->nvert + 2;
		*spares = sp;
		loops->nvert = *nvert;
		loops->base  = p0;

		np->prev->tried = 0;
		np->tried = 0;

		if (srprint)
		{
		    printf("sequence after adding nested loop:");
		    for (j = 0, maxj = *nvert, p = p0; j < maxj; j++, p = p->next)
			 printf("%d\n", p->id);
		}

		*first = np->prev;
		*end   = nq->next;

		return (TRUE);
	    
	    }
	}
    }

    return FALSE;
}

static void LoopSign (Loop *lp)
{
    BVertex	*p;
    BVertex	*p0, *p1, *p2;
    BVertex	*tmp;
    Point2	pos;
    Point2	ll;			/* the lower left */
    Point2	ur;
    int		psign;

    /*
     * Find the "leftmost bottom-most" vertex for use as our starting midpoint.
     * While were at it, find the rest of the bounding box too.
     */

    p = lp->base;

    for (p1 = p2 = p, ll = ur = p->pos, tmp = p->next;
	 tmp != p;
	 tmp = tmp->next)
    {
	pos = tmp->pos;

	if (pos.u < ll.u)
	{
	    ll.u = pos.u;
	    p1   = tmp;
	}
	else if (pos.u == ll.u)
	{
	    if (pos.v < p1->pos.v)
		p1 = tmp;
	}
	else if (pos.u > ur.u)
	    ur.u = pos.u;
	
	if (pos.v < ll.v)
	    ll.v = pos.v;
	else if (pos.v > ur.v)
	    ur.v = pos.v;
    }

    /*
     * Walk around the polygon until we either get a non-zero sign for the
     * determinant or until we've come all the way around.  If we do get
     * a sign then that is the characteristic direction for the polygon
     * and we can return it immediately.  Otherwise all of the points in
     * the polygon are colinear.
     */

    tmp = p1;
    p0  = p1->prev;
    p2  = p1->next;

    do
    {
	psign = DetSign (p0->pos, p1->pos, p2->pos);
	if (psign)
	{
	    lp->sign     = psign;
	    lp->box.ll   = ll;
	    lp->box.ur   = ur;
	    return;
	}

	p0 = p1;
	p1 = p2;
	p2 = p2->next;
    } while (p2 != tmp);

    lp->sign = 0;
    ll.u = ll.v = 0.0;
    lp->box.ll = ll;
    lp->box.ur = ll;
}