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

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

/******************************************************************************
* Evaluate the given tensor product Bspline surface at a given point, by      *
* extracting an isoparamteric curve along u and evaluating v in it.	      *
*									      *
*		u -->							      *
*     +----------------------+						      *
*     |P0		 Pi-1|						      *
*   V |Pi		P2i-1|	Parametric space orientation - control mesh.  *
*    ||			     |						      *
*    v|Pn-i		 Pn-1|						      *
*     +----------------------+						      *
*									      *
******************************************************************************/
CagdRType *BspSrfEvalAtParam(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, v))
	FATAL_ERROR(CAGD_ERR_V_NOT_IN_SRF);

    Pt = BspCrvEvalAtParam(IsoCrv, v);

    CagdCrvFree(IsoCrv);

    return Pt;
}

/******************************************************************************
* Extract an isoline curve out of the given tensor product Bspline surface.   *
* Operations should prefer the CONST_U_DIR, in which the extraction is	      *
* somewhat faster if it is possible.					      *
******************************************************************************/
CagdCrvStruct *BspSrfCrvFromSrf(CagdSrfStruct *Srf, CagdRType t,
							CagdSrfDirType dir)
{
    CagdCrvStruct *Crv;
    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, t))
		FATAL_ERROR(CAGD_ERR_U_NOT_IN_SRF);
	    Crv = BspCrvNew(CrvLen = Srf -> VLength,
			    Srf -> VOrder, Srf -> PType);
	    GEN_COPY(Crv -> KnotVector, Srf -> VKnotVector,
			sizeof(CagdRType) * (Srf -> VLength + 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, t);
		    SrfP += CAGD_NEXT_V(Srf);
		}
	    }
	    break;
	case CAGD_CONST_V_DIR:
	    if (!BspKnotParamInDomain(Srf -> VKnotVector, Srf -> VLength,
							  Srf -> VOrder, t))
		FATAL_ERROR(CAGD_ERR_V_NOT_IN_SRF);
	    Crv = BspCrvNew(CrvLen = Srf -> ULength,
			    Srf -> UOrder, Srf -> PType);
	    GEN_COPY(Crv -> KnotVector, Srf -> UKnotVector,
			sizeof(CagdRType) * (Srf -> ULength + 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, t);
		    SrfP += CAGD_NEXT_U(Srf);
		}
	    }
	    break;
	default:
	    FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV);
	    break;
    }
    return Crv;
}

/******************************************************************************
* Extract an isoline curve out of the given mesh row/col.		      *
* The provided (zero based) Index specifies which row/col Index to extract.   *
******************************************************************************/
CagdCrvStruct *BspSrfCrvFromMesh(CagdSrfStruct *Srf, int Index,
							CagdSrfDirType Dir)
{
    CagdCrvStruct *Crv;
    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 + 1 > Srf -> ULength)
		FATAL_ERROR(CAGD_ERR_INDEX_NOT_IN_MESH);

	    Crv = BspCrvNew(CrvLen = Srf -> VLength,
			    Srf -> VOrder, Srf -> PType);
	    GEN_COPY(Crv -> KnotVector, Srf -> VKnotVector,
			sizeof(CagdRType) * (Srf -> VLength + 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 + 1 > Srf -> VLength)
		FATAL_ERROR(CAGD_ERR_INDEX_NOT_IN_MESH);

	    Crv = BspCrvNew(CrvLen = Srf -> ULength,
			    Srf -> UOrder, Srf -> PType);
	    GEN_COPY(Crv -> KnotVector, Srf -> UKnotVector,
			sizeof(CagdRType) * (Srf -> ULength + 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:
	    FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV);
	    break;
    }
    return Crv;
}