/*****************************************************************************
*   "Irit" - the 3d (not only polygonal) solid modeller.		     *
*									     *
* Written by:  Gershon Elber				Ver 0.2, Mar. 1990   *
******************************************************************************
*   Module to handle the high level Boolean operations. The other modules    *
* should only call this module to perform Boolean operations. All the        *
* operations are none-destructives, meaning the given data is not modified.  *
*   Note all the polygons of the two given objects must be convex, and the   *
* returned object will also have only convex polygons!			     *
*****************************************************************************/

/* #define DEBUG2	         If defined, defines some printing routines. */

#include <stdio.h>
#include <ctype.h>
#include <math.h>
#include <string.h>
#include <signal.h>
#include "program.h"
#include "allocate.h"
#include "attribut.h"
#include "booleang.h"
#include "booleanl.h"
#include "convex.h"
#include "poly_cln.h"
#include "ctrl-brk.h"
#include "geomat3d.h"
#include "objects.h"
#include "windows.h"

int BooleanOutputInterCurve = FALSE;	/* Kind of output from Boolean oper. */

static jmp_buf LclLongJumpBuffer;	     /* Used in fatal Boolean error. */
static int FatalErrorType,		     /* Type of fatal Boolean error. */
	   BooleanOperation;	       /* One of BooleanOR, BooleanAND, etc. */

static IPObjectStruct *BooleanCombineThreeObjs(IPObjectStruct *PObj1,
					       IPObjectStruct *PObj2,
					       IPObjectStruct *PObj3);
static IPObjectStruct *BooleanCoplanar(IPObjectStruct *PObj1,
				       IPObjectStruct *PObj2, int BoolOper);
static IPObjectStruct *VerifyBooleanInput(IPObjectStruct *PObj1,
					  IPObjectStruct *PObj2, int Oper);
static IPPolygonStruct *ComputeRotatedPolys(IPPolygonStruct *Pl, int CopyPl,
							     MatrixType RotMat);

static void BooleanFPE(int Type);
static void SetBooleanOutputKind(void);

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

/*****************************************************************************
*   Verify input for Booleans. Ret. NULL if OK, otherwise an obj to return.  *
*****************************************************************************/
static IPObjectStruct *VerifyBooleanInput(IPObjectStruct *PObj1,
					  IPObjectStruct *PObj2, int Oper)
{
    IPObjectStruct *PObj;
    IPPolygonStruct *Pl;

    BooleanOperation = Oper;

    SetBooleanOutputKind();

    if (!IP_IS_POLY_OBJ(PObj1) || (PObj2 != NULL && !IP_IS_POLY_OBJ(PObj2)))
	IritFatalError("Boolean: operation on non polygonal object(s).");

    signal(SIGFPE, BooleanFPE);		 /* Will trap floating point errors. */

    switch (Oper) {
	case BOOL_OPER_OR:
	    if (IP_IS_POLYLINE_OBJ(PObj1) && IP_IS_POLYLINE_OBJ(PObj2)) {
		if (PObj1 -> U.Pl == NULL)
		    PObj = CopyObject(NULL, PObj2, FALSE);
		else {
		    PObj = CopyObject(NULL, PObj1, FALSE);
		    Pl = IritPrsrGetLastPoly(PObj -> U.Pl);
		    Pl -> Pnext = CopyPolygonList(PObj2 -> U.Pl);
		}
	    	return PObj;
	    }
	case BOOL_OPER_AND:
	case BOOL_OPER_SUB:
	case BOOL_OPER_CUT:
	case BOOL_OPER_MERGE:
	case BOOL_OPER_NEG:
    	    if (IP_IS_POLYLINE_OBJ(PObj1) ||
		(PObj2 != NULL && IP_IS_POLYLINE_OBJ(PObj2))) {
    	    	WndwInputWindowPutStr(
		"Boolean: illegal operation on mixed polygon/line geometric object(s).");
		PObj = GenPolyObject("", NULL, NULL);
		return PObj;
    	    }

	    if (Oper != BOOL_OPER_NEG) {
	    	ConvexPolyObject(PObj1);/* Make sure all polygons are convex.*/
	    	ConvexPolyObject(PObj2);
	    }

	    return NULL;
	default:
    	    IritFatalError("Boolean: undefined Boolean operation.");
    	    return NULL;			     /* Make warning silent. */
    }
}

