/******************************************************************************
* Trim2Ply.c - Converts a trimmed surface into polygons.                      *
*******************************************************************************
* (C) Gershon Elber, Technion, Israel Institute of Technology                 *
*******************************************************************************
* Written by Gershon Elber, June. 95.					      *
******************************************************************************/

#include "trim_loc.h"
#include "allocate.h"
#include "iritprsr.h"
#include "geomat3d.h"
#include "convex.h"

#define COLINEAR_EPSILON 1e-6

static CagdBType GlblMessageReported;
static CagdPolygonStruct
    *GlblPolys = NULL;

static void TrimSrf2PolygonsAux(TrimSrfStruct *TrimSrf,
				int LogFineNess,
				CagdBType ComputeNormals,
				CagdBType ComputeUV);
static void ConvertTrimmingLoopToPolys(CagdSrfStruct *Srf,
				       TrimCrvStruct *UVTCrvs,
				       CagdBType ComputeNormals,
				       CagdBType ComputeUV);

/*****************************************************************************
* DESCRIPTION:                                                               M
* Routine to convert a single trimmed surface to set of triangles            M
* approximating it. 							     M
*   FineNess is a fineness control on result and the larger it is more       M
* triangles may result.							     M
*   A value of 10 is a good starting value.				     M
* NULL is returned in case of an error, otherwise list of CagdPolygonStruct. M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf:              To approximate into triangles.                         M
*   FineNess:         Control on accuracy, the higher the finer.             M
*   ComputeNormals:   If TRUE, normal information is also computed.          M
*   ComputeUV:        If TRUE, UV values are stored and returned as well.    M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdPolygonStruct *:   A list of polygons with optional normal and/or    M
*                         UV parametric information.                         M
*                         NULL is returned in case of an error.              M
*                                                                            *
* KEYWORDS:                                                                  M
*   TrimSrf2Polygons, polygonization, surface approximation                  M
*****************************************************************************/
CagdPolygonStruct *TrimSrf2Polygons(TrimSrfStruct *TrimSrf,
				    int FineNess,
				    CagdBType ComputeNormals,
				    CagdBType ComputeUV)
{
    int i,
	LogFineNess = 0;

    for (i = 1; i < FineNess; i *= 2, LogFineNess++);

    GlblMessageReported = FALSE;
    GlblPolys = NULL;

    TrimSrf2PolygonsAux(TrimSrf, LogFineNess, ComputeNormals, ComputeUV);

    return GlblPolys;
}

/*****************************************************************************
* DESCRIPTION:                                                               *
*   Auxiliary function of TrimSrf2Polygons which does the recursive	     *
* traversal. Only parameter different is LogFineNess which sets the depth of *
* the recursion. Polygons are placed on the global list GlblPolys.           *
*                                                                            *
* PARAMETERS:                                                                *
*   Srf:              To approximate into triangles.                         *
*   FineNess:         Control on accuracy, the higher the finer.             *
*   ComputeNormals:   If TRUE, normal information is also computed.          *
*   ComputeUV:        If TRUE, UV values are stored and returned as well.    *
*                                                                            *
* RETURN VALUE:                                                              *
*   void                                                                     *
*****************************************************************************/
static void TrimSrf2PolygonsAux(TrimSrfStruct *TrimSrf,
				int LogFineNess,
				CagdBType ComputeNormals,
				CagdBType ComputeUV)
{
    TrimSrfStruct *SubdivTSrfs;
    CagdSrfStruct *Srf;
    CagdRType t, UMin, UMax, VMin, VMax;
    
    TrimSrfDomain(TrimSrf, &UMin, &UMax, &VMin, &VMax);

#ifdef DEBUG_PRINT_TRIMSRF_DOMAINS
    fprintf(stderr, "Trimming surface domain [%g  %g], [%g  %g]\n",
	    UMin, UMax, VMin, VMax);
#endif /* DEBUG_PRINT_TRIMSRF_DOMAINS */

    if (LogFineNess < 1) {
	/* Convert into polygons. */
	if (TrimSrf -> TrimCrvList &&
	    TrimSrf -> TrimCrvList -> Pnext &&
	    !GlblMessageReported) {
	    fprintf(stderr,
		    "More than one trimming curve at the lowest level.\n");
	    GlblMessageReported = TRUE;
	}

	ConvertTrimmingLoopToPolys(TrimSrf -> Srf, TrimSrf -> TrimCrvList,
				   ComputeNormals, ComputeUV);
    }
    else {
	Srf = TrimSrf -> Srf;

	if ((LogFineNess & 0x01) == 1) {
	    /* Subdivide in the U direction. */
	    if (CAGD_IS_BEZIER_SRF(Srf) || Srf -> UOrder == Srf -> ULength)
		t = (UMin + UMax) / 2;
	    else
	        t = Srf -> UKnotVector[(Srf -> ULength + Srf -> UOrder) / 2];

	    SubdivTSrfs = TrimSrfSubdivAtParam(TrimSrf, t, CAGD_CONST_U_DIR);
	}
	else {
	    /* Subdivide in the V direction.*/
	    if (CAGD_IS_BEZIER_SRF(Srf) || Srf -> VOrder == Srf -> VLength)
	        t = (VMin + VMax) / 2;
	    else
	        t = Srf -> VKnotVector[(Srf -> VLength + Srf -> VOrder) / 2];

	    SubdivTSrfs = TrimSrfSubdivAtParam(TrimSrf, t, CAGD_CONST_V_DIR);
	}

	TrimSrf2PolygonsAux(SubdivTSrfs, LogFineNess - 1,
			    ComputeNormals, ComputeUV);
	if (SubdivTSrfs -> Pnext != NULL)
	    TrimSrf2PolygonsAux(SubdivTSrfs -> Pnext, LogFineNess - 1,
				ComputeNormals, ComputeUV);
    }
}

