/*****************************************************************************
*   "Irit" - the 3d polygonal solid modeller.				     *
*									     *
* Written by:  Gershon Elber				Ver 0.2, Mar. 1990   *
******************************************************************************
*   Module to handle the low level boolean operations. The routines in this  *
* module should be accessed from "bool-hi.c" module only.		     *
*   Note the polygons of the two given objects may not be convex, and must   *
* be subdivided into convex ones in the boolean upper level routines (module *
* bool-hi.c). All the routines of this module expects objects with convex    *
* polygons only although they might return objects with non convex polygons, *
* but marked as so (via polygons CONVEX tags - see Irit.h)!		     *
*   Because Bool-low.c module was too big, it was subdivided to two:	     *
* Bool1Low.c - mainly handles the intersection polyline between the oper.    *
* Bool2Low.c - mainly handles the polygon extraction from operands given the *
*	       polyline of intersection and the adjacencies (see ADJACNCY.C) *
* Note we said mainly has routines CAN call one another!		     *
*****************************************************************************/

/* #define DEBUG	    If defined, return intersection POLYLINE object. */
/* #define DEBUG2	         If defined, defines some printing routines. */
/* #define DEBUG3	         Print messages to entry/exit from routines. */

#ifdef __MSDOS__
#include <alloc.h>
#include <dos.h>
#endif /* __MSDOS__ */

#include <stdio.h>
#include <ctype.h>
#include <math.h>
#include <string.h>
#include "program.h"
#include "allocatg.h"
#include "convexg.h"
#include "objects.h"
#include "booleang.h"
#include "booleanl.h"
#include "geomat3d.h"

static int TestAinB(VertexStruct *V, PolygonStruct *Pl, int AinB);
static struct PolygonStruct *ClipOpenLoopFromPoly(PolygonStruct *Pl,
					InterSegListStruct *OLoop, int AinB);
static void ClosedLoopsToPolys(InterSegListStruct *PClosed, PolygonStruct *Pl);
static VertexStruct *CutPolygonAtRay(PolygonStruct *Pl, PointType Pt);
static CombineClosedLoops(PolygonStruct *Pl, InterSegListStruct *PClosed,
								int AinB);
static void PushAdjOnStack(PolygonStruct *Pl, PolygonStruct *AdjStack[],
							 int *StackPointer);
static struct PolygonStruct *ChainPolyLists(PolygonStruct *Pl1,
							PolygonStruct *Pl2);

#ifdef DEBUG2
static void PrintVrtxList(VertexStruct *V);
static void PrintInterList(InterSegmentStruct *PInt);
#endif /* DEBUG2 */

/*****************************************************************************
*   Test on which side of polygon Pl, given Point V is, and according to     *
* the requiremnet (AinB) returns TRUE/FALSE.				     *
*   If test on V fails, we tries the next one V -> Pnext until success...    *
*****************************************************************************/
static int TestAinB(VertexStruct *V, PolygonStruct *Pl, int AinB)
{
    int In;
    RealType Distance;
    VertexStruct *VHead = V;

    do {
	Distance = Pl -> Plane[0] * V -> Pt[0] +
		   Pl -> Plane[1] * V -> Pt[1] +
		   Pl -> Plane[2] * V -> Pt[2] + Pl -> Plane[3];
	In = Distance > 0.0;
	V = V -> Pnext;
    }
    while (ABS(Distance) < EPSILON && V != NULL && V != VHead);

    return (In && AinB) || (!In && !AinB);   /* I wish I has logical XOR ... */
}

/*****************************************************************************
*   Convert an inter loop into an open vertex list.			     *
*****************************************************************************/
struct VertexStruct *InterLoopToVrtxList(InterSegmentStruct *PIHead)
{
    VertexStruct *VHead, *V;

    VHead = V = AllocVertex(0, 0, NULL, NULL);

    PT_COPY(VHead -> Pt, PIHead -> PtSeg[0]);

    while (PIHead != NULL) {
	V -> Pnext = AllocVertex(0, 0, NULL, NULL);
	V = V -> Pnext;
	PT_COPY(V -> Pt, PIHead -> PtSeg[1]);

	PIHead = PIHead -> Pnext;
    }

    V -> Pnext = NULL;

    return VHead;
}