/*****************************************************************************
*   Perform a Boolean OR between two objects.				     *
*****************************************************************************/
IPObjectStruct *BooleanOR(IPObjectStruct *PObj1, IPObjectStruct *PObj2)
{
    IPObjectStruct *PObj;
    IPPolygonStruct *Pl;

    if ((PObj = VerifyBooleanInput(PObj1, PObj2, BOOL_OPER_OR)) != NULL)
	return PObj;
    else {
	if (setjmp(LclLongJumpBuffer) == 0) { /* Its the setjmp itself call! */
	    signal(SIGFPE, BooleanFPE);	 /* Will trap floating point errors. */
	    if (BooleanOutputInterCurve)
		PObj = BooleanLow1Out2(PObj1, PObj2);/* Ret intersection crv.*/
	    else
		PObj = BooleanCombineThreeObjs(BooleanLow1Out2(PObj1, PObj2),
					       BooleanLow1Out2(PObj2, PObj1),
					       BooleanCoplanar(PObj1, PObj2,
							       BOOL_OPER_OR));

	}
	else {
	    /* We gain control from fatal error long jump - usually we should*/
	    /* return empty object, but if error is no intersection between  */
	    /* the two objects, we assume they have no common volume and     */
	    /* return a new object consists of the concat. of all polygons!  */
	    if (FatalErrorType != FTL_BOOL_NO_INTER) {
		PObj = GenPolyObject("", NULL, NULL);/* Return empty object. */
	    }
	    else {
		if (PObj1 -> U.Pl == NULL)
		    PObj = CopyObject(NULL, PObj2, FALSE);
		else {
		    PObj = CopyObject(NULL, PObj1, FALSE);/* Copy Obj1 polys.*/
		    Pl = PObj -> U.Pl;
		    while (Pl -> Pnext) Pl = Pl -> Pnext;
		    Pl -> Pnext = CopyPolygonList(PObj2 -> U.Pl);/*Obj2 poly.*/
		}
	    }
	}
    }

    signal(SIGFPE, DefaultFPEHandler);	            /* Default FPE trapping. */

    AttrSetObjectColor(PObj, BooleanOutputInterCurve ?
			     GlblICrvColor :
			     GlblBoolColor);   /* Default bool object color. */

    return PObj;
}

/*****************************************************************************
*   Perform a Boolean AND between two objects.				     *
*****************************************************************************/
IPObjectStruct *BooleanAND(IPObjectStruct *PObj1, IPObjectStruct *PObj2)
{
    IPObjectStruct *PObj;

    if ((PObj = VerifyBooleanInput(PObj1, PObj2, BOOL_OPER_AND)) != NULL)
	return PObj;
    else {
	if (setjmp(LclLongJumpBuffer) == 0) { /* Its the setjmp itself call! */
	    signal(SIGFPE, BooleanFPE);	 /* Will trap floating point errors. */
	    if (BooleanOutputInterCurve)
		PObj = BooleanLow1In2(PObj1, PObj2);/* Ret intersection crv. */
	    else
		PObj = BooleanCombineThreeObjs(BooleanLow1In2(PObj1, PObj2),
					       BooleanLow1In2(PObj2, PObj1),
					       BooleanCoplanar(PObj1, PObj2,
							       BOOL_OPER_AND));

	}
	else {/* We gain control from fatal error long jump - ret empty obj. */
	    PObj = GenPolyObject("", NULL, NULL);
	}
    }

    signal(SIGFPE, DefaultFPEHandler);		    /* Default FPE trapping. */

    AttrSetObjectColor(PObj, BooleanOutputInterCurve ?
			     GlblICrvColor :
			     GlblBoolColor);   /* Default bool object color. */

    return PObj;
}

