/*****************************************************************************
*   Routines to	prepare objects according to view file matrix:		     *
*									     *
* Written by:  Gershon Elber				Ver 1.0, Jan. 1989   *
*****************************************************************************/

#include <math.h>
#include <stdio.h>
#include "program.h"
#include "genmat.h"

/* #define DEBUG1			      Print edge/hash table content. */

#define SAME_VERTEX(V1, V2)	(APX_EQ(V1 -> Coord[0], V2 -> Coord[0]) && \
				 APX_EQ(V1 -> Coord[1], V2 -> Coord[1]) && \
				 APX_EQ(V1 -> Coord[2], V2 -> Coord[2]))

#ifdef DEBUG1
static void PrintEdgeContent(EdgeStruct *PEdge);
static void DrawEdgeHashTable(void);
static void DrawEdge(EdgeStruct *PEdge);
#endif /* DEBUG1 */

static int MinYLevel, MaxYLevel, CrntYLevel, PrintYLevel;

static void PrepareAllObjects(IPObjectStruct *PObjects);
static void PrepareOneObject(IPObjectStruct *PObject);
static int PrepareOnePolygon(IPPolygonStruct *PPolygon, int IsPolygon);
static void UpdateBBoxPolygon(IPPolygonStruct *PPolygon);
static int UpdateEqnPolygon(IPPolygonStruct *PPolygon, int *ReversePoly);
static IPVertexStruct *ReverseLinList(IPVertexStruct *VList);
static void GenEdgesFromPoly(IPPolygonStruct *PPolygon, int IsPolygon);
static void InsertEdgeToHashTbl1(EdgeStruct *PEdge);
static void IntersectAllEdges(void);
static void InsertEdgeToHashTbl2(EdgeStruct *PEdge);
static int IntersectEdgeList(EdgeStruct *PEdge, EdgeStruct *PEList,
								int TestYMin);
static int IntersectEdgeEdge(EdgeStruct *PEdge1, EdgeStruct *PEdge2,
			EdgeStruct **PEdgeNew1, EdgeStruct **PEdgeNew2);
static void PrintPolyContent(IPPolygonStruct *PPoly);

/*****************************************************************************
* Routine to prepare to draw NumOfObjects given in Objects from		     *
* FileDescription FD according to view matrix Mat. If NumOfObjects == 0	then *
* all the objects defined by the data sturcture	are drawn.		     *
* If GlblNumEdge != 0 then only GlblNumEdge first edges of each polygons are *
* tested for visibility	(usefull in case in input polygons has known	     *
* repitition edges sequence which is redundent).			     *
*****************************************************************************/
void PrepareViewData(IPObjectStruct *PObjects)
{
    int	i;

    IritCPUTime(TRUE);

    for	(i = 0; i < EDGE_HASH_TABLE_SIZE; i++)
	EdgeHashTable[i] = NULL;
    for	(i = 0; i < POLY_HASH_TABLE_SIZE; i++)
	PolyHashTable[i] = NULL;

    if (!GlblQuiet)
	fprintf(stderr, "\nPass 2, Edges =      ");

    PrepareAllObjects(PObjects);

    if (!GlblQuiet)
	fprintf(stderr, ",  %6.2lf seconds.", IritCPUTime(FALSE));

    IntersectAllEdges();       /* Break edges to visibily uniform sub-edges. */
}

/*****************************************************************************
* Scan all objects.							     *
*****************************************************************************/
static void PrepareAllObjects(IPObjectStruct *PObjects)
{
    while (PObjects) {
	if (IP_IS_POLY_OBJ(PObjects))
	    PrepareOneObject(PObjects);
	PObjects = PObjects -> Pnext;
    }
}

/*****************************************************************************
* Routine to prepare one object PObject.				     *
*****************************************************************************/
static void PrepareOneObject(IPObjectStruct *PObject)
{
    int Level;
    IPPolygonStruct *Ptemp,
	*PList = PObject -> U.Pl;

    while (PList) {
	Ptemp = PList -> Pnext;

	if (!IritPrsrIsConvexPolygon(PList)) {
	    static int Printed = FALSE;

	    if (!Printed) {
		fprintf(stderr,
			"\nWARNING: Non convex polygon(s) might be in data (see CONVEX in IRIT),\n\t\t\t\toutput can be wrong as the result!\n     ");
		Printed = TRUE;
	    }
	}

	if (PrepareOnePolygon(PList, IP_IS_POLYGON_OBJ(PObject))) {
	    /* And add polygon into polygon hash table sorted by Ymin: */
	    Level = (int) ((PList -> BBox[0][1] + 1.0) * POLY_HASH_TABLE_SIZE2);
	    Level = BOUND(Level, 0, POLY_HASH_TABLE_SIZE1); /* Be 100% safe. */
	    PList -> Pnext = PolyHashTable[Level];
	    PolyHashTable[Level] = PList;   /* Concat it to poly hash table. */
	}

	PList =	Ptemp;
    }
}

