/******************************************************************************
* Trim_sub.c - subdivision of trimmed surfaces.                               *
*******************************************************************************
* (C) Gershon Elber, Technion, Israel Institute of Technology                 *
*******************************************************************************
* Written by Gershon Elber, June 95.					      *
******************************************************************************/

#include "trim_loc.h"

#define SUBDIV_TCRV_EPSILON	1e-3
#define PERTURB_SUBDIV		3.01060123456789e-6

static void ClipTrimmingCurves(CagdCrvStruct *Crvs, CagdSrfStruct *Srf);
static void SplitTrimmingCurve(CagdCrvStruct *Crv,
			       CagdPtStruct *Pts,
			       int Axis,
			       CagdRType t,
			       CagdCrvStruct **UVCrv1,
			       CagdCrvStruct **UVCrv2);
static void MergeTrimmingCurves(CagdCrvStruct **Crvs, int OAxis);
static void MergeExactEndPoints(CagdCrvStruct **Crvs);
static CagdCrvStruct *GetMinimalCurve(CagdCrvStruct **Crvs,
				      int OAxis,
				      int *StartPt);

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given a trimmed surface - subdivides it into two sub-surfaces at given     M
* parametric value t in the given direction Dir.                             M
*    Returns pointer to a list of two trimmed surfaces, at most. It can very M
* well may happen that the subdivided surface is completely trimmed out and  M
* hence nothing is returned for it.					     M
*                                                                            *
* PARAMETERS:                                                                M
*   TrimSrf:  To subdivide at the prescibed parameter value t.               M
*   t:        The parameter to subdivide the curve Crv at.                   M
*   Dir:      Direction of subdivision. Either U or V.                       M
*                                                                            *
* RETURN VALUE:                                                              M
*   TrimSrfStruct *:  The subdivided surfaces. Usually two, but can have     M
*		      only one, if other is totally trimmed away.            M
*                                                                            *
* KEYWORDS:                                                                  M
*   TrimSrfSubdivAtParam, subdivision                                        M
*****************************************************************************/
TrimSrfStruct *TrimSrfSubdivAtParam(TrimSrfStruct *TrimSrf,
				    CagdRType t,
				    CagdSrfDirType Dir)
{
    int Axis, OAxis;
    CagdRType UMin, UMax, VMin, VMax;
    CagdCrvStruct *Crv,
	*UVCrv1All = NULL,
	*UVCrv2All = NULL;
    CagdSrfStruct *Srf, *DividedSrf;
    TrimCrvStruct *TrimCrvTmp,
	*TrimCrv1 = NULL,
	*TrimCrv2 = NULL,
	*TrimCrv = TrimSrf -> TrimCrvList;
    TrimSrfStruct *TSrf1, *TSrf2;
    CagdLType Line;

    /* Coerce Bezier surfaces to Bspline so we do not have to scale the */
    /* trimming curves to the unit domain all the time.		        */
    Srf = TrimSrf -> Srf;
    if (CAGD_IS_BEZIER_SRF(Srf))
	Srf = CnvrtBezier2BsplineSrf(Srf);
    DividedSrf = CagdSrfSubdivAtParam(Srf, t, Dir);
    if (Srf != TrimSrf -> Srf)
	CagdSrfFree(Srf);

    /* Compute a line to subdivide along and represent it as a curve. */
    TrimSrfDomain(TrimSrf, &UMin, &UMax, &VMin, &VMax);
    switch (Dir) {
	case CAGD_CONST_U_DIR:
	    Axis = 0;
	    OAxis = 1;
	    Line[0] = 1;
	    Line[1] = 0;
	    Line[2] = PERTURB_SUBDIV - t;
	    break;
	case CAGD_CONST_V_DIR:
	    Axis = 1;
	    OAxis = 0;
	    Line[0] = 0;
	    Line[1] = 1;
	    Line[2] = PERTURB_SUBDIV - t;
	    break;
	default:
	    TRIM_FATAL_ERROR(TRIM_ERR_DIR_NOT_CONST_UV);
	    return NULL;
    }

    /* Compute the hard part - subdividing the trimming curves into two. */
    for( ; TrimCrv != NULL; TrimCrv = TrimCrv -> Pnext) {
	TrimCrvSegStruct
	    *TrimCrvSeg = TrimCrv -> TrimCrvSegList;

	for ( ; TrimCrvSeg != NULL; TrimCrvSeg = TrimCrvSeg -> Pnext) {
	    CagdCrvStruct
		*UVCrv = TrimCrvSeg -> UVCrv;
	    CagdPtStruct
		*Pts = SymbLclDistCrvLine(UVCrv, Line, SUBDIV_TCRV_EPSILON,
					  TRUE, FALSE);
	    CagdUVType UV;

	    if (Pts == NULL) {/* No intersection - classify to one of halfs. */
		CagdCoerceToE2(UV, UVCrv -> Points, 0, UVCrv -> PType);
		if (UV[Axis] > t) {
		    /* Goes to second half. */
		    TrimCrvTmp = TrimCrvNew(TrimCrvSegCopy(TrimCrvSeg));
		    TrimCrvTmp -> Pnext = TrimCrv2;
		    TrimCrv2 = TrimCrvTmp;
		}
		else {
		    /* Goes to first half. */
		    TrimCrvTmp = TrimCrvNew(TrimCrvSegCopy(TrimCrvSeg));
		    TrimCrvTmp -> Pnext = TrimCrv1;
		    TrimCrv1 = TrimCrvTmp;
		}
	    }
	    else { /* Needs to split the curve at the intersections. */
		CagdCrvStruct *UVCrv1, *UVCrv2;

		SplitTrimmingCurve(UVCrv, Pts, Axis, t, &UVCrv1, &UVCrv2);
		if (UVCrv1 != NULL) {
		    Crv = CagdListLast(UVCrv1);
		    Crv -> Pnext = UVCrv1All;
		    UVCrv1All = UVCrv1;
		}
		if (UVCrv2 != NULL) {
		    Crv = CagdListLast(UVCrv2);
		    Crv -> Pnext = UVCrv2All;
		    UVCrv2All = UVCrv2;
		}

		CagdPtFreeList(Pts);
	    }
	}
    }

    /* Merge all the resulting segments into closed loops and form segments. */
    MergeTrimmingCurves(&UVCrv1All, OAxis);
    MergeTrimmingCurves(&UVCrv2All, OAxis);
    ClipTrimmingCurves(UVCrv1All, DividedSrf);
    ClipTrimmingCurves(UVCrv2All, DividedSrf -> Pnext);
	
    for (Crv = UVCrv1All; Crv != NULL; ) {
	CagdCrvStruct
	    *NextCrv = Crv -> Pnext;

	TrimCrvTmp = TrimCrvNew(TrimCrvSegNew(Crv, NULL));
	Crv -> Pnext = NULL;

	TrimCrvTmp -> Pnext = TrimCrv1;
	TrimCrv1 = TrimCrvTmp;

	Crv = NextCrv;
    }
    for (Crv = UVCrv2All; Crv != NULL; ) {
	CagdCrvStruct
	    *NextCrv = Crv -> Pnext;

	TrimCrvTmp = TrimCrvNew(TrimCrvSegNew(Crv, NULL));
	Crv -> Pnext = NULL;

	TrimCrvTmp -> Pnext = TrimCrv2;
	TrimCrv2 = TrimCrvTmp;

	Crv = NextCrv;
    }

    /* Construct the (upto) two new trimmed surfaces and return them. */
    if (TrimCrv1 == NULL) {
	TSrf1 = TrimSrfNew(DividedSrf -> Pnext, TrimCrv2, TRUE);
	CagdSrfFree(DividedSrf);
    }
    else if (TrimCrv2 == NULL) {
	TSrf1 = TrimSrfNew(DividedSrf, TrimCrv1, TRUE),
	CagdSrfFree(DividedSrf -> Pnext);
	DividedSrf -> Pnext = NULL;
    }
    else {
	TSrf1 = TrimSrfNew(DividedSrf, TrimCrv1, TRUE),
	TSrf2 = TrimSrfNew(DividedSrf -> Pnext, TrimCrv2, TRUE);
	DividedSrf -> Pnext = NULL;
	TSrf1 -> Pnext = TSrf2;
    }

    return TSrf1;
}

