/******************************************************************************
* TrimCntr.c - Converts a set of contours in parametric space to trimmed srfs *
*******************************************************************************
* (C) Gershon Elber, Technion, Israel Institute of Technology                 *
*******************************************************************************
* Written by Gershon Elber, May 96.					      *
******************************************************************************/

#include "irit_sm.h"
#include "trim_lib.h"
#include "user_lib.h"
#include "allocate.h"
#include "iritprsr.h"
#include "geomat3d.h"
#include "convex.h"

#define CONTOUR_EPS	1e-3

#define VRTX_BOUNDARY_TAG	0x10
#define IS_BOUNDARY_VRTX(Vrtx)	((Vrtx) -> Tags & VRTX_BOUNDARY_TAG)
#define SET_BOUNDARY_VRTX(Vrtx)	((Vrtx) -> Tags |= VRTX_BOUNDARY_TAG)
#define RST_BOUNDARY_VRTX(Vrtx)	((Vrtx) -> Tags &= ~VRTX_BOUNDARY_TAG)

#define ASSIGN_VERTEX_COORDS(V, X, Y, Z) { (V) -> Coord[0] = (X); \
					   (V) -> Coord[1] = (Y); \
					   (V) -> Coord[2] = (Z); \
					   SET_BOUNDARY_VRTX(V); }

static CagdCrvStruct *TrimPolylines2LinBsplineCrvs(IPPolygonStruct *Polys,
						   CagdRType UMin,
						   CagdRType UMax,
						   CagdRType VMin,
						   CagdRType VMax);
static IPPolygonStruct *ClosedCntrsFromOpenCntrs(IPPolygonStruct *OpenCntrs,
						 IPPolygonStruct *Domain,
						 CagdRType UMin,
						 CagdRType UMax,
						 CagdRType VMin,
						 CagdRType VMax);
static IPPolygonStruct *SplitDomainByCntr(IPPolygonStruct *Domain,
					  IPPolygonStruct *Cntr,
					  CagdRType UMin,
					  CagdRType UMax,
					  CagdRType VMin,
					  CagdRType VMax);
static CagdRType VertexWeight(IPVertexStruct *V,
			      CagdRType UMin,
			      CagdRType UMax,
			      CagdRType VMin,
			      CagdRType VMax);

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

