/******************************************************************************
* FFCnvHul.c - computation of convex hull of a freeform curve.		      *
*******************************************************************************
* (C) Gershon Elber, Technion, Israel Institute of Technology                 *
*******************************************************************************
* Written by Gershon Elber, January 96  				      *
******************************************************************************/

#include "symb_loc.h"
#include "user_lib.h"
#include "iritprsr.h"
#include "allocate.h"

#define PLANE_Z_EPS		1e-6
#define CNTR_DIAGONAL_EPS	1e-3

static PlaneType
    GlblPlane = { 1.0, 0.0, 0.0, PLANE_Z_EPS }; /* Used to contour surfaces. */

static CagdCrvStruct *ProcessContours(CagdCrvStruct *Crv,
				      IPPolygonStruct *Cntrs);
static void InsertDomainIntoList(CagdPtStruct **PtList,
				 CagdPtStruct *PtPrev,
				 CagdRType TMin,
				 CagdRType TMax);

/*****************************************************************************
* DESCRIPTION:                                                               M
*   Computes the convex hull of a C1 freeform planar curve, in the           M
* XY plane.  The convex hull is computed by symbolically isolating the non   M
* negative set (in t) of:						     M
*                                                                            M
*		C'(t) x (C(r) - C(t)) >= 0,   forall r.			     V
*                                                                            M
* Note the above equation yields a scalar value since C(t) is planar. The    M
* resulting set in t contains all the subdomain in C(t) that is on the       M
* convex hull of C(t).  Connecting these pieces with straight lines yeilds   M
* the final convex hull curve.						     M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:       To compute its convex hull.                                   M
*   FineNess:  Of numeric search for the zero set (for surface subdivision). M
*	       A positive value (10 is a good start).			     M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:   A curve representing the convex hull of Crv.          M
*                                                                            *
* SEE ALSO:                                                                  M
*   SymbCrvPtTangents, SymbCrvDiameter                                       M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbCrvCnvxHull                                                          M
*****************************************************************************/
CagdCrvStruct *SymbCrvCnvxHull(CagdCrvStruct *Crv, RealType FineNess)
{
    int IsClosed = CagdIsClosedCrv(Crv);
    IPPolygonStruct *Cntrs;
    CagdCrvStruct *DCrv, *CHCrv, *CrvDomainsInCH;
    CagdSrfStruct *Srf, *DSrf, *SrfT, *SrfR, *Srf1, *Srf2,
	*SrfW, *SrfX, *SrfY, *SrfZ, *DSrfW, *DSrfX, *DSrfY, *DSrfZ;
    CagdPtStruct Pt;

    /* Make sure the given curve is open end conditioned curve. */
    if (CAGD_IS_BEZIER_CRV(Crv))
	Crv = CnvrtBezier2BsplineCrv(Crv);
    else
        Crv = CagdCrvCopy(Crv);

    if (CAGD_IS_PERIODIC_CRV(Crv)) {
        CagdCrvStruct
	    *TCrv = CnvrtPeriodic2FloatCrv(Crv);

	CagdCrvFree(Crv);
	Crv = TCrv;
    }

    if (CAGD_IS_BSPLINE_CRV(Crv) && !BspCrvHasOpenEC(Crv)) {
	CagdCrvStruct
	    *TCrv = BspCrvOpenEnd(Crv);

	CagdCrvFree(Crv);
	Crv = TCrv;
    }

    /* Make sure the domain is zero to one. */
    BspKnotAffineTrans(Crv -> KnotVector,
		       Crv -> Length + Crv -> Order,
		       -Crv -> KnotVector[0],
		       1.0 / (Crv -> KnotVector[Crv -> Length +
						Crv -> Order - 1] -
			      Crv -> KnotVector[0]));

    DCrv = CagdCrvDerive(Crv);
    DSrf = CagdPromoteCrvToSrf(DCrv, CAGD_CONST_U_DIR);
    CagdCrvFree(DCrv);
    SymbSrfSplitScalar(DSrf, &DSrfW, &DSrfX, &DSrfY, &DSrfZ);
    CagdSrfFree(DSrf);
    if (DSrfW != NULL)
	CagdSrfFree(DSrfW);
    if (DSrfZ != NULL)
	CagdSrfFree(DSrfZ);

    SrfT = CagdPromoteCrvToSrf(Crv, CAGD_CONST_U_DIR);
    SrfR = CagdPromoteCrvToSrf(Crv, CAGD_CONST_V_DIR);

    Srf1 = SymbSrfSub(SrfR, SrfT);
    CagdSrfFree(SrfR);
    CagdSrfFree(SrfT);
    SymbSrfSplitScalar(Srf1, &SrfW, &SrfX, &SrfY, &SrfZ);
    CagdSrfFree(Srf1);
    if (SrfW != NULL)
	CagdSrfFree(SrfW);
    if (SrfZ != NULL)
	CagdSrfFree(SrfZ);

    Srf1 = SymbSrfMult(SrfX, DSrfY);
    CagdSrfFree(SrfX);
    CagdSrfFree(DSrfY);
    Srf2 = SymbSrfMult(SrfY, DSrfX);
    CagdSrfFree(SrfY);
    CagdSrfFree(DSrfX);
    Srf = SymbSrfSub(Srf1, Srf2);
    CagdSrfFree(Srf1);
    CagdSrfFree(Srf2);

#ifdef DUMP_CH_SYMBOLIC_FUNC
    IritPrsrStdoutObject(GenSRFObject(Srf));        /* Also a memory leak... */
#endif /* DUMP_CH_SYMBOLIC_FUNC */

    Cntrs = UserCntrSrfWithPlane(Srf, GlblPlane, FineNess);

    CagdSrfFree(Srf);

#ifdef DUMP_CH_CONTOURS
    IritPrsrStdoutObject(GenPOLYObject(Cntrs));     /* Also a memory leak... */
#endif /* DUMP_CH_CONTOURS */

    CrvDomainsInCH = ProcessContours(Crv, Cntrs);
    IPFreePolygonList(Cntrs);

    CHCrv = CagdMergeCrvList(CrvDomainsInCH, TRUE);
    CagdCrvFreeList(CrvDomainsInCH);

    /* Make sure the end points are indeed closed. */
    if (IsClosed) {
	CagdCrvStruct *TCrv;

	CagdCoerceToE3(Pt.Pt, Crv -> Points, 0, Crv -> PType);
	TCrv = CagdMergeCrvPt(CHCrv, &Pt);
	CagdCrvFree(CHCrv);
	CHCrv = TCrv;
    }

    CagdCrvFree(Crv);

    return CHCrv;
}