/*****************************************************************************
* DESCRIPTION:                                                               *
*   Clip the given set of curve to remain inside the parameteric space       *
* of surface Srf by moving them inside if found outside.                     *
*                                                                            *
* PARAMETERS:                                                                *
*   Crvs:     To clip into the parametric space of Srf.                      *
*   Srf:      that should be using Crvs as trimming curves.                  *
*                                                                            *
* RETURN VALUE:                                                              *
*   void                                                                     *
*****************************************************************************/
static void ClipTrimmingCurves(CagdCrvStruct *Crvs, CagdSrfStruct *Srf)
{
    CagdRType UMin, UMax, VMin, VMax;
    CagdCrvStruct *Crv;

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

    for (Crv = Crvs; Crv != NULL; Crv = Crv -> Pnext) {
	int i;
	CagdRType
	    **Points = Crv -> Points;
	CagdPointType
	    PType = Crv -> PType;
	CagdBType
	    IsRational = CAGD_IS_RATIONAL_PT(PType);

	for (i = 0; i < Crv -> Length; i++) {
	    CagdPType Pt;

	    CagdCoerceToE2(Pt, Points, i, PType);

	    if (Pt[0] < UMin)
		Points[1][i] = IsRational ? UMin * Points[0][i] : UMin;
	    if (Pt[0] > UMax)
		Points[1][i] = IsRational ? UMax * Points[0][i] : UMax;
	    if (Pt[1] < VMin)
		Points[2][i] = IsRational ? VMin * Points[0][i] : VMin;
	    if (Pt[1] > VMax)
		Points[2][i] = IsRational ? VMax * Points[0][i] : VMax;
	}
    }
}