/*****************************************************************************
* DESCRIPTION:                                                               M
*   Creats a set of trimmed surfaces as defined by the given set of contours M
* than can contain either closed or open contours.  Open contours must       M
* terminate at the boundary of the parameteric domain of the surface. Closed M
* contours must be completely contained in the parametric domain with last   M
* point equals first.                                                        M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf:        To trim into pieces.                                         M
*   Cntrs:      Polylines to use as seperating edges.                        M
*                                                                            *
* RETURN VALUE:                                                              M
*   TrimSrfStruct *:   List of trimmed surface pieces.                       M
*                                                                            *
* SEE ALSO:                                                                  M
*                                                                            M
*                                                                            *
* KEYWORDS:                                                                  M
*   TrimSrfsFromContours                                                     M
*****************************************************************************/
TrimSrfStruct *TrimSrfsFromContours(CagdSrfStruct *Srf,
				    IPPolygonStruct *Cntrs)
{
    CagdRType UMin, UMax, VMin, VMax;
    IPPolygonStruct *Cntr, *SrfDomain,
	*OpenCntrs = NULL,
	*ClosedCntrs = NULL;
    IPVertexStruct *VTmp;
    TrimSrfStruct
	*RetTrimSrfs = NULL;

    CagdSrfDomain(Srf, &UMin, &UMax, &VMin, &VMax);

    Cntrs = CopyPolygonList(Cntrs);            /* Keep the original intact. */

    /* Split the contours into closed and open loops. */
    while (Cntrs) {
	IPPolygonStruct
	    *Cntr = Cntrs;
	IPVertexStruct
	    *V = Cntr -> PVertex,
	    *VLast = IritPrsrGetLastVrtx(V);

	Cntrs = Cntrs -> Pnext;
	Cntr -> Pnext = NULL;

	if (PT_APX_EQ_EPS(V -> Coord, VLast -> Coord, CONTOUR_EPS)) {
	    LIST_PUSH(Cntr, ClosedCntrs);
	}
	else {
	    LIST_PUSH(Cntr, OpenCntrs);
	}
    }

    /* Create a domain polygon of entire surface as top level. */
    VTmp = IPAllocVertex(0, NULL, NULL);
    ASSIGN_VERTEX_COORDS(VTmp, UMin, VMin, 0.0);
    VTmp = IPAllocVertex(0, NULL, VTmp);
    ASSIGN_VERTEX_COORDS(VTmp, UMax, VMin, 0.0);
    VTmp = IPAllocVertex(0, NULL, VTmp);
    ASSIGN_VERTEX_COORDS(VTmp, UMax, VMax, 0.0);
    VTmp = IPAllocVertex(0, NULL, VTmp);
    ASSIGN_VERTEX_COORDS(VTmp, UMin, VMax, 0.0);
    VTmp = IPAllocVertex(0, NULL, VTmp);
    ASSIGN_VERTEX_COORDS(VTmp, UMin, VMin, 0.0);
    SrfDomain = IPAllocPolygon(0, IritPrsrReverseVrtxList2(VTmp), NULL);

    /* Split open contours into closed regions with the aid of the boundary. */
    if (OpenCntrs != NULL) {

	/* Mark all vertices in the open contours as non boundary. */
	for (Cntr = OpenCntrs; Cntr != NULL; Cntr = Cntr -> Pnext) {
	    IPVertexStruct *V;

	    for (V = Cntr -> PVertex; V != NULL; V = V -> Pnext)
	        RST_BOUNDARY_VRTX(V);	    
	}
	OpenCntrs = ClosedCntrsFromOpenCntrs(OpenCntrs, SrfDomain,
					     UMin, UMax, VMin, VMax);
	IPFreePolygon(SrfDomain);

#	ifdef DEBUG_OPEN_CONTOURS
	{
	    static int
		j = -5;
	    int i = 0;

	    while (OpenCntrs) { /* Distructive!!! */
		IPObjectStruct *PObj, *PObj2;
		MatrixType Mat;

		Cntr = OpenCntrs;
		OpenCntrs = OpenCntrs -> Pnext;
		Cntr -> Pnext = NULL;

		PObj = GenPOLYObject(Cntr);
		IP_SET_POLYLINE_OBJ(PObj);

		MatGenMatTrans(j, 0.0, i++, Mat);
		PObj2 = GMTransformObject(PObj, Mat);

		IritPrsrStdoutObject(PObj2);
		IPFreeObject(PObj);
		IPFreeObject(PObj2);
	    }

	    j += 5.0;
	}
#	endif /* DEBUG_OPEN_CONTOURS */
    }
    else {
	/* Place the entire domain as one open contour. */
	OpenCntrs = SrfDomain;
    }

    /* We now have only open contours. OpenCntrs contains a set of disjoint */
    /* loops formed with the boundary.  However, ClosedCntrs contains	    */
    /* contours that must be classified into the contours in OpenCntrs as a */
    /* whole hierarchy.							    */

    while (OpenCntrs) {
	IPPolygonStruct *ClosedCntr,
	    *PrevCntr = NULL,
	    *ClosedCntrsOfOpenCntr = NULL,
	    *OpenCntr = OpenCntrs;
	TrimSrfStruct *TrimSrf;

	OpenCntrs = OpenCntrs -> Pnext;
	OpenCntr -> Pnext = NULL;

	for (ClosedCntr = ClosedCntrs; ClosedCntr != NULL; ) {
	    if (CGPolygonRayInter(OpenCntr,
				  ClosedCntr -> PVertex -> Coord, 0) & 0x01) {
		IPPolygonStruct
		    *NextCntr = ClosedCntr -> Pnext;

		/* We have odd number of intersections - this closed loop is */
		/* to be classified inside the open loop, OpenCntr.	     */
		if (PrevCntr == NULL)
		    ClosedCntrs = NextCntr;
		else
		    PrevCntr -> Pnext = NextCntr;

		LIST_PUSH(ClosedCntr, ClosedCntrsOfOpenCntr);

		ClosedCntr = NextCntr;
	    }
	    else {
		PrevCntr = ClosedCntr;
		ClosedCntr = ClosedCntr -> Pnext;
	    }
	}

	/* Create a trimmed surface using OpenCntr as outerloop and all the */
	/* other ClosedCntrsOfOpenCntr loops inside it.			    */
	OpenCntr -> Pnext = ClosedCntrsOfOpenCntr;
	TrimSrf = TrimSrfNew2(CagdSrfCopy(Srf),
			      TrimPolylines2LinBsplineCrvs(OpenCntr,
						       UMin, UMax, VMin, VMax),
			      TRUE);
	IPFreePolygonList(OpenCntr);
	LIST_PUSH(TrimSrf, RetTrimSrfs);
    }

    if (ClosedCntrs != NULL) {
	/* We should have classified all ClosedCntrs at this time. */
	TrimFatalError(TRIM_ERR_INCONSISTENT_CNTRS);
    }

    return RetTrimSrfs;
}

