/*****************************************************************************
*   "Irit" - the 3d (not only 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 IPPolygonStruct def.)!    *
*   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 DEBUG2	         If defined, defines some printing routines. */
/* #define DEBUG3	         Print messages to entry/exit from routines. */

#include <stdio.h>
#include <ctype.h>
#include <math.h>
#include <string.h>
#include "irit_sm.h"
#include "allocate.h"
#include "bool_loc.h"
#include "geomat3d.h"

static int
    GlblObjsIntersects,
    GlblPolySortAxis = 0;

static IPObjectStruct *BooleanLowGenInOut(IPObjectStruct *PObj1, int InOut);
static void BooleanLowInterAll(IPObjectStruct *PObj1, IPObjectStruct *PObj2);
static void BooleanLowInterOne(IPPolygonStruct *Pl1, IPPolygonStruct *Pl2);
static InterSegmentStruct *InterSegmentPoly(IPPolygonStruct *Pl,
					    IPPolygonStruct *SegPl,
					    PointType Segment[2]);
static void SwapPointInterList(InterSegmentStruct *PSeg);
static void DeleteSegInterList(InterSegmentStruct *PSeg,
			       InterSegmentStruct **PSegList);
static InterSegmentStruct *FindMatchInterList(PointType Pt,
					      InterSegmentStruct **PSegList);
static void SortOpenReverseLoop(SortOpenStruct *PSHead);
static RealType SortOpenInsertOne(SortOpenStruct **PSHead,
				  SortOpenStruct *PSTemp,
				  PointType Pt,
				  IPVertexStruct *V,
				  IPPolygonStruct *Pl);
static IPObjectStruct *PolylineFromInterSeg(IPObjectStruct *PObj);
#if defined(ultrix) && defined(mips)
static int BooleanPrepObjectCmpr(VoidPtr VPoly1, VoidPtr VPoly2);
#else
static int BooleanPrepObjectCmpr(const VoidPtr VPoly1, const VoidPtr VPoly2);
#endif /* ultrix && mips (no const support) */
static void BooleanPrepObject(IPObjectStruct *PObj);

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

/*****************************************************************************
* DESCRIPTION:                                                               M
*   Finds the part of PObj1 which is out of PObj2:			     M
*                                                                            *
* PARAMETERS:                                                                M
*   PObj1:    First object of Boolean operation.                             M
*   PObj2:    Second object of Boolean operation.                            M
*                                                                            *
* RETURN VALUE:                                                              M
*   IPObjectStruct: Result of one out two.                                   M
*                                                                            *
* KEYWORDS:                                                                  M
*   BooleanLow1Out2, Booleans                                                M
*****************************************************************************/
IPObjectStruct *BooleanLow1Out2(IPObjectStruct *PObj1, IPObjectStruct *PObj2)
{
#ifdef DEBUG3
    printf("Enter BooleanLow1Out2\n");
#endif /* DEBUG3 */

    BoolGenAdjacencies(PObj1);

    /* Find all intersections of PObj1 polygons with PObj2 polygons. */
    GlblObjsIntersects = FALSE;
    BooleanLowInterAll(PObj1, PObj2);

    /* Generate all the polygons in PObj1 which are out of PObj2. */
    if (!GlblObjsIntersects) {
	FatalBooleanError(FTL_BOOL_NO_INTER);
    }

    return BooleanLowGenInOut(PObj1, FALSE);
}

/*****************************************************************************
* DESCRIPTION:                                                               M
*   Finds the part of PObj1 which is in of PObj2:			     M
*                                                                            *
* PARAMETERS:                                                                M
*   PObj1:    First object of Boolean operation.                             M
*   PObj2:    Second object of Boolean operation.                            M
*                                                                            *
* RETURN VALUE:                                                              M
*   IPObjectStruct: Result of one in two.                                    M
*                                                                            *
* KEYWORDS:                                                                  M
*   BooleanLow1In2, Booleans                                                 M
*****************************************************************************/
IPObjectStruct *BooleanLow1In2(IPObjectStruct *PObj1, IPObjectStruct *PObj2)
{
#ifdef DEBUG3
    printf("Enter BooleanLow1In2\n");
#endif /* DEBUG3 */

    BoolGenAdjacencies(PObj1);

    /* Find all intersections of PObj1 polygons with PObj2 polygons. */
    GlblObjsIntersects = FALSE;
    BooleanLowInterAll(PObj1, PObj2);

    /* Generate all the polygons in PObj1 which are in PObj2. */
    if (!GlblObjsIntersects) {
	FatalBooleanError(FTL_BOOL_NO_INTER);
    }
    return BooleanLowGenInOut(PObj1, TRUE);
}

/*****************************************************************************
* DESCRIPTION:                                                               *
*   Scans the InterSegmentList of each polygon and decides using Inter.      *
* list, if it is IN relative to the other object. Note that for polygons     *
* that do not intersect at all, we use the polygon adjacencies to decide     *
* if they are IN or not.						     *
*                                                                            *
* PARAMETERS:                                                                *
*   PObj:      To decide if it is in out out.                                *
*   InOut:     What are we looking for? in or out.                           *
*                                                                            *
* RETURN VALUE:                                                              *
*   IPObjectStruct:   Resulting set of polygons that are in the requested    *
*                     half space.                                            *
*****************************************************************************/
static IPObjectStruct *BooleanLowGenInOut(IPObjectStruct *PObj, int InOut)
{
    if (BoolOutputInterCurve) {
	/* Return a polyline object - extract the InterSegment list of each */
	/* polygon into open/closed polyline loops object.		    */
	return PolylineFromInterSeg(PObj);
    }
    else {
	if (BoolHandleCoplanarPoly) {
	    IPObjectStruct
		*PObjNew = BoolExtractPolygons(PObj, InOut);
	    IPPolygonStruct
		*Pl = PObjNew -> U.Pl;

	    /* Purge all polygons that were tagged being coplanar. */
	    while (Pl != NULL && IS_COPLANAR_POLY(Pl)) {
	        PObjNew -> U.Pl = PObjNew -> U.Pl -> Pnext;
	        IPFreePolygon(Pl);
		Pl = PObjNew -> U.Pl;
	    }
	    if (Pl != NULL) {
		while (Pl -> Pnext != NULL) {
		    if (IS_COPLANAR_POLY(Pl -> Pnext)) {
			IPPolygonStruct
			    *PlTemp = Pl -> Pnext;

			Pl -> Pnext = PlTemp -> Pnext;
			IPFreePolygon(PlTemp);
	            }
	            else
		        Pl = Pl -> Pnext;
		}
	    }

	    return PObjNew;
	}
	else
	    return BoolExtractPolygons(PObj, InOut);
    }
}