/*****************************************************************************
* Routine to prepare one polygon PPolygon.				     *
* Returns TRUE iff this object is a valid POLYGON (not a POLYLINE!).	     *
*****************************************************************************/
static int PrepareOnePolygon(IPPolygonStruct *PPolygon, int IsPolygon)
{
    int	i, Reverse;
    RealType CpCoord[3];
    IPVertexStruct
	*PList = PPolygon -> PVertex;

    /* Make a circluar polygon list inlt a linear list by breaking the end.  */
    while (PList != NULL && PList -> Pnext != PPolygon -> PVertex)
	PList = PList -> Pnext;
    if (PList)
	PList -> Pnext = NULL;
    PList = PPolygon -> PVertex;

    while (PList) {
	/* Convert the coordinate to screen space (in RealType pres.). */
	MatMultVecby4by4(CpCoord, PList -> Coord, GlblViewMat);
	for (i = 0; i < 3; i++)
	    PList -> Coord[i] = CpCoord[i];

	PList =	PList -> Pnext;
    }
    if (IsPolygon) {
	/* Find plane equation of poly, and let know if need to reverse. */
	if (!UpdateEqnPolygon(PPolygon, &Reverse))
	    return FALSE;
	UpdateBBoxPolygon(PPolygon);/* Find X, Y extrem in screen space. */
	GenEdgesFromPoly(PPolygon, TRUE);	  /* Generate all its edges. */
	if (Reverse)
    	    PPolygon -> PVertex = ReverseLinList(PPolygon -> PVertex);
	return TRUE;
    }
    else {
	GenEdgesFromPoly(PPolygon, FALSE);	  /* Generate all its edges. */
	return FALSE;
    }
}

/*****************************************************************************
* Routine to update polygon boundary box in screen space:		     *
* Note this routine is called after the	polygons was checked for validity -  *
* all the list of objects was found to be vertices only.		     *
*****************************************************************************/
static void UpdateBBoxPolygon(IPPolygonStruct *PPolygon)
{
    RealType *Coord, Xmin, Xmax, Ymin, Ymax; /* Bounding box of the polygon. */
    IPVertexStruct
	*PList = PPolygon -> PVertex;

    Xmin = Xmax	= PList -> Coord[0];
    Ymin = Ymax	= PList	-> Coord[1];
    PList = PList -> Pnext;
    while (PList) {
	Coord = PList -> Coord;
	if (Coord[0] > Xmax)
	    Xmax = Coord[0];
	if (Coord[0] < Xmin)
	    Xmin = Coord[0];
	if (Coord[1] > Ymax)
	    Ymax = Coord[1];
	if (Coord[1] < Ymin)
	    Ymin = Coord[1];
	PList =	PList -> Pnext;
    }
    PPolygon ->	BBox[0][0] = Xmin;
    PPolygon -> BBox[1][0] = Xmax;
    PPolygon ->	BBox[0][1] = Ymin;
    PPolygon -> BBox[1][1] = Ymax;
}

