/****************************************************************************** * SBspEval.c - Bspline surfaces handling routines - evaluation routines. * ******************************************************************************* * (C) Gershon Elber, Technion, Israel Institute of Technology * ******************************************************************************* * Written by Gershon Elber, Mar. 90. * ******************************************************************************/ #include #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 * * * SEE ALSO: M * CagdSrfEval, BzrSrfEvalAtParam, BspSrfEvalAtParam2, TrimSrfEval M * * * KEYWORDS: M * BspSrfEvalAtParam, evaluation, Bsplines M *****************************************************************************/ CagdRType *BspSrfEvalAtParam(CagdSrfStruct *Srf, CagdRType u, CagdRType v) { static CagdCrvStruct *IsoSubCrv = NULL; CagdSrfEvalCashStruct *SrfEvalCash; 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 - IRIT_EPS || u > UMax + IRIT_EPS) CAGD_FATAL_ERROR(CAGD_ERR_U_NOT_IN_SRF); if (v < VMin - IRIT_EPS || v > VMax + IRIT_EPS) CAGD_FATAL_ERROR(CAGD_ERR_V_NOT_IN_SRF); if (u > UMax - IRIT_UEPS * 2) u = UMax - IRIT_UEPS * 2; if (v > VMax - IRIT_UEPS * 2) v = VMax - IRIT_UEPS * 2; if (u < UMin) u = UMin; if (v < VMin) v = VMin; /* Create the cashed data structure, if surface has none. */ if (Srf -> PAux == NULL) { SrfEvalCash = (CagdSrfEvalCashStruct *) IritMalloc(sizeof(CagdSrfEvalCashStruct)); Srf -> PAux = (VoidPtr) SrfEvalCash; UIndexFirst = BspKnotLastIndexLE(Srf -> UKnotVector, (CAGD_SRF_UPT_LST_LEN(Srf) + UOrder), u) - (UOrder - 1); VIndexFirst = BspCrvCoxDeBoorIndexFirst(Srf -> VKnotVector, VOrder, CAGD_SRF_VPT_LST_LEN(Srf), v); IsoSubCrv = SrfEvalCash -> IsoSubCrv = BspCrvNew(UOrder, UOrder, Srf -> PType); CAGD_GEN_COPY(IsoSubCrv -> KnotVector, &Srf -> UKnotVector[UIndexFirst], sizeof(CagdRType) * UOrder * 2); } else { SrfEvalCash = (CagdSrfEvalCashStruct *) Srf -> PAux; UIndexFirst = BspKnotLastIndexLE(Srf -> UKnotVector, (CAGD_SRF_UPT_LST_LEN(Srf) + UOrder), u) - (UOrder - 1); VIndexFirst = BspCrvCoxDeBoorIndexFirst(Srf -> VKnotVector, VOrder, CAGD_SRF_VPT_LST_LEN(Srf), v); IsoSubCrv = SrfEvalCash -> IsoSubCrv; } /* Make sure the cashed curve is in the proper location. */ CAGD_GEN_COPY(IsoSubCrv -> KnotVector, &Srf -> UKnotVector[UIndexFirst], sizeof(CagdRType) * UOrder * 2); /* Make sure the cashed V basis functions are at the proper location. */ VBasisFunc = BspCrvCoxDeBoorBasis(Srf -> VKnotVector, VOrder, CAGD_SRF_VPT_LST_LEN(Srf), v, &VIndexFirst); 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 * This function is the same as BspSrfEvalAtParam above. Cleaner, but much M * less efficient. 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 * * * 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 * * * SEE ALSO: M * CagdSrfEval, BzrSrfEvalAtParam, BspSrfEvalAtParam, TrimSrfEval M * * * KEYWORDS: M * BspSrfEvalAtParam2, evaluation, Bsplines 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 * * * SEE ALSO: M * CagdCrvFromSrf, BzrSrfCrvFromSrf, CagdCrvFromMesh, BzrSrfCrvFromMesh, M * BspSrfCrvFromMesh 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 * * * SEE ALSO: M * CagdCrvFromSrf, BzrSrfCrvFromSrf, BspSrfCrvFromSrf, M * CagdCrvFromMesh, BzrSrfCrvFromMesh 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; }