/*****************************************************************************
* DESCRIPTION:                                                               *
* Returns a list of linear Bspline curves constructed from given polylines.  *
* Adjust the trimming curves so that the importance of curves on the surface *
* boundary will get much more importance, parametrically.		     *
*                                                                            *
* PARAMETERS:                                                                *
*   Polys:    To convert to linear bspline curves.                           *
*   UMin, UMax, VMin, VMax:   Domain of surface whose Polys are to serve as  *
*			      trimming curves.		                     *
*                                                                            *
* RETURN VALUE:                                                              *
*   CagdCrvStruct *:  Linear Bspline curves representing Poly.               *
*****************************************************************************/
static CagdCrvStruct *TrimPolylines2LinBsplineCrvs(IPPolygonStruct *Polys,
						   CagdRType UMin,
						   CagdRType UMax,
						   CagdRType VMin,
						   CagdRType VMax)
{
    CagdCrvStruct *Crv,
	*Crvs = UserPolylines2LinBsplineCrvs(Polys);

    for (Crv = Crvs; Crv != NULL; Crv = Crv -> Pnext) {
	int i;
	CagdRType
	    *PtsU = Crv -> Points[1],
	    *PtsV = Crv -> Points[2],
	    *KV = &Crv -> KnotVector[2],
	    ParamOffset = 0.0;

	for (i = Crv -> Length; i > 1; i--, PtsU++, PtsV++) {
	    if ((APX_EQ(PtsU[0], UMin) ||
		 APX_EQ(PtsU[0], UMax) ||
		 APX_EQ(PtsV[0], VMin) ||
		 APX_EQ(PtsV[0], VMax)) &&
		(APX_EQ(PtsU[1], UMin) ||
		 APX_EQ(PtsU[1], UMax) ||
		 APX_EQ(PtsV[1], VMin) ||
		 APX_EQ(PtsV[1], VMax)))
		ParamOffset += 1.0;
	    *KV++ += ParamOffset;
	}
    }

    return Crvs;
}