/*****************************************************************************
*   Perform a Boolean SUBSTRACT between two objects (PObj1 - PObj2).	     *
*****************************************************************************/
IPObjectStruct *BooleanSUB(IPObjectStruct *PObj1, IPObjectStruct *PObj2)
{
    IPObjectStruct *PObj, *PTemp, *PTempRev;

    if ((PObj = VerifyBooleanInput(PObj1, PObj2, BOOL_OPER_SUB)) != NULL)
	return PObj;
    else {
	if (setjmp(LclLongJumpBuffer) == 0) { /* Its the setjmp itself call! */
	    signal(SIGFPE, BooleanFPE);	 /* Will trap floating point errors. */
	    /* The 1 in 2 must be reversed (the inside/outside orientation): */
	    if (BooleanOutputInterCurve) {
		PObj = BooleanLow1In2(PObj2, PObj1);/* Ret intersection crv. */
	    }
	    else {
		PTemp = BooleanLow1In2(PObj2, PObj1);
		PTempRev = BooleanNEG(PTemp);
		IPFreeObject(PTemp);

		PObj = BooleanCombineThreeObjs(BooleanLow1Out2(PObj1, PObj2),
					       PTempRev,
					       BooleanCoplanar(PObj1, PObj2,
							       BOOL_OPER_SUB));
	    }
	}
	else {/* We gain control from fatal error long jump - ret empty obj. */
	    PObj = GenPolyObject("", NULL, NULL);
	}
    }

    signal(SIGFPE, DefaultFPEHandler);	            /* Default FPE trapping. */

    AttrSetObjectColor(PObj, BooleanOutputInterCurve ?
			     GlblICrvColor :
			     GlblBoolColor);   /* Default bool object color. */

    return PObj;
}

/*****************************************************************************
*   Perform a Boolean CUT between two objects (PObj1 / PObj2).		     *
*****************************************************************************/
IPObjectStruct *BooleanCUT(IPObjectStruct *PObj1, IPObjectStruct *PObj2)
{
    IPObjectStruct *PObj;

    if ((PObj = VerifyBooleanInput(PObj1, PObj2, BOOL_OPER_CUT)) != NULL)
	return PObj;
    else {
	if (setjmp(LclLongJumpBuffer) == 0) { /* Its the setjmp itself call! */
	    signal(SIGFPE, BooleanFPE);	 /* Will trap floating point errors. */
	    /* The 1 in 2 must be reversed (the inside/outside orientation): */
	    if (BooleanOutputInterCurve) {
		PObj = BooleanLow1In2(PObj2, PObj1);/* Ret intersection crv. */
	    }
	    else {
		PObj = BooleanLow1Out2(PObj1, PObj2);
	    }
	}
	else {/* We gain control from fatal error long jump - ret empty obj. */
	    PObj = GenPolyObject("", NULL, NULL);
	}
    }

    signal(SIGFPE, DefaultFPEHandler);	            /* Default FPE trapping. */

    AttrSetObjectColor(PObj, BooleanOutputInterCurve ?
			     GlblICrvColor :
			     GlblBoolColor);   /* Default bool object color. */

    return PObj;
}

/*****************************************************************************
*   Perform a Boolean MERGE between two objects.			     *
*****************************************************************************/
IPObjectStruct *BooleanMERGE(IPObjectStruct *PObj1, IPObjectStruct *PObj2)
{
    IPObjectStruct *PObj;
    IPPolygonStruct *Pl;

    if ((PObj = VerifyBooleanInput(PObj1, PObj2, BOOL_OPER_MERGE)) != NULL)
	return PObj;
    else {
	if (PObj1 -> U.Pl == NULL)
    	    PObj = CopyObject(NULL, PObj2, FALSE);
    	else {
    	    PObj = CopyObject(NULL, PObj1, FALSE);       /* Copy Obj1 polys. */
	    Pl = PObj -> U.Pl;
	    while (Pl -> Pnext) Pl = Pl -> Pnext;
	    Pl -> Pnext = CopyPolygonList(PObj2 -> U.Pl);     /* Obj2 polys. */
	}
    }

    AttrSetObjectColor(PObj, GlblBoolColor);   /* Default bool object color. */

    return PObj;
}