/*****************************************************************************
* DESCRIPTION:                                                               *
*   Creates a polyline object out of the intersection list of the polygons.  *
*                                                                            *
* PARAMETERS:                                                                *
*   PObj:     Polygon for which intersection loops are to be formed.         *
*                                                                            *
* RETURN VALUE:                                                              *
*   IPObjectStruct:  Resulting loops.                                        *
*****************************************************************************/
static IPObjectStruct *PolylineFromInterSeg(IPObjectStruct *PObj)
{
    IPObjectStruct *PObjRet;
    IPPolygonStruct *PlTemp, *PlObj,
	*PlHead = NULL;
    InterSegmentStruct *PInter, *PInterNext;
    InterSegListStruct *PClosed, *POpen, *PClosedNext, *POpenNext;
    IPVertexStruct *V;

    PlObj = PObj -> U.Pl;
    while (PlObj != NULL) {
	BoolLoopsFromInterList(PlObj, &PClosed, &POpen);
	while (POpen != NULL) {
	    /* Make one polyline from each loop in list: */
	    PInter = POpen -> PISeg;
	    POpenNext = POpen -> Pnext;

	    PlTemp = IPAllocPolygon(0, 0, NULL, NULL);
	    PlTemp -> PVertex = V = IPAllocVertex(0, 0, NULL, NULL);
	    PT_COPY(V -> Coord, PInter -> PtSeg[0]);
	    while (PInter) {
		PInterNext = PInter -> Pnext;

		V -> Pnext = IPAllocVertex(0, 0, NULL, NULL); V = V -> Pnext;
		PT_COPY(V -> Coord, PInter -> PtSeg[1]);

		IritFree((VoidPtr) PInter);
		PInter = PInterNext;
	    }
	    PlTemp -> Pnext = PlHead;
	    PlHead = PlTemp;

	    IritFree((VoidPtr) POpen);
	    POpen = POpenNext;
	}
	while (PClosed != NULL) {
	    /* Make one polyline from each loop in list: */
	    PInter = PClosed -> PISeg;
	    PClosedNext = PClosed -> Pnext;

	    PlTemp = IPAllocPolygon(0, 0, NULL, NULL);
	    PlTemp -> PVertex = V = IPAllocVertex(0, 0, NULL, NULL);
	    PT_COPY(V -> Coord, PInter -> PtSeg[0]);
	    while (PInter) {
		V -> Pnext = IPAllocVertex(0, 0, NULL, NULL); V = V -> Pnext;
		PT_COPY(V -> Coord, PInter -> PtSeg[1]);
		PInter = PInter -> Pnext;
	    }
	    PlTemp -> Pnext = PlHead;
	    PlHead = PlTemp;

	    IritFree((VoidPtr) PClosed);
	    PClosed = PClosedNext;
	}
	PlObj = PlObj -> Pnext;
    }

    PObjRet = GenPolyObject("", PlHead, NULL);
    IP_SET_POLYLINE_OBJ(PObjRet);	    /* Mark it as polyline object. */
    return PObjRet;
}

/*****************************************************************************
* DESCRIPTION:                                                               *
*   Routine to compare two polygon's bbox value for sorting purposes.        *
*                                                                            *
* PARAMETERS:                                                                *
*   VPoly1, VPoly2:  Two pointers to polygons.                               *
*                                                                            *
* RETURN VALUE:                                                              *
*   int:   >0, 0, or <0 as the relation between the two polygons.            *
*****************************************************************************/
#if defined(ultrix) && defined(mips)
static int BooleanPrepObjectCmpr(VoidPtr VPoly1, VoidPtr VPoly2)
#else
static int BooleanPrepObjectCmpr(const VoidPtr VPoly1, const VoidPtr VPoly2)
#endif /* ultrix && mips (no const support) */
{
    RealType
	Diff = (*((IPPolygonStruct **) VPoly1)) -> BBox[0][GlblPolySortAxis] -
	       (*((IPPolygonStruct **) VPoly2)) -> BBox[0][GlblPolySortAxis];

    return SIGN(Diff);
}