/*****************************************************************************
* DESCRIPTION:                                                               *
*   Converts the UV trimming curve UVCrv into piecewise linear, convert the  *
* piecewise linear loop into triangles and constructs a list of polygons.    *
*                                                                            *
*                                                                            *
* PARAMETERS:                                                                *
*   Srf:             That UVCrv is its parametric space trimming curve.      *
*   UVTCrvs:         Curves in UV space of Srf, that also trims Srf.         *
*   ComputeNormals:  If TRUE, normals are to be computed as well.            *
*   ComputeUV:       If TRUE, UV values are to be saved as well.             *
*                                                                            *
* RETURN VALUE:                                                              *
*   void								     *
*****************************************************************************/
static void ConvertTrimmingLoopToPolys(CagdSrfStruct *Srf,
				       TrimCrvStruct *UVTCrvs,
				       CagdBType ComputeNormals,
				       CagdBType ComputeUV)
{
    int i,
	OldCircVal = IritPrsrSetPolyListCirc(TRUE);
    IPVertexStruct *V, *VLast,
	*VHead = NULL;
    IPPolygonStruct *Pl, *ConvexPls;
    CagdCrvStruct
	*UVCrv = UVTCrvs -> TrimCrvSegList -> UVCrv;
    CagdPolylineStruct
	*PolyLine = SymbCrv2Polyline(UVCrv, UVCrv -> Length * 2,
				     FALSE, TRUE);
    CagdPolylnStruct
	*PolyLn = PolyLine -> Polyline;

    if (UVTCrvs -> TrimCrvSegList -> Pnext != NULL) {
	fprintf(stderr, "A trimming loop using several curves is not supported for polygonization\n");
    }

    /* Convert the PolyLine into IPPolygonStruct format. */
    for (i = 0; i < PolyLine -> Length; i++) {
        VHead = IPAllocVertex(0, 0, NULL, VHead);

	VHead -> Normal[0] = VHead -> Normal[1] = 0.0;
	VHead -> Normal[2] = 1.0;

	VHead -> Coord[0] = PolyLn[i].Pt[0];
	VHead -> Coord[1] = PolyLn[i].Pt[1];
	VHead -> Coord[2] = 0.0;
    }

    /* Filter out collinear vertices and make into a circular list. */
    for (V = VHead; V -> Pnext -> Pnext != NULL; ) {
	CagdVType VCross, V12, V23;
	CagdRType l;

	PT_SUB(V12, V -> Coord, V -> Pnext -> Coord);
	PT_SUB(V23, V -> Pnext -> Coord, V -> Pnext -> Pnext -> Coord);

	GMVecCrossProd(VCross, V12, V23);

	l = PT_LENGTH(VCross);

	if (APX_EQ_EPS(l, 0.0, COLINEAR_EPSILON)) {
	    IPVertexStruct
		*VTmp = V -> Pnext;

	    V -> Pnext = VTmp -> Pnext;

	    IPFreeVertex(VTmp);
	}
	else
	    V = V -> Pnext;
    }
    IPFreeVertex(V -> Pnext);
    V -> Pnext = VHead;

    Pl = IPAllocPolygon(0, 0, VHead, NULL);
    IritPrsrReverseVrtxList(Pl);
    IP_RST_CONVEX_POLY(Pl);
    IritPrsrUpdatePolyPlane(Pl);

#ifdef DEBUG_TRIM_POLY
    fprintf(stderr, "Trimming polygon:\n");
    V = Pl -> PVertex;
    do {
	fprintf(stderr, "\t%10.7g %10.7g\n",
		V -> Coord[0], V -> Coord[1]);
	V = V -> Pnext;
    }
    while (V != Pl -> PVertex && V != NULL);
#endif /* DEBUG_TRIM_POLY */

    if (ConvexPolygon(Pl)) {
	ConvexPls = Pl;
    }
    else {
	ConvexPls = SplitNonConvexPoly(Pl);
	IPFreePolygon(Pl);
    }

    IritPrsrSetPolyListCirc(OldCircVal);

    for (Pl = ConvexPls; Pl != NULL; Pl = Pl -> Pnext) {
	CagdRType *R;
	CagdPType VHeadPt, VLastPt;
	CagdVType VHeadNrml, VLastNrml;

	VHead = Pl -> PVertex;
	VLast = VHead -> Pnext;

	R = CagdSrfEval(Srf, VHead -> Coord[0], VHead -> Coord[1]);
	CagdCoerceToE3(VHeadPt, &R, -1, Srf -> PType);
	if (ComputeNormals) {
	    CagdVecStruct
	      *Nrml = CagdSrfNormal(Srf, VHead -> Coord[0], VHead -> Coord[1]);

	    PT_COPY(VHeadNrml, Nrml -> Vec);
	}

	R = CagdSrfEval(Srf, VLast -> Coord[0], VLast -> Coord[1]);
	CagdCoerceToE3(VLastPt, &R, -1, Srf -> PType);
	if (ComputeNormals) {
	    CagdVecStruct
	        *Nrml = CagdSrfNormal(Srf,
				      VLast -> Coord[0], VLast -> Coord[1]);

	    PT_COPY(VLastNrml, Nrml -> Vec);
	}

	for (V = VLast -> Pnext;
	     V != NULL && V != VHead;
	     VLast = V, V = V -> Pnext) {
	    CagdPolygonStruct
	        *CagdPl = CagdPolygonNew();

	    /* Form a triangle between VHead, V, and V -> Pnext. */

	    PT_COPY(CagdPl -> Polygon[0].Pt, VHeadPt);
	    PT_COPY(CagdPl -> Polygon[1].Pt, VLastPt);
	    R = CagdSrfEval(Srf, V -> Coord[0], V -> Coord[1]);
	    CagdCoerceToE3(CagdPl -> Polygon[2].Pt, &R, -1, Srf -> PType);

	    if (ComputeNormals) {
		CagdVecStruct
		    *Nrml = CagdSrfNormal(Srf, V -> Coord[0], V -> Coord[1]);

		PT_COPY(CagdPl -> Polygon[0].Nrml, VHeadNrml);
		PT_COPY(CagdPl -> Polygon[1].Nrml, VLastNrml);
		PT_COPY(CagdPl -> Polygon[2].Nrml, Nrml -> Vec);
	    }
	    if (ComputeUV) {
		CagdPl -> Polygon[0].UV[0] = VHead -> Coord[0];
		CagdPl -> Polygon[0].UV[1] = VHead -> Coord[1];
		CagdPl -> Polygon[1].UV[0] = VLast -> Coord[0];
		CagdPl -> Polygon[1].UV[1] = VLast -> Coord[1];
		CagdPl -> Polygon[2].UV[0] = V -> Coord[0];
		CagdPl -> Polygon[2].UV[1] = V -> Coord[1];
	    }

	    PT_COPY(VLastPt, CagdPl -> Polygon[2].Pt);
	    PT_COPY(VLastNrml, CagdPl -> Polygon[2].Nrml);
	    
	    LIST_PUSH(CagdPl, GlblPolys);
	}
    }

    IPFreePolygonList(ConvexPls);
}