/*****************************************************************************
*   Clip an open loop from a polygon:					     *
* 1. Clip the section (S) of the polygon between the two loop end points and *
*    replace it by the loop itself.					     *
* 2. If the polygon formed from S and the loop should be in the output       *
*    (as tested by AinB) return that polygon. Otherwise return NULL.	     *
*   The open loop itself (excluding the header OLoop) is freed.		     *
* Note it is assumed (ordered by the sorting routines above) that the open   *
* loop starts from the second end back to first:			     *
*									     *
*			       L1-----------------------L2		     *
*				|			|		     *
*				|L0			|L3		     *
*	*---------------*-------+----*-------------*----+-----------*	     *
*	P0		P1	     P2		   P3		    P0	     *
*****************************************************************************/
static struct PolygonStruct *ClipOpenLoopFromPoly(PolygonStruct *Pl,
					InterSegListStruct *OLoop, int AinB)
{
    int GenClipPoly;	    /* TRUE if needs to form polygon from S & OLoop. */
    PointType Pt;
    VertexStruct *VStart,  *VEnd, *VEnd1,     /* Corresponds to L0 and L3... */
		 *ClipPoly,		/* The clipped element from polygon. */
		 *PLoop, *PRevLoop,	  /* The loop itself as vertex list. */
		 *PLoopEnd, *PLoopEnd1, *Ptemp1, *Ptemp2;
    InterSegmentStruct *PISeg, *PItemp;
    PolygonStruct *ClipPl;

#ifdef DEBUG3
    printf("Enter ClipOpenLoopFromPoly\n");
#endif /* DEBUG3 */

    PISeg = OLoop -> PISeg;
    VStart = PISeg -> V[0];
    while (PISeg -> Pnext != NULL) PISeg = PISeg -> Pnext;
    VEnd = PISeg -> V[1];
    if (VStart == NULL || VEnd == NULL)
	FatalError("ClipOpenLoopFromPoly: None open loop\n");
    VEnd1 = VEnd;		   /* Find the pointer thats points on VEnd. */
    while (VEnd1 -> Pnext != VEnd) VEnd1 = VEnd1 -> Pnext;

    PLoop = InterLoopToVrtxList(OLoop -> PISeg);/* Make vertex list out of it*/
    PLoopEnd = PLoop;			  /* Prepare pointer on last vertex. */
    while (PLoopEnd -> Pnext != NULL) PLoopEnd = PLoopEnd -> Pnext;
    PLoopEnd1 = PLoop;
    while (PLoopEnd1 -> Pnext != PLoopEnd) PLoopEnd1 = PLoopEnd1 -> Pnext;

    if (VStart != VEnd) {
	/* Now we test if we need to create a polygon formed from S & open   */
	/* loop by testing on which side of the polygon that caused	     */
	/* intersection L0L1, point P2 is, and compare with requirement AinB.*/
	GenClipPoly = TestAinB(VStart -> Pnext, OLoop -> PISeg -> Pl, AinB);

	/* Keep the clipped VertexList P2, P3 & substitute with open loop:   */
	/* Note we must keep vertex VEnd in the polygon as another InterSeg. */
	/* may point on it, so we replace InterSeg point (L3) by (P3) and    */
	/* leave VEnd intact.						     */
	Ptemp1 = VEnd -> Pnext;		     /* Save that pointer temporary. */
	VEnd -> Pnext = NULL;		   /* Close the clipped vertex list. */

	PT_COPY(Pt, VEnd -> Pt);
	PT_COPY(VEnd -> Pt, PLoopEnd -> Pt);
	PT_COPY(PLoopEnd -> Pt, Pt);
	VEnd1 -> Pnext = PLoopEnd;
	PLoopEnd1 -> Pnext = VEnd;
	PLoopEnd -> Count = VEnd -> Count;
	PLoopEnd -> Tags = VEnd -> Tags;

	Ptemp2 = VEnd;
	VEnd = PLoopEnd;
	PLoopEnd = Ptemp2;

	ClipPoly = VStart -> Pnext;

	/* New ClipPoly is isolated (Open loop of P2, P3 only). Save	     */
	/* reversed list of open loop if we need to form an S/OLoop polygon, */
	/* otherwise free ClipPoly. Chain the OLoop instead of S.	     */
	if (GenClipPoly) PRevLoop = GenReverseVrtxList(PLoop);
	else MyFree((char *) ClipPoly, VERTEX_TYPE);
	VStart -> Pnext = PLoop;	    /* Chain the OLoop instead of S. */
	PLoopEnd -> Pnext = Ptemp1;
    }
    else { /* VStart == VEnd */
	/* Now we test if we need to create a polygon formed from S & open   */
	/* loop by testing on which side of the polygon that caused	     */
	/* intersection L0L1, point L3 is, and compare with requirement AinB.*/
	GenClipPoly = TestAinB(PLoopEnd, OLoop -> PISeg -> Pl, AinB);

	/* In this case the clipped part is empty, so its simpler: */
	ClipPoly = NULL;

	/* Save reversed list of open loop if we need to form an S/OLoop     */
	/* polygon. Chain the OLoop instead of S.			     */
	if (GenClipPoly) PRevLoop = GenReverseVrtxList(PLoop);

	PLoopEnd -> Pnext = VEnd -> Pnext;
	VStart -> Pnext = PLoop;	    /* Chain the OLoop instead of S. */
    }

    /* Time to free the InterSegment list pointed by OLoop: */
    PISeg = OLoop -> PISeg;
    while (PISeg != NULL) {
	PItemp = PISeg;
	PISeg = PISeg -> Pnext;
	MyFree((char *) PItemp, OTHER_TYPE);
    }
    OLoop -> PISeg = NULL;			/* To be on the safe side... */

    /* There is a chance that Pl -> V will point on vertex in the clipped    */
    /* part so we update it to point on VStart, which for sure is in polygon.*/
    Pl -> V = VStart;
    if (GenClipPoly) {	   /* Generate the polygon from ClipPoly & PRevLoop: */
	PLoopEnd = PRevLoop;
	while (PLoopEnd -> Pnext != NULL) PLoopEnd = PLoopEnd -> Pnext;

	if (ClipPoly == NULL) {
	    PLoopEnd -> Pnext = PRevLoop;  /* Close that loop and return it. */
	    ClipPl = AllocPolygon(0, 0, PRevLoop, NULL);
	    PLANE_COPY(ClipPl -> Plane, Pl -> Plane);
	    RST_CONVEX_POLY(ClipPl);		   /* May be not convex now. */
	    return ClipPl;
	}

	PLoopEnd -> Pnext = ClipPoly;
	PLoopEnd -> Tags = VStart -> Tags;
	PLoopEnd -> Count = VStart -> Count;
	Ptemp1 = ClipPoly;
	while (Ptemp1 -> Pnext != NULL) Ptemp1 = Ptemp1 -> Pnext;
	Ptemp1 -> Pnext = PRevLoop;

	ClipPl = AllocPolygon(0, 0, ClipPoly, NULL);
	PLANE_COPY(ClipPl -> Plane, Pl -> Plane);
	RST_CONVEX_POLY(ClipPl);		   /* May be not convex now. */
	return ClipPl;
    }
    else return NULL;
}