/*****************************************************************************
* Routine to update plane equation of the given	polygon:		     *
*   It is assumed that at list 3 points in polygon do exists, and pick the   *
* tuple that has biggest length for maximum accuracy.			     *
*   Note we IGNORE PLANE if was in data file.				     *
*   In addition a test is made if all polygon vertices are ordered such that *
* the cross product of each 3 consecutive vertices (projected to Z=0 plane)  *
* is allways positive. Note the	polygon	must be	convex,	so result might	be   *
* all positive or all negative.	In the later case the order is reversed	     *
*****************************************************************************/
static int UpdateEqnPolygon(IPPolygonStruct *PPolygon, int *ReversePoly)
{
    static int
	PolygonCount = 0;
    int	i;
    RealType V1[3], V2[3], *Coord, *CoordNext, *CoordNextNext,
	Plane[3], Len, MaxPlane[3],
	MaxLen = 0.0;
    IPVertexStruct
	*PList = PPolygon -> PVertex;

    PolygonCount++;

    *ReversePoly = FALSE;

    do {	/* Search for 3 consequtive non-colinear point from polygon: */
	Coord = PList -> Coord;
	CoordNext = PList -> Pnext -> Coord;
	CoordNextNext = PList -> Pnext -> Pnext -> Coord;
	for (i = 0; i < 3; i++) {   /* Prepare two vectors on polygon plane. */
	    V1[i] = Coord[i] - CoordNext[i];
	    V2[i] = CoordNext[i] - CoordNextNext[i];
	}

	/* Find plane normal by a cross product of the two vectors on plane: */
	Plane[0] = V1[1] * V2[2] - V1[2] * V2[1];
	Plane[1] = V1[2] * V2[0] - V1[0] * V2[2];
	Plane[2] = V1[0] * V2[1] - V1[1] * V2[0];

	/* Find vector Len. - we are looking for the biggest: */
	Len = sqrt(SQR(Plane[0]) + SQR(Plane[1]) + SQR(Plane[2]));
	if (Len > MaxLen) {
	    for (i = 0; i < 3; i++)
		MaxPlane[i] = Plane[i];
	    MaxLen = Len;
	}
	PList = PList -> Pnext;				  /* Try next tuple. */
    } while (PList -> Pnext -> Pnext != NULL);

    if (ABS(MaxLen) < SQR(EPSILON)) { /* Fail to find 3 non-colinear points. */
	if (GlblMore) {
	    fprintf(stderr,
	        "\nError: Invalid polygon (%d) found in file (zero edge length/colinear vertices):\n",
		PolygonCount);
	    PrintPolyContent(PPolygon);
	}
	return FALSE;
    }

    for (i = 0; i < 3; i++)
	PPolygon -> Plane[i] = MaxPlane[i] / MaxLen;

    /* Make sure the Z component of the	plane is positive: */
    if (PPolygon -> Plane[2] < 0.0) {
	for (i = 0; i < 3; i++)
	    PPolygon -> Plane[i] = (-PPolygon -> Plane[i]);
	*ReversePoly = TRUE;
    }
    else
        if (GlblBackFacing)
	    return FALSE;

    PPolygon ->	Plane[3] =
	(- Coord[0] * PPolygon -> Plane[0]
	 - Coord[1] * PPolygon -> Plane[1]
	 - Coord[2] * PPolygon -> Plane[2]);

    return TRUE;
}

/*****************************************************************************
* Routine to evaluate the cross	product	of 3 points projected to Z = 0 plane *
* and return the sign of the result (Only Z component).			     *
*****************************************************************************/
int CrossProd(RealType Pt1[3], RealType Pt2[3], RealType Pt3[3])
{
    RealType Zout;

    /* U = Pt2 - Pt1,  V = Pt3 - Pt2,		Zoutput	= Ux * Vy - Uy * Vx. */
    Zout = (Pt2[0] - Pt1[0]) /*	Ux */  * (Pt3[1] - Pt2[1]) /* Vy */  -
	   (Pt2[1] - Pt1[1]) /*	Uy */  * (Pt3[0] - Pt2[0]) /* Vx */;
    if (APX_EQ(Zout, 0.0))
	return 0;
    if (Zout < 0.0)
	return -1;
    else
	return 1;
}

/*****************************************************************************
* Routine to reverse linear list PList - return	pointer	to the reversed	list *
* Although this routine is basically generic, it should be called on the     *
* VertexStruct only as it updates their Internal flags as well.		     *
*****************************************************************************/
static IPVertexStruct *ReverseLinList(IPVertexStruct *VList)
{
    int i;
    IPVertexStruct *VLtemp,
	*VLreverse = NULL;

    while (VList != NULL) {
	VLtemp = VList -> Pnext;/* Save pointer to next element in old list. */

	VList -> Pnext = VLreverse; /* Add old element in front of new list. */
	VLreverse = VList;

	VList =	VLtemp;		    /* Continue to next element in old list. */
    }

    VLtemp = VLreverse;
    i = VLtemp -> Tags;
    while (VLtemp != NULL) {
	if (VLtemp -> Pnext != NULL) {
	    VLtemp -> Tags = VLtemp -> Pnext -> Tags;
	}
	else {
	    VLtemp -> Tags = i;
	}

	VLtemp = VLtemp -> Pnext;
    }

    return VLreverse;
}