/*****************************************************************************
* DESCRIPTION:                                                               *
*   Computes the u domain of the given Cntrs that signals the domain that is *
* not in the CH. the rest of the domain is hence, the the CH.  Then, the     *
* proper sub-domains that are in CH are extracted from Crv and returned as a *
* list.					                                     *
*                                                                            *
* PARAMETERS:                                                                *
*   Crv:      To extract the domain of it on the CH.                         *
*   Cntrs:    the (almost) zero set of the symbolic functional. The contours *
*	      has the Height (Z), U, and V as the XYZ in this order.	     *
*                                                                            *
* RETURN VALUE:                                                              *
*   CagdCrvStruct *:  A list of sub curves of Crv that are on the CH.        *
*****************************************************************************/
static CagdCrvStruct *ProcessContours(CagdCrvStruct *Crv,
				      IPPolygonStruct *Cntrs)
{
    IPPolygonStruct *Cntr;
    CagdCrvStruct
	*CrvList = NULL;
    CagdPtStruct *Pt,
	*InCHDomain = CagdPtNew();

    /* We start with the entire zero to one domain as in the CH. */
    InCHDomain -> Pt[0] = 0.0;
    InCHDomain -> Pt[1] = 1.0;

    for (Cntr = Cntrs; Cntr != NULL; Cntr = Cntr -> Pnext) {
	IPVertexStruct *V;
	CagdRType
	    TMin = IRIT_INFNTY,
	    TMax = -IRIT_INFNTY;

	for (V = Cntr -> PVertex; V != NULL; V = V -> Pnext) {
	    /* Consider only points that are off the diagonal. */
	    if (!APX_EQ_EPS(V -> Coord[1], V -> Coord[2], CNTR_DIAGONAL_EPS)) {
		if (TMin > V -> Coord[1])
		    TMin = V -> Coord[1];
		if (TMax < V -> Coord[1])
		    TMax = V -> Coord[1];
	    }
	}

	if (TMin < TMax) {
#ifdef DUMP_CH_CNTR_DOMAINS
	    fprintf(stderr, "TMin = %f, TMax = %f\n", TMin, TMax);
#endif /* DUMP_CH_CNTR_DOMAINS */

	    InsertDomainIntoList(&InCHDomain, NULL, TMin, TMax);
	}
    }

    for (Pt = InCHDomain; Pt != NULL; Pt = Pt -> Pnext) {
#ifdef DUMP_CH_CNTR_DOMAINS
	fprintf(stderr, "Valid domain at [%f : %f]\n",
		Pt -> Pt[0], Pt -> Pt[1]);
#endif /* DUMP_CH_CNTR_DOMAINS */

	if (!APX_EQ(Pt -> Pt[0], Pt -> Pt[1])) {
	    CagdCrvStruct
	        *TCrv = CagdCrvRegionFromCrv(Crv, Pt -> Pt[0], Pt -> Pt[1]);

#ifdef DUMP_CH_CRV_DOMAINS
	    IritPrsrStdoutObject(GenCRVObject(TCrv));/* Also a memory leak...*/
#endif /* DUMP_CH_CRV_DOMAINS */

	    TCrv -> Pnext = CrvList;
	    CrvList = TCrv;
	}
    }

    return CagdListReverse(CrvList);
}