/*****************************************************************************
* Find the intersection of the ray fired from Pt to +X direction with the    *
* given polygon. Note Pt MUST be in the polygon. Two vertices equal to	     *
* ray/polygon intersection point are added to polygon vertex list, and a     *
* pointer to the first one is also returned. This routine is exclusively     *
* used by the CombineClosedLoops below.					     *
* The polygon is NOT assumed to be convex and we look for the minimum X      *
* intersection. The polygon might not be convex as a result of combinning    *
* some other closed loop before the current one...			     *
*****************************************************************************/
static VertexStruct *CutPolygonAtRay(PolygonStruct *Pl, PointType Pt)
{
    int OnVertex;
    RealType MinX = INFINITY, X;
    VertexStruct *V, *Vnext, *VMinX = NULL;

    V = Pl -> V;
    do {
	Vnext = V -> Pnext;

	/* A polygon edge might intersect the ray iff one of the following:  */
	/* 1. The first vertex is exactly on the ray Y level. (if this is    */
	/*    true for the second edge, it will be detected next iteration). */
	/* 2. The first vertex is below ray Y level, and second is above.    */
	/* 3. The first vertex is above ray Y level, and second is below.    */
	if (APX_EQ(V -> Pt[1], Pt[1])) {		    /* Case 1 above. */
	    if (MinX > V -> Pt[0] && Pt[0] < V -> Pt[0]) {
		OnVertex = TRUE;
		MinX = V -> Pt[0];
		VMinX = V;
	    }
	}
	else
	if ((V -> Pt[1] < Pt[1] && Vnext -> Pt[1] > Pt[1]) ||/* Case 2 above.*/
	    (V -> Pt[1] > Pt[1] && Vnext -> Pt[1] < Pt[1])) {/* Case 3 above.*/
	    X = ((Vnext -> Pt[1] - Pt[1]) * V -> Pt[0] +
		 (Pt[1] - V -> Pt[1]) * Vnext -> Pt[0]) /
		(Vnext -> Pt[1] - V -> Pt[1]);
	    if (MinX > X && Pt[0] < X) {
		OnVertex = FALSE;
		MinX = X;
		VMinX = V;
	    }

	}

	V = Vnext;
    }
    while (V != NULL && V != Pl -> V);

    if ((V = VMinX) == NULL)
	FatalError("CutPolygonAtRay: fail to find intersection");

    /* Now that we have the intersection point - create two new vertices     */
    /* (one if OnVertex), insert them (it) after V and return the first.     */
    if (OnVertex) {
	V -> Pnext = AllocVertex(V -> Count, V -> Tags, NULL, V -> Pnext);
	PT_COPY(V -> Pnext -> Pt, V -> Pt);
	V -> Tags = V -> Count = 0;
    }
    else {
	V -> Pnext = AllocVertex(V -> Count, V -> Tags, NULL, V -> Pnext);
	Vnext = V -> Pnext;
	Vnext -> Pt[0] = MinX;		 /* X - as intersection point found. */
	Vnext -> Pt[1] = Pt[1];		      /* Y - must be as ray Y level. */
	Vnext -> Pt[2] = V -> Pt[2];    /* Z - all polygon has same Z value. */

	V -> Pnext = AllocVertex(0, 0, NULL, V -> Pnext);
	V = V -> Pnext;
	PT_COPY(V -> Pt, Vnext -> Pt);
    }

    return V;
}

/*****************************************************************************
*   Convert the given closed loop list to polygons, and return them. the     *
* original given polygon vertex list is freed, and the first loop is subst.  *
* instead.								     *
*****************************************************************************/
static void ClosedLoopsToPolys(InterSegListStruct *PClosed, PolygonStruct *Pl)
{
    int LoopNum = 0;
    VertexStruct *V, *VHead;
    InterSegmentStruct *PISeg, *PItemp;
    InterSegListStruct *PClosedTemp;

    MyFree((char *) Pl -> V, VERTEX_TYPE);
    Pl -> Pnext = NULL;
    SET_INOUTPUT_POLY(Pl);		   /* Mark the polygon as in output. */

    while (PClosed != NULL) {
	/* Convert the InterList to vertex list and free the first: */
	V = VHead = InterLoopToVrtxList(PClosed -> PISeg);
	if (V -> Pnext == NULL || V -> Pnext -> Pnext == NULL)
	    FatalError("ClosedLoopsToPolys: Closed loop with less than 3 vertices");

	PISeg = PClosed -> PISeg;      /* Time to free the InterSegmentList: */
	while (PISeg != NULL) {
	    PItemp = PISeg;
	    PISeg = PISeg -> Pnext;
	    MyFree((char *) PItemp, OTHER_TYPE);
	}

	while (V -> Pnext -> Pnext != NULL) V = V -> Pnext; /* Go to last pt */
	MyFree((char *) V -> Pnext, VERTEX_TYPE);/* and free - same as first.*/
	V -> Pnext = VHead;		       /* Make vertex list circular. */

	if (LoopNum++) {
	    Pl -> Pnext = AllocPolygon(0, 0, VHead, Pl -> Pnext);
	    PLANE_COPY(Pl -> Pnext -> Plane, Pl -> Plane);
	}
	else {
	    Pl -> V = VHead;
	}

	PClosedTemp = PClosed;
	PClosed = PClosed -> Pnext;
	MyFree((char *) PClosedTemp, OTHER_TYPE);
    }
}