/*****************************************************************************
* Routine to generate all the edges from the given polygon in screen space.  *
* The polygon must be valid - only vertices in its list.		     *
* Edges	are inserted to	an edge	hash table of EDGE_HASH_TABLE_SIZE entries.  *
* If global variable GlblNumEdge != 0 then only the first PGlblNumEdge edges *
* are generated.							     *
* If edge is INTERNAL it is marked as so.				     *
* If this is polyline, the last edge is NOT generated.			     *
*****************************************************************************/
static void GenEdgesFromPoly(IPPolygonStruct *PPolygon, int IsPolygon)
{
    int	CountEdges = GlblNumEdge;
    EdgeStruct *PEdge;
    IPVertexStruct
	*PList = PPolygon -> PVertex;

    if (!PList || !PList -> Pnext)
	return;					 /* If less than 2 vertices. */

    while (PList -> Pnext) {
	PEdge =	(EdgeStruct *) IritMalloc(sizeof(EdgeStruct));
	PEdge -> Vertex[0] = PList;
	PEdge -> Vertex[1] = PList -> Pnext;
	PEdge -> Internal = IP_IS_INTERNAL_VRTX(PList);
	PEdge -> Pnext = NULL;
	InsertEdgeToHashTbl1(PEdge);

	if (!--CountEdges)
	    return;

	PList =	PList -> Pnext;
    }
    /* Close the contour to first vertex in list (if not polyline): */
    if (IsPolygon) {
	PEdge = (EdgeStruct *) IritMalloc(sizeof(EdgeStruct));
	PEdge -> Vertex[0] = PList;
	PEdge -> Vertex[1] = PPolygon -> PVertex;
	PEdge -> Internal = IP_IS_INTERNAL_VRTX(PList);
	PEdge -> Pnext = NULL;
	InsertEdgeToHashTbl1(PEdge);
    }
}

/*****************************************************************************
* Routine to insert new	edge to	edge hash table	structure sorted (hashed) by *
* the edge Y min value.	The edge is tested for duplicated entry	(if	     *
* interior edge	- entered twice	and ignored in this case.		     *
* Also the edge	is updated such	that Ymin will be Vertex[0], Ymax Vertex[1]. *
*****************************************************************************/
static void InsertEdgeToHashTbl1(EdgeStruct *PEdge)
{
    int	Level;
    RealType Xmin, Xmax, Ymin, Ymax;
    IPVertexStruct *PVertex;
    EdgeStruct *PEtemp;

    if (PEdge -> Vertex[0] -> Coord[1] > PEdge -> Vertex[1] -> Coord[1]) {
	PVertex	= PEdge	-> Vertex[0];
	PEdge -> Vertex[0] = PEdge -> Vertex[1];
	PEdge -> Vertex[1] = PVertex;
    }
    Xmin = MIN(PEdge -> Vertex[0] -> Coord[0],
	       PEdge -> Vertex[1] -> Coord[0]);
    Xmax = MAX(PEdge -> Vertex[0] -> Coord[0],
	       PEdge -> Vertex[1] -> Coord[0]);
    Ymin = PEdge -> Vertex[0] -> Coord[1];
    Ymax = PEdge -> Vertex[1] -> Coord[1];

    if (GlblClipScreen && ((Ymin > 1.0) || (Ymax < -1.0) ||
			   (Xmin > 1.0) || (Xmax < -1.0)))
	IritFree((VoidPtr) PEdge);			   /* Out of screen. */
    else {
	/* Normalize [-1..1] to	[0..EDGE_HASH_TABLE_SIZE]: */
	Level =	(int) ((Ymin + 1.0) * EDGE_HASH_TABLE_SIZE2);
	Level =	BOUND(Level, 0, EDGE_HASH_TABLE_SIZE1);	 /* To be 100% safe. */

	/* Look	for duplicate entry - it must have the same two	vertices: */
	PEtemp = EdgeHashTable[Level];
	while (PEtemp) {
	    /* Test to see if same edge	by comparing vertices pointers.	*/
	    if ((SAME_VERTEX(PEdge -> Vertex[0], PEtemp -> Vertex[0]) &&
		 SAME_VERTEX(PEdge -> Vertex[1], PEtemp -> Vertex[1])) ||
		(SAME_VERTEX(PEdge -> Vertex[0], PEtemp -> Vertex[1]) &&
		 SAME_VERTEX(PEdge -> Vertex[1], PEtemp -> Vertex[0]))) {
		IritFree((VoidPtr) PEdge);
		return;					/* Ignore new entry! */
	    }
	    PEtemp = PEtemp -> Pnext;
	}

	if (!GlblQuiet)
	    fprintf(stderr, "\b\b\b\b\b%5d", ++EdgeCount);
	PEdge -> Pnext = EdgeHashTable[Level];	     /* Concat to main list. */
	EdgeHashTable[Level] = PEdge;
    }
}

