/******************************************************************************
* Bsp2Poly.c - Bezier to polygon/polylines conversion routines.		      *
*******************************************************************************
* Written by Gershon Elber, Mar. 90.					      *
******************************************************************************/

#include "cagd_loc.h"

static CagdPolygonStruct *BspC1Srf2Polygons(CagdSrfStruct *Srf, int FineNess,
			 CagdBType ComputeNormals, CagdBType FourPerFlat);

/*****************************************************************************
* Routine to convert a single bspline surface to set of triangles	     *
* approximating it. FineNess is a finess control on result and the bigger it *
* is more triangles may result. a value of 10 is a good start value.	     *
* NULL is returned in case of an error, otherwise list of CagdPolygonStruct. *
*   This routine looks for C1 discontinuities in the surface and split it    *
* into C1 continuous patches to invoke BspC1Srf2Polygons to gen. polygons.   *
*****************************************************************************/
CagdPolygonStruct *BspSrf2Polygons(CagdSrfStruct *Srf, int FineNess,
			 CagdBType ComputeNormals, CagdBType FourPerFlat)
{
    CagdRType u, v;
    int UOrder = Srf -> UOrder,
	VOrder = Srf -> VOrder,
	ULength = Srf -> ULength,
	VLength = Srf -> VLength;
    CagdBType
	HasUDiscont = BspKnotC1Discont(Srf -> UKnotVector, UOrder, ULength, &u),
	HasVDiscont = BspKnotC1Discont(Srf -> VKnotVector, VOrder, VLength, &v);
    CagdPolygonStruct *Poly;

    if (HasUDiscont || HasVDiscont) {
	CagdSrfStruct
	    *Srf1 = HasUDiscont ? BspSrfSubdivAtParam(Srf, u, CAGD_CONST_U_DIR)
				: BspSrfSubdivAtParam(Srf, v, CAGD_CONST_V_DIR),
	    *Srf2 = Srf1 -> Pnext;
	CagdPolygonStruct
	    *Poly1 = BspSrf2Polygons(Srf1, FineNess, ComputeNormals, FourPerFlat),
	    *Poly2 = BspSrf2Polygons(Srf2, FineNess, ComputeNormals, FourPerFlat);

	CagdSrfFreeList(Srf1);

	/* Chain the two lists together: */
	for (Poly = Poly1; Poly -> Pnext != NULL; Poly = Poly -> Pnext);
	Poly -> Pnext = Poly2;
	Poly = Poly1;
    }
    else
	Poly = BspC1Srf2Polygons(Srf, FineNess, ComputeNormals, FourPerFlat);

    return Poly;
}

