/*****************************************************************************
*   "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 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 "adjacncy.h"
#include "allocate.h"
#include "booleang.h"
#include "booleanl.h"
#include "geomat3d.h"
#include "objects.h"

static int ObjectsIntersect;

static ObjectStruct *BooleanLowGenInOut(ObjectStruct *PObj1, int InOut);
static void BooleanLowInterAll(ObjectStruct *PObj1, ObjectStruct *PObj2);
static void BooleanLowInterOne(PolygonStruct *Pl1, PolygonStruct *Pl2);
static InterSegmentStruct *InterSegmentPoly(PolygonStruct *Pl,
				PolygonStruct *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, VertexStruct *V, PolygonStruct *Pl);
static ObjectStruct *PolylineFromInterSeg(ObjectStruct *PObj);

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

/*****************************************************************************
*   Routine to find the part of PObj1 which is out of PObj2:		     *
*****************************************************************************/
ObjectStruct *BooleanLow1Out2(ObjectStruct *PObj1, ObjectStruct *PObj2)
{
#ifdef DEBUG3
    printf("Enter BooleanLow1Out2\n");
#endif /* DEBUG3 */

    GenAdjacencies(PObj1);

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

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

/*****************************************************************************
*   Routine to find the part of PObj1 which is in PObj2:		     *
*****************************************************************************/
ObjectStruct *BooleanLow1In2(ObjectStruct *PObj1, ObjectStruct *PObj2)
{
#ifdef DEBUG3
    printf("Enter BooleanLow1In2\n");
#endif /* DEBUG3 */

    GenAdjacencies(PObj1);

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

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

/*****************************************************************************
*   Scan the InterSegmentList of each polygon and decide using Intersection  *
* list, if it is IN relative to the other object. Note that for polygons     *
* that does not intersect at all, we use the polygon adjacencies to decide   *
* if they are IN or not.						     *
*****************************************************************************/
static ObjectStruct *BooleanLowGenInOut(ObjectStruct *PObj, int InOut)
{
    if (BooleanOutputInterCurve) {
	/* Return a polyline object - extract the InterSegment list of each  */
	/* polygon into open/closed polyline loops object.		     */
	return PolylineFromInterSeg(PObj);
    }
    else
	return ExtractPolygons(PObj, InOut);
}

/*****************************************************************************
*   Create a polyline object out of the intersection list of the polygons.   *
*****************************************************************************/
static ObjectStruct *PolylineFromInterSeg(ObjectStruct *PObj)
{
    ObjectStruct *PObjRet;
    PolygonStruct *PlTemp, *PlHead = NULL, *PlObj;
    InterSegmentStruct *PInter, *PInterNext;
    InterSegListStruct *PClosed, *POpen, *PClosedNext, *POpenNext;
    VertexStruct *V;

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

	    PlTemp = AllocPolygon(0, 0, NULL, NULL);
	    PlTemp -> V = V = AllocVertex(0, 0, NULL, NULL);
	    PT_COPY(V -> Pt, PInter -> PtSeg[0]);
	    while (PInter) {
		PInterNext = PInter -> Pnext;

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

		MyFree((char *) PInter, ALLOC_OTHER);
		PInter = PInterNext;
	    }
	    PlTemp -> Pnext = PlHead;
	    PlHead = PlTemp;

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

	    PlTemp = AllocPolygon(0, 0, NULL, NULL);
	    PlTemp -> V = V = AllocVertex(0, 0, NULL, NULL);
	    PT_COPY(V -> Pt, PInter -> PtSeg[0]);
	    while (PInter) {
		V -> Pnext = AllocVertex(0, 0, NULL, NULL); V = V -> Pnext;
		PT_COPY(V -> Pt, PInter -> PtSeg[1]);
		PInter = PInter -> Pnext;
	    }
	    PlTemp -> Pnext = PlHead;
	    PlHead = PlTemp;

	    MyFree((char *) PClosed, ALLOC_OTHER);
	    PClosed = PClosedNext;
	}
	PlObj = PlObj -> Pnext;
    }

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

/*****************************************************************************
*   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  *
* PolygonStruct). 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		     *
*****************************************************************************/
static void BooleanLowInterAll(ObjectStruct *PObj1, ObjectStruct *PObj2)
{
    PolygonStruct *Pl1, *Pl2;

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

    Pl1 = PObj1 -> U.Pl.P;
    while (Pl1 != NULL) {
	Pl1 -> PAux = NULL;	   /* Empty InterSegment list to start with: */

	Pl2 = PObj2 -> U.Pl.P;
	while (Pl2 != NULL) {
	    BooleanLowInterOne(Pl1, Pl2);
	    Pl2 = Pl2 -> Pnext;
	}
	ObjectsIntersect |= (Pl1 -> PAux != NULL);   /* If any intersection. */

	Pl1 = Pl1 -> Pnext;
    }

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

/*****************************************************************************
*   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, only one segment can result from   *
* such intersection.							     *
*   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 PolygonStruct).						     *
*****************************************************************************/
static void BooleanLowInterOne(PolygonStruct *Pl1, PolygonStruct *Pl2)
{
    int NumOfInter = 0;
    RealType TInter[2],
	*Plane = Pl1 -> Plane;			       /* For faster access. */
    PointType Inter[2], Dir;
    VertexStruct *Vnext, *V = Pl2 -> V;
    InterSegmentStruct *PSeg, *PLSeg;

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

    TestBooleanCtrlBrk();

    do {
	Vnext = V -> Pnext;
	PT_SUB(Dir, Vnext -> Pt, V -> Pt);
	if (CGPointFromLinePlane01(V -> Pt, Dir, Plane, Inter[NumOfInter],
							&TInter[NumOfInter]) &&
	    ((NumOfInter == 0) || (NumOfInter == 1 &&
				   !PT_EQ(Inter[0], Inter[1]))))
	    NumOfInter++;
	if (NumOfInter == 2) break;	   /* Cannt have more intersections. */

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

    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]) {
		    MyFree((char *) PSeg, ALLOC_OTHER);	       /* Ignore it! */
		    return;
		}
		PLSeg = (InterSegmentStruct *) Pl1 -> PAux;
		while (PLSeg != NULL) {
		    if ((PT_EQ(PSeg -> PtSeg[0], PLSeg -> PtSeg[0]) &&
			 PT_EQ(PSeg -> PtSeg[1], PLSeg -> PtSeg[1])) ||
			(PT_EQ(PSeg -> PtSeg[0], PLSeg -> PtSeg[1]) &&
			 PT_EQ(PSeg -> PtSeg[1], PLSeg -> PtSeg[0]))) {
			MyFree((char *) PSeg, ALLOC_OTHER);     /* Ignore it! */
			return;
		    }
		    PLSeg = PLSeg -> Pnext;
		}

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

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

/*****************************************************************************
*   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.      *
*****************************************************************************/
static InterSegmentStruct *InterSegmentPoly(PolygonStruct *Pl,
				PolygonStruct *SegPl, PointType Segment[2])
{
    int i, NumOfInter = 0, Reverse, Res;
    RealType TInter[2], temp, Min, Max, *PtSeg0, *PtSeg1;
    PointType Dir, Inter[2], SegDir, Pt1, Pt2;
    VertexStruct *VInter[2], *V = Pl -> V, *Vnext;
    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 -> Pt, V -> Pt);
	/* 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 -> Pt, Dir, Pt1, &TInter[NumOfInter], Pt2, &temp)) != 0 &&
	    (temp > 0.0 || APX_EQ(temp, 0.0)) &&
	    (temp < 1.0 || APX_EQ(temp, 1.0)) && PT_EQ(Pt1, Pt2) &&
	    (NumOfInter == 0 ||
	     (NumOfInter == 1 && !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 && CGDistPointPoint(Pt1, Pt2) < EPSILON) {
		/* Wipe out adjacency of this vertex as we dont want to      */
		/* propagate through this one nothing - its on in/out border.*/
		VertexStruct *Vtemp;

		if (V -> PAdj == NULL) return NULL;

		Vtemp = V -> PAdj -> V;
		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 -> V);
		V -> PAdj = NULL;		/* And wipe out ours also... */
		return NULL;
	    }
	}

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

    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 || APX_EQ(Min, 1.0) ||
		Max <= 0.0 || APX_EQ(Max, 0.0)) return NULL;

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

	    /* Handle the Min end point: */
	    if (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 (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]) < EPSILON)
		    PSeg -> PtSeg[0][i] = 0.0;
		if (ABS(PSeg -> PtSeg[1][i]) < EPSILON)
		    PSeg -> PtSeg[1][i] = 0.0;
	    }
	    if (PT_EQ(PSeg -> PtSeg[0], PSeg -> PtSeg[1])) {
		MyFree((char *) PSeg, ALLOC_OTHER);
		return NULL;
	    }
	    return PSeg;
    }
    return NULL;				    /* Makes warning silent. */
}