/*****************************************************************************
* Routine to collect all edges in hash table into one big list and intersect *
* them among themselves. In any	intersection both edges	are broken into	two. *
* The resulting	edges are inserted back	into the hash table:		     *
*****************************************************************************/
static void IntersectAllEdges(void)
{
    RealType Ymin;
    int	i, Level;
    EdgeStruct *PEmain = NULL, *PEtemp;

    IritCPUTime(TRUE);

    EdgeCount =	0;
    MinYLevel =	EDGE_HASH_TABLE_SIZE; /* Set "clip" levels in table entries. */
    MaxYLevel =	0;

    /* Clear the hash table and	collect	all edges into one big list: */
    for	(i = EDGE_HASH_TABLE_SIZE1; i >= 0; i--)
	if ((PEtemp = EdgeHashTable[i]) != NULL) {
	    while (PEtemp -> Pnext) PEtemp = PEtemp -> Pnext;
	    PEtemp -> Pnext = PEmain;
	    PEmain = EdgeHashTable[i];
	    EdgeHashTable[i] = NULL;
	    if (i > MaxYLevel)
		MaxYLevel = i;
	    if (i < MinYLevel)
		MinYLevel = i;
	}

    PrintYLevel	= CrntYLevel = 0;	 /* Have to start from some place... */
    if (!GlblQuiet)
	fprintf(stderr, "\nPass 3, Level [%5d] =      ", MaxYLevel);

    while (PEmain) { /* Insert back after intersecting with all other edges. */
	PEtemp = PEmain	-> Pnext;      /* As PEmain->Pnext might be changed. */
	InsertEdgeToHashTbl2(PEmain);
	PEmain = PEtemp;

	/* Now test to see if we can update current y level: */
	if (CrntYLevel < MaxYLevel) {
	    Ymin = MIN(PEmain -> Vertex[0] -> Coord[1],
		       PEmain -> Vertex[1] -> Coord[1]);
	    /* Normalize [-1..1] to [0..EDGE_HASH_TABLE_SIZE]: */
	    Level = (int) ((Ymin + 1.0) * EDGE_HASH_TABLE_SIZE2);
	    Level = BOUND(Level, 0, EDGE_HASH_TABLE_SIZE1); /* Be 100% safe. */
	    if (Level > CrntYLevel)
		CrntYLevel = Level;
	}
    }
    if (!GlblQuiet)
	fprintf(stderr, ",  %6.2lf seconds.", IritCPUTime(FALSE));
}

/*****************************************************************************
* Routine to insert old	edge to	edge hash table	structure sorted (hashed) by *
* the edge Y min value.	The edge is tested for intersections with other	     *
* edges	allready in structure and both edges are broken	if found one.	     *
*****************************************************************************/
static void InsertEdgeToHashTbl2(EdgeStruct *PEdge)
{
    int	i, Level, UpperLevel,
	FoundIntersection = FALSE;
    RealType Ymin, Ymax;
    EdgeStruct *PEtemp;

    Ymin = PEdge -> Vertex[0] -> Coord[1];
    Ymax = PEdge -> Vertex[1] -> Coord[1];

    /* Normalize [-1..1] to [0..EDGE_HASH_TABLE_SIZE]: */
    Level = (int) ((Ymin + 1.0)	* EDGE_HASH_TABLE_SIZE2);
    Level = BOUND(Level, 0, EDGE_HASH_TABLE_SIZE1);      /* To be 100% safe. */
    UpperLevel = 1 + (int) ((Ymax + 1.0) * EDGE_HASH_TABLE_SIZE2);
    UpperLevel = BOUND(UpperLevel, 0, EDGE_HASH_TABLE_SIZE1);

    if (CrntYLevel > PrintYLevel) {
	PrintYLevel = CrntYLevel;
	if (!GlblQuiet)
	    fprintf(stderr, "\b\b\b\b\b%5d", PrintYLevel);
    }

    /* Test for	intersections while we find intersections... */
    for	(i=MinYLevel; i<=UpperLevel; i++)
	if (EdgeHashTable[i])
	     if ((FoundIntersection =
		  IntersectEdgeList(PEdge, EdgeHashTable[i],
				    i == MinYLevel)) != 0)
		break;

    if (FoundIntersection) {	   /* Call recursively with the edge pieces: */
	while (PEdge) {
	    PEtemp = PEdge -> Pnext;  /* As Pedge->Pnext might point to new. */
	    InsertEdgeToHashTbl2(PEdge);   /* Place after the recursive ins. */
	    PEdge = PEtemp;
	}
    }
    else {		      /* Its a single edge - insert it in its place: */
	EdgeCount++;
	PEdge -> Pnext = EdgeHashTable[Level];	     /* Concat to main list. */
	EdgeHashTable[Level] = PEdge;
    }
}