/*****************************************************************************
*   Perform a Boolean NEGATION of given object PObj.			     *
*  Negation is simply reversing the direction of the plane equation of each  *
* polygon - the simplest Boolean operation...				     *
*****************************************************************************/
IPObjectStruct *BooleanNEG(IPObjectStruct *PObj)
{
    int i;
    IPObjectStruct *PTemp;
    IPPolygonStruct *Pl;
    IPVertexStruct *V;

    if ((PTemp = VerifyBooleanInput(PObj, NULL, BOOL_OPER_NEG)) != NULL)
	return PTemp;
    else {
	PTemp = CopyObject(NULL, PObj, FALSE); /* Make fresh copy of object. */

	/* Scans all polygons and reverse plane equation and their vetrex    */
	/* list (cross prod. of consecutive edges must be in normal dir.).   */
	Pl = PTemp -> U.Pl;
	while (Pl != NULL) {
	    for (i = 0; i < 4; i++) Pl -> Plane[i] = (-Pl -> Plane[i]);
	    IP_RST_CONVEX_POLY(Pl);

	    /* Invert vertices normals as well. */
	    V = Pl -> PVertex;
	    do {
		PT_SCALE(V -> Normal, -1.0);
		V = V -> Pnext;
	    }
	    while (V != NULL && V != Pl -> PVertex);

	    Pl = Pl -> Pnext;
	}
    }

    signal(SIGFPE, DefaultFPEHandler);		    /* Default FPE trapping. */

    AttrSetObjectColor(PTemp, BooleanOutputInterCurve ?
			      GlblICrvColor :
			      GlblBoolColor);  /* Default bool object color. */

    return PTemp;
}

/*****************************************************************************
*  Combining three geometric objects, by simply concat. their polygon lists: *
*  Any object may be NULL in which only the other two are merged.	     *
*****************************************************************************/
static IPObjectStruct *BooleanCombineThreeObjs(IPObjectStruct *PObj1,
					       IPObjectStruct *PObj2,
					       IPObjectStruct *PObj3)
{
    IPPolygonStruct *Pl;

    if (PObj1 == NULL) {
	PObj1 = PObj2;
	PObj2 = PObj3;
	PObj3 = NULL;
    }
    if (PObj2 == NULL) {
	PObj2 = PObj3;
	PObj3 = NULL;
    }

    if (PObj1 != NULL)
	CleanUpPolygonList(&PObj1 -> U.Pl);
    if (PObj2 != NULL)
	CleanUpPolygonList(&PObj2 -> U.Pl);
    if (PObj3 != NULL)
	CleanUpPolygonList(&PObj3 -> U.Pl);

    if (PObj2 != NULL) {
	/* Concat the polygons of PObj2 after the polygons of PObj1. */
	Pl = PObj1 -> U.Pl;
	while (Pl -> Pnext != NULL)
	    Pl = Pl -> Pnext;
	Pl -> Pnext = PObj2 -> U.Pl;  /* Concat. the polygons into one list. */
	PObj2 -> U.Pl = NULL;	    /* And release the second object header. */
	IPFreeObject(PObj2);

        if (PObj3 != NULL) {
	    /* Concat the polygons of PObj3 after the polygons of PObj1/2. */
	    while (Pl -> Pnext != NULL)
		Pl = Pl -> Pnext;
	    Pl -> Pnext = PObj3 -> U.Pl;  /* Concat. polygons into one list. */
	    PObj3 -> U.Pl = NULL;   /* And release the second object header. */
	    IPFreeObject(PObj3);
        }
    }

    return PObj1;
}