/*****************************************************************************
* DESCRIPTION:                                                               *
*   Inserts one new invalid domain from TMin to TMax to the list of valid    *
* domains PtList.			                                     *
*                                                                            *
* PARAMETERS:                                                                *
*   PtList:      Current list of valid domains, to update in place.          *
*   PtPrev:	 Previous item in list. Internally only. NULL from outside.  *
*   TMin, TMax:  New invalid domain to insert.                               *
*                                                                            *
* RETURN VALUE:                                                              *
*   void                                                                     *
*****************************************************************************/
static void InsertDomainIntoList(CagdPtStruct **PtList,
				 CagdPtStruct *PtPrev,
				 CagdRType TMin,
				 CagdRType TMax)
{
    CagdPtStruct *Pt;

    if (TMax - TMin < IRIT_EPS)
	return;

    for (Pt = *PtList;
	 Pt != NULL && TMin > Pt -> Pt[1];
	 PtPrev = Pt, Pt = Pt -> Pnext);

    if (Pt != NULL && TMax > Pt -> Pt[0]) {
	if (TMin <= Pt -> Pt[0] && TMax >= Pt -> Pt[1]) {
	    /* This domain becomes invalid: */
	    if (PtPrev != NULL) {
		PtPrev -> Pnext = Pt -> Pnext;
		CagdPtFree(Pt);
		Pt = PtPrev -> Pnext;
		InsertDomainIntoList(&Pt, PtPrev, TMin, TMax);
	    }
	    else {
		*PtList = (*PtList) -> Pnext;
		CagdPtFree(Pt);
		InsertDomainIntoList(PtList, NULL, TMin, TMax);
	    }
	}
	else if (TMin >= Pt -> Pt[0] && TMax <= Pt -> Pt[1]) {
	    CagdPtStruct
		*PtTemp = CagdPtNew();

	    /* New domain is completely contained in this domain. */
	    PtTemp -> Pnext = Pt -> Pnext;
	    Pt -> Pnext = PtTemp;

	    PtTemp -> Pt[0] = TMax;
	    PtTemp -> Pt[1] = Pt -> Pt[1];
	    PtTemp -> Pt[2] = 0.0;

	    Pt -> Pt[1] = TMin;
	}
	else if (TMin <= Pt -> Pt[0] && TMax < Pt -> Pt[1]) {
	    Pt -> Pt[0] = TMax;
	}
	else if (TMin > Pt -> Pt[0] && TMax >= Pt -> Pt[1]) {
	    Pt -> Pt[1] = TMin;
	    /* We also might affect the objects following us. */
	    InsertDomainIntoList(&Pt -> Pnext, Pt, TMin, TMax);
	}
    }
}