/*****************************************************************************
* DESCRIPTION:                                                               *
*   Sorts, in place, the open contours into order so we can clip the         *
* rectangular surface domain's boundary with.                                *
*                                                                            *
* PARAMETERS:                                                                *
*   OpenCntrs:                To sort along the boundary.                    *
*   Domain:		      A polyline representing the current domain,    *
*			      starting from the entire surface's domain.     *
*   UMin, UMax, VMin, VMax:   Domain of surface.		             *
*                                                                            *
* RETURN VALUE:                                                              *
*   IPPolygonStruct *:        Sorted, in place, list of contours.            *
*****************************************************************************/
static IPPolygonStruct *ClosedCntrsFromOpenCntrs(IPPolygonStruct *OpenCntrs,
						 IPPolygonStruct *Domain,
						 CagdRType UMin,
						 CagdRType UMax,
						 CagdRType VMin,
						 CagdRType VMax)
{
    IPPolygonStruct
	*InsideCntrs = NULL,
	*OutsideCntrs = NULL,
	*SubdivCntr = OpenCntrs;
    IPVertexStruct
	*VFirst = SubdivCntr -> PVertex,
	*VLast = IritPrsrGetLastVrtx(VFirst);
    CagdRType
        VFirstWeight = VertexWeight(VFirst, UMin, UMax, VMin, VMax),
        VLastWeight = VertexWeight(VLast, UMin, UMax, VMin, VMax);

    OpenCntrs = OpenCntrs -> Pnext;
    SubdivCntr -> Pnext = NULL;
    if (VFirstWeight > VLastWeight)
	SWAP(CagdRType, VFirstWeight, VLastWeight);

    if (OpenCntrs != NULL) {
	IPPolygonStruct *LastCntr, *Cntrs1, *Cntrs2,
	    *InsideDomain, *OutsideDomain;

	/* Split all other open contours to contours inside weights' range */
	/* assigned to the end points of SplitCntr and outside that range. */
	while (OpenCntrs) {
	    IPPolygonStruct
	        *Cntr = OpenCntrs;
	    IPVertexStruct
	        *V2First = Cntr -> PVertex,
	        *V2Last = IritPrsrGetLastVrtx(V2First);
	    CagdRType
	        V2FirstWeight = VertexWeight(V2First, UMin, UMax, VMin, VMax),
	        V2LastWeight = VertexWeight(V2Last, UMin, UMax, VMin, VMax);

	    OpenCntrs = OpenCntrs -> Pnext;
	    Cntr -> Pnext = NULL;

	    if (V2FirstWeight > VFirstWeight && V2FirstWeight < VLastWeight &&
		V2LastWeight > VFirstWeight && V2LastWeight < VLastWeight) {
		LIST_PUSH(Cntr, InsideCntrs);
	    }
	    else if ((V2FirstWeight < VFirstWeight ||
		      V2FirstWeight > VLastWeight) &&
		     (V2LastWeight < VFirstWeight ||
		      V2LastWeight > VLastWeight)) {
		LIST_PUSH(Cntr, OutsideCntrs);
	    }
	    else {
		TrimFatalError(TRIM_ERR_INCONSISTENT_CNTRS);
	    }
	}

	/* Split Domain into two using SubdivCntr. */
#	ifdef DEBUG_PRINT_VERTEX_LIST_SPLIT
	    fprintf(stderr, "Original:\n");
	    PrintVrtxList(Domain -> PVertex);
	    fprintf(stderr, "Subdiv:\n");
	    PrintVrtxList(SubdivCntr -> PVertex);
#	endif /* DEBUG_PRINT_VERTEX_LIST_SPLIT */

	InsideDomain = SplitDomainByCntr(Domain, SubdivCntr,
					 UMin, UMax, VMin, VMax);
	OutsideDomain = InsideDomain -> Pnext;
	InsideDomain -> Pnext = NULL;

#	ifdef DEBUG_PRINT_VERTEX_LIST_SPLIT
	    fprintf(stderr, "Splitted:\n");
	    PrintVrtxList(InsideDomain -> PVertex);
	    PrintVrtxList(OutsideDomain -> PVertex);
#	endif /* DEBUG_PRINT_VERTEX_LIST_SPLIT */

	/* Do the two recursive calls, merge results, and return. */
	if (InsideCntrs == NULL)
	    Cntrs1 = InsideDomain;
	else {
	    Cntrs1 = ClosedCntrsFromOpenCntrs(InsideCntrs, InsideDomain,
					      UMin, UMax, VMin, VMax);
	    IPFreePolygon(InsideDomain);
	}

	if (OutsideCntrs == NULL)
	    Cntrs2 = OutsideDomain;
	else {
	    Cntrs2 = ClosedCntrsFromOpenCntrs(OutsideCntrs, OutsideDomain,
					      UMin, UMax, VMin, VMax);
	    IPFreePolygon(OutsideDomain);
	}

	LastCntr = IritPrsrGetLastPoly(Cntrs1);
	LastCntr -> Pnext = Cntrs2;
	return Cntrs1;
    }
    else {
	IPPolygonStruct *Cntrs;

	/* We have no other contour - split and return the two domains. */

#	ifdef DEBUG_PRINT_VERTEX_LIST_SPLIT
	    PrintVrtxList(Domain -> PVertex);
	    fprintf(stderr, "Subdiv:\n");
	    PrintVrtxList(SubdivCntr -> PVertex);
#	endif /* DEBUG_PRINT_VERTEX_LIST_SPLIT */

	Cntrs = SplitDomainByCntr(Domain, SubdivCntr, UMin, UMax, VMin, VMax);

#	ifdef DEBUG_PRINT_VERTEX_LIST_SPLIT
	    fprintf(stderr, "Splitted:\n");
	    PrintVrtxList(OutsideCntrs -> PVertex);
	    PrintVrtxList(OutsideCntrs -> Pnext -> PVertex);
#	endif /* DEBUG_PRINT_VERTEX_LIST_SPLIT */

	return Cntrs;
    }
}