/*****************************************************************************
*   This routines cuts the given polygon into its interior closed loops by   *
* adding an edge from the polygon boundary to each of its closed loops.      *
*   For example:							     *
*									     *
* +-----------------------+      +-----------------------+		     *
* |                       |	 |                       |		     *
* |  / \         / \      |	 |  / \________ / \      |		     *
* |  \ /        |   |     |	 |  \ /        |   |_____|		     *
* |      _       \ /      | -->  |      _       \ /      |		     *
* |     / \_              | -->  |     / \_              |		     *
* |    |    |             |	 |    |    |_____________|		     *
* |     \__/              |	 |     \__/              |		     *
* |                       |      |                       |		     *
* +-----------------------+      +-----------------------+		     *
*									     *
*   Algorithm:								     *
* 1. Transform the polygon and the closed loops to a plane parallel to XY    *
*    plane (Z = Const plane).						     *
* 2. For each loop find its MaxX while keeping the information on the        *
*    vertex on that extremum.						     *
* 3. For the loop with the biggest MaxX:				     *
*    3.1. Use that extremum vertex (which must be non concave corner) to     *
*         test if loop is in the reverse direction the polygon itself is,    *
*         and reverse it if not.					     *
*    3.2. Fire a ray from the extremum vertex, to the +X direction outside   *
*         of the loop till it intersect the polygon, break the polygon at    *
*         that point and add two edges to beginning of loop from breaking    *
*         point and from end of loop to breaking point (beginning/end point  *
*	  of loop is the extremum vertex point).			     *
* 4. Repeat step 3, with all loops.					     *
* 5. Transfrom the new polygon back (using the inverse matrix of step 1)     *
*    to its original orientation.					     *
*									     *
* Returns TRUE iff the original polygon boundary is in output.		     *
*									     *
*****************************************************************************/
static int CombineClosedLoops(PolygonStruct *Pl, InterSegListStruct *PClosed,
								int AinB)
{
    RealType MaxX;
    PointType V1, V2, Normal, PlNormal;
    MatrixType RotMat;
    VertexStruct *V, *Vnext, *Vprev, *VHead, *VMaxX, VStruct;
    InterSegListStruct *PClosedTemp, *PClosedMaxX, *PClosedLast,
	*PClosedMaxXLast;
    InterSegmentStruct *PISeg, *PItemp, *PISegMaxX;

    Pl -> Pnext = NULL;		  /* Make sure this polygon is disconnected. */

    /* Stage 1 - Transform the vertices to a plane parallel to XY plane: */
    GenRotateMatrix(RotMat, Pl -> Plane);
    V = VHead = Pl -> V;
    do {				    /* Transform the polygon itself. */
	MatMultVecby4by4(V -> Pt, V -> Pt, RotMat);
	V = V -> Pnext;
    }
    while (V != NULL && V != VHead);

    PClosedTemp = PClosed;
    while (PClosedTemp != NULL) {	    /* And transform the loops also. */
	PISeg = PClosed -> PISeg;
	while (PISeg != NULL) {
	    MatMultVecby4by4(PISeg -> PtSeg[0], PISeg -> PtSeg[0], RotMat);
	    MatMultVecby4by4(PISeg -> PtSeg[1], PISeg -> PtSeg[1], RotMat);

	    PISeg = PISeg -> Pnext;
	}

	PClosedTemp = PClosedTemp -> Pnext;
    }
    if (!MatInverseMatrix(RotMat, RotMat))	 /* Find the inverse matrix. */
	FatalError("CombineClosedLoops: Inverse matrix does not exits");

    /* Evalaute the normal to the polygon (which must be convex!). Note we   */
    /* cannt simply use the Plane normal as the polygon was transformed.     */
    V = Pl -> V;
    do {
	Vnext = V -> Pnext;
	PT_SUB(V1, Vnext -> Pt, V -> Pt);
	PT_SUB(V2, Vnext -> Pnext -> Pt, Vnext -> Pt);
	VecCrossProd(PlNormal, V1, V2);
	V = Vnext;
    }
    while (PT_LENGTH(PlNormal) < EPSILON);

    PClosedTemp = PClosed;
    while (PClosedTemp != NULL) {
	/* Stage 2 - find MaxX extremum value of given loop, test the loop   */
	/* direction, and reverse it if wrong. Note we convert the loop      */
	/* given as InterSegListStruct list into VertexStruct list first.    */
	PISegMaxX = PISeg = PClosedTemp -> PISeg;
	MaxX = PISeg -> PtSeg[0][0];		  /* Get first vertex X val. */
	PISeg = PISeg -> Pnext;
	while (PISeg)
	{
	    if (PISeg -> PtSeg[0][0] > MaxX) {
		MaxX = PISeg -> PtSeg[0][0];
		PISegMaxX = PISeg;
	    }
	    PISeg = PISeg -> Pnext;
	}
	PClosedTemp -> PISegMaxX = PISegMaxX;

	PClosedTemp = PClosedTemp -> Pnext;
    }

    /* Stage 3/4 - find the loop with biggest MaxX and combine it with the   */
    /* polygon itself. Do it for all closed loops, in list:		     */
    PClosedTemp = PClosed;
    while (PClosed != NULL) {
	/* Find the one with maximum MaxX, and delete it from PClosed list.  */
	PClosedLast = PClosedMaxX = PClosedTemp = PClosed;
	MaxX = PClosedMaxX -> PISegMaxX -> PtSeg[0][0];
	while (PClosedTemp != NULL)
	{
	    if (PClosedTemp -> PISegMaxX -> PtSeg[0][0] > MaxX) {
		PClosedMaxX = PClosedTemp;
		PClosedMaxXLast = PClosedLast;
	    }
	    PClosedLast = PClosedTemp;
	    PClosedTemp = PClosedTemp -> Pnext;
	}
	if (PClosedMaxX == PClosed) PClosed = PClosed -> Pnext; /* Del. from */
	else PClosedMaxXLast -> Pnext = PClosedMaxX -> Pnext;/* PClosed list.*/

	/* Create a vertex list from the loop: */
	V = VHead = InterLoopToVrtxList(PClosedMaxX -> PISeg);
	if (V -> Pnext == NULL || V -> Pnext -> Pnext == NULL)
	    FatalError("CombineClosedLoops: Closed loop with less than 3 vertices");

	V = VHead;
	while (V -> Pnext -> Pnext != NULL) V = V -> Pnext; /* Go to last pt */
	MyFree((char *) V -> Pnext, VERTEX_TYPE);/* and free - same as first.*/
	V -> Pnext = VHead;		       /* Make vertex list circular. */

	PISegMaxX = PClosedMaxX -> PISegMaxX;

	/* Now test if the vertex list should be reversed. Find a vertex     */
	/* which is below MaxX but its successor is on MaxX. The 3 vertices  */
	/* V, V -> Pnext (on MaxX), V -> Pnext -> Pnext, must form convex    */
	/* corner which we use to test if the loop needs to be reversed:     */
	while (APX_EQ(V -> Pt[0], MaxX)) V = V -> Pnext;
	while (!APX_EQ(V -> Pnext -> Pt[0], MaxX)) V = V -> Pnext;
	VMaxX = V -> Pnext;
	PT_COPY(VStruct.Pt, V -> Pt);  /* Prepare in point in REAL position. */
	MatMultVecby4by4(VStruct.Pt, VStruct.Pt, RotMat);
	if (TestAinB(&VStruct, PISegMaxX -> Pl, AinB)) {
	    /* The Inside of the object is actually the loop itself. In that */
	    /* case we simply return all the loops converted into polygon.   */
	    /* This case is simple...					     */
	    MyFree((char *) VHead, VERTEX_TYPE);   /* Free loop vertex list. */
	    PClosedMaxX -> Pnext = PClosed;/* Put back first loop into list. */
	    PClosedTemp = PClosed = PClosedMaxX;
	    while (PClosedTemp != NULL) {	/* Transform the loops back. */
		PISeg = PClosed -> PISeg;
		while (PISeg != NULL) {
		    MatMultVecby4by4(PISeg -> PtSeg[0], PISeg -> PtSeg[0], RotMat);
		    MatMultVecby4by4(PISeg -> PtSeg[1], PISeg -> PtSeg[1], RotMat);

		    PISeg = PISeg -> Pnext;
		}
		PClosedTemp = PClosedTemp -> Pnext;
	    }

	    ClosedLoopsToPolys(PClosedMaxX, Pl);/* And convert them to polys.*/
	    return FALSE;	   /* Boundary is NOT part of object result. */
	}
	PT_SUB(V1, VMaxX -> Pt, V -> Pt);
	PT_SUB(V2, VMaxX -> Pnext -> Pt, VMaxX -> Pt);
	VecCrossProd(Normal, V1, V2);
	if (DOT_PROD(Normal, PlNormal) > 0) {	    /* Need to reverse list. */
	    Vprev = VHead;
	    V = Vprev -> Pnext;
	    Vnext = V -> Pnext;
	    do {
		V -> Pnext = Vprev;/* Reverse to point on prev instead next. */
		Vprev = V;
		V = Vnext;
		Vnext = V -> Pnext;
	    }
	    while (Vprev != VHead);
	}

	PISeg = PClosedMaxX -> PISeg;  /* Time to free the InterSegmentList: */
	while (PISeg != NULL) {
	    PItemp = PISeg;
	    PISeg = PISeg -> Pnext;
	    MyFree((char *) PItemp, OTHER_TYPE);
	}
	MyFree((char *) PClosedMaxX, OTHER_TYPE);

	/* O.k. lets fire a ray from VMaxX to +X direction and see where it  */
	/* intersects the polygon. The routine CutPolygonAtRay will add two  */
	/* vertices at the ray intersection into polygon vertex list and     */
	/* return a pointer to first one, so we can chain our loop directly  */
	/* between these two new vertices.				     */
	V = CutPolygonAtRay(Pl, VMaxX -> Pt);
	Vnext = V -> Pnext;
	/* Introduce a copy of VMaxX and successor to VMaxX: */
	VMaxX -> Pnext = AllocVertex(VMaxX -> Count, VMaxX -> Tags, NULL,
								VMaxX -> Pnext);
	PT_COPY(VMaxX -> Pnext -> Pt, VMaxX -> Pt);
	V -> Pnext = VMaxX -> Pnext;
	SET_INTERNAL_EDGE(V);
	VMaxX -> Pnext = Vnext;
	SET_INTERNAL_EDGE(VMaxX);
    }

    /* Stage 5 - Time to rotate polygon back to its original position. */
    V = VHead = Pl -> V;
    do {
	MatMultVecby4by4(V -> Pt, V -> Pt, RotMat);
	V = V -> Pnext;
    }
    while (V != NULL && V != VHead);

    SET_INOUTPUT_POLY(Pl);		   /* Mark the polygon as in output. */
    RST_CONVEX_POLY(Pl);	     /* This polygon is not convex any more. */

    return TRUE;
}