/*****************************************************************************
* DESCRIPTION:                                                               M
*   Given a freeform curve, compute its diameter as a function.	             M
* If the curve is a convex, probably as a result of a convex hull	     M
* computation of an original curve, the matching will be one to one.	     M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:       A curve to process its diameter function.                     M
*   FineNess:  Of numeric search for the zero set (for surface subdivision). M
*	       A positive value (10 is a good start).			     M
*                                                                            *
* RETURN VALUE:                                                              M
*   IPPolygonStruct *:  Contours of the matched parallel lines on Crv. Each  M
*		vertex will hold two parameter values on Crv.		     M
*                                                                            *
* SEE ALSO:                                                                  M
*   SymbCrvCnvxHull, SymbCrvPtTangents, SymbCrvDiameterMinMax                M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbCrvDiameter                                                          M
*****************************************************************************/
struct IPPolygonStruct *SymbCrvDiameter(CagdCrvStruct *Crv, RealType FineNess)
{
    IPPolygonStruct *Cntrs, *Cntr, *PrevC;
    CagdCrvStruct
	*DCrv = CagdCrvDerive(Crv);
    CagdSrfStruct *Srf, *Srf1, *Srf2, *SrfT,*SrfR,
	*SrfRW, *SrfRX, *SrfRY, *SrfRZ, *SrfTW, *SrfTX, *SrfTY, *SrfTZ;

    /* Make sure the domain is zero to one. */
    BspKnotAffineTrans(DCrv -> KnotVector,
		       DCrv -> Length + DCrv -> Order,
		       -DCrv -> KnotVector[0],
		       1.0 / (DCrv -> KnotVector[DCrv -> Length +
						 DCrv -> Order - 1] -
			      DCrv -> KnotVector[0]));

    SrfT = CagdPromoteCrvToSrf(DCrv, CAGD_CONST_U_DIR),
    SrfR = CagdPromoteCrvToSrf(DCrv, CAGD_CONST_V_DIR);
    CagdCrvFree(DCrv);

    SymbSrfSplitScalar(SrfT, &SrfTW, &SrfTX, &SrfTY, &SrfTZ);
    CagdSrfFree(SrfT);
    if (SrfTW != NULL)
	CagdSrfFree(SrfTW);
    if (SrfTZ != NULL)
	CagdSrfFree(SrfTZ);

    SymbSrfSplitScalar(SrfR, &SrfRW, &SrfRX, &SrfRY, &SrfRZ);
    CagdSrfFree(SrfR);
    if (SrfRW != NULL)
	CagdSrfFree(SrfRW);
    if (SrfRZ != NULL)
	CagdSrfFree(SrfRZ);

    Srf1 = SymbSrfMult(SrfTX, SrfRY);
    CagdSrfFree(SrfTX);
    CagdSrfFree(SrfRY);
    Srf2 = SymbSrfMult(SrfTY, SrfRX);
    CagdSrfFree(SrfTY);
    CagdSrfFree(SrfRX);
    Srf = SymbSrfSub(Srf1, Srf2);
    CagdSrfFree(Srf1);
    CagdSrfFree(Srf2);

#ifdef DUMP_DIAMETER_SYMBOLIC_FUNC
    IritPrsrStdoutObject(GenSRFObject(Srf));        /* Also a memory leak... */
#endif /* DUMP_DIAMETER_SYMBOLIC_FUNC */

    Cntrs = UserCntrSrfWithPlane(Srf, GlblPlane, FineNess);

    CagdSrfFree(Srf);