/*****************************************************************************
*  If required (I.e. Object "COPLANAR" is set to TRUE) search all coplanar   *
* polygons in PObj1/2 and invoke the necessary two dimenstional Boolean.     *
*****************************************************************************/
static IPObjectStruct *BooleanCoplanar(IPObjectStruct *PObj1,
				       IPObjectStruct *PObj2, int BoolOper)
{
    int HandleCoplanar;
    MatrixType RotMat;
    IPObjectStruct *PObj,
	*PCoplanar = GetObject("COPLANAR");
    IPPolygonStruct *Pl, *PlTmp, *Pl1, *Pl2, *Pl1XY, *Pl2XY, *Pl1XYR, *Pl2XYR,
	*PlOut = NULL;

    if (PCoplanar == NULL || !IP_IS_NUM_OBJ(PCoplanar)) {
	WndwInputWindowPutStr("No numeric object name COPLANAR is defined");
	HandleCoplanar = FALSE;
    }
    else
	HandleCoplanar = (int) (PCoplanar -> U.R);

    if (!HandleCoplanar)
	return NULL;

    for ( Pl1 = PObj1 -> U.Pl; Pl1 != NULL; Pl1 = Pl1 -> Pnext) {
        for ( Pl2 = PObj2 -> U.Pl; Pl2 != NULL; Pl2 = Pl2 -> Pnext) {
	    RealType
		*Plane1 = Pl1 -> Plane,
		*Plane2 = Pl2 -> Plane;
	    int Shared = APX_EQ(Plane1[0], Plane2[0]) &&
		         APX_EQ(Plane1[1], Plane2[1]) &&
		         APX_EQ(Plane1[2], Plane2[2]) &&
		         APX_EQ(Plane1[3], Plane2[3]),
		AntiShared = APX_EQ(Plane1[0], -Plane2[0]) &&
			     APX_EQ(Plane1[1], -Plane2[1]) &&
			     APX_EQ(Plane1[2], -Plane2[2]) &&
			     APX_EQ(Plane1[3], -Plane2[3]);

	    if (!Shared && !AntiShared) {
		/* The two polygons are not coplanar, do not intersect. */
		continue;
	    }

	    if (Shared) {
		switch (BoolOper) {
		    case BOOL_OPER_AND:
		    case BOOL_OPER_OR:
		    case BOOL_OPER_SUB:
			GenRotateMatrix(RotMat, Pl1 -> Plane);
			Pl1XY = ComputeRotatedPolys(Pl1, TRUE, RotMat);
			Pl2XY = ComputeRotatedPolys(Pl2, TRUE, RotMat);
			/* Find the inverse matrix. */
			if (!MatInverseMatrix(RotMat, RotMat))
			    IritFatalError("BooleanCoplanar: Inverse matrix does not exists");

			if ((Pl = Boolean2D(Pl1XY, Pl2XY, BoolOper)) != NULL) {
			    Pl = ComputeRotatedPolys(Pl, FALSE, RotMat);
			    for (PlTmp = Pl;
				 PlTmp -> Pnext != NULL;
				 PlTmp = PlTmp -> Pnext);
			    PlTmp -> Pnext = PlOut;
			    PlOut = Pl;
			}

			IPFreePolygonList(Pl1XY);
			IPFreePolygonList(Pl2XY);
			break;

		    default:
		        IritFatalError("BooleanCoplanar: Unsupported Boolean operation BooleanCoplanar");
		        return NULL;
		}
	    }
	    else if (AntiShared) {
		switch (BoolOper) {
		    case BOOL_OPER_AND:
		    case BOOL_OPER_SUB:
			WndwInputWindowPutStr("Antishared coplanar polygons are ignored.");
			break;

		    case BOOL_OPER_OR:
			GenRotateMatrix(RotMat, Pl1 -> Plane);
			Pl1XY = ComputeRotatedPolys(Pl1, TRUE, RotMat);
			Pl1XYR = CopyPolygonList(Pl1XY);
			IritPrsrReverseVrtxList(Pl1XYR);
			Pl2XY = ComputeRotatedPolys(Pl2, TRUE, RotMat);
			Pl2XYR = CopyPolygonList(Pl2XY);
			IritPrsrReverseVrtxList(Pl2XYR);

			if (!MatInverseMatrix(RotMat, RotMat))
			    IritFatalError("BooleanCoplanar: Inverse matrix does not exists");

			if ((Pl = Boolean2D(Pl1XY, Pl2XYR, BOOL_OPER_SUB)) != NULL) {
			    Pl = ComputeRotatedPolys(Pl, FALSE, RotMat);
			    for (PlTmp = Pl;
				 PlTmp -> Pnext != NULL;
				 PlTmp = PlTmp -> Pnext);
			    PlTmp -> Pnext = PlOut;
			    PlOut = Pl;
			}
			if ((Pl = Boolean2D(Pl2XY, Pl1XYR, BOOL_OPER_SUB)) != NULL) {
			    Pl = ComputeRotatedPolys(Pl, FALSE, RotMat);
			    for (PlTmp = Pl;
				 PlTmp -> Pnext != NULL;
				 PlTmp = PlTmp -> Pnext);
			    PlTmp -> Pnext = PlOut;
			    PlOut = Pl;
			}

			IPFreePolygonList(Pl1XY);
			IPFreePolygonList(Pl2XY);
			IPFreePolygonList(Pl1XYR);
			IPFreePolygonList(Pl2XYR);
			break;

		    default:
		        IritFatalError("BooleanCoplanar: Unsupported Boolean operation BooleanCoplanar");
		        return NULL;
		}
	    }
	}
    }

    PObj = IPAllocObject("", IP_OBJ_POLY, NULL);
    PObj -> U.Pl = PlOut;
    return PObj;
}