/*****************************************************************************
* DESCRIPTION:                                                               *
*   Splits into two the given Domain that represents the current trimmed     *
* surface domain, by the given open contour.                                 *
*                                                                            *
* PARAMETERS:                                                                *
*   Domain:      A closed polyline of the current trimmed surface domain.    *
*   Cntr:        An open polyline that has its end points on the boundary of *
*		 Domain.						     *
*   UMin, UMax, VMin, VMax:   Domain of surface.		             *
*                                                                            *
* RETURN VALUE:                                                              *
*   IPPolygonStruct *:   Two sub domains, results of splitting Domain,       *
*		inside sub-domain first, followed by outside sub-domain.     *
*****************************************************************************/
static IPPolygonStruct *SplitDomainByCntr(IPPolygonStruct *Domain,
					  IPPolygonStruct *Cntr,
					  CagdRType UMin,
					  CagdRType UMax,
					  CagdRType VMin,
					  CagdRType VMax)
{
    CagdBType
	ReversedSubdivCntr = FALSE;
    IPVertexStruct *VDomain, *VDmnStart, *VDmnEnd;
    IPVertexStruct
	*VFirst = Cntr -> PVertex,
	*VLast = IritPrsrGetLastVrtx(VFirst);
    CagdRType
        VFirstWeight = VertexWeight(VFirst, UMin, UMax, VMin, VMax),
        VLastWeight = VertexWeight(VLast, UMin, UMax, VMin, VMax);

    if (VFirstWeight > VLastWeight) {
	ReversedSubdivCntr = TRUE;
	SWAP(CagdRType, VFirstWeight, VLastWeight);
    }

    /* Split Domain into two using SubdivCntr. */
    VDomain = CopyVertexList(Domain -> PVertex);
    for (VDmnStart = VDomain;
	 VDmnStart -> Pnext &&
	 (!IS_BOUNDARY_VRTX(VDmnStart -> Pnext) ||
	  VertexWeight(VDmnStart -> Pnext, UMin, UMax, VMin, VMax) <
								VFirstWeight);
	 VDmnStart = VDmnStart -> Pnext);
    for (VDmnEnd = VDmnStart;
	 VDmnEnd -> Pnext != NULL &&
	 (!IS_BOUNDARY_VRTX(VDmnEnd -> Pnext) ||
	  VertexWeight(VDmnEnd -> Pnext, UMin, UMax, VMin, VMax) <
								VLastWeight);
	 VDmnEnd = VDmnEnd -> Pnext);

    if (VDmnStart -> Pnext == NULL ||
	VDmnEnd -> Pnext == NULL ||
	!IS_BOUNDARY_VRTX(VDmnStart) ||
	!IS_BOUNDARY_VRTX(VDmnStart -> Pnext) ||
	!IS_BOUNDARY_VRTX(VDmnEnd) ||
	!IS_BOUNDARY_VRTX(VDmnEnd -> Pnext)) {
	TrimFatalError(TRIM_ERR_INCONSISTENT_CNTRS);
	return NULL;
    }
    else {
	IPPolygonStruct *InsideDomain, *OutsideDomain;

	IPVertexStruct *VLastCopy, *VDmnNext,
	    *VFirstCopy = CopyVertexList(VFirst);

	VFirstCopy = IritPrsrReverseVrtxList2(VFirstCopy);
	VLastCopy = IritPrsrGetLastVrtx(VFirstCopy);

	SET_BOUNDARY_VRTX(VFirst);
	SET_BOUNDARY_VRTX(VLast);
	SET_BOUNDARY_VRTX(VFirstCopy);
	SET_BOUNDARY_VRTX(VLastCopy);

	if (ReversedSubdivCntr) {
	    VDmnNext = VDmnEnd -> Pnext;

	    /* Chain the portion of Domain that is inside. */ 
	    if (VDmnStart != VDmnEnd) {
		VLast -> Pnext = VDmnStart -> Pnext;
		VLast = VDmnEnd;
		VDmnEnd -> Pnext = NULL;
	    }

	    /* Duplicate the last point as first. */
	    VFirst = IPAllocVertex(0, NULL, VFirst);
	    PT_COPY(VFirst -> Coord, VLast -> Coord);
	    SET_BOUNDARY_VRTX(VFirst);
	    InsideDomain = IPAllocPolygon(0, VFirst, NULL);

	    VLastCopy -> Pnext = VDmnNext;
	    VDmnStart -> Pnext = VFirstCopy;
	    OutsideDomain = IPAllocPolygon(0, VDomain, NULL);
	}
	else {
	    VDmnNext = VDmnEnd -> Pnext;

	    /* Chain the portion of Domain that is inside. */ 
	    if (VDmnStart != VDmnEnd) {
		VLastCopy -> Pnext = VDmnStart -> Pnext;
		VLastCopy = VDmnEnd;
		VDmnEnd -> Pnext = NULL;
	    }

	    /* Duplicate the last point as first. */
	    VFirstCopy = IPAllocVertex(0, NULL, VFirstCopy);
	    PT_COPY(VFirstCopy -> Coord, VLastCopy -> Coord);
	    SET_BOUNDARY_VRTX(VFirstCopy);
	    InsideDomain = IPAllocPolygon(0, VFirstCopy, NULL);

	    VLast -> Pnext = VDmnNext;
	    VDmnStart -> Pnext = VFirst;
	    OutsideDomain = IPAllocPolygon(0, VDomain, NULL);
	}

	InsideDomain -> Pnext = OutsideDomain;
	return InsideDomain;
    }
}