/*****************************************************************************
* DESCRIPTION:                                                               *
*   Splits one trimming curve at the t value in axis Axis. The split pieces  *
* are placed in either UVCrv1 (< t) or UVCrv2 (>t) while being merged into   *
* new loops.								     *
*                                                                            *
* PARAMETERS:                                                                *
*   UVCrv:   Curve to split at t in axis Axis.                               *
*   Pts:     Parameter values at which to split the UVCrv curve.             *
*   Axis:    Axis that is considered (either 0 for U or 1 for V).            *
*   t:       Parameter value of surface at which the trimming curve is split.*
*   UVCrv1:  To hold the <t side of the trimming curve.                      *
*   UVCrv2:  To hold the >t side of the trimming curve.                      *
*                                                                            *
* RETURN VALUE:                                                              *
*   void                                                                     *
*****************************************************************************/
static void SplitTrimmingCurve(CagdCrvStruct *UVCrv,
			       CagdPtStruct *Pts,
			       int Axis,
			       CagdRType t,
			       CagdCrvStruct **UVCrv1,
			       CagdCrvStruct **UVCrv2)
{
    CagdRType *R, TMin, TMax;
    CagdPtStruct *Pt;
    CagdUVType UV;

    *UVCrv1 = *UVCrv2 = NULL;
    UVCrv = CagdCrvCopy(UVCrv);

    for (Pt = Pts; Pt != NULL; Pt = Pt -> Pnext) {
	CagdCrvStruct
	    *Crv1 = CagdCrvSubdivAtParam(UVCrv, Pt -> Pt[0]),
	    *Crv2 = Crv1 -> Pnext;

	Crv1 -> Pnext = NULL;

	/* Make sure intersection point is exactly on the line. */
	if (CAGD_IS_RATIONAL_CRV(UVCrv)) {
	    Crv1 -> Points[Axis + 1][Crv1 -> Length - 1] =
		t * Crv1 -> Points[0][Crv1 -> Length - 1];
	    Crv2 -> Points[Axis + 1][0] = t * Crv2 -> Points[0][0];
	}
	else {
	    Crv1 -> Points[Axis + 1][Crv1 -> Length - 1] = t;
	    Crv2 -> Points[Axis + 1][0] = t;
	}

	CagdCrvDomain(Crv1, &TMin, &TMax);
	R = CagdCrvEval(Crv1, (TMin + TMax) / 2.0);
	CagdCoerceToE2(UV, &R, -1, Crv1 -> PType);
	if (UV[Axis] < t) {
	    Crv1 -> Pnext = *UVCrv1;
	    *UVCrv1 = Crv1;
	}
	else {
	    Crv1 -> Pnext = *UVCrv2;
	    *UVCrv2 = Crv1;
	}

	CagdCrvFree(UVCrv);
	UVCrv = Crv2;
    }

    /* Classify the last segment. */
    CagdCrvDomain(UVCrv, &TMin, &TMax);
    R = CagdCrvEval(UVCrv, (TMin + TMax) / 2.0);
    CagdCoerceToE2(UV, &R, -1, UVCrv -> PType);
    if (UV[Axis] < t) {
	UVCrv -> Pnext = *UVCrv1;
	*UVCrv1 = UVCrv;
    }
    else {
	UVCrv -> Pnext = *UVCrv2;
	*UVCrv2 = UVCrv;
    }

    /* In case the opening of the closed loop did not occur on subdivision */
    /* value, which is very likely, we must hook together the opening.     */
    MergeExactEndPoints(UVCrv1);
    MergeExactEndPoints(UVCrv2);
}

