/******************************************************************************
* SBspEval.c - Bspline surfaces handling routines - evaluation routines.      *
*******************************************************************************
* Written by Gershon Elber, Mar. 90.					      *
******************************************************************************/

#include <string.h>
#include "cagd_loc.h"

/*****************************************************************************
* DESCRIPTION:                                                               M
* Evaluates the given tensor product Bspline surface at a given point, by    M
* extracting an isoparamteric curve along u from the surface and evaluating  M
* the curve at parameter v.                                                  M
*									     M
*		u -->							     V
*     +----------------------+						     V
*     |P0		 Pi-1|						     V
*   V |Pi		P2i-1|	Parametric space orientation - control mesh. V
*    ||			     |						     V
*    v|Pn-i		 Pn-1|						     V
*     +----------------------+						     V
*                                                                            *
* PARAMETERS:                                                                M
*   Srf:       Surface to evaluate at the given (u, v) location.             M
*   u, v:      Location where to evaluate the surface.                       M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdRType *:  A vector holding all the coefficients of all components    M
*                 of curve Crv's point type. If for example the curve's      M
*                 point type is P2, the W, X, and Y will be saved in the     M
*                 first three locations of the returned vector. The first    M
*                 location (index 0) of the returned vector is reserved for  M
*                 the rational coefficient W and XYZ always starts at second M
*                 location of the returned vector (index 1).                 M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspSrfEvalAtParam, evaluation, Bsplines                                  M
*****************************************************************************/
CagdRType *BspSrfEvalAtParam(CagdSrfStruct *Srf, CagdRType u, CagdRType v)
{
    static CagdCrvStruct
	*IsoSubCrv = NULL;
    CagdRType *Pt, *VBasisFunc, UMin, UMax, VMin, VMax;
    CagdBType
	IsNotRational = !CAGD_IS_RATIONAL_SRF(Srf);
    int k, UIndexFirst, VIndexFirst,
	UOrder = Srf -> UOrder,
	VOrder = Srf -> VOrder,
        ULength = Srf -> ULength,
        VLength = Srf -> VLength,
	MaxCoord = CAGD_NUM_OF_PT_COORD(Srf -> PType);

    CagdSrfDomain(Srf, &UMin, &UMax, &VMin, &VMax);
    if (u < UMin - EPSILON || u > UMax + EPSILON)
	CAGD_FATAL_ERROR(CAGD_ERR_U_NOT_IN_SRF);
    if (v < VMin - EPSILON || v > VMax + EPSILON)
	CAGD_FATAL_ERROR(CAGD_ERR_V_NOT_IN_SRF);
    if (u > UMax - IRIT_EPSILON * 2)
	u = UMax - IRIT_EPSILON * 2;
    if (v > VMax - IRIT_EPSILON * 2)
	v = VMax - IRIT_EPSILON * 2;
    if (u < UMin)
        u = UMin;
    if (v < VMin)
        v = VMin;

    UIndexFirst = BspKnotLastIndexLE(Srf -> UKnotVector,
				     (CAGD_SRF_UPT_LST_LEN(Srf) + UOrder), u) -
								  (UOrder - 1);
    VBasisFunc = BspCrvCoxDeBoorBasis(Srf -> VKnotVector, VOrder,
				      VLength +
					 (Srf -> VPeriodic ? VOrder - 1 : 0),
				      v, &VIndexFirst);

    if (IsoSubCrv != NULL &&
	(Srf -> PType != IsoSubCrv -> PType ||
	 UOrder != IsoSubCrv -> Order)) {
	/* The cached curve is not the proper type - release it. */
	CagdCrvFree(IsoSubCrv);
	IsoSubCrv = NULL;
    }
    if (IsoSubCrv == NULL) {
        IsoSubCrv = BspPeriodicCrvNew(UOrder, UOrder, FALSE, Srf -> PType);
    }
    CAGD_GEN_COPY(IsoSubCrv -> KnotVector,
		  &Srf -> UKnotVector[UIndexFirst],
		  sizeof(CagdRType) * UOrder * 2);

    for (k = 0; k < UOrder; k++) {
	int l;

	for (l = IsNotRational; l <= MaxCoord; l++) {
	    int i,
		VIndexFirstTmp = VIndexFirst;
	    CagdRType
	        *SrfP = &Srf -> Points[l][CAGD_MESH_UV(Srf,
						  UIndexFirst,
						  VIndexFirstTmp)],
	        *CrvP = &IsoSubCrv -> Points[l][k];

	    *CrvP = 0.0;
	    for (i = 0; i < VOrder; i++) {
		*CrvP += VBasisFunc[i] * *SrfP;
		if (++VIndexFirstTmp >= VLength) {
		    VIndexFirstTmp -= VLength;
		    SrfP -= VLength * CAGD_NEXT_V(Srf);
		}
		else
		    SrfP += CAGD_NEXT_V(Srf);
	    }
	}

	if (++UIndexFirst >= ULength)
	    UIndexFirst -= ULength;
    }

    Pt = BspCrvEvalAtParam(IsoSubCrv, u);

    return Pt;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Same as BspSrfEvalAtParam above. Cleaner, but much less efficient.	     M
*****************************************************************************/
CagdRType *BspSrfEvalAtParam2(CagdSrfStruct *Srf, CagdRType u, CagdRType v)
{
    CagdRType *Pt;
    CagdCrvStruct
	*IsoCrv = BspSrfCrvFromSrf(Srf, u, CAGD_CONST_U_DIR);

    if (!BspKnotParamInDomain(IsoCrv -> KnotVector, IsoCrv -> Length,
			      IsoCrv -> Order, IsoCrv -> Periodic, v))
	CAGD_FATAL_ERROR(CAGD_ERR_V_NOT_IN_SRF);

    Pt = BspCrvEvalAtParam(IsoCrv, v);

    CagdCrvFree(IsoCrv);

    return Pt;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Extracts an isoparametric curve out of the given tensor product Bspline    M
* surface in direction Dir at the parameter value of t.                      M
*   Operations should prefer the CONST_U_DIR, in which the extraction is     M
* somewhat faster if that is possible.					     M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf:       To extract an isoparametric ocurve from.                      M
*   t:         Parameter value of extracted isoparametric curve.             M
*   dir:       Direction of the isocurve on the surface. Either U or V.      M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:  An isoparametric curve of Srf. This curve inherits the M
*		      order and continuity of surface Srf in direction Dir.  M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspSrfCrvFromSrf, isoparametric curves, curve from surface               M
*****************************************************************************/
CagdCrvStruct *BspSrfCrvFromSrf(CagdSrfStruct *Srf,
				CagdRType t,
				CagdSrfDirType dir)
{
    CagdCrvStruct
	*Crv = NULL;
    CagdBType
	IsNotRational = !CAGD_IS_RATIONAL_SRF(Srf);
    int i, j, CrvLen,
	MaxCoord = CAGD_NUM_OF_PT_COORD(Srf -> PType);
    CagdRType *CrvP, *SrfP;

    switch (dir) {
	case CAGD_CONST_U_DIR:
	    if (!BspKnotParamInDomain(Srf -> UKnotVector, Srf -> ULength,
				      Srf -> UOrder, Srf -> UPeriodic, t))
		CAGD_FATAL_ERROR(CAGD_ERR_U_NOT_IN_SRF);
	    Crv = BspPeriodicCrvNew(CrvLen = Srf -> VLength, Srf -> VOrder,
				    Srf -> VPeriodic, Srf -> PType);
	    CAGD_GEN_COPY(Crv -> KnotVector, Srf -> VKnotVector,
			  sizeof(CagdRType) *
			      (CAGD_SRF_VPT_LST_LEN(Srf) + Srf -> VOrder));

	    for (i = IsNotRational; i <= MaxCoord; i++) {
		CrvP = Crv -> Points[i];
		SrfP = Srf -> Points[i];
		for (j = 0; j < CrvLen; j++) {
		    *CrvP++ = BspCrvEvalVecAtParam(SrfP, CAGD_NEXT_U(Srf),
					Srf -> UKnotVector, Srf -> UOrder,
					Srf -> ULength, Srf -> UPeriodic, t);
			SrfP += CAGD_NEXT_V(Srf);
		}
	    }
	    break;
	case CAGD_CONST_V_DIR:
	    if (!BspKnotParamInDomain(Srf -> VKnotVector, Srf -> VLength,
				      Srf -> VOrder, Srf -> VPeriodic, t))
		CAGD_FATAL_ERROR(CAGD_ERR_V_NOT_IN_SRF);
	    Crv = BspPeriodicCrvNew(CrvLen = Srf -> ULength, Srf -> UOrder,
				    Srf -> UPeriodic, Srf -> PType);
	    CAGD_GEN_COPY(Crv -> KnotVector, Srf -> UKnotVector,
			  sizeof(CagdRType) *
			      (CAGD_SRF_UPT_LST_LEN(Srf) + Srf -> UOrder));

	    for (i = IsNotRational; i <= MaxCoord; i++) {
		CrvP = Crv -> Points[i];
		SrfP = Srf -> Points[i];
		for (j = 0; j < CrvLen; j++) {
		    *CrvP++ = BspCrvEvalVecAtParam(SrfP, CAGD_NEXT_V(Srf),
					Srf -> VKnotVector, Srf -> VOrder,
					Srf -> VLength, Srf -> VPeriodic, t);
		    SrfP += CAGD_NEXT_U(Srf);
		}
	    }
	    break;
	default:
	    CAGD_FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV);
	    break;
    }
    return Crv;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Extracts a curve from the mesh of a tensor product Bspline surface Srf in  M
* direction Dir at index Index.                                              M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf:       To extract a curve from.  		                     M
*   Index:     Index along the mesh of Srf to extract the curve from.        M
*   Dir:       Direction of extracted curve. Either U or V.		     M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:   A curve from Srf. This curve inherit the order and    M
*                      continuity of surface Srf in direction Dir. However,  M
*                      thiscurve is not on surface Srf, in general.          M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspSrfCrvFromMesh, isoparametric curves, curve from mesh                 M
*****************************************************************************/
CagdCrvStruct *BspSrfCrvFromMesh(CagdSrfStruct *Srf,
				 int Index,
				 CagdSrfDirType Dir)
{
    CagdCrvStruct
	*Crv = NULL;
    CagdBType
	IsNotRational = !CAGD_IS_RATIONAL_SRF(Srf);
    int i, j, CrvLen,
	MaxCoord = CAGD_NUM_OF_PT_COORD(Srf -> PType);
    CagdRType *CrvP, *SrfP;

    switch (Dir) {
	case CAGD_CONST_U_DIR:
	    if (Index >= CAGD_SRF_UPT_LST_LEN(Srf) || Index < 0)
		CAGD_FATAL_ERROR(CAGD_ERR_INDEX_NOT_IN_MESH);

	    Index %= Srf -> ULength;
	    Crv = BspPeriodicCrvNew(CrvLen = Srf -> VLength, Srf -> VOrder,
				    Srf -> VPeriodic, Srf -> PType);
	    CAGD_GEN_COPY(Crv -> KnotVector, Srf -> VKnotVector,
			  sizeof(CagdRType) *
			      (CAGD_SRF_VPT_LST_LEN(Srf) + Srf -> VOrder));

	    for (i = IsNotRational; i <= MaxCoord; i++) {
		CrvP = Crv -> Points[i];
		SrfP = Srf -> Points[i] + Index * CAGD_NEXT_U(Srf);
		for (j = 0; j < CrvLen; j++) {
		    *CrvP++ = *SrfP;
		    SrfP += CAGD_NEXT_V(Srf);
		}
	    }
	    break;
	case CAGD_CONST_V_DIR:
	    if (Index >= CAGD_SRF_VPT_LST_LEN(Srf) || Index < 0)
		CAGD_FATAL_ERROR(CAGD_ERR_INDEX_NOT_IN_MESH);

	    Index %= Srf -> VLength;
	    Crv = BspPeriodicCrvNew(CrvLen = Srf -> ULength, Srf -> UOrder,
				    Srf -> UPeriodic, Srf -> PType);
	    CAGD_GEN_COPY(Crv -> KnotVector, Srf -> UKnotVector,
			  sizeof(CagdRType) *
			      (CAGD_SRF_UPT_LST_LEN(Srf) + Srf -> UOrder));

	    for (i = IsNotRational; i <= MaxCoord; i++) {
		CrvP = Crv -> Points[i];
		SrfP = Srf -> Points[i] + Index * CAGD_NEXT_V(Srf);
		for (j = 0; j < CrvLen; j++) {
		    *CrvP++ = *SrfP;
		    SrfP += CAGD_NEXT_U(Srf);
		}
	    }
	    break;
	default:
	    CAGD_FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV);
	    break;
    }
    return Crv;
}