/*****************************************************************************
* DESCRIPTION:                                                               *
*   Computes a unique weights according to the location where the vertex     *
* meets the boundary:							     *
*                 VMin is (0-1), UMax is (1-2), VMax is (2-3), UMin is (3-4) *
*                                                                            *
* PARAMETERS:                                                                *
*   V:			      To assign a weight.                            *
*   UMin, UMax, VMin, VMax:   Domain of surface.		             *
*                                                                            *
* RETURN VALUE:                                                              *
*   CagdRType:   Weight of V                                                 *
*****************************************************************************/
static CagdRType VertexWeight(IPVertexStruct *V,
			      CagdRType UMin,
			      CagdRType UMax,
			      CagdRType VMin,
			      CagdRType VMax)
{
    if (APX_EQ_EPS(V -> Coord[0], UMin, CONTOUR_EPS)) {
	return 3.0 + (V -> Coord[1] - VMax) / (VMin - VMax);
    }
    else if (APX_EQ_EPS(V -> Coord[0], UMax, CONTOUR_EPS)) {
	return 1.0 + (V -> Coord[1] - VMin) / (VMax - VMin);
    }
    else if (APX_EQ_EPS(V -> Coord[1], VMin, CONTOUR_EPS)) {
	return 0.0 + (V -> Coord[0] - UMin) / (UMax - UMin);
    }
    else if (APX_EQ_EPS(V -> Coord[1], VMax, CONTOUR_EPS)) {
	return 2.0 + (V -> Coord[0] - UMax) / (UMin - UMax);
    }
    else {
	TrimFatalError(TRIM_ERR_INCONSISTENT_CNTRS);
	return IRIT_INFNTY;
    }    
}

#ifdef DEBUG_PRINT_VERTEX_LIST

/*****************************************************************************
* 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;

    fprintf(stderr, "\nPolyline:\n");

    do {
	fprintf(stderr, "    %12f %12f %12f\n",
	       V -> Coord[0], V -> Coord[1], V -> Coord[2]);
	V = V -> Pnext;
    }
    while (V!= NULL && V != VHead);
}
#endif /* DEBUG_PRINT_VERTEX_LIST */