/*****************************************************************************
* DESCRIPTION:                                                               *
*   Merges the curves into loops by adding segments along the OAxis which is *
* where they were split as part of the trimmed surface sunbdivision process. *
*                                                                            *
* PARAMETERS:                                                                *
*   Crvs:    To merge together, in place.                                    *
*   OAxis:   Other axis where Crvs were split along (0 or 1 for U or V).     *
*                                                                            *
* RETURN VALUE:                                                              *
*   void                                                                     *
*****************************************************************************/
static void MergeTrimmingCurves(CagdCrvStruct **Crvs, int OAxis)
{
    CagdCrvStruct
	*NewCrvs = NULL;

    while (*Crvs != NULL) {
	CagdBType StartPt1, StartPt2;
	CagdCrvStruct *Crv,
	    *MinCrv1 = GetMinimalCurve(Crvs, OAxis, &StartPt1),
	    *MinCrv2 = GetMinimalCurve(Crvs, OAxis, &StartPt2);

	if (MinCrv2 == NULL) {
	    CagdPtStruct Pt;

	    /* We have only one curve. Make MinCrv1 closed. */
	    CagdCoerceToE3(Pt.Pt, MinCrv1 -> Points, 0, MinCrv1 -> PType);
	    Crv = CagdMergeCrvPt(MinCrv1, &Pt);
	    CagdCrvFree(MinCrv1);
	    Crv -> Pnext = NewCrvs;
	    NewCrvs = Crv;
	}
	else {
	    CagdCrvStruct *TCrv;
	    CagdRType
		Crv1Max = MinCrv1 -> Points[OAxis + 1]
					   [StartPt1 ? MinCrv1 -> Length - 1
						     : 0],
		Crv2Min = MinCrv2 -> Points[OAxis + 1]
					   [StartPt2 ? 0
						     : MinCrv2 -> Length - 1];

	    if (Crv1Max > Crv2Min) {
		/* Marge MinCrv1 and MinCrv2 and push result back in. */
		if (StartPt1 && StartPt2) {
		    TCrv = CagdCrvReverse(MinCrv1);
		    Crv = CagdMergeCrvCrv(TCrv, MinCrv2, TRUE);
		    CagdCrvFree(TCrv);
		}
		else if (!StartPt1 && !StartPt2) {
		    TCrv = CagdCrvReverse(MinCrv2);
		    Crv = CagdMergeCrvCrv(MinCrv1, TCrv, TRUE);
		    CagdCrvFree(TCrv);
		}
		else if (StartPt1) {
		    Crv = CagdMergeCrvCrv(MinCrv2, MinCrv1, TRUE);
		}
		else {
		    Crv = CagdMergeCrvCrv(MinCrv1, MinCrv2, TRUE);
		}

		CagdCrvFree(MinCrv1);
		CagdCrvFree(MinCrv2);
		Crv -> Pnext = *Crvs;
		*Crvs = Crv;
	    }
	    else {
		CagdPtStruct Pt;

		/* Marge MinCrv1 with itself and push MinCrv2 back in. */
		CagdCoerceToE3(Pt.Pt, MinCrv1 -> Points, 0, MinCrv1 -> PType);
		Crv = CagdMergeCrvPt(MinCrv1, &Pt);
		CagdCrvFree(MinCrv1);
		Crv -> Pnext = NewCrvs;
		NewCrvs = Crv;

		MinCrv2 -> Pnext = *Crvs;
		*Crvs = MinCrv2;
	    }
	}
    }

    *Crvs = NewCrvs;
}