/*****************************************************************************
* DESCRIPTION:                                                               *
*   Routine to compute BBox for all polygons in provided object PObj. Also,  *
* the polygons are sorted in the list with according to their minimal BBox   *
* value in GlblPolySortAxis axis.					     *
*                                                                            *
* PARAMETERS:                                                                *
*   PObj:       To prepare.                                                  *
*                                                                            *
* RETURN VALUE:                                                              *
*   void                                                                     *
*****************************************************************************/
static void BooleanPrepObject(IPObjectStruct *PObj)
{
    int i,
	PolyCount = 0;
    IPPolygonStruct **PlArray, **PlTemp,
	*Pl = PObj -> U.Pl;

    /* Count how many polygons this object have and make sure all have bbox. */
    while (Pl != NULL)
    {
	PolyCount++;

	if (BoolHandleCoplanarPoly)
	    RST_COPLANAR_POLY(Pl);

	if (!IP_HAS_BBOX_POLY(Pl)) {			  /* Compute it now. */
	    RealType
		*BBoxMin = Pl -> BBox[0],
		*BBoxMax = Pl -> BBox[1];
	    IPVertexStruct
		*V = Pl -> PVertex;

	    BBoxMin[0] = BBoxMin[1] = BBoxMin[2] = INFINITY;
	    BBoxMax[0] = BBoxMax[1] = BBoxMax[2] = -INFINITY;

	    do {
	    	for (i = 0; i < 3; i++) {
	    	    if (BBoxMin[i] > V -> Coord[i])
			BBoxMin[i] = V -> Coord[i];
	    	    if (BBoxMax[i] < V -> Coord[i])
			BBoxMax[i] = V -> Coord[i];
	    	}

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

	    IP_SET_BBOX_POLY(Pl);
	}
	    
    	Pl = Pl -> Pnext;
    }

    /* Sort the polygons with according to minimal BBox in GlblPolySortAxis. */
    PlArray = (IPPolygonStruct **) IritMalloc(sizeof(IPPolygonStruct *) *
							PolyCount);
    for (Pl = PObj -> U.Pl, PlTemp = PlArray; Pl != NULL; Pl = Pl -> Pnext)
	*PlTemp++ = Pl;

    qsort(PlArray, PolyCount, sizeof(IPPolygonStruct *), BooleanPrepObjectCmpr);

    for (i = 1; i < PolyCount; i++)
	PlArray[i - 1] -> Pnext = PlArray[i];
    PlArray[PolyCount - 1] -> Pnext = NULL;
    PObj -> U.Pl = PlArray[0];
    IritFree((VoidPtr) PlArray);
}

/*****************************************************************************
* DESCRIPTION:                                                               M
*  Routine to set polygonal sorting axis.	  			     M
*                                                                            *
* PARAMETERS:                                                                M
*   PolySortAxis:    Sorting axis. Either 0(x), 1 (y), or 2 (z).             M
*                                                                            *
* RETURN VALUE:                                                              M
*   void                                                                     M
*                                                                            *
* KEYWORDS:                                                                  M
*   BoolSetPolySortAxis, Booleans                                            M
*****************************************************************************/
void BoolSetPolySortAxis(int PolySortAxis)
{
    GlblPolySortAxis = PolySortAxis;
}

/*****************************************************************************
* DESCRIPTION:                                                               *
*   Routine to find all the intersections between all PObj1 polygons with    *
* PObj2 polygons. The intersections are saved as a list of segments (struct  *
* InterSegmentStruct) in each of PObj1 polygons using the PAux pointer (see  *
* IPPolygonStruct). Note PObj2 is not modified at all, and in PObj1, only    *
* PAux of each polygon is set to the segment list, or NULL if none.	     *
*                                                                            *
* PARAMETERS:                                                                *
*   PObj1:        First object to compute intersection for.                  *
*   PObj2:        Second object to compute intersection for.                 *
*                                                                            *
* RETURN VALUE:                                                              *
*   void                                                                     *
*****************************************************************************/
static void BooleanLowInterAll(IPObjectStruct *PObj1, IPObjectStruct *PObj2)
{
    IPPolygonStruct *Pl1, *Pl2, *Pl2Start;

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

    /* Sort polygons and compute BBox for them if none exists. */
    BooleanPrepObject(PObj1);
    BooleanPrepObject(PObj2);
    Pl2Start = PObj2 -> U.Pl;
    
    Pl1 = PObj1 -> U.Pl;
    while (Pl1 != NULL) {
	Pl1 -> PAux = NULL;	   /* Empty InterSegment list to start with: */

	/* Intersect Pl1 against the Pl2 list until end of Pl2 list or       */
	/* until the minimum of Pl2 polygon is larger than Pl1 maximum.      */
	Pl2 = Pl2Start;
	while (Pl2 != NULL &&
	       Pl1 -> BBox[1][GlblPolySortAxis] >=
					Pl2 -> BBox[0][GlblPolySortAxis]) {
	    if (Pl1 -> BBox[1][0] < Pl2 -> BBox[0][0] ||
		Pl1 -> BBox[1][1] < Pl2 -> BBox[0][1] ||
		Pl1 -> BBox[1][2] < Pl2 -> BBox[0][2] ||
		Pl2 -> BBox[1][0] < Pl1 -> BBox[0][0] ||
		Pl2 -> BBox[1][1] < Pl1 -> BBox[0][1] ||
		Pl2 -> BBox[1][2] < Pl1 -> BBox[0][2]) {
		/* The Bounding boxes do not intersect, skip polygons. */
	    }
	    else
		BooleanLowInterOne(Pl1, Pl2);

	    /* If Pl2 maximum is smaller than Pl1 minimum there is no reason */
	    /* to consider this Pl2 polygon any more as it cannot intersect  */
	    /* any more polygons in the Pl1 list.                            */
	    if (Pl2 -> BBox[1][GlblPolySortAxis] <
                          Pl2 -> BBox[1][GlblPolySortAxis] &&
		Pl2Start == Pl2)
		Pl2Start = Pl2 -> Pnext;
                                                        
	    Pl2 = Pl2 -> Pnext;
	}

	if (Pl1 -> PAux != NULL)		     /* If any intersection. */
	    GlblObjsIntersects = TRUE;

	Pl1 = Pl1 -> Pnext;
    }

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

/*****************************************************************************
* DESCRIPTION:                                                               *
*   Routine to intersect polygon Pl1, with polygon Pl2. If found common      *
* intersection, that segment will be added to the InterSegmentStruct list    *
* saved in Pl1 PAux list.						     *
*   Note that as the two polygons convex, at most one line segment can       *
* result from such intersection (of two non coplanar polygons).		     *
*   Algorithm: intersect all Pl2 edges with Pl1 plane. If found that         *
* (exactly) two vertices (one segment) of Pl2 do intersect Pl1 plane then:   *
* Perform clipping of the segment against Pl1. If result is not empty, add   *
* the result segment to Pl1 InterSegmentStruct list (saved at PAux of	     *
* polygon - see IPPolygonStruct).					     *
*                                                                            *
* PARAMETERS:                                                                *
*   Pl1:       First polygon to compute intersection for.                    *
*   Pl2:       Second polygon to compute intersection for.                   *
*                                                                            *
* RETURN VALUE:                                                              *
*   void                                                                     *
*****************************************************************************/
static void BooleanLowInterOne(IPPolygonStruct *Pl1, IPPolygonStruct *Pl2)
{
    int NumOfInter = 0;
    RealType TInter[2],
	*Plane1 = Pl1 -> Plane,			       /* For faster access. */
	*Plane2 = Pl2 -> Plane;
    PointType Inter[2], Dir;
    IPVertexStruct *Vnext,
	*V = Pl2 -> PVertex;
    InterSegmentStruct *PSeg, *PLSeg;

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

    if ((BOOL_APX_EQ(Plane1[0], Plane2[0]) &&
	 BOOL_APX_EQ(Plane1[1], Plane2[1]) &&
	 BOOL_APX_EQ(Plane1[2], Plane2[2]) &&
	 BOOL_APX_EQ(Plane1[3], Plane2[3])) ||
	(BOOL_APX_EQ(Plane1[0], -Plane2[0]) &&
	 BOOL_APX_EQ(Plane1[1], -Plane2[1]) &&
	 BOOL_APX_EQ(Plane1[2], -Plane2[2]) &&
	 BOOL_APX_EQ(Plane1[3], -Plane2[3]))) {
	/* The two polygons are coplanar, do not attempt to intersect. */
	if (BoolHandleCoplanarPoly) {
	    SET_COPLANAR_POLY(Pl1);
	    SET_COPLANAR_POLY(Pl2);
	}
	else {
	    IritWarningError("Boolean: coplanar polygons detected. Enable COPLANAR state.");
	}
    }
    else {
	do {
	    Vnext = V -> Pnext;
	    PT_SUB(Dir, Vnext -> Coord, V -> Coord);
	    if (CGPointFromLinePlane01(V -> Coord, Dir, Plane1,
				 Inter[NumOfInter], &TInter[NumOfInter]) &&
	        ((NumOfInter == 0) || (NumOfInter == 1 &&
				       !BOOL_PT_APX_EQ(Inter[0], Inter[1]))))
	        NumOfInter++;
	    if (NumOfInter == 2)
		break;			   /* Cannt have more intersections. */

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

    switch (NumOfInter) {
	case 0:
	    break;
	case 1:
	    /* One intersection is possible if only one vertex of Pl2 is in  */
	    /* the plane of Pl1, all other vertices are on one side of plane.*/
	    break;
	case 2:
	    /* Clip the segment against the polygon and insert if not empty: */
	    if ((PSeg = InterSegmentPoly(Pl1, Pl2, Inter)) != NULL) {
		/* insert that segment to list of Pl1. Note however that the */
		/* intersection may be exactly on 2 other polygons boundary, */
		/* And therefore creates the same intersection edge TWICE!   */
		/* Another possiblity is on same case, the other polygon     */
		/* will have that inter. edge on its edge, and its ignored.  */
		/* We therefore test for duplicates and ignore edge if so.   */
		if (PSeg -> V[0] != NULL && PSeg -> V[0] == PSeg -> V[1]) {
		    IritFree((VoidPtr) PSeg);		       /* Ignore it! */
		    return;
		}
		PLSeg = (InterSegmentStruct *) Pl1 -> PAux;
		while (PLSeg != NULL) {
		    if ((BOOL_PT_APX_EQ(PSeg -> PtSeg[0], PLSeg -> PtSeg[0]) &&
			 BOOL_PT_APX_EQ(PSeg -> PtSeg[1], PLSeg -> PtSeg[1])) ||
			(BOOL_PT_APX_EQ(PSeg -> PtSeg[0], PLSeg -> PtSeg[1]) &&
			 BOOL_PT_APX_EQ(PSeg -> PtSeg[1], PLSeg -> PtSeg[0]))) {
			IritFree((VoidPtr) PSeg);	       /* Ignore it! */
			return;
		    }
		    PLSeg = PLSeg -> Pnext;
		}

		PSeg -> Pnext = (InterSegmentStruct *) Pl1 -> PAux;
		Pl1 -> PAux = (VoidPtr) PSeg;
	    }
	    break;
    }

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

/*****************************************************************************
* DESCRIPTION:                                                               *
*   Intersects the given segment (given as two end points), with the given   *
* polygon (which must be convex). Upto two intersections are possible, as    *
* again, the polygon is convex. Note Segment polygon is given as SegPl.      *
*                                                                            *
* PARAMETERS:                                                                *
*   Pl:           Full complete polygon.                                     *
*   SegPl:        Origin of Segment.                                         *
*   Segment:      A single linear segment as two points.                     *
*                                                                            *
* RETURN VALUE:                                                              *
*   InterSegmentStruct:      The (upto) two intersections information.       *
*****************************************************************************/
static InterSegmentStruct *InterSegmentPoly(IPPolygonStruct *Pl,
					    IPPolygonStruct *SegPl,
					    PointType Segment[2])
{
    int i, Reverse, Res,
	NumOfInter = 0;
    RealType TInter[3], temp, Min, Max, *PtSeg0, *PtSeg1;
    PointType Dir, Inter[2], SegDir, Pt1, Pt2;
    IPVertexStruct *VInter[2], *Vnext,
	*V = Pl -> PVertex;
    InterSegmentStruct *PSeg;

    /* Find the segment direction vector: */
    PT_SUB(SegDir, Segment[1], Segment[0]);

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

    do {
	Vnext = V -> Pnext;
	PT_SUB(Dir, Vnext -> Coord, V -> Coord);
	/* Find the intersection of the segment with all the polygon edges: */
	/* Note the t parameter value of the edge (temp) must be in [0..1]. */
	if ((Res = CG2PointsFromLineLine(Segment[0], SegDir,
		V -> Coord, Dir, Pt1, &TInter[NumOfInter], Pt2, &temp)) != 0 &&
	    (temp > 0.0 || BOOL_APX_EQ(temp, 0.0)) &&
	    (temp < 1.0 || BOOL_APX_EQ(temp, 1.0)) &&
	    BOOL_PT_APX_EQ(Pt1, Pt2) &&
	    (NumOfInter == 0 ||
	     (NumOfInter == 1 && !BOOL_APX_EQ(TInter[0], TInter[1])))) {
	    PT_COPY(Inter[NumOfInter], Pt1);
	    VInter[NumOfInter++] = V;/* Save pointer to intersected edge. */
	}
	else {
	    /* If Res == 0 its parallel to edge. If distance is zero then    */
	    /* line is on the edge line itself - quit from this one!         */
	    if (!Res && CGDistPointLine(Segment[0], V -> Coord, Dir) <
								BOOL_EPSILON) {
		/* Wipe out adjacency of this vertex if not shared any more  */
		IPVertexStruct *Vtemp;

		if (V -> PAdj == NULL || BoolVerifyShareEdge(V, V -> PAdj))
		    return NULL;

		Vtemp = V -> PAdj -> PVertex;
		do {/* Find the edge on the other polygon to wipe out first. */
		    if (Vtemp -> PAdj == Pl) {
			Vtemp -> PAdj = NULL;
			break;
		    }
		    Vtemp = Vtemp -> Pnext;
		}
		while (Vtemp != NULL && Vtemp != V -> PAdj -> PVertex);
		V -> PAdj = NULL;		/* And wipe out ours also... */
		return NULL;
	    }
	}

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

    switch (NumOfInter) {
	case 0:
	    return NULL;
	case 1:
	    /* One intersection is possible if segment intersects one vertex */
	    /* of polygon and all other vertices are on same side of segment.*/
	    return NULL;
	case 2:
	    /* In order the segment to really intersect the polygon, it must */
	    /* at list part of t in [0..1] in the polygon. Test it:	     */
	    Min = MIN(TInter[0], TInter[1]);
	    Max = MAX(TInter[0], TInter[1]);
	    Reverse = TInter[0] > TInter[1];
	    if (Min >= 1.0 || BOOL_APX_EQ(Min, 1.0) ||
		Max <= 0.0 || BOOL_APX_EQ(Max, 0.0))
		return NULL;

	    PSeg = (InterSegmentStruct *)
		IritMalloc(sizeof(InterSegmentStruct));
	    PSeg -> Pl = SegPl;     /* Pointer to other (intersect) polygon. */

	    /* Handle the Min end point: */
	    if (BOOL_APX_EQ(Min, 0.0)) {
		PtSeg0 = Segment[0];
		PSeg -> V[0] = (Reverse ? VInter[1] : VInter[0]);
	    }
	    else if (Min < 0.0) {
		PtSeg0 = Segment[0];
		PSeg -> V[0] = NULL;			 /* End is internal. */
	    }
	    else { /* Min > 0.0 */
		PtSeg0 = (Reverse ? Inter[1] : Inter[0]);
		PSeg -> V[0] = (Reverse ? VInter[1] : VInter[0]);
	    }

	    /* Handle the Max end point: */
	    if (BOOL_APX_EQ(Max, 1.0)) {
		PtSeg1 = Segment[1];
		PSeg -> V[1] = (Reverse ? VInter[0] : VInter[1]);
	    }
	    else if (Max > 1.0) {
		PtSeg1 = Segment[1];
		PSeg -> V[1] = NULL;			 /* End is internal. */
	    }
	    else { /* Max < 1.0 */
		PtSeg1 = (Reverse ? Inter[0] : Inter[1]);
		PSeg -> V[1] = (Reverse ? VInter[0] : VInter[1]);
	    }

	    PT_COPY(PSeg -> PtSeg[0], PtSeg0); /* The two segment end point. */
            PT_COPY(PSeg -> PtSeg[1], PtSeg1);

	    for (i = 0; i < 3; i++) {		 /* Make zeros look nicer... */
		if (ABS(PSeg -> PtSeg[0][i]) < BOOL_EPSILON)
		    PSeg -> PtSeg[0][i] = 0.0;
		if (ABS(PSeg -> PtSeg[1][i]) < BOOL_EPSILON)
		    PSeg -> PtSeg[1][i] = 0.0;
	    }
	    if (BOOL_PT_APX_EQ(PSeg -> PtSeg[0], PSeg -> PtSeg[1])) {
		IritFree((VoidPtr) PSeg);
		return NULL;
	    }
	    return PSeg;
    }
    return NULL;				    /* Makes warning silent. */
}

/*****************************************************************************
* DESCRIPTION:                                                               M
*   Given a polygon with the intersection list, creates the polylines        M
* loop(s) out of it, which can be one of the two:			     M
* 1. Closed loop - all the intersections create a loop in one polygon.	     M
* 2. Open polyline - if the intersections cross the polygon boundary. In     M
*    this case the two end point of the polyline, must lay on polygon	     M
*    boundary.								     M
*									     M
* In both cases, the polyline will be as follows:			     M
* First point at first list element at PtSeg[0] (see InterSegmentStruct).    M
* Second point at first list element at PtSeg[1] (see InterSegmentStruct).   M
* Point i at list element (i-1) at PtSeg[0] (PtSeg[1] is not used!).         M
* In the closed loop case the last point is equal to first.		     M
* Both cases returns NULL terminated list.				     M
*                                                                            *
* PARAMETERS:                                                                M
*   Pl:         Polygon with intersection information in its PAux slot.      M
*   PClosed:    To be updated with the closed loops found in Pl.             M
*   POpen:      To be updated with the open loops found in Pl.               M
*                                                                            *
* RETURN VALUE:                                                              M
*   void                                                                     M
*                                                                            *
* KEYWORDS:                                                                  M
*   BoolLoopsFromInterList                                                   M
*****************************************************************************/
void BoolLoopsFromInterList(IPPolygonStruct *Pl,
			    InterSegListStruct **PClosed,
			    InterSegListStruct **POpen)
{
    InterSegmentStruct *PSeg, *PSegHead, *PSegTemp, *PSegNewTemp;
    InterSegListStruct *PSLTemp;

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

    *PClosed = (*POpen) = NULL;

    if ((PSeg = (InterSegmentStruct *) Pl -> PAux) == NULL) {
	return;
    }
    else {
	/* Do intersect - find if there are closed loops and/or open ones:   */
	Pl -> PAux = NULL;
	while (TRUE) {
	    /* First, we look for open loops - scans linearly (maybe should  */
	    /* be done more efficiently) for segment on boundary and start   */
	    /* build chain from it (up to the other end, on the boundary).   */
	    PSegHead = PSeg;
	    while (PSegHead) {
		if (PSegHead -> V[0] != NULL || PSegHead -> V[1] != NULL) {
		    /* Found one - make it in correct order, del. from list: */
		    if (PSegHead -> V[0] == NULL)
			SwapPointInterList(PSegHead);
		    DeleteSegInterList(PSegHead, &PSeg);
		    break;
		}
		else
		    PSegHead = PSegHead -> Pnext;
	    }
	    if (PSegHead == NULL)
		break;			    /* No more open segments here... */

	    PSegTemp = PSegHead;
	    while (PSegTemp -> V[1] == NULL) {
		/* Search for matching to the second boundary end: */
		PSegNewTemp = FindMatchInterList(PSegTemp -> PtSeg[1], &PSeg);
		PSegTemp -> Pnext = PSegNewTemp;
		PSegTemp = PSegNewTemp;
	    }
	    PSegTemp -> Pnext = NULL;
	    PSLTemp = (InterSegListStruct *)
		IritMalloc(sizeof(InterSegListStruct));
	    PSLTemp -> PISeg = PSegHead;
	    PSLTemp -> PISegMaxX = NULL;
	    PSLTemp -> Pnext = *POpen;
	    *POpen = PSLTemp;
	}

	while (TRUE) {
	    /* Now, we look for closed loops - pick one segment and search   */
	    /* for matching until you close the loop.			     */
	    PSegHead = PSeg;
	    if (PSegHead == NULL)
		break;			  /* No more closed segments here... */
	    DeleteSegInterList(PSegHead, &PSeg);

	    PSegTemp = PSegHead;
	    while (!BOOL_PT_APX_EQ(PSegTemp -> PtSeg[1],
				   PSegHead -> PtSeg[0])) {
		/* Search for matching until we back at first point: */
		PSegNewTemp = FindMatchInterList(PSegTemp -> PtSeg[1], &PSeg);
		PSegTemp -> Pnext = PSegNewTemp;
		PSegTemp = PSegNewTemp;
	    }
	    PSegTemp -> Pnext = NULL;
	    PSLTemp = (InterSegListStruct *)
		IritMalloc(sizeof(InterSegListStruct));
	    PSLTemp -> PISeg = PSegHead;
	    PSLTemp -> PISegMaxX = NULL;
	    PSLTemp -> Pnext = *PClosed;
	    *PClosed = PSLTemp;
	}
    }

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

/*****************************************************************************
* DESCRIPTION:                                                               *
* Swaps the two points in an InterSegmentStruct (modifies PtSeg & V entries) *
*                                                                            *
* PARAMETERS:                                                                *
*   PSeg:   To swap the two points in it.                                    *
*                                                                            *
* RETURN VALUE:                                                              *
*   void                                                                     *
*****************************************************************************/
static void SwapPointInterList(InterSegmentStruct *PSeg)
{
    PointType Pt;
    IPVertexStruct *V;

    PT_COPY(Pt,		      PSeg -> PtSeg[0]);
    PT_COPY(PSeg -> PtSeg[0], PSeg -> PtSeg[1]);
    PT_COPY(PSeg -> PtSeg[1], Pt);

    V		 = PSeg -> V[0];
    PSeg -> V[0] = PSeg -> V[1];
    PSeg -> V[1] = V;
}

/*****************************************************************************
* DESCRIPTION:                                                               *
*   Deletes one InterSegment PSeg, from InterSegmentList PSegList.	     *
*                                                                            *
* PARAMETERS:                                                                *
*   PSeg:         To be deleted from the list (not delete itself).           *
*   PSegList:     List containing PSeg. Updated in place.                    *
*                                                                            *
* RETURN VALUE:                                                              *
*   void                                                                     *
*****************************************************************************/
static void DeleteSegInterList(InterSegmentStruct *PSeg,
			       InterSegmentStruct **PSegList)
{
    InterSegmentStruct *PSegTemp;

    if (*PSegList == PSeg) { /* Its the first one - simply advance list ptr: */
	*PSegList = (*PSegList) -> Pnext;
    }
    else {
	PSegTemp = (*PSegList);
	while (PSegTemp -> Pnext != NULL && PSegTemp -> Pnext != PSeg)
	    PSegTemp = PSegTemp -> Pnext;
	if (PSegTemp -> Pnext == PSeg)
	    PSegTemp -> Pnext = PSegTemp -> Pnext -> Pnext;
	else
	    IritFatalError("DeleteSegInterList: element to delete not found");
    }
}

/*****************************************************************************
* DESCRIPTION:                                                               *
*   Searches for matching point, in the given inter segment list. Returns    *
* the pointer to that element after swapping its points if needed (the match *
* must be with point 0 of new segment returned), and deletes it from list.   *
*                                                                            *
* PARAMETERS:                                                                *
*   Pt:          To find a match for.                                        *
*   PSegList:    To search for a match for Pt at.                            *
*                                                                            *
* RETURN VALUE:                                                              *
*   InterSegmentStruct:   Matching edge with mtching point to Pt, if any.    *
*****************************************************************************/
static InterSegmentStruct *FindMatchInterList(PointType Pt,
					      InterSegmentStruct **PSegList)
{
    InterSegmentStruct
	*PSegTemp = (*PSegList);

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

    while (PSegTemp != NULL) {
	if (BOOL_PT_APX_EQ(Pt, PSegTemp -> PtSeg[0])) {
	    DeleteSegInterList(PSegTemp, PSegList);
	    return PSegTemp;
	}
	if (BOOL_PT_APX_EQ(Pt, PSegTemp -> PtSeg[1])) {
	    DeleteSegInterList(PSegTemp, PSegList);
	    SwapPointInterList(PSegTemp);
	    return PSegTemp;
	}
	PSegTemp = PSegTemp -> Pnext;
    }

    IritFatalError("FindMatchInterList: No match found - Empty Object Result");
    return NULL;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
*   Sorts the open loops of given polygon to an order that can be used in    M
* subdividing into sub polygons later (see comment of BoolExtractPolygons).  M
*   This order is such that each loops will have no other loop between its   M
* end points, if we walk along the polygon in the (linked list direction)    M
* perimeter from one end to the other, before it. For example:		     M
*	         -----------------------------L3			     V
*		|  ---------------L1  -----L2 |          --------L4   --L5   V
*		| |               |  |     |  |         |        |   |  |    V
*	  P0 ------ P1 ------- P2 ----- P3 -------- P4 ------ P5 -------- P0 V
* In this case, any order such that L1, L2 are before L3 will do. Obviously  M
* this is not a total order, and they are few correct ways to sort it.	     M
* Algorithm:								     M
*   For each open loop, for each of its two end, evaluate a RealType key for M
* the end point P between segment P(i) .. P(i+1) to be i + t, where:	     M
* t is the ratio  (P - P(i)) / (P(i+1) - P(i)) . This maps the all perimeter M
* of the polygon onto 0..N-1, where N is number of vertices of that polygon. M
*   Sort the keys, and while they are keys in data sturcture, search and     M
* remove a consecutive pair of keys assosiated with same loop, and output it.M
*   Note that each open loop point sequence is tested to be such that it     M
* starts on the first point (first and second along vertex list) on polygon  M
* perimeter, and the sequence end is on the second point, and the sequence   M
* is reversed if not so. This order will make the replacement of the         M
* perimeter from first to second points by the open loop much easier.	     M
*   This may be real problem if there are two intersection points almost     M
* identical - floating point errors may cause it to loop forever. We use     M
* some reordering heuristics in this case, and return fatal error if fail!   M
*                                                                            *
* PARAMETERS:                                                                M
*   Pl:       To sort the loops for.                                         M
*   POpen:    The set of open loops. Updated in place.                       M
*                                                                            *
* RETURN VALUE:                                                              M
*   void                                                                     M
*                                                                            *
* KEYWORDS:                                                                  M
*   BoolSortOpenInterList                                                    M
*****************************************************************************/
void BoolSortOpenInterList(IPPolygonStruct *Pl, InterSegListStruct **POpen)
{
    int Found = TRUE,
	ReOrderFirst = FALSE,
	NumOfReOrders = 0;
    RealType Key1, Key2;
    InterSegmentStruct *PSeg;
    InterSegListStruct *PLSeg, *PLNext,
	*PResTemp = NULL,
	*PResHead = NULL;
    SortOpenStruct *PSTemp, *PSTemp1, *PSTemp2,
	*PSHead = NULL;

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

    PLSeg = (*POpen);
    while (PLSeg != NULL) {    /* Evaluate the two end keys and insert them: */
        PSeg = PLSeg -> PISeg;
	PLNext = PLSeg -> Pnext;
	PLSeg -> Pnext = NULL;

	PSTemp1 = (SortOpenStruct *) IritMalloc(sizeof(SortOpenStruct));
	PSTemp1 -> PLSeg = PLSeg;
	/* Insert PSTemp1 into PLHead list according to position of Pt in Pl.*/
	Key1 = SortOpenInsertOne(&PSHead, PSTemp1, PSeg -> PtSeg[0],
						   PSeg -> V[0], Pl);

	while (PSeg -> Pnext != NULL) PSeg = PSeg -> Pnext; /* Go other end. */
	PSTemp2 = (SortOpenStruct *) IritMalloc(sizeof(SortOpenStruct));
	PSTemp2 -> PLSeg = PLSeg;
	/* Insert PSTemp2 into PLHead list according to position of Pt in Pl.*/
	Key2 = SortOpenInsertOne(&PSHead, PSTemp2, PSeg -> PtSeg[1],
						   PSeg -> V[1], Pl);

	if (Key1 > Key2)	      /* Reverse the open loop points order: */
	    SortOpenReverseLoop(PSTemp1);

	PLSeg = PLNext;
    }

    while (PSHead != NULL) {	/* Search for consecutive pair of same loop. */
	if (NumOfReOrders++ > 10)
	    IritFatalError("BoolSortOpenInterList: fail to sort intersection list");
	if (Found)
	    NumOfReOrders = 0;

	Found = FALSE;
	PSTemp = PSHead;
	if (PSTemp -> PLSeg == PSTemp -> Pnext -> PLSeg) { /* First is pair! */
	    if (PResHead == NULL)
		PResHead = PResTemp = PSTemp -> PLSeg;
	    else {
		PResTemp -> Pnext = PSTemp -> PLSeg;
		PResTemp = PSTemp -> PLSeg;
	    }
	    PSHead = PSHead -> Pnext -> Pnext;	     /* Skip the first pair. */
	    IritFree((VoidPtr) PSTemp -> Pnext);
	    IritFree((VoidPtr) PSTemp);
	    Found = TRUE;
	    continue;
	}
	/* If we are here, first pair is not of same loop - search on: */
	while (PSTemp -> Pnext -> Pnext != NULL) {
	    if (PSTemp -> Pnext -> PLSeg == PSTemp -> Pnext -> Pnext -> PLSeg) {
		if (PResHead == NULL)
		    PResHead = PResTemp = PSTemp -> Pnext -> PLSeg;
		else {
		    PResTemp -> Pnext = PSTemp -> Pnext -> PLSeg;
		    PResTemp = PSTemp -> Pnext -> PLSeg;
		}
		PSTemp2 = PSTemp -> Pnext;
		PSTemp -> Pnext = PSTemp -> Pnext -> Pnext -> Pnext;
		IritFree((VoidPtr) PSTemp2 -> Pnext);
		IritFree((VoidPtr) PSTemp2);
		Found = TRUE;
		break;
	    }
	    PSTemp = PSTemp -> Pnext;
	}
	/* The only way we might found nothing is in floating point round */
	/* off error - two curve ends has almost the same Key...	  */
	/* Note, obviously, that there are at list 4 entries in list.     */
	if (!Found) {
	    if (!ReOrderFirst &&
		BOOL_APX_EQ(PSHead -> Pnext -> Key, PSHead -> Key)) {
		ReOrderFirst = TRUE;
		PSTemp1 = PSHead -> Pnext;
		PSHead -> Pnext = PSTemp1 -> Pnext;
		PSTemp1 -> Pnext = PSHead;
		PSHead = PSTemp1;
		continue;
	    }
	    else
		ReOrderFirst = FALSE;
	    PSTemp = PSHead;
	    while (PSTemp -> Pnext -> Pnext != NULL) {
		if (BOOL_APX_EQ(PSTemp -> Pnext -> Key,
			   PSTemp -> Pnext -> Pnext -> Key)) {
		    PSTemp1 = PSTemp -> Pnext;
		    PSTemp2 = PSTemp1 -> Pnext;
		    PSTemp1 -> Pnext = PSTemp2 -> Pnext;
		    PSTemp2 -> Pnext = PSTemp1;
		    PSTemp -> Pnext = PSTemp2;
		    break;
		}
		PSTemp = PSTemp -> Pnext;
	    }
	}
    }

    *POpen = PResHead;

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

/*****************************************************************************
* DESCRIPTION:                                                               *
*   Reverses the order of the open loop pointed by PSHead.		     *
*                                                                            *
* PARAMETERS:                                                                *
*   PSHead:    Head of open loop to be reversed, in place.                   *
*                                                                            *
* RETURN VALUE:                                                              *
*   void                                                                     *
*****************************************************************************/
static void SortOpenReverseLoop(SortOpenStruct *PSHead)
{
    InterSegmentStruct *PIHead, *PITemp,
	*PINewHead = NULL;

    PIHead = PSHead -> PLSeg -> PISeg;

    while (PIHead != NULL) {
	PITemp = PIHead;
	PIHead = PIHead -> Pnext;
	SwapPointInterList(PITemp);
	PITemp -> Pnext = PINewHead;
	PINewHead = PITemp;
    }

    PSHead -> PLSeg -> PISeg = PINewHead;
}

/*****************************************************************************
* DESCRIPTION:                                                               *
*   Inserts new loop PSTemp, defines key by Pt and V (V defines the edge     *
* and Pt is the points on it), into (decreasingly) ordered list PSHead.	     *
*                                                                            *
* PARAMETERS:                                                                *
*   PSHead:      Ordered list of loops.                                      *
*   PSTemp:      Loop to be inserted.                                        *
*   Pt, V:       Defining the order's key.                                   *
*   Pl:          The polygon that is responsible for all the mess.           *
*                                                                            *
* RETURN VALUE:                                                              *
*   RealType:    The key that the new loop was inserted according to.        *
*****************************************************************************/
static RealType SortOpenInsertOne(SortOpenStruct **PSHead,
				  SortOpenStruct *PSTemp,
				  PointType Pt,
				  IPVertexStruct *V,
				  IPPolygonStruct *Pl)
{
    int i = 0;
    RealType Key;
    PointType Pt1, Pt2;
    IPVertexStruct
	*VTemp = Pl -> PVertex;
    SortOpenStruct *PSTail;

    PSTemp -> Pnext = NULL;

    while (VTemp != V && VTemp != NULL) {
	i++;
	VTemp = VTemp -> Pnext;
    }
    if (VTemp == NULL)
	IritFatalError("SortOpenInsertOne: fail to find vertex");

    PT_SUB(Pt1, V -> Pnext -> Coord, V -> Coord);
    PT_SUB(Pt2, Pt, V -> Coord);
    Key = PSTemp -> Key = i + PT_LENGTH(Pt2) / PT_LENGTH(Pt1);

    /* Now insert PSTemp into the ordered list: */
    if (*PSHead == NULL) {
	*PSHead = PSTemp;
	return Key;
    }
    if (PSTemp -> Key > (*PSHead) -> Key) {		 /* Insert as first? */
	PSTemp -> Pnext = (*PSHead);
	*PSHead = PSTemp;
	return Key;
    }
    PSTail = (*PSHead);
    while (PSTail -> Pnext != NULL && PSTemp -> Key < PSTail -> Pnext -> Key)
	PSTail = PSTail -> Pnext;
    PSTemp -> Pnext = PSTail -> Pnext;
    PSTail -> Pnext = PSTemp;

    return Key;
}

#ifdef DEBUG2

/*****************************************************************************
* DESCRIPTION:                                                               *
*   Print the content of the given vertex list, to standard output.	     *
*                                                                            *
* PARAMETERS:                                                                *
*   V:       The vertex list to be printed.                                  *
*                                                                            *
* RETURN VALUE:                                                              *
*   void                                                                     *
*****************************************************************************/
static void PrintVrtxList(IPVertexStruct *V)
{
    IPVertexStruct
	*VHead = V;

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

/*****************************************************************************
* DESCRIPTION:                                                               *
*   Print the content of the given InterSegment list, to standard output.    *
*                                                                            *
* PARAMETERS:                                                                *
*   PInt:       The InterSegment list to be printed.                         *
*                                                                            *
* RETURN VALUE:                                                              *
*   void                                                                     *
*****************************************************************************/
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],
	    PInt->V[0],
	    PInt->PtSeg[1][0], PInt->PtSeg[1][1], PInt->PtSeg[1][2],
	    PInt->V[1]);
	if (PInt -> Pnext == NULL)
	    printf("Exit vertex pointer %p\n", PInt -> V[1]);

	PInt = PInt -> Pnext;
    }
}

#endif /* DEBUG2 */