/*****************************************************************************
* Routine to scan all edges in list and	intersect everything against	     *
* the given edge. intersected edges are	broken into two	parts each. The	edge *
* is updated to	a list of 2 pieces, and	the list edge is broken	and inserted *
* to the hash table (one piece in same entry as	it has the same	Ymin).	     *
* Note this routine returns TRUE after the first intersection found - no     *
* test is made for ALL intersections if	more than one exists.		     *
* A test is made if MinYLevel can be updated if	TestYMin == TRUE.	     *
*****************************************************************************/
static int IntersectEdgeList(EdgeStruct *PEdge, EdgeStruct *PEList,
								int TestYMin)
{
    int	Level,
	UpdateYMin = TRUE;
    RealType Ymin, Ymax;
    EdgeStruct *PEdgeNew, *PEListNew;

    if (!PEdge || !PEList)
	return FALSE;				       /* NULL entry - quit. */

    while (PEList) {
	if (IntersectEdgeEdge(PEdge, PEList, &PEdgeNew,	&PEListNew)) {
	    PEdge -> Pnext = PEdgeNew;
	    /* PEListNew can be	inserted to the	hash table with	no check as  */
	    /* its cannt intersect anything - it is part of checked edge!    */
	    if (PEListNew) {
		Ymin = PEListNew -> Vertex[0] -> Coord[1];
		/* Normalize [-1..1] to	[0..EDGE_HASH_TABLE_SIZE]: */
		Level =	(int) ((Ymin + 1.0) * EDGE_HASH_TABLE_SIZE2);
		Level =	BOUND(Level, 0, EDGE_HASH_TABLE_SIZE1);
		EdgeCount++;
		PEListNew -> Pnext = EdgeHashTable[Level];
		EdgeHashTable[Level] = PEListNew;
	    }
	    return TRUE;
	}
	if (TestYMin &&	UpdateYMin) {
	    Ymax = PEList -> Vertex[1] -> Coord[1];
	    /* Normalize [-1..1] to [0..EDGE_HASH_TABLE_SIZE]: */
	    Level = (int) ((Ymax + 1.0)	* EDGE_HASH_TABLE_SIZE2);
	    Level = BOUND(Level, 0, EDGE_HASH_TABLE_SIZE1);
	    if (Level >= CrntYLevel)
		UpdateYMin = FALSE;
	}
	PEList = PEList	-> Pnext;
    }
    if (TestYMin && UpdateYMin)			/* No need to test any more. */
	do
	    MinYLevel++;
	while (MinYLevel < CrntYLevel && !EdgeHashTable[MinYLevel]);

    return FALSE;
}