/*****************************************************************************
*   Push on the adjacency stack, all adjacent polygons which are complete    *
* (no intersection) and are adjacent to complete edges (originated from	     *
* input polygons) of the given polygon.					     *
*****************************************************************************/
static void PushAdjOnStack(PolygonStruct *Pl, PolygonStruct *AdjStack[],
							 int *StackPointer)
{
    struct VertexStruct *V = Pl -> V;

    do {
	/* If you wants to propagate iff both vertices are original then     */
	/* uncomment the next line.					     */
	if (IS_ORIGINAL_VRTX(V) /* && IS_ORIGINAL_VRTX(V -> Pnext) */
	    && V -> PAdj != NULL)
	    AdjStack[++*StackPointer] = V -> PAdj;
	if (*StackPointer >= ADJ_STACK_SIZE) {
	    FatalError("Adjacency stack overflow, object too big\n");
	}
	V = V -> Pnext;
    }
    while (V != Pl -> V);
}

/*****************************************************************************
*   Chain two Polygon lists into one. For fast processing it is prefered the *
* first one to be to shorter. Returns pointer to chained list.		     *
*****************************************************************************/
static struct PolygonStruct *ChainPolyLists(PolygonStruct *Pl1,
							PolygonStruct *Pl2)
{
    PolygonStruct *Ptemp;

    if (Pl1 == NULL) return Pl2;
    else
    if (Pl2 == NULL) return Pl1;
    else {
	Ptemp = Pl1;
	while (Ptemp -> Pnext != NULL) Ptemp = Ptemp -> Pnext;
	Ptemp -> Pnext = Pl2;
	return Pl1;
    }
}