    /* Eliminate points below or on the diagonal. */
    for (Cntr = Cntrs; Cntr != NULL; Cntr = Cntr -> Pnext) {
	IPVertexStruct *V, *PrevV;
	CagdRType TMin, TMax, Dt;

	CagdCrvDomain(Crv, &TMin, &TMax);
	Dt = TMax - TMin;

	for (V = Cntr -> PVertex, PrevV = NULL; V != NULL; ) {
	    /* Consider only points that are off the diagonal. */
	    if (APX_EQ_EPS(V -> Coord[1], V -> Coord[2], CNTR_DIAGONAL_EPS) ||
		V -> Coord[1] < V -> Coord[2]) {
		/* Delete this vertex. */
		if (PrevV != NULL) {
		    V = V -> Pnext;
		    IPFreeVertex(PrevV -> Pnext);
		    PrevV -> Pnext = V;
		}
		else {
		    Cntr -> PVertex = V -> Pnext;
		    IPFreeVertex(V);
		    V = Cntr -> PVertex;
		}
	    }
	    else {
		V -> Coord[0] = TMin + Dt * V -> Coord[1];
		V -> Coord[0] = BOUND(V -> Coord[0], TMin, TMax);
		V -> Coord[1] = TMin + Dt * V -> Coord[2];
		V -> Coord[1] = BOUND(V -> Coord[1], TMin, TMax);
		V -> Coord[2] = 0.0;
		PrevV = V;
		V = V -> Pnext;
	    }
	}
    }

    /* Eliminate empty contours. */
    for (Cntr = Cntrs, PrevC = NULL; Cntr != NULL; ) {
	if (Cntr -> PVertex == NULL) {
	    if (PrevC != NULL) {
		Cntr = Cntr -> Pnext;
		IPFreePolygon(PrevC -> Pnext);
		PrevC -> Pnext = Cntr;
	    }
	    else {
		Cntrs = Cntrs -> Pnext;
		IPFreePolygon(Cntr);
		Cntr = Cntrs;
	    }
	}
	else {
	    PrevC = Cntr;
	    Cntr = Cntr -> Pnext;
	}
    }

#ifdef DUMP_DIAMETER_CONTOURS
    IritPrsrStdoutObject(GenPOLYObject(Cntrs));     /* Also a memory leak... */
#endif /* DUMP_DIAMETER_CONTOURS */

    return Cntrs;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
*   Computes the maximum or minimum diameter out of diameter matched list    M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:   To compute its diameter.                                          M
*   Cntrs: Output of SymbCrvDiameter - the matched paraller tangents.        M
*   Min:   TRUE of minimum diameter, FALSE for maximum diameter.	     M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdPype *:  Two parameter values on Crv of tangent lines extreme value. M
*		 Returns an address to a statically allocated point.	     M
*                                                                            *
* SEE ALSO:                                                                  M
*   SymbCrvDiameter                                                          M
*                                                                            *
* KEYWORDS:                                                                  M
*   SymbCrvDiameterMinMax                                                    M
*****************************************************************************/
CagdRType *SymbCrvDiameterMinMax(CagdCrvStruct *Crv,
				 struct IPPolygonStruct *Cntrs,
				 int Min)
{
    static CagdRType Params[2];
    IPPolygonStruct *Cntr;
    CagdRType
	ExtremeDist = Min ? IRIT_INFNTY : -IRIT_INFNTY;

    for (Cntr = Cntrs; Cntr != NULL; Cntr = Cntr -> Pnext) {
	IPVertexStruct *V;

	for (V = Cntr -> PVertex; V != NULL; V = V -> Pnext) {
	    CagdRType Dist, *R;
	    CagdPType Pt1, Pt2;

	    R = CagdCrvEval(Crv, V -> Coord[0]);
	    CagdCoerceToE3(Pt1, &R, -1, Crv -> PType);
	    R = CagdCrvEval(Crv, V -> Coord[1]);
	    CagdCoerceToE3(Pt2, &R, -1, Crv -> PType);

	    Dist = PT_PT_DIST(Pt1, Pt2);
	    if (Min) {
		if (Dist < ExtremeDist) {
		    ExtremeDist = Dist;
		    PT_COPY(Params, V -> Coord);
		}
	    }
	    else {
		if (Dist > ExtremeDist) {
		    ExtremeDist = Dist;
		    PT_COPY(Params, V -> Coord);
		}
	    }
	}
    }

    return Params;
}