/*****************************************************************************
* Routine to test if two edges intersects. If they do, it brakes the bottom  *
* edge into two pieces, leaving the lower part (with the same Ymin) in	     *
* original struct and allocated and updates new struct with upper edge part. *
* Returns TRUE if found intersection, FALSE otherwise.			     *
* Note the intersection	is tested in the XY axes (Z is ignored!).	     *
*****************************************************************************/
static int IntersectEdgeEdge(EdgeStruct *PEdge1, EdgeStruct *PEdge2,
			     EdgeStruct **PEdgeNew1, EdgeStruct **PEdgeNew2)
{
    int	i, OneInter1, OneInter2;
    RealType Xmin1, Xmax1, Ymin1, Ymax1, Xmin2, Xmax2, Ymin2, Ymax2,
	  a1, b11, b12, a2, b21, b22, det, t1, t2, Z1, Z2;
    /* To speed	up the intensive access	of the coordinates: */
    RealType
	*Crd10 = PEdge1 -> Vertex[0] -> Coord,
	*Crd11 = PEdge1 -> Vertex[1] -> Coord,
	*Crd20 = PEdge2 -> Vertex[0] -> Coord,
	*Crd21 = PEdge2 -> Vertex[1] -> Coord;

    Xmin1 = MIN(Crd10[0], Crd11[0]);
    Xmax1 = MAX(Crd10[0], Crd11[0]);
    Ymin1 = Crd10[1];
    Ymax1 = Crd11[1];

    Xmin2 = MIN(Crd20[0], Crd21[0]);
    Xmax2 = MAX(Crd20[0], Crd21[0]);
    Ymin2 = Crd20[1];
    Ymax2 = Crd21[1];
    if ((Xmin1 > Xmax2)	|| (Xmax1 < Xmin2) ||/* Test if out of Boundary Box. */
	(Ymin1 > Ymax2)	|| (Ymax1 < Ymin2))
	return FALSE;

    /* Let the line equations of the two edges be defined as:		     */
    /* L1 = p11	+ t1 * (pt12 - pt11) , t1 = [0..1]			     */
    /* L2 = p21	+ t2 * (pt22 - pt21) , t2 = [0..1]			     */
    /* at intersection point (if any) we have:				     */
    /* pt11 + t1 * (pt12 - pt11) == pt21 + t2 *	(pt22 -	pt21)  for x, y	     */
    /* or two equations	(for x, y) with two unknown (t1, t2) to solve:	     */
    /* a1 = b11	* t1 + b12 * t2		from x				     */
    /* a2 = b21	* t1 + b22 * t2		from y				     */
    /* and we have interesection if both t1, t2	in the range [0..1]	     */
    a1 =  Crd10[0] - Crd20[0];
    b11	= Crd10[0] - Crd11[0];
    b12	= Crd21[0] - Crd20[0];
    a2 =  Crd10[1] - Crd20[1];
    b21	= Crd10[1] - Crd11[1];
    b22	= Crd21[1] - Crd20[1];

    /* If the detereminant is zero, the	two lines are parellel - no inter. */
    if (APX_EQ((det = b11 * b22 - b21 * b12), 0.0))
	return FALSE;

    t1 = (a1 * b22 - a2	* b12) / det;
    t2 = (b11 *	a2 - b21 * a1) / det;

    /* Test if intersection is happening in one	edge END - in that case	*/
    /* we break	only the second	edge into two parts.			*/
    OneInter1 =	((t1 < 1.0) && (t1 > 0.0) &&
		 !(APX_EQ(t1, 0.0) || APX_EQ(t1, 1.0)) &&
		  (APX_EQ(t2, 0.0) || APX_EQ(t2, 1.0)));
    OneInter2 =	((t2 < 1.0) && (t2 > 0.0) &&
		 !(APX_EQ(t2, 0.0) || APX_EQ(t2, 1.0)) &&
		  (APX_EQ(t1, 0.0) || APX_EQ(t1, 1.0)));

    /* If out of 0..1 range in one of edges - no intersection: */
    if ((!(OneInter1 ||	OneInter2)) &&
	((t1 >=	1.0) ||	(t1 <= 0.0) || (t2 >= 1.0) || (t2 <= 0.0) ||
	 APX_EQ(t1, 0.0) || APX_EQ(t1, 1.0) ||
	 APX_EQ(t2, 0.0) || APX_EQ(t2, 1.0))) return FALSE;

    /* If we are here, we have intersection - find the bottom edge and split */
    /* it - allocated new edge struct and update to new upper (in Y) part.   */
    Z1 = Crd10[2] * (1.0 - t1) + Crd11[2] * t1;
    Z2 = Crd20[2] * (1.0 - t2) + Crd21[2] * t2;
    if (!OneInter2 && Z1 < Z2) {
	*PEdgeNew1 = (EdgeStruct *) IritMalloc(sizeof(EdgeStruct));
	(*PEdgeNew1) -> Pnext = NULL;
	(*PEdgeNew1) -> Internal = PEdge1 -> Internal;
	(*PEdgeNew1) ->	Vertex[0] = IPAllocVertex(0, 0, NULL, NULL);
	for (i = 0; i < 2; i++)
	    (*PEdgeNew1) -> Vertex[0] -> Coord[i] =
		Crd10[i] * (1.0	- t1) +	Crd11[i] * t1;
	(*PEdgeNew1) -> Vertex[0] -> Coord[2] = Z1;
	/* Now update the second vertex	of both	PEdge1 & PEdgeNew1:	   */
	/* Note	we assume Vertex[0] -> Coord[1]	< Vertex[1] -> Coord[1]	as */
	/* all input edges are sorted this way when entered to hash table. */
	(*PEdgeNew1) ->	Vertex[1] = PEdge1 -> Vertex[1];
	PEdge1 -> Vertex[1] = (*PEdgeNew1) -> Vertex[0];

    }
    else
	*PEdgeNew1 = NULL;

    if (!OneInter1 && Z2 < Z1) {
	*PEdgeNew2 = (EdgeStruct *) IritMalloc(sizeof(EdgeStruct));
	(*PEdgeNew2) -> Pnext = NULL;
	(*PEdgeNew2) -> Internal = PEdge2 -> Internal;
	(*PEdgeNew2) ->	Vertex[0] = IPAllocVertex(0, 0, NULL, NULL);
	for (i = 0; i < 2; i++)
	    (*PEdgeNew2) -> Vertex[0] -> Coord[i] =
		Crd20[i] * (1.0	- t2) +	Crd21[i] * t2;
	(*PEdgeNew2) -> Vertex[0] -> Coord[2] = Z2;
	/* Now update the second vertex	of both	PEdge2 & PEdgeNew2: */
	(*PEdgeNew2) ->	Vertex[1] = PEdge2 -> Vertex[1];
	PEdge2 -> Vertex[1] = (*PEdgeNew2) -> Vertex[0];
    }
    else
	*PEdgeNew2 = NULL;

    return (*PEdgeNew1 != NULL) || (*PEdgeNew2 != NULL);
}