/*****************************************************************************
*   This routine coordinates all the extraction of the polygons from the     *
* intersecting lists. Does it in the following steps:			     *
* 1. Mark all polygons with no intersection at all as complete polygons.     *
*    (this is cause that this polygon will be totally in or out, according   *
*    to inter-polygon adjacencies propagation...)			     *
*    Also Mark them as undefined (if in output or undefined) yet.	     *
*    Uses PolygonStruct Tags to save these bits.			     *
* 2. 2.1. Convert the unordered segment list of each polygon to closed loops *
*         (if create a hole in polygon) or open (if crosses its boundary).   *
*    2.2. Order the open loops along the perimeter of the polygon (note      *
*	  these loops cannt intersect. For example (5 vertices polygon):     *
*	         -----------------------------L3			     *
*		|  ---------------L1  -----L2 |          --------L4   --L5   *
*		| |               |  |     |  |         |        |   |  |    *
*	  P0 ------ P1 ------- P2 ----- P3 -------- P4 ------ P5 -------- P0 *
*         Note L1, L2 are enclosed in L3 loop, and that the order is circular*
*    2.3. "Break" the polygon at each open loop that has no enclosed loops   *
*	  in it. For example we can start at L1, L2, L4, L5 and then L3.     *
*	  "Break" means - replace the vertex chain between the two loop end  *
*	  points, with the loops itself. Depends upon the relation required  *
*	  we may need to output a new polygon form from the deleted chain    *
*	  and that loop. In addition we may form a new polygon from last     *
*	  loop and was was left from the original polygon		     *
*	  For each formed polygon, for each complete edge of it (i.e. edge   *
*	  which was originally in the polygon) test the adjacent polygon     *
*	  if it is complete (as marked in 1.) and if in or undefined (marked *
*	  undefined in 1.) is still undefined:				     *
*	  2.3.1. set it to be in.					     *
*	  2.3.2. push it on adjacency stack.				     *
*    2.4. For each closed loop - find in which polygon (from polygons	     *
*	  created in 2.3.) it is enclosed, and decompose it.		     *
* 3. While adjacencies stack not empty do:				     *
*    3.1. pop top polygon from stack and output it.			     *
*    3.2. For each of its edges (which obviousely must be complete edges)    *
*	  if adjacent polygon is complete and undefined:		     *
*	  3.3.1. set it to be in.					     *
*	  3.3.2. push it on adjacency stack.				     *
*    3.3  go back to 3.							     *
*									     *
*   The above algorithm defines in as in output, but dont be confused with   *
* the required inter-object AinB (or AoutB if FALSE), which used to	     *
* determine which side of the trimming loop should be output.		     *
*   Note this routine may return non-convex polygons (but marked as so) even *
* though the input for the booleans must be convex polygons only!	     *
*   In order to keep the given object unchanged, a whole new copy off the    *
* polygon list is made. The polygons of the list that are not in the output  *
* are freed: a global list of all polygons (pointers) is used to scan them   *
* in the end and free the unused ones (list PolysPtr).			     *
*****************************************************************************/
ObjectStruct *ExtractPolygons(ObjectStruct *PObj, int AinB)
{
    int NumOfPolys = 0, StackPointer = -1;
    struct PolygonStruct *AdjStack[ADJ_STACK_SIZE], **PolysPtr,
		*Pl, *PlHead, *PlNext, *OutputPl = NULL, *SplPl, *NewPl;
    struct VertexStruct *V, *Vnext;
    struct InterSegListStruct *PClosed, *POpen, *Ptemp;

#ifdef DEBUG3
    printf("Enter ExtractPolygons\n");
#endif /* DEBUG3 */

    TestBooleanCtrlBrk();
    
    /* Stage 1. - mark all polygons as needed: */
    PlHead = Pl = PObj -> U.Pl;
    /* Gen. a copy of Pl, so we can modify the original polygon list: */
    PObj -> U.Pl = CopyPolygonList(Pl);

    while (Pl != NULL) {
	NumOfPolys++;	     /* Count number of polygons in original object. */
	if (Pl -> PAux != NULL) {     /* The intersection list is not empty! */
	    RST_COMPLETE_POLY(Pl);	 /* Mark it as non complete polygon. */
	    V = Pl -> V;
	    do {
		SET_ORIGINAL_VRTX(V); /* Mark vertices from original object. */
		V = V -> Pnext;
	    }
	    while (V != Pl -> V);
	}
	else {
	    SET_COMPLETE_POLY(Pl);	     /* Mark it as complete polygon. */
	    RST_INOUTPUT_POLY(Pl);	 /* And as undefined (if in output). */
	    RST_ADJPUSHED_POLY(Pl);	 /* And as not pushed on adj. stack. */
	}
	Pl = Pl -> Pnext;
    }

    /* Stage 2. - scan the polygons and subdivide the intersecting ones: */
    Pl = PlHead;

    /* We will save pointers to ALL polygons in list so we could free in the */
    /* end the ones that are not in the output list, and therefore unused.   */
    PolysPtr = (PolygonStruct **)
		MyMalloc(sizeof(PolygonStruct *) * NumOfPolys, OTHER_TYPE);
    NumOfPolys = 0;

    while (Pl != NULL) {
	TestBooleanCtrlBrk();

	PolysPtr[NumOfPolys++] = Pl;	    /* Save pointer to this polygon. */

	PlNext = Pl -> Pnext;

	if (!IS_COMPLETE_POLY(Pl)) {/* They are intersections with this one. */
	    /* Convert the intersecting segments into open/closed loops. */
	    LoopsFromInterList(Pl, &PClosed, &POpen);

	    if (PClosed != NULL && POpen != NULL)
		FatalError("ExtractPolygons: Polygon with both open & closed loops - not supported\n");

	    if (POpen != NULL) {
		/* Sort the Open loops into an order we can handle... */
		SortOpenInterList(Pl, &POpen);
		SplPl = NewPl = NULL;	 /* Keep the list of split polygons. */

		while (POpen != NULL) {	/* Clip the open loops from polygon: */
		    /* Note ClipOpenLoopFromPoly also frees the InterSegment */
		    /* list pointed by POpen (but not POpen itself). 	     */
		    NewPl = ClipOpenLoopFromPoly(Pl, POpen, AinB);
		    if (NewPl != NULL) {   /* If loop clipped a new polygon, */
			NewPl -> Pnext = SplPl;/* add to split polygon list. */
			SplPl = NewPl;
			/* And push adj. polygons of complete edges on stack.*/
			PushAdjOnStack(NewPl, AdjStack, &StackPointer);
		    }
		    Ptemp = POpen;
		    POpen = POpen -> Pnext;
		    MyFree((char *) Ptemp, OTHER_TYPE);
		}
		/* Last clip generated nothing (NewPl == NULL) so part that  */
		/* left from Pl (PlCopy) is IN output list! Add this poly:   */
		if (NewPl == NULL) {
		    Pl -> Pnext = SplPl; /* And chain what was left from it. */
		    SplPl = Pl;
		    /* And push adjacent polygons of complete edges on stack.*/
		    PushAdjOnStack(Pl, AdjStack, &StackPointer);
		    SET_INOUTPUT_POLY(Pl);/* So we wouldnt free that in end. */
		    RST_CONVEX_POLY(Pl);	   /* May be not convex now. */
		}
		OutputPl = ChainPolyLists(SplPl, OutputPl);
	    }
	    else { /* PClosed != NULL */
		/* Make a "cut" from the loop(s)  +-------+      +---+---+   */
		/* to boundary if possible, and   |       |      |       |   */
		/* converting Pl to a nonconvex   |  / \  |  ->  |  / \__|   */
		/* polygon, that has an edge (the |  \ /  |  ->  |  \ /  |   */
		/* cut) which is shared twice in  |       |      |       |   */
		/* the same polygon		  +-------+      +-------+   */
		if (CombineClosedLoops(Pl, PClosed, AinB)) {
		    /* If returned with TRUE - the polygon boundary is in    */
		    /* output, so add all its neighbours to adj. stack.      */
		    PushAdjOnStack(Pl, AdjStack, &StackPointer);
		}

		OutputPl = ChainPolyLists(Pl, OutputPl);
	    }
	}
	Pl = PlNext;
    }

    /* Stage 3. - handling adjacencies and propagate them in polygons:	     */
    /* Pop off the elements from the stack, and propagate them using their   */
    /* adjacencies.							     */
    while (StackPointer >= 0) {
	Pl = AdjStack[StackPointer--];			 /* Pop top element. */
	if (!IS_COMPLETE_POLY(Pl) ||	    /* Ignore non complete polygons. */
	     IS_INOUTPUT_POLY(Pl)) continue;	      /* If already handled. */

	SET_INOUTPUT_POLY(Pl);	  /* Mark this one as handled for next time. */

	V = Pl -> V;	     /* Push all adjacent ones that not handled yet. */
	do {
	    if (V -> PAdj &&
		IS_COMPLETE_POLY(V -> PAdj) &&
		!IS_INOUTPUT_POLY(V -> PAdj) &&
		!IS_ADJPUSHED_POLY(V -> PAdj)) {
		SET_ADJPUSHED_POLY(V -> PAdj);
		AdjStack[++StackPointer] = V -> PAdj;   /* Push it on stack. */
		if (StackPointer >= ADJ_STACK_SIZE)
		    FatalError("Adjacency stack overflow, object too big\n");
	    }
	    V = V -> Pnext;
	}
	while (V != Pl -> V);

	Pl -> Pnext = OutputPl;		   /* And chain it into output list. */
	OutputPl = Pl;
    }

    /* Free all polygons which are not in the output list: */
    while (--NumOfPolys >= 0) {
	if (!IS_INOUTPUT_POLY(PolysPtr[NumOfPolys])) {
	    PolysPtr[NumOfPolys] -> Pnext = NULL; /* Free only this polygon. */
	    MyFree((char *) (PolysPtr[NumOfPolys]), POLYGON_TYPE);
	}
    }
    MyFree((char *) PolysPtr, OTHER_TYPE);

    /* Another floating point kludge: a polygon may have zero length edge so */
    /* search for those and remove them - someone may die because of one...  */
    Pl = OutputPl;
    while (Pl != NULL) {
	V = Pl -> V;
	do {
	    Vnext = V -> Pnext;
	    if (PT_EQ(V -> Pt, Vnext -> Pt)) {
		/* Ahh - we got you. Simply skip Vnext vertex and free it: */
		V -> Pnext = Vnext -> Pnext;
		/* Update polygon vertex pointer if point on freed vertex: */
		if (Pl -> V == Vnext) Pl -> V = Vnext -> Pnext;
		Vnext -> Pnext = NULL;		      /* Dont free too much! */
		MyFree((char *) Vnext, VERTEX_TYPE);
		Vnext = V -> Pnext;
	    }
	    V = Vnext;
	}
	while (V != Pl -> V && V != NULL);

	Pl = Pl -> Pnext;
    }

#ifdef DEBUG3
    printf("Exit ExtractPolygons\n");
#endif /* DEBUG3 */

    return GenGeomObject("", OutputPl, NULL);	 /* Return resulting object. */
}