/*****************************************************************************
*   Given a polygon with the intersection list, create the polylines loop(s) *
* out of it, which can be one of the two:				     *
* 1. Closed loop - all the intersection create a loop in one polygon.	     *
* 2. Open polyline - if the intersection crosses the polygon boundary. In    *
*    this case the two end point of the polyline, must lay on polygon	     *
*    boundary.								     *
* In both cases, the polyline will be as follows:			     *
* First point at first list element at PtSeg[0] (see InterSegmentStruct).    *
* Second point at first list element at PtSeg[1] (see InterSegmentStruct).   *
* Point i at list element (i-1) at PtSeg[0] (PtSeg[1] is not used!).         *
* In the closed loop case the last point is equal to first.		     *
* Both cases returns NULL terminated list.				     *
*****************************************************************************/
void LoopsFromInterList(PolygonStruct *Pl,
		InterSegListStruct **PClosed, InterSegListStruct **POpen)
{
    InterSegmentStruct *PSeg, *PSegHead, *PSegTemp, *PSegNewTemp;
    InterSegListStruct *PSLTemp;

#ifdef DEBUG3
    printf("Enter LoopsFromInterList\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 *)
		MyMalloc(sizeof(InterSegListStruct), ALLOC_OTHER);
	    PSLTemp -> PISeg = PSegHead;
	    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 (!PT_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 *)
		MyMalloc(sizeof(InterSegListStruct), ALLOC_OTHER);
	    PSLTemp -> PISeg = PSegHead;
	    PSLTemp -> Pnext = (*PClosed);
	    *PClosed = PSLTemp;
	}
    }

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