/*****************************************************************************
*   Routine to optionally copy (if CpolyOnePl) a single polygon and rotate   *
* according to the rotation matrix provided. If, however, CopyOnePl is False *
* all polygons in list are converted.					     *
*****************************************************************************/
static IPPolygonStruct *ComputeRotatedPolys(IPPolygonStruct *Pl, int CopyOnePl,
							     MatrixType RotMat)
{
    IPVertexStruct *V, *VHead;
    IPPolygonStruct *PlNext, *PlTemp,
	*PlOut = NULL;

    while (Pl != NULL) {
	PlNext = Pl -> Pnext;

	if (CopyOnePl) {
	    PlTemp = IPAllocPolygon(Pl -> Count, Pl -> Tags,
				    CopyVertexList(Pl -> PVertex), NULL);
	    PLANE_COPY(PlTemp ->Plane, Pl -> Plane);
	    Pl = PlTemp;
	}

	V = VHead = Pl -> PVertex;
	do {				    /* Transform the polygon itself. */
	    PointType PTemp;

	    PT_ADD(PTemp, V -> Coord, V -> Normal);

	    MatMultVecby4by4(V -> Coord, V -> Coord, RotMat);

	    MatMultVecby4by4(PTemp, PTemp, RotMat);   /* Update normal. */
	    PT_SUB(V -> Normal, PTemp, V -> Coord);
	    PT_NORMALIZE(V -> Normal);

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

	Pl -> Pnext = PlOut;
	PlOut = Pl;

	if (CopyOnePl)
	    break;
	else
	    Pl = PlNext;
    }

    return PlOut;
}

/*****************************************************************************
*   Routine that is called from the floating point package in case of fatal  *
* floating point error. Print error message, and quit the Boolean module.    *
*****************************************************************************/
static void BooleanFPE(int Type)
{
    char Line[LINE_LEN];

    sprintf(Line, "Floating point error %d.", Type);
    WndwInputWindowPutStr(Line);

    FatalErrorType = FTL_BOOL_FPE;

    longjmp(LclLongJumpBuffer, 1);
}

/*****************************************************************************
*   Routine to select the output kind needed. Output of a Boolean operator   *
* can be the real Boolean result model, or the intersection curves only.     *
*****************************************************************************/
static void SetBooleanOutputKind(void)
{
    IPObjectStruct *PObj = GetObject("INTERCRV");

    if (PObj == NULL || !IP_IS_NUM_OBJ(PObj)) {
	WndwInputWindowPutStr("No numeric object name INTERCRV is defined");
	BooleanOutputInterCurve = FALSE;
    }
    else
	BooleanOutputInterCurve = !APX_EQ(PObj -> U.R, 0.0);
}

/*****************************************************************************
*   Routine that is called by the bool-low module in fatal cases errors.     *
* Will print error message and long jump using LclLongJumpBuffer.	     *
*****************************************************************************/
void FatalBooleanError(int ErrorType)
{
    char s[LINE_LEN_LONG];

    FatalErrorType = ErrorType;

    switch (ErrorType) {
	case FTL_BOOL_NO_INTER:
	    /* If operation is union (OR), then if no intersection we assume */
	    /* the objects have no common volume and simply combine them!    */
	    if (BooleanOperation != BOOL_OPER_OR) {
		sprintf(s, "Boolean: %s",
			"Objects do not intersect - Empty object result");
		WndwInputWindowPutStr(s);
	    }
	    break;
	default:
	    sprintf(s, "Boolean: Undefined Fatal Error (%d !?)", ErrorType);
	    WndwInputWindowPutStr(s);
	    break;
    }

    longjmp(LclLongJumpBuffer, 1);
}

#ifdef DEBUG2

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

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

#endif /* DEBUG2 */