/*****************************************************************************
* Routine to print the content of a given edge:				     *
*****************************************************************************/
static void PrintPolyContent(IPPolygonStruct *PPoly)
{
    IPVertexStruct
	*PList = PPoly -> PVertex;

    while (PList) {
	fprintf(stderr, "   %12f %12f %12f\n",
	    PList -> Coord[0],
	    PList -> Coord[1],
	    PList -> Coord[2]);
	PList =	PList -> Pnext;
    }
}

#ifdef DEBUG1

/*****************************************************************************
* Routine to print the content of a given edge:				     *
*****************************************************************************/
static void PrintEdgeContent(EdgeStruct *PEdge)
{
    fprintf(stderr, "   %11f %11f %11f : %11f %11f %11f\n",
	PEdge -> Vertex[0] -> Coord[0],
	PEdge -> Vertex[0] -> Coord[1],
	PEdge -> Vertex[0] -> Coord[2],
	PEdge -> Vertex[1] -> Coord[0],
	PEdge -> Vertex[1] -> Coord[1],
	PEdge -> Vertex[1] -> Coord[2]);
}

/*****************************************************************************
* Routine to draw all the segments in the EdgeHashTable:		     *
*****************************************************************************/
static void DrawEdgeHashTable(void)
{
    int	i;
    EdgeStruct *PEtemp;

    for	(i=0; i<EDGE_HASH_TABLE_SIZE; i++) {
	PEtemp = EdgeHashTable[i];
	while(PEtemp) {
	    DrawEdge(PEtemp);
	    PEtemp = PEtemp -> Pnext;
	}
    }
}

/*****************************************************************************
* Routine to draw an edge:						     *
*****************************************************************************/
static void DrawEdge(EdgeStruct *PEdge)
{
    printf("    [POINTLIST 2\n\t[%lf %lf %lf]\n\t[%lf %lf %lf]\n    ]\n",
	   PEdge->Vertex[0]->Coord[0],
	   PEdge->Vertex[0]->Coord[1],
	   PEdge->Vertex[0]->Coord[2],
	   PEdge->Vertex[1]->Coord[0],
	   PEdge->Vertex[1]->Coord[1],
	   PEdge->Vertex[1]->Coord[2]);
    printf("    [POLYLINE 2\n\t[%lf %lf %lf]\n\t[%lf %lf %lf]\n    ]\n",
	   PEdge->Vertex[0]->Coord[0],
	   PEdge->Vertex[0]->Coord[1],
	   PEdge->Vertex[0]->Coord[2],
	   PEdge->Vertex[1]->Coord[0],
	   PEdge->Vertex[1]->Coord[1],
	   PEdge->Vertex[1]->Coord[2]);
}

#endif /* DEBUG1 */