#ifdef DEBUG2

/*****************************************************************************
*   Print the content of the given polygon, to standard output.		     *
*****************************************************************************/
static void PrintVrtxList(VertexStruct *V)
{
    struct VertexStruct *VHead = V;

    do {
	printf("    %12lf %12lf %12lf", V -> Pt[0], V -> Pt[1], V -> Pt[2]);
	if (IS_INTERNAL_EDGE(V))
	     printf(" (Internal)\n");
	else printf("\n");
	V = V -> Pnext;
    }
    while (V!= NULL && V != VHead);
}

/*****************************************************************************
*   Print the content of the given InterSegment list, to standard output.    *
*****************************************************************************/
static void PrintInterList(InterSegmentStruct *PInt)
{
    printf("INTER SEGMENT LIST:\n");

    if (PInt) printf("Entry vertex pointer %p\n", PInt -> V[0]);
    while (PInt) {
	printf("%9lg %9lg %9lg (%04x), %9lg %9lg %9lg (%04x)\n",
	    PInt->PtSeg[0][0], PInt->PtSeg[0][1], PInt->PtSeg[0][2],
	    FP_SEG(PInt->V[0]),
	    PInt->PtSeg[1][0], PInt->PtSeg[1][1], PInt->PtSeg[1][2],
	    FP_SEG(PInt->V[1]));
	if (PInt -> Pnext == NULL)
	    printf("Exit vertex pointer %p\n", PInt -> V[1]);

	PInt = PInt -> Pnext;
    }
}

#endif /* DEBUG2 */