/*****************************************************************************
* DESCRIPTION:                                                               *
*   Merges, in place, curves by detecting identical end points. 	     *
*                                                                            *
* PARAMETERS:                                                                *
*   Crvs:    To merge together, in place.                                    *
*                                                                            *
* RETURN VALUE:                                                              *
*   void                                                                     *
*****************************************************************************/
static void MergeExactEndPoints(CagdCrvStruct **Crvs)
{
    CagdCrvStruct *Crv2,
	*NewCrvs = NULL;

    while (*Crvs) {
	CagdPType Pt1Start, Pt1End;
	CagdCrvStruct *Crv = NULL,
	    *Crv1 = *Crvs;

	*Crvs = (*Crvs) -> Pnext;
	Crv1 -> Pnext = NULL;
	CagdCoerceToE3(Pt1Start, Crv1 -> Points, 0, Crv1 -> PType);
	CagdCoerceToE3(Pt1End, Crv1 -> Points,
		       Crv1 -> Length - 1, Crv1 -> PType);

	if (AttrGetIntAttrib(Crv1 -> Attr, "_used") == TRUE) {
	    CagdCrvFree(Crv1);
	    continue;
	}

	for (Crv2 = *Crvs; Crv2 != NULL; Crv2 = Crv2 -> Pnext) {
	    CagdPType Pt2Start, Pt2End;

	    if (AttrGetIntAttrib(Crv2 -> Attr, "_used") == TRUE)
	        continue;

	    CagdCoerceToE3(Pt2Start, Crv2 -> Points, 0, Crv2 -> PType);
	    CagdCoerceToE3(Pt2End, Crv2 -> Points,
			   Crv2 -> Length - 1, Crv2 -> PType);

	    if (PT_APX_EQ(Pt1Start, Pt2Start)) {
		CagdCrvStruct
		    *TCrv = CagdCrvReverse(Crv1);

		Crv = CagdMergeCrvCrv(TCrv, Crv2, TRUE);
		CagdCrvFree(Crv1);
		CagdCrvFree(TCrv);
		Crv1 = Crv;
	    }
	    else if (PT_APX_EQ(Pt1Start, Pt2End)) {
		Crv = CagdMergeCrvCrv(Crv2, Crv1, TRUE);
		CagdCrvFree(Crv1);
		Crv1 = Crv;
	    }
	    else if (PT_APX_EQ(Pt1End, Pt2Start)) {
		Crv = CagdMergeCrvCrv(Crv1, Crv2, TRUE);
		CagdCrvFree(Crv1);
		Crv1 = Crv;
	    }
	    else if (PT_APX_EQ(Pt1End, Pt2End)) {
		CagdCrvStruct
		    *TCrv = CagdCrvReverse(Crv2);

		Crv = CagdMergeCrvCrv(Crv1, TCrv, TRUE);
		CagdCrvFree(Crv1);
		CagdCrvFree(TCrv);
		Crv1 = Crv;
	    }

	    if (Crv != NULL) {			   /* We merged two curves. */
		AttrSetIntAttrib(&Crv2 -> Attr, "_used", TRUE);
		Crv2 = *Crvs;
		Crv = NULL;
	    }
	}

	Crv1 -> Pnext = NewCrvs;
	NewCrvs = Crv1;
    }

    *Crvs = NewCrvs;
}

/*****************************************************************************
* DESCRIPTION:                                                               *
*   Finds the curve with the minimal entry and remove from list.             *
*                                                                            *
* PARAMETERS:                                                                *
*   Crvs:     List of curves to search.                                      *
*   OAxis:    Axis to check(0 or 1 for U or V) .                             *
*   StartPt:  If this is starting point (TRUE) or end point (FALSE).         *
*                                                                            *
* RETURN VALUE:                                                              *
*   CagdCrvStruct *:   Curve with minimal end point                          *
*****************************************************************************/
static CagdCrvStruct *GetMinimalCurve(CagdCrvStruct **Crvs,
				      int OAxis,
				      int *StartPt)
{
    CagdRType
	MinVal = INFINITY;
    CagdCrvStruct *Crv,
        *MinCrv = NULL;

    /* Compute the minimal location. */
    for (Crv = *Crvs; Crv != NULL; Crv = Crv -> Pnext) {
	CagdRType
	    *Pts = Crv -> Points[OAxis + 1];

	if (Pts[0] < MinVal) {
	    MinVal = Pts[0];
	    MinCrv = Crv;
	    *StartPt = TRUE;
	}
	if (Pts[Crv -> Length - 1] < MinVal) {
	    MinVal = Pts[Crv -> Length - 1];
	    MinCrv = Crv;
	    *StartPt = FALSE;
	}
    }

    if (MinCrv == NULL)
	return NULL;

    /* Found the minimal location - remove from list. */
    if (MinCrv == *Crvs) {
	*Crvs = (*Crvs) -> Pnext;
    }
    else {
	for (Crv = *Crvs; Crv -> Pnext != MinCrv; Crv = Crv -> Pnext);
	Crv -> Pnext = MinCrv -> Pnext;
    }
    MinCrv -> Pnext = NULL;

    return MinCrv;
}