/*****************************************************************************
* Routine to convert a single C1 continuouse bspline srf to set of triangles *
* approximating it. FineNess is a finess control on result and the bigger it *
* is more triangles may result. a value of 10 is a good start value.	     *
* NULL is returned in case of an error, otherwise list of CagdPolygonStruct. *
*****************************************************************************/
static CagdPolygonStruct *BspC1Srf2Polygons(CagdSrfStruct *Srf, int FineNess,
			 CagdBType ComputeNormals, CagdBType FourPerFlat)
{
    int i, j, FineNessU1, FineNessV1, FineNessU, FineNessV, BaseIndex;
    CagdRType u, v, UMin, UMax, VMin, VMax, *Pt;
    CagdPointType
	PType = Srf -> PType;
    CagdPtStruct PtCenter, *Pt1, *Pt2, *Pt3, *Pt4, *PtMesh, *PtMeshPtr;
    CagdVecStruct NlCenter, *Nl1, *Nl2, *Nl3, *Nl4, *PtNrml;
    CagdCrvStruct *Crv;
    CagdPolygonStruct *Poly,
	*PolyHead = NULL;

    if (!CAGD_IS_BSPLINE_SRF(Srf)) return NULL;

    /* Simple heuristic to estimate how many samples to compute. */
    FineNessU = Srf -> ULength * FineNess / 10;
    FineNessV = Srf -> VLength * FineNess / 10;


    if (FineNessU < 2) FineNessU = 2;
    if (FineNessV < 2) FineNessV = 2;
    switch (_CagdLin2Poly) {
	case CAGD_REG_POLY_PER_LIN:
	    break;
    	case CAGD_ONE_POLY_PER_LIN:
	    if (Srf -> UOrder == 2) FineNessU = 2;
	    if (Srf -> VOrder == 2) FineNessV = 2;
	    break;
    	case CAGD_ONE_POLY_PER_COLIN:
	    break;
    }
    FineNessU1 = FineNessU - 1;
    FineNessV1 = FineNessV - 1;

    /* Current to surface property such as curvature is used as subdivison   */
    /* criterion and the surface is subdivided, equally spaced in parametric */
    /* space, using FineNess as number of subdivisions per axis.	     */

    /* Allocate a mesh to hold all vertices so common vertices need not be   */
    /* Evaluated twice, and evaluate the surface at these mesh points.	     */
    PtMeshPtr = PtMesh = (CagdPtStruct *) CagdMalloc(FineNessU * FineNessV *
							sizeof(CagdPtStruct));

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

    for (i = 0; i < FineNessU; i++) {
	u = UMin + (UMax - UMin) * i / ((CagdRType) FineNessU1);
	Crv = BspSrfCrvFromSrf(Srf, u, CAGD_CONST_U_DIR);

	for (j = 0; j < FineNessV; j++) {
	    v = VMin + (VMax - VMin) * j / ((CagdRType) FineNessV1);
	    Pt = BspCrvEvalAtParam(Crv, v);
	    CagdCoerceToE3(PtMeshPtr -> Pt, &Pt, -1, PType);
	    PtMeshPtr++;
	}

	CagdCrvFree(Crv);
    }

    if (ComputeNormals) PtNrml = BspSrfMeshNormals(Srf, FineNessU, FineNessV);

    /* Now that we have the mesh, create the polygons. */
    for (i = 0; i < FineNessU1; i++)
	for (j = 0; j < FineNessV1; j++) {
	    BaseIndex = i * FineNessV + j;
	    Pt1 = &PtMesh[BaseIndex];        /* Cache the four flat corners. */
	    Pt2 = &PtMesh[BaseIndex + 1];
	    Pt3 = &PtMesh[BaseIndex + FineNessV + 1];
	    Pt4 = &PtMesh[BaseIndex + FineNessV];
	    if (ComputeNormals) {
		Nl1 = &PtNrml[BaseIndex];
		Nl2 = &PtNrml[BaseIndex + 1];
		Nl3 = &PtNrml[BaseIndex + FineNessV + 1];
		Nl4 = &PtNrml[BaseIndex + FineNessV];
	    }

	    if (FourPerFlat) {  /* Eval middle point and create 4 triangles. */
		CAGD_COPY_POINT(PtCenter, *Pt1);
		CAGD_ADD_POINT(PtCenter, *Pt2);
		CAGD_ADD_POINT(PtCenter, *Pt3);
		CAGD_ADD_POINT(PtCenter, *Pt4);
		CAGD_MULT_POINT(PtCenter, 0.25);

		if (ComputeNormals) {
		    /* Average the four normals to find the middle one. */
		    CAGD_COPY_VECTOR(NlCenter, *Nl1);
		    CAGD_ADD_VECTOR(NlCenter, *Nl2);
		    CAGD_ADD_VECTOR(NlCenter, *Nl3);
		    CAGD_ADD_VECTOR(NlCenter, *Nl4);
		    CAGD_NORMALIZE_VECTOR(NlCenter);
		}

		Poly = _CagdMakePolygon(ComputeNormals, Pt1, Pt2, &PtCenter,
							Nl1, Nl2, &NlCenter);
		CAGD_LIST_PUSH(Poly, PolyHead);
		Poly = _CagdMakePolygon(ComputeNormals, Pt2, Pt3, &PtCenter,
						        Nl2, Nl3, &NlCenter);
		CAGD_LIST_PUSH(Poly, PolyHead);
		Poly = _CagdMakePolygon(ComputeNormals, Pt3, Pt4, &PtCenter,
						        Nl3, Nl4, &NlCenter);
		CAGD_LIST_PUSH(Poly, PolyHead);
		Poly = _CagdMakePolygon(ComputeNormals, Pt4, Pt1, &PtCenter,
						        Nl4, Nl1, &NlCenter);
		CAGD_LIST_PUSH(Poly, PolyHead);
	    }
	    else {			   /* Only two along the diagonal... */
		Poly = _CagdMakePolygon(ComputeNormals, Pt1, Pt2, Pt3,
						        Nl1, Nl2, Nl3);
		CAGD_LIST_PUSH(Poly, PolyHead);
		Poly = _CagdMakePolygon(ComputeNormals, Pt3, Pt4, Pt1,
						        Nl3, Nl4, Nl1);
		CAGD_LIST_PUSH(Poly, PolyHead);
	    }
	}

    CagdFree((VoidPtr) PtMesh);
    if (ComputeNormals) CagdFree((VoidPtr) PtNrml);

    return PolyHead;
}

