/****************************************************************************** * BspCoxDB.c - Bspline evaluation using Cox - de Boor recursive algorithm. * ******************************************************************************* * (C) Gershon Elber, Technion, Israel Institute of Technology * ******************************************************************************* * Written by Gershon Elber, Aug. 90. * ******************************************************************************/ #include #include #include #include "cagd_loc.h" #define COX_DB_IRIT_EPS 1e-20 #define COX_DB_APX_EQ(x, y) (ABS((x) - (y)) < COX_DB_IRIT_EPS) /***************************************************************************** * DESCRIPTION: M * Returns a pointer to a static data, holding the value of the curve at M * the prescribed parametric location t. M * Uses the recursive Cox de-Boor algorithm, to evaluate the spline, which M * is not very efficient if many evaluations of the same curve are necessary M * Use knot insertion when multiple evaluations are to be performed. M * * * PARAMETERS: M * Crv: To evaluate at the given parametric location t. M * t: The parameter value at which the curve Crv is to be evaluated. 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 * BspCrvEvalCoxDeBoor, evaluation, Bsplines M *****************************************************************************/ CagdRType *BspCrvEvalCoxDeBoor(CagdCrvStruct *Crv, CagdRType t) { static CagdRType Pt[CAGD_MAX_PT_COORD]; CagdBType IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv); CagdRType *p, *BasisFunc; int i, j, l, IndexFirst, k = Crv -> Order, Length = Crv -> Length, MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType); BasisFunc = BspCrvCoxDeBoorBasis(Crv -> KnotVector, k, CAGD_CRV_PT_LST_LEN(Crv), t, &IndexFirst); /* And finally multiply the basis functions with the control polygon. */ for (i = IsNotRational; i <= MaxCoord; i++) { Pt[i] = 0; p = Crv -> Points[i]; for (j = IndexFirst, l = 0; l < k; ) Pt[i] += p[j++ % Length] * BasisFunc[l++]; } return Pt; } /***************************************************************************** * DESCRIPTION: M * Returns a pointer to a vector of size Order, holding values of the non M * zero basis functions of a given curve at given parametric location t. M * This vector SHOULD NOT BE FREED. Although it is dynamically allocated, M * the returned pointer does not point to the beginning of this memory and it M * it be maintained by this routine (i.e. it might be freed next time this M * routine is being called). M * IndexFirst returns the index of first non zero basis function for the M * given parameter value t. M * Uses the recursive Cox de Boor algorithm, to evaluate the Bspline basis M * functions. M * Algorithm: M * Use the following recursion relation with B(i,0) == 1. M * M * t - t(i) t(i+k) - t V * B(i,k) = --------------- B(i,k-1) + --------------- B(i+1,k-1) V * t(i+k-1) - t(i) t(i+k) - t(i+1) V * M * Starting with constant Bspline (k == 0) only one basis function is non M * zero and is equal to one. This is the constant Bspline spanning interval M * t(i)...t(i+1) such that t(i) <= t < t(i+1). We then raise this constant M * Bspline to the prescribed Order and find in this process all the basis M * functions that are non zero in t for order Order. Sound simple hah!? M * * * PARAMETERS: M * KnotVector: To evaluate the Bspline Basis functions for this space. M * Order: Of the geometry. M * Len: Number of control points in the geometry. The length of M * KnotVector is equal to Len + Order. M * t: At which the Bspline basis functions are to be evaluated. M * IndexFirst: Index of the first Bspline basis function that might be M * non zero. M * * * RETURN VALUE: M * CagdRType *: A vector of length Order thats holds the values of the M * Bspline basis functions for the given t. A Bspline of M * order Order might have at most Order non zero basis M * functions that will hence start at IndexFirst and upto M * (*IndexFirst + Order - 1). M * * * KEYWORDS: M * BspCrvCoxDeBoorBasis, evaluation, Bsplines M *****************************************************************************/ CagdRType *BspCrvCoxDeBoorBasis(CagdRType *KnotVector, int Order, int Len, CagdRType t, int *IndexFirst) { static int BSize = 0; static CagdRType *B = NULL; CagdRType s1, s2, *BasisFunc; int i, l, Index, KVLen = Order + Len; if (!BspKnotParamInDomain(KnotVector, Len, Order, FALSE, t)) CAGD_FATAL_ERROR(CAGD_ERR_T_NOT_IN_CRV); if (t == KnotVector[Len]) t -= IRIT_UEPS; Index = BspKnotLastIndexLE(KnotVector, KVLen, t); /* Starting the recursion from constant splines - one spline is non */ /* zero and is equal to one. This is the spline that starts at Index. */ /* As we are going to reference index-1 we increment the buffer by one */ /* and save 0.0 at index-1. We then initialize the constant spline */ /* values - all are zero but the one from t(i) to t(i+1). */ if (BSize < 1 + Order) { if (B != NULL) IritFree((VoidPtr) B); B = (CagdRType *) IritMalloc(sizeof(CagdRType) * (1 + Order)); BSize = 1 + Order; } B[0] = 0.0; BasisFunc = &B[1]; if (Index >= KVLen - 1) { /* We are at the end of the parametric domain and this is open */ /* end condition - simply return last point. */ for (i = 0; i < Order; i++) BasisFunc[i] = (CagdRType) (i == Order - 1); *IndexFirst = Len - Order; return BasisFunc; } else for (i = 0; i < Order; i++) BasisFunc[i] = (CagdRType) (i == 0); /* Here is the tricky part. we raise these constant splines to the */ /* required order of the curve Crv for the given parameter t. There are */ /* at most order non zero function at param. value t. These functions */ /* start at Index-order+1 up to Index (order functions overwhole). */ for (i = 2; i <= Order; i++) { /* Goes through all orders... */ for (l = i - 1; l >= 0; l--) { /* And all basis funcs. of order i. */ s1 = (KnotVector[Index + l] - KnotVector[Index + l - i + 1]); s1 = COX_DB_APX_EQ(s1, 0.0) ? 0.0 : (t - KnotVector[Index + l - i + 1]) / s1; s2 = (KnotVector[Index + l + 1] - KnotVector[Index + l - i + 2]); s2 = COX_DB_APX_EQ(s2, 0.0) ? 0.0 : (KnotVector[Index + l + 1] - t) / s2; BasisFunc[l] = s1 * BasisFunc[l - 1] + s2 * BasisFunc[l]; } } *IndexFirst = Index - Order + 1; return BasisFunc; } /***************************************************************************** * DESCRIPTION: M * Computes the index of the first non zero basis function as returned by the M * BspCrvCoxDeBoorBasis function. M * * * PARAMETERS: M * KnotVector: To evaluate the Bspline Basis functions for this space. M * Order: Of the geometry. M * Len: Number of control points in the geometry. The length of M * KnotVector is equal to Len + Order. M * t: At which the Bspline basis functions are to be evaluated. M * * * RETURN VALUE: M * int: The index. M * * * KEYWORDS: M * BspCrvCoxDeBoorIndexFIrst, evaluation, Bsplines M *****************************************************************************/ int BspCrvCoxDeBoorIndexFirst(CagdRType *KnotVector, int Order, int Len, CagdRType t) { int Index, KVLen = Order + Len; Index = BspKnotLastIndexLE(KnotVector, KVLen, t); if (Index >= KVLen - 1) return Len - Order; else return Index - Order + 1; }