/*****************************************************************************
* Swap the two points in the InterSegmentStruct (modifies PtSeg & V entries) *
*****************************************************************************/
static void SwapPointInterList(InterSegmentStruct *PSeg)
{
    PointType Pt;
    VertexStruct *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;
}

/*****************************************************************************
* Delete one InterSegment PSeg, from InterSegmentList PSegList:		     *
*****************************************************************************/
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
	    FatalError("DeleteSegInterList: element to delete not found\n");
    }
}

/*****************************************************************************
* Search 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 deleting it from list   *
*****************************************************************************/
static InterSegmentStruct *FindMatchInterList(PointType Pt,
						InterSegmentStruct **PSegList)
{
    InterSegmentStruct *PSegTemp = (*PSegList);

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

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

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

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

#ifdef DEBUG3
    printf("Enter SortOpenInterList\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 *) MyMalloc(sizeof(SortOpenStruct),
					      ALLOC_OTHER);
	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 *) MyMalloc(sizeof(SortOpenStruct),
					      ALLOC_OTHER);
	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)
	    FatalError("SortOpenInterList: fail to sort intersection list\n");
	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. */
	    MyFree((char *) PSTemp -> Pnext, ALLOC_OTHER);
	    MyFree((char *) PSTemp, ALLOC_OTHER);
	    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;
		MyFree((char *) PSTemp2 -> Pnext, ALLOC_OTHER);
		MyFree((char *) PSTemp2, ALLOC_OTHER);
		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 &&
		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 (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 SortOpenInterList\n");
#endif /* DEBUG3 */
}

/*****************************************************************************
*   Reverse the order of the open loop pointed by PSHead.		     *
*****************************************************************************/
static void SortOpenReverseLoop(SortOpenStruct *PSHead)
{
    InterSegmentStruct *PINewHead = NULL, *PIHead, *PITemp;

    PIHead = PSHead -> PLSeg -> PISeg;

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

    PSHead -> PLSeg -> PISeg = PINewHead;
}

/*****************************************************************************
*   Insert new loop PSTemp, defines key by Pt and V (V defines the vertex    *
* and P is the points on it), into (decreasing) ordered list PSHead.	     *
*****************************************************************************/
static RealType SortOpenInsertOne(SortOpenStruct **PSHead,
       SortOpenStruct *PSTemp, PointType Pt, VertexStruct *V, PolygonStruct *Pl)
{
    int i = 0;
    RealType Key;
    PointType Pt1, Pt2;
    VertexStruct *VTemp = Pl -> V;
    SortOpenStruct *PSTail;

    PSTemp -> Pnext = NULL;

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

    PT_SUB(Pt1, V -> Pnext -> Pt, V -> Pt);
    PT_SUB(Pt2, Pt, V -> Pt);
    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;
}

/*****************************************************************************
*   Create a new vertex list, in reverse order. The original is not modified.*
*****************************************************************************/
VertexStruct *GenReverseVrtxList(VertexStruct *VIn)
{
    VertexStruct *VOutHead, *VOut;

    VOutHead = AllocVertex(0, 0, NULL, NULL);

    PT_COPY(VOutHead -> Pt, VIn -> Pt);
    VIn = VIn -> Pnext;

    while (VIn != NULL) {
	VOut = AllocVertex(0, 0, NULL, NULL);
	PT_COPY(VOut -> Pt, VIn -> Pt);
	PT_COPY(VOut -> Normal, VIn -> Normal);

	VOut -> Pnext = VOutHead;	     /* Chain them in reverse order. */
	VOutHead = VOut;

	VIn = VIn -> Pnext;
    }

    return VOutHead;
}

#ifdef DEBUG2

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

    do {
	printf("    %12lf %12lf %12lf\n", V -> Pt[0], V -> Pt[1], V -> Pt[2]);
	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],
#ifdef __MSDOS__
	    FP_SEG(PInt->V[0]),
#else
	    PInt->V[0],
#endif /* __MSDOS__ */
	    PInt->PtSeg[1][0], PInt->PtSeg[1][1], PInt->PtSeg[1][2],
#ifdef __MSDOS__
	    FP_SEG(PInt->V[1]));
#else
	    PInt->V[1]);
#endif /* __MSDOS__ */
	if (PInt -> Pnext == NULL)
	    printf("Exit vertex pointer %p\n", PInt -> V[1]);

	PInt = PInt -> Pnext;
    }
}

#endif /* DEBUG2 */