/*****************************************************************************
* Routine to convert a single bspline surface to NumOfIsolines polylines list*
* in each param. direction with SamplesPerCurve in each isoparametric curve. *
* Polyline are always E3 of CagdPolylineStruct type.			     *
* Iso parametric curves are sampled equally spaced in parametric space.	     *
* NULL is returned in case of an error, otherwise list of CagdPolylineStruct.*
* Attempt is made to extract isolines along C1 discontinuities first.	     *
*****************************************************************************/
CagdPolylineStruct *BspSrf2Polylines(CagdSrfStruct *Srf, int NumOfIsocurves,
							 int SamplesPerCurve)
{
    int i, NumC1Disconts;
    CagdRType u, v, UMin, UMax, VMin, VMax, *C1Disconts, *IsoParams;
    CagdCrvStruct *Crv;
    CagdPolylineStruct *Poly,
	*PolyList = NULL;

    if (!CAGD_IS_BSPLINE_SRF(Srf)) return NULL;

    /* Make sure requested format is something reasonable. */
    if (SamplesPerCurve < 1) SamplesPerCurve = 1;
    if (SamplesPerCurve > CAGD_MAX_BEZIER_CACHE_ORDER)
	SamplesPerCurve = CAGD_MAX_BEZIER_CACHE_ORDER;
    if (NumOfIsocurves < 2) NumOfIsocurves = 2;

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

    /* Compute discontinuities along the u axis and use that to determine    */
    /* where to extract isolines along u.				     */
    C1Disconts = BspKnotAllC1Discont(Srf -> UKnotVector, Srf -> UOrder,
				     Srf -> ULength, &NumC1Disconts);
    IsoParams = BspKnotParamValues(UMin, UMax, NumOfIsocurves, C1Disconts,
								NumC1Disconts);
    for (i = 0; i < NumOfIsocurves; i++) {
	u = IsoParams[i];

	Crv = BspSrfCrvFromSrf(Srf, u, CAGD_CONST_U_DIR);
	Poly = BspCrv2Polyline(Crv, SamplesPerCurve);
	Poly -> Pnext = PolyList;
	PolyList = Poly;
	CagdCrvFree(Crv);
    }
    CagdFree((VoidPtr) IsoParams);

    /* Compute discontinuities along the v axis and use that to determine    */
    /* where to extract isolines along v.				     */
    C1Disconts = BspKnotAllC1Discont(Srf -> VKnotVector, Srf -> VOrder,
				     Srf -> VLength, &NumC1Disconts);
    IsoParams = BspKnotParamValues(VMin, VMax, NumOfIsocurves, C1Disconts,
								NumC1Disconts);
    for (i = 0; i < NumOfIsocurves; i++) {
	v = IsoParams[i];

	Crv = BspSrfCrvFromSrf(Srf, v, CAGD_CONST_V_DIR);
	Poly = BspCrv2Polyline(Crv, SamplesPerCurve);
	Poly -> Pnext = PolyList;
	PolyList = Poly;
	CagdCrvFree(Crv);
    }
    CagdFree((VoidPtr) IsoParams);

    return PolyList;
}

/*****************************************************************************
* Routine to approx. a single bspline curve as a polyline with		     *
* 2^SamplesPerCurve samples. Polyline is always E3 CagdPolylineStruct type.  *
* Curve is sampled equally spaced in parametric space, unless the curve is   *
* linear in which the control polygon is simply being copied.		     *
* NULL is returned in case of an error, otherwise CagdPolylineStruct.	     *
*****************************************************************************/
CagdPolylineStruct *BspCrv2Polyline(CagdCrvStruct *Crv, int SamplesPerCurve)
{
    int i, j, n,
	IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv),
	MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType);
    CagdRType *Polyline[CAGD_MAX_PT_SIZE], Scaler;
    CagdPtStruct *NewPolyline;
    CagdPolylineStruct *P;

    if (!CAGD_IS_BSPLINE_CRV(Crv)) return NULL;

    /* Make sure requested format is something reasonable. */
    if (SamplesPerCurve < 1) SamplesPerCurve = 1;
    if (SamplesPerCurve > CAGD_MAX_BEZIER_CACHE_ORDER)
	SamplesPerCurve = CAGD_MAX_BEZIER_CACHE_ORDER;
    if (Crv -> Order == 2 && (1 << SamplesPerCurve) < Crv -> Length) {
    	/* Make sure 2^SamplesPerCurve can hold the entire control polygon. */
    	for (i = 1, SamplesPerCurve = 0;
	     i < Crv -> Length;
	     i <<= 1, SamplesPerCurve++);
    }

    n = 1 << SamplesPerCurve;

    P = CagdPolylineNew(n);
    NewPolyline = P -> Polyline;

    /* Allocate temporary memory to hold evaluated curve. */
    for (i = 0; i < CAGD_MAX_PT_SIZE; i++)
	Polyline[i] = (CagdRType *) CagdMalloc(sizeof(CagdRType) * n);

    if (MaxCoord > 3) MaxCoord = 3;

    n = P -> Length = BspCrvEvalToPolyline(Crv, SamplesPerCurve, Polyline);

    for (i = n - 1; i >= 0; i--) {	          /* Convert to E3 polyline. */
	if (IsNotRational)
	    Scaler = 1.0;
	else
	    Scaler = Polyline[0][i];

	for (j = 0; j < MaxCoord; j++)
	   NewPolyline[i].Pt[j] = Polyline[j+1][i] / Scaler;
	for (j = MaxCoord; j < 3; j++)
	   NewPolyline[i].Pt[j] = 0.0;
    }
    for (i = 0; i < CAGD_MAX_PT_SIZE; i++) CagdFree((VoidPtr) Polyline[i]);

    return P;
}