/****************************************************************************** * Bsp_knot.c - Various bspline routines to handle knot vectors. * ******************************************************************************* * (C) Gershon Elber, Technion, Israel Institute of Technology * ******************************************************************************* * Written by Gershon Elber, Aug. 90. * ******************************************************************************/ #include #include #include #include "cagd_loc.h" #define KNOT_IS_THE_SAME 1e-10 /***************************************************************************** * DESCRIPTION: M * Returns TRUE iff the given curve has no interior knot open end KV. M * * * PARAMETERS: M * Crv: To check for KV that mimics Bezier polynomial curve. M * * * RETURN VALUE: M * CagdBType: TRUE if same as Bezier curve, FALSE otherwise. M * * * KEYWORDS: M * BspCrvHasBezierKV, conversion M *****************************************************************************/ CagdBType BspCrvHasBezierKV(CagdCrvStruct *Crv) { return BspKnotHasBezierKV(Crv -> KnotVector, Crv-> Length, Crv -> Order); } /***************************************************************************** * DESCRIPTION: M * Returns TRUE iff the given surface has no interior knot open end KVs. M * * * PARAMETERS: M * Srf: To check for KVs that mimics Bezier polynomial surface. M * * * RETURN VALUE: M * CagdBType: TRUE if same as Bezier curve, FALSE otherwise. M * * * KEYWORDS: M * BspSrfHasBezierKVs, conversion M *****************************************************************************/ CagdBType BspSrfHasBezierKVs(CagdSrfStruct *Srf) { return BspKnotHasBezierKV(Srf -> UKnotVector, Srf-> ULength, Srf -> UOrder) && BspKnotHasBezierKV(Srf -> VKnotVector, Srf-> VLength, Srf -> VOrder); } /***************************************************************************** * DESCRIPTION: M * Returns TRUE iff the given knot vector has no interior knots and it has M * an open end cond. M * * * PARAMETERS: M * KnotVector: To check for open end and no interior knots conditions. M * Len: Of knot vector. M * Order: Of curve/surface the exploits this knot vector. M * * * RETURN VALUE: M * CagdBType: TRUE if has open end conditions and no interior knots, M * FALSE otherwise. M * * * KEYWORDS: M * BspKnotHasBezierKV, knot vectors, conversion M *****************************************************************************/ CagdBType BspKnotHasBezierKV(CagdRType *KnotVector, int Len, int Order) { return Len == Order * 2 && BspKnotHasOpenEC(KnotVector, Len, Order); } /***************************************************************************** * DESCRIPTION: M * Returns TRUE iff the given Bspline curve has open end coditions. M * * * PARAMETERS: M * Crv: To check for open end conditions. M * * * RETURN VALUE: M * CagdBType: TRUE, if curve has open end conditions, FALSE otherwise. M * * * KEYWORDS: M * BspCrvHasOpenEC, open end conditions M *****************************************************************************/ CagdBType BspCrvHasOpenEC(CagdCrvStruct *Crv) { return BspKnotHasOpenEC(Crv -> KnotVector, Crv -> Length, Crv -> Order); } /***************************************************************************** * DESCRIPTION: M * Returns TRUE iff the given Bspline surface has open end coditions. M * * * PARAMETERS: M * Srf: To check for open end conditions. M * * * RETURN VALUE: M * CagdBType: TRUE, if surface has open end conditions, FALSE otherwise. M * * * KEYWORDS: M * BspSrfHasOpenEC, open end conditions M *****************************************************************************/ CagdBType BspSrfHasOpenEC(CagdSrfStruct *Srf) { return BspKnotHasOpenEC(Srf -> UKnotVector, Srf -> ULength, Srf -> UOrder) && BspKnotHasOpenEC(Srf -> VKnotVector, Srf -> VLength, Srf -> VOrder); } /***************************************************************************** * DESCRIPTION: M * Returns TRUE iff the given Bspline surface has open end coditions in the M * specified direction. M * * * PARAMETERS: M * Srf: To check for open end conditions. M * Dir: Either the U or the V parametric direction. M * * * RETURN VALUE: M * CagdBType: TRUE, if surface has open end conditions, FALSE otherwise. M * * * KEYWORDS: M * BspSrfHasOpenECDir, open end conditions M *****************************************************************************/ CagdBType BspSrfHasOpenECDir(CagdSrfStruct *Srf, CagdSrfDirType Dir) { switch (Dir) { case CAGD_CONST_U_DIR: return BspKnotHasOpenEC(Srf -> UKnotVector, Srf -> ULength, Srf -> UOrder); case CAGD_CONST_V_DIR: return BspKnotHasOpenEC(Srf -> VKnotVector, Srf -> VLength, Srf -> VOrder); default: CAGD_FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV); return FALSE; } } /***************************************************************************** * DESCRIPTION: M * Returns TRUE iff the given knot vector has open end conditions. M * * * PARAMETERS: M * KnotVector: To check for open end condition. M * Len: Of knot vector. M * Order: Of curve/surface the exploits this knot vector. M * * * RETURN VALUE: M * CagdBType: TRUE if KV has open end conditions. M * * * KEYWORDS: M * BspKnotHasOpenEC, knot vectors, open end conditions M *****************************************************************************/ CagdBType BspKnotHasOpenEC(CagdRType *KnotVector, int Len, int Order) { int i, LastIndex = Len + Order - 1; for (i = 1; i < Order; i++) if (!APX_EQ_EPS(KnotVector[0], KnotVector[i], IRIT_UEPS)) return FALSE; for (i = LastIndex - 1; i >= Len; i--) if (!APX_EQ_EPS(KnotVector[LastIndex], KnotVector[i], IRIT_UEPS)) return FALSE; return TRUE; } /***************************************************************************** * DESCRIPTION: M * Returns TRUE iff t is in the parametric domain as define by the knot M * vector KnotVector, its length Len, and the order Order. M * * * PARAMETERS: M * KnotVector: To verify t is indeed in. M * Len: Length of curve/surface using KnotVector. This is NOT the M * length of KnotVector which is equal to (Length + Order). M * Order: Order of curve/surface using KnotVector. M * Periodic: TRUE if this KnotVector is periodic. M * t: Parametric value to verify. M * * * RETURN VALUE: M * CagdBType: TRUE, if t is contained in the parametric domain, FALSE M * otherwise. M * * * KEYWORDS: M * BspKnotParamInDomain, parametric domain, knot vectors M *****************************************************************************/ CagdBType BspKnotParamInDomain(CagdRType *KnotVector, int Len, int Order, CagdBType Periodic, CagdRType t) { CagdRType r1 = KnotVector[Order - 1], r2 = KnotVector[Len + (Periodic ? Order - 1 : 0)]; return (r1 < t || APX_EQ_EPS(r1, t, IRIT_UEPS)) && (r2 > t || APX_EQ_EPS(r2, t, IRIT_UEPS)); } /***************************************************************************** * DESCRIPTION: M * Returns the index of the last knot which is less than or equal to t in the M * given knot vector KnotVector of length Len. APX_EQ_EPS is used in equality.M * Parameter t is assumed to be in the parametric domain for the knot M * vector. M * * * PARAMETERS: M * KnotVector: To search for a knot with the LE relation to t. M * Len: Of knot vector. This is not the length of the M * curve/surface using this KnotVector. M * t: The parameter value to search for the LE relation. M * * * RETURN VALUE: M * int: Index of last knot in KnotVector that is LE t, or -1 if M * t is below the first knot. M * * * KEYWORDS: M * BspKnotLastIndexLE, knot vectors M *****************************************************************************/ int BspKnotLastIndexLE(CagdRType *KnotVector, int Len, CagdRType t) { int i; for (i = 0; i < Len && (*KnotVector <= t || APX_EQ_EPS(*KnotVector, t, IRIT_UEPS)); i++, KnotVector++); return i - 1; } /***************************************************************************** * DESCRIPTION: M * Returns the index of the last knot which is less t in the given knot M * vector KnotVector of length Len. APX_EQ_EPS is used for equality. M * Parameter t is assumed to be in the parametric domain for the knot M * vector. M * * * PARAMETERS: M * KnotVector: To search for a knot with the L relation to t. M * Len: Of knot vector. This is not the length of the M * curve/surface using this KnotVector. M * t: The parameter value to search for the L relation. M * * * RETURN VALUE: M * int: Index of last knot in KnotVector that is less than t or M * -1 if t is below the first knot. M * * * KEYWORDS: M * BspKnotLastIndexL, knot vectors M *****************************************************************************/ int BspKnotLastIndexL(CagdRType *KnotVector, int Len, CagdRType t) { int i; for (i = 0; i < Len && *KnotVector < t && !APX_EQ_EPS(*KnotVector, t, IRIT_UEPS); i++, KnotVector++); return i - 1; } /***************************************************************************** * DESCRIPTION: M * Returns the index of the first knot which is greater than t in the given M * knot vector KnotVector of length Len. APX_EQ_EPS is used for equality. M * Parameter t is assumed to be in the parametric domain for the knot M * vector. M * * * PARAMETERS: M * KnotVector: To search for a knot with the G relation to t. M * Len: Of knot vector. This is not the length of the M * curve/surface using this KnotVector. M * t: The parameter value to search for the G relation. M * * * RETURN VALUE: M * int: Index of first knot in KnotVector that is greater than t M * or Len if t is greater than last knot in KnotVector. M * * * KEYWORDS: M * BspKnotFirstIndexG, knot vectors M *****************************************************************************/ int BspKnotFirstIndexG(CagdRType *KnotVector, int Len, CagdRType t) { int i; CagdRType *KV; for (i = Len - 1, KV = &KnotVector[i]; i >= 0 && *KV > t && !APX_EQ_EPS(*KV, t, IRIT_UEPS); i--, KV--); return i + 1; } /***************************************************************************** * DESCRIPTION: M * Returns a uniform periodic knot vector for Len Control points and order M * Order. The actual length of the KV is Len + Order + Order - 1. M * If KnotVector is NULL it is being allocated dynamically. M * * * PARAMETERS: M * Len: Of control polygon/mesh of curve/surface that are to use M * this knot vector. M * Order: Of the curve/surface that are to use this knot vector. M * KnotVector: If new knot vector is to be saved here, otherwise a new M * space is allocated. M * * * RETURN VALUE: M * CagdRType *: The created uniform periodic knot vector, either newly M * allocated or given in Knotvector and just initialized. M * * * KEYWORDS: M * BspKnotUniformPeriodic, knot vectors, periodic end conditions, end M * conditions M *****************************************************************************/ CagdRType *BspKnotUniformPeriodic(int Len, int Order, CagdRType *KnotVector) { int i, KVLen = Len + Order + Order - 1; CagdRType *KV; if (KnotVector == NULL) KV = KnotVector = (CagdRType *) IritMalloc(sizeof(CagdRType) * KVLen); else KV = KnotVector; for (i = 0; i < KVLen; i++) *KV++ = (CagdRType) i; return KnotVector; } /***************************************************************************** * DESCRIPTION: M * Returns a uniform floating knot vector for Len Control points and order M * Order. The actual length of the KV is Len + Order. M * If KnotVector is NULL it is being allocated dynamically. M * * * PARAMETERS: M * Len: Of control polygon/mesh of curve/surface that are to use M * this knot vector. M * Order: Of the curve/surface that are to use this knot vector. M * KnotVector: If new knot vector is to be saved here, otherwise a new M * space is allocated. M * * * RETURN VALUE: M * CagdRType *: The created uniform floating knot vector, either newly M * allocated or given in Knotvector and just initialized. M * * * KEYWORDS: M * BspKnotUniformFloat, knot vectors, floating end conditions, end M * conditions M *****************************************************************************/ CagdRType *BspKnotUniformFloat(int Len, int Order, CagdRType *KnotVector) { int i; CagdRType *KV; if (KnotVector == NULL) KV = KnotVector = (CagdRType *) IritMalloc(sizeof(CagdRType) * (Len + Order)); else KV = KnotVector; for (i = 0; i < Len + Order; i++) *KV++ = (CagdRType) i; return KnotVector; } /***************************************************************************** * DESCRIPTION: M * Returns a uniform open knot vector for Len Control points and order M * Order. The actual length of the KV is Len + Order. M * If KnotVector is NULL it is being allocated dynamically. M * * * PARAMETERS: M * Len: Of control polygon/mesh of curve/surface that are to use M * this knot vector. M * Order: Of the curve/surface that are to use this knot vector. M * KnotVector: If new knot vector is to be saved here, otherwise a new M * space is allocated. M * * * RETURN VALUE: M * CagdRType *: The created uniform open knot vector, either newly M * allocated or given in Knotvector and just initialized. M * * * KEYWORDS: M * BspKnotUniformOpen, knot vectors, open end conditions, end M * conditions M *****************************************************************************/ CagdRType *BspKnotUniformOpen(int Len, int Order, CagdRType *KnotVector) { int i, j; CagdRType *KV; if (KnotVector == NULL) KV = KnotVector = (CagdRType *) IritMalloc(sizeof(CagdRType) * (Len + Order)); else KV = KnotVector; for (i = 0; i < Order; i++) *KV++ = 0.0; for (i = 1, j = Len - Order; i <= j; ) *KV++ = (CagdRType) i++; for (j = 0; j < Order; j++) *KV++ = (CagdRType) i; return KnotVector; } /***************************************************************************** * DESCRIPTION: M * Applies a scale transformation to the given knot vector. Note scale M * transformation on the knot vector does not change the Bspline curve. M * Knot vector is scaled by Scale as KV[i] = KV[i] * Scale. M * Scaling is taken place in place. M * * * PARAMETERS: M * KnotVector: To affinely transform. M * Len: Of knot vector. This is not the length of the curve or M * surface using this knot vector. M * Scale: Amount to scale the knot vector. M * * * RETURN VALUE: M * void M * * * SEE ALSO: M * BspKnotAffineTrans, BspKnotAffineTrans2 M * * * KEYWORDS: M * BspKnotScale, knot vectors, affine transformation M *****************************************************************************/ void BspKnotScale(CagdRType *KnotVector, int Len, CagdRType Scale) { int i; for (i = 0; i < Len; i++) KnotVector[i] = KnotVector[i] * Scale; } /***************************************************************************** * DESCRIPTION: M * Applies an affine transformation to the given knot vector. Note affine M * transformation on the knot vector does not change the Bspline curve. M * Knot vector is translated by Translate amount and scaled by Scale as M * KV[i] = (KV[i] - KV[0]) * Scale + (KV0 + Translate). V * All transformation as taken place in place. M * * * PARAMETERS: M * KnotVector: To affinely transform. M * Len: Of knot vector. This is not the length of the curve or M * surface using this knot vector. M * Translate: Amount to translate the knot vector. M * Scale: Amount to scale the knot vector. M * * * RETURN VALUE: M * void M * * * SEE ALSO: M * BspKnotScale, BspKnotAffineTrans2 M * * * KEYWORDS: M * BspKnotAffineTrans, knot vectors, affine transformation M *****************************************************************************/ void BspKnotAffineTrans(CagdRType *KnotVector, int Len, CagdRType Translate, CagdRType Scale) { int i; CagdRType KV0 = KnotVector[0]; for (i = 0; i < Len; i++) KnotVector[i] = (KnotVector[i] - KV0) * Scale + KV0 + Translate; } /***************************************************************************** * DESCRIPTION: M * Applies an affine transformation to the given knot vector. Note affine M * transformation on the knot vector does not change the Bspline curve. M * Knot vector is translated and scaled so as to span the domain from M * MinVal top MaxVal. This works for open end condition curves only. M * KV[i] = (KV[i] - KV[0]) * Scale + MinVal, V * where Scale = (MaxVal - MinVal) / (KV[Len - 1] - KV[0]). M * All transformation as taken place in place. M * * * PARAMETERS: M * KnotVector: To affinely transform. M * Len: Of knot vector. This is not the length of the curve or M * surface using this knot vector. M * MinVal, MaxVal: New parametric domain of knot vector. M * * * RETURN VALUE: M * void M * * * SEE ALSO: M * BspKnotScale, BspKnotAffineTrans M * * * KEYWORDS: M * BspKnotAffineTrans2, knot vectors, affine transformation M *****************************************************************************/ void BspKnotAffineTrans2(CagdRType *KnotVector, int Len, CagdRType MinVal, CagdRType MaxVal) { int i; CagdRType KV0 = KnotVector[0], KVn = KnotVector[Len - 1], Scale = (MaxVal - MinVal) / (KVn - KV0); for (i = 0; i < Len; i++) KnotVector[i] = (KnotVector[i] - KV0) * Scale + MinVal; } /***************************************************************************** * DESCRIPTION: M * Creates an identical copy of a given knot vector KnotVector of length Len. M * * * PARAMETERS: M * KnotVector: Knot vector to duplicate M * Len: Length of knot vector. This is not the length of the curve M * or surface using this knot vector. M * * * RETURN VALUE: M * CagdRType *: The duplicated knot vector. M * * * KEYWORDS: M * BspKnotCopy, allocation, knot vectors M *****************************************************************************/ CagdRType *BspKnotCopy(CagdRType *KnotVector, int Len) { CagdRType *NewKnotVector = (CagdRType *) IritMalloc(sizeof(CagdRType) * Len); CAGD_GEN_COPY(NewKnotVector, KnotVector, sizeof(CagdRType) * Len); return NewKnotVector; } /***************************************************************************** * DESCRIPTION: M * Reverse a knot vector of length Len. Reversing of knot vector keeps the M * knots monotonically non-decreasing as well as the parametric domain. Only M * the spaces between the knots are being flipped. For example the knot M * vector: M * [0 0 0 0 1 2 2 6 6 6 6] V * is reversed to be: M * [0 0 0 0 4 4 5 6 6 6 6] V * * * PARAMETERS: M * KnotVector: Knot vector to be reversed. M * Len: Length of knot vector. This is not the length of the curve M * or surface using this knot vector. M * * * RETURN VALUE: M * CagdRType *: The reversed knot vector. M * * * KEYWORDS: M * BspKnotReverse, knot vectors, reverse M *****************************************************************************/ CagdRType *BspKnotReverse(CagdRType *KnotVector, int Len) { int i; CagdRType *NewKnotVector = (CagdRType *) IritMalloc(sizeof(CagdRType) * Len), t = KnotVector[Len - 1] + KnotVector[0]; for (i = 0; i < Len; i++) NewKnotVector[i] = t - KnotVector[Len - i - 1]; return NewKnotVector; } /***************************************************************************** * DESCRIPTION: M * Returns a knot vector that contains all the knots in KnotVector1 that are M * not in KnotVector2. M * NewLen is set to new KnotVector length. M * * * PARAMETERS: M * KnotVector1: First knot vector. M * Len1: Length of KnotVector1. This is not the length of the M * curve or surface using this knot vector. M * KnotVector2: Second knot vector. M * Len1: Length of KnotVector2. This is not the length of the M * curve or surface using this knot vector. M * NewLen: To save the size of the knot vector that contains the M * computed subset of KnotVector1 / KnotVector2. M * * * RETURN VALUE: M * CagdRType *: The subset of knot in KnotVector1 that are not in M * KnotVector2 (KnotVector1 / KnotVector2). M * * * KEYWORDS: M * BspKnotSubtrTwo, knot vectors, compatibility, refinement M *****************************************************************************/ CagdRType *BspKnotSubtrTwo(CagdRType *KnotVector1, int Len1, CagdRType *KnotVector2, int Len2, int *NewLen) { int i = 0, j = 0; CagdRType *NewKnotVector = (CagdRType *) IritMalloc(sizeof(CagdRType) * Len1), *t = NewKnotVector; *NewLen = 0; while (i < Len1 && j < Len2) { if (APX_EQ_EPS(KnotVector1[i], KnotVector2[j], IRIT_UEPS)) { i++; j++; } else if (KnotVector1[i] > KnotVector2[j]) { j++; } else { /* Current KnotVector1 value is less than current KnotVector2: */ *t++ = KnotVector1[i++]; (*NewLen)++; } } return NewKnotVector; } /***************************************************************************** * DESCRIPTION: M * Merges two knot vector KnotVector1 and KnotVector2 of length Len1 and Len2 M * respectively into one. If Mult is not zero then knot multiplicity is M * tested not to be larger than Mult value. M * NewLen is set to new KnotVector length. M * * * PARAMETERS: M * KnotVector1: First knot vector. M * Len1: Length of KnotVector1. This is not the length of the M * curve or surface using this knot vector. M * KnotVector2: Second knot vector. M * Len2: Length of KnotVector2. This is not the length of the M * curve or surface using this knot vector. M * Mult: Maximum multiplicity to allow in merged knot vector. M * NewLen: To save the size of the knot vector that contains the M * merged knot vectors. M * * * RETURN VALUE: M * CagdRType *: The merged knot verctor (KnotVector1 U KnotVector2). M * * * KEYWORDS: M * BspKnotMergeTwo, knot vectors, compatibility, refinement M *****************************************************************************/ CagdRType *BspKnotMergeTwo(CagdRType *KnotVector1, int Len1, CagdRType *KnotVector2, int Len2, int Mult, int *NewLen) { int i = 0, j = 0, k = 0, m = 0, Len = Len1 + Len2; CagdRType t, *NewKnotVector = (CagdRType *) IritMalloc(sizeof(CagdRType) * Len); if (Mult == 0) Mult = Len + 1; while (i < Len1 && j < Len2) { if (KnotVector1[i] < KnotVector2[j]) { /* Use value from KnotVector1: */ t = KnotVector1[i++]; } else { /* Use value from KnotVector2: */ t = KnotVector2[j++]; } if (k == 0 || (k > 0 && !APX_EQ_EPS(NewKnotVector[k - 1], t, IRIT_UEPS))) { NewKnotVector[k++] = t; m = 1; } else if (m < Mult) { NewKnotVector[k++] = t; m++; } } while (i < Len1) NewKnotVector[k++] = KnotVector1[i++]; while (j < Len2) NewKnotVector[k++] = KnotVector1[j++]; /* It should be noted that k <= Len so there is a chance some of the new */ /* KnotVector space will not be used (it is not memory leak!). If you */ /* really care that much - copy it to a new allocated vector of size k */ /* and return it while freeing the original of size Len. */ *NewLen = k; return NewKnotVector; } /***************************************************************************** * DESCRIPTION: M * Merges two knot vector KnotVector1 and KnotVector2 of length Len1 and Len2 M * respectively into one, from geometries of orders Order1 and Order2. M * Merged knot vector is for order ResOrder so that the resulting curve can M * represent the discontinuities in both geometries. M * Assumes both knot vectors are open end spanning the same domain. M * * * PARAMETERS: M * KnotVector1: First knot vector. M * Len1: Length of KnotVector1. This is not the length of the M * curve or surface using this knot vector. M * Order1: Order of first knot vector's geometry. M * KnotVector2: Second knot vector. M * Len2: Length of KnotVector2. This is not the length of the M * curve or surface using this knot vector. M * Order2: Order of secon knot vector's geometry. M * ResOrder: Expected order of geometry that will use the merged M * knot vector. M * NewLen: To save the size of the knot vector that contains the M * merged knot vectors. M * * * RETURN VALUE: M * CagdRType *: The merged knot verctor (KnotVector1 U KnotVector2). M * * * KEYWORDS: M * BspKnotContinuityMergeTwo, knot vectors, compatibility, refinement M *****************************************************************************/ CagdRType *BspKnotContinuityMergeTwo(CagdRType *KnotVector1, int Len1, int Order1, CagdRType *KnotVector2, int Len2, int Order2, int ResOrder, int *NewLen) { int l1, l2, m, cont, i = 0, j = 0, k = 0, Len = (Len1 + Len2) * (1 + ResOrder - MAX(Order1, Order2)); CagdRType t, *NewKnotVector = (CagdRType *) IritMalloc(sizeof(CagdRType) * Len); while (i < Len1 && j < Len2) { if (APX_EQ_EPS(KnotVector1[i], KnotVector2[j], IRIT_UEPS)) { /* Compute multiplicity and continuity. */ for (l1 = 1; i + l1 < Len1 && APX_EQ_EPS(KnotVector1[i], KnotVector1[i + l1], IRIT_UEPS); l1++); for (l2 = 1; j + l2 < Len2 && APX_EQ_EPS(KnotVector2[j], KnotVector2[j + l2], IRIT_UEPS); l2++); cont = MIN(Order1 - l1, Order2 - l2) - 1; /* Use value from KnotVector1: */ t = KnotVector1[i]; i += l1; j += l2; } else if (KnotVector1[i] < KnotVector2[j]) { /* Compute multiplicity and continuity. */ for (l1 = 1; i + l1 < Len1 && APX_EQ_EPS(KnotVector1[i], KnotVector1[i + l1], IRIT_UEPS); l1++); cont = Order1 - l1 - 1; /* Use value from KnotVector1: */ t = KnotVector1[i]; i += l1; } else { /* Compute multiplicity and continuity. */ for (l2 = 1; j + l2 < Len2 && APX_EQ_EPS(KnotVector2[j], KnotVector2[j + l2], IRIT_UEPS); l2++); cont = Order2 - l2 - 1; /* Use value from KnotVector2: */ t = KnotVector2[j++]; j += l2; } for (m = 0; m < ResOrder - cont - 1; m++) NewKnotVector[k++] = t; } /* It should be noted that k <= Len so there is a chance some of the new */ /* KnotVector space will not be used (it is not memory leak!). If you */ /* really care that much - copy it to a new allocated vector of size k */ /* and return it while freeing the original of size Len. */ *NewLen = k; return NewKnotVector; } /***************************************************************************** * DESCRIPTION: M * Creates a new knot vector from the given KnotVector that averages Ave M * consequetive knots. M * Resulting vector will have (Len - Ave + 1) elements. M * * * PARAMETERS: M * KnotVector: To average out. M * Len: Length of KnotVector. This is not the length of the M * curve or surface using this knot vector. M * Ave: How many knots to average each time. M * * * RETURN VALUE: M * CagdRType *: The averaged knot vector of length (Len - Ave + 1). M * * * SEE ALSO: M * BspKnotNodes M * * * KEYWORDS: M * BspKnotAverage, knot vectors, node values M *****************************************************************************/ CagdRType *BspKnotAverage(CagdRType *KnotVector, int Len, int Ave) { int i, AveLen = Len - Ave + 1; CagdRType *AveVector = (CagdRType *) IritMalloc(sizeof(CagdRType) * AveLen); if (Ave > Len || Ave < 1) CAGD_FATAL_ERROR(CAGD_ERR_OUT_OF_RANGE); /* Compute the first average value. */ for (AveVector[0] = 0.0, i = 0; i < Ave; i++) AveVector[0] += KnotVector[i]; for (i = 1; i < AveLen; i++) AveVector[i] = AveVector[i - 1] + KnotVector[i + Ave - 1] - KnotVector[i - 1]; for (i = 0; i < AveLen; i++) AveVector[i] /= Ave; return AveVector; } /****************************************************************************** * DESCRIPTION: M * Creates a new vector with the given KnotVector Node values. The given M * knot vector is assumed to have open end conditions. M * The nodes are the approximated parametric value associated with the each M * control point. Therefore for a knot vector of length Len and order Order M * there are Len - Order control points and therefore nodes. M * Nodes are defined as (k = Order - 1 or degree): M * M * i+k V * - First Node N(i = 0) V * \ V * / KnotVector(j) Last Node N(i = Len - k - 2) V * - V * j=i+1 V * N(i) = ----------------- V * k V * M * PARAMETERS: M * KnotVector: To average out as nodes. M * Len: Length of KnotVector. This is not the length of the M * curve or surface using this knot vector. M * Order: Of curve or surface that exploits this knot vector. M * * * RETURN VALUE: M * CagdRType *: The nodes computed for the given knot vector. M * * * SEE ALSO: M * BspKnotAverage M * * * KEYWORDS: M * BspKnotNodes, knot vectors, node values M *****************************************************************************/ CagdRType *BspKnotNodes(CagdRType *KnotVector, int Len, int Order) { int i, j, k = MAX(Order - 1, 1), NodeLen = Len - Order; CagdRType TMin = KnotVector[k], TMax = KnotVector[NodeLen], *NodeVector = (CagdRType *) IritMalloc(sizeof(CagdRType) * NodeLen); for (i = 0; i < NodeLen; i++) { NodeVector[i] = 0.0; for (j = i + 1; j <= i + k; j++) NodeVector[i] += BOUND(KnotVector[j], TMin, TMax); NodeVector[i] /= k; /* This is due to round off errors. We cannot allow even a small */ /* fraction value below TMin or above TMax. */ NodeVector[i] = BOUND(NodeVector[i], TMin, TMax); } return NodeVector; } /***************************************************************************** * DESCRIPTION: M * Returns the nodes of a freeform curve. M * * * PARAMETERS: M * Crv: To compute node values for. M * * * RETURN VALUE: M * CagdRType *: Node values of the given curve. M * * * KEYWORDS: M * CagdCrvNodes, node values M *****************************************************************************/ CagdRType *CagdCrvNodes(CagdCrvStruct *Crv) { int i, Order = Crv -> Order, Length = CAGD_CRV_PT_LST_LEN(Crv); CagdRType *NodeVector; switch (Crv -> GType) { case CAGD_CBEZIER_TYPE: NodeVector = (CagdRType *) IritMalloc(sizeof(CagdRType) * Order); for (i = 0; i < Order; i++) NodeVector[i] = i / ((RealType) (Order - 1)); break; case CAGD_CBSPLINE_TYPE: NodeVector = BspKnotNodes(Crv -> KnotVector, Length + Order, Order); break; default: NodeVector = NULL; break; } return NodeVector; } /***************************************************************************** * DESCRIPTION: M * Returns the nodes of a freeform surface. M * * * PARAMETERS: M * Srf: To compute node values for. M * Dir: Either the U or the V parametric direction. M * * * RETURN VALUE: M * CagdRType *: Node values of the given surface and given parametric M * direction. M * * * KEYWORDS: M * CagdSrfNodes, node values M *****************************************************************************/ CagdRType *CagdSrfNodes(CagdSrfStruct *Srf, CagdSrfDirType Dir) { int i, Order = Dir == CAGD_CONST_U_DIR ? Srf -> UOrder : Srf -> VOrder, Length = Dir == CAGD_CONST_U_DIR ? Srf -> ULength : Srf -> VLength; CagdRType *NodeVector; if (Dir != CAGD_CONST_U_DIR && Dir != CAGD_CONST_V_DIR) { CAGD_FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV); return NULL; } switch (Srf -> GType) { case CAGD_SBEZIER_TYPE: NodeVector = (CagdRType *) IritMalloc(sizeof(CagdRType) * Order); for (i = 0; i < Order; i++) NodeVector[i] = i / ((RealType) (Order - 1)); break; case CAGD_SBSPLINE_TYPE: NodeVector = BspKnotNodes(Dir == CAGD_CONST_U_DIR ? Srf -> UKnotVector : Srf -> VKnotVector, Length + Order, Order); break; default: NodeVector = NULL; break; } return NodeVector; } /***************************************************************************** * DESCRIPTION: M * Finds the parameter value with the largest coefficient of the curve using M * nodes values to estimate the coefficients' parameters. M * * * PARAMETERS: M * Crv: To compute the parameter node value of the largest M * coefficient. M * Axis: Which axis should we search for maximal coefficient? M * 1 for X, 2 for Y, etc. M * MaxVal: The coefficient itself will be place herein. M * * * RETURN VALUE: M * CagdRType: The node parameter value of the detected maximal M * coefficient. M * * * KEYWORDS: M * BspCrvMaxCoefParam, extremum M *****************************************************************************/ CagdRType BspCrvMaxCoefParam(CagdCrvStruct *Crv, int Axis, CagdRType *MaxVal) { int i, Index = 0, Length = Crv -> Length, Order = Crv -> Order; CagdRType *Points = Crv -> Points[Axis], R = *Points, *Nodes = BspKnotNodes(Crv -> KnotVector, Length + Order, Order); for (i = 0; i < Length; i++, Points++) { if (*Points > R) { R = *Points; Index = i; } } *MaxVal = R; R = Nodes[Index]; IritFree(Nodes); return R; } /***************************************************************************** * DESCRIPTION: M * Finds the parameter value with the largest coefficient of the surface M * using nodes values to estimate the coefficients' parameters. M * Returns a pointer to a static array of two elements holding U and V. M * * * PARAMETERS: M * Srf: To compute the parameter node value of the largest M * coefficient. M * Axis: Which axis should we search for maximal coefficient? M * 1 for X, 2 for Y, etc. M * MaxVal: The coefficient itself will be place herein. M * * * RETURN VALUE: M * CagdRType *: The node UV parameter values of the detected maximal M * coefficient. M * * * KEYWORDS: M * BspSrfMaxCoefParam, extremum M *****************************************************************************/ CagdRType *BspSrfMaxCoefParam(CagdSrfStruct *Srf, int Axis, CagdRType *MaxVal) { static CagdRType UV[2]; int i, UIndex = 0, VIndex = 0, ULength = Srf -> ULength, VLength = Srf -> VLength, UOrder = Srf -> UOrder, VOrder = Srf -> VOrder; CagdRType *Points = Srf -> Points[Axis], R = *Points, *UNodes = BspKnotNodes(Srf -> UKnotVector, ULength + UOrder, UOrder), *VNodes = BspKnotNodes(Srf -> VKnotVector, VLength + VOrder, VOrder); for (i = 0; i < ULength * VLength; i++, Points++) { if (*Points > R) { R = *Points; UIndex = i % ULength; VIndex = i / ULength; } } *MaxVal = R; UV[0] = UNodes[UIndex]; UV[1] = VNodes[VIndex]; IritFree(UNodes); IritFree(VNodes); return UV; } /***************************************************************************** * DESCRIPTION: M * Creates a new vector with t inserted as a new knot. Attempt is made to M * make sure t is in the knot vector domain. M * No test is made for the current multiplicity of knot t in KnotVector. M * This function only constructs a refined knot vector and does not M * compute the actual refined coefficients. M * * * PARAMETERS: M * KnotVector: To insert new knot t in. M * Order: Of geometry that exploits KnotVector. M * Len: Length of curve/surface using KnotVector. This is NOT the M * length of KnotVector which is equal to (Length + Order). M * t: The new knot t to insert. M * * * RETURN VALUE: M * CagdRType *: A new knot vector larger by one than KnotVector that M * contains t. M * * * KEYWORDS: M * BspKnotInsertOne, knot vectors, knot insertion, refinement M *****************************************************************************/ CagdRType *BspKnotInsertOne(CagdRType *KnotVector, int Order, int Len, CagdRType t) { int Mult = BspKnotFindMult(KnotVector, Order, Len, t) + 1; return BspKnotInsertMult(KnotVector, Order, &Len, t, Mult); } /***************************************************************************** * DESCRIPTION: M * Inserts Mult knots with value t into the knot vector KnotVector. M * Attempt is made to make sure t in knot vector domain. M * If a knot equal to t (up to APX_EQ) already exists with multiplicity i M * only (Mult - i) knot are being inserted into the new knot vector. M * Len is updated to the resulting knot vector. M * It is possible to DELETE a knot using this routine by specifying M * multiplicity less then current multiplicity! M * This function only constructs a refined knot vector and does not M * compute the actual refined coefficients. M * * * PARAMETERS: M * KnotVector: To insert new knot t in. M * Order: Of geometry that exploits KnotVector. M * Len: Length of curve/surface using KnotVector. This is NOT the M * length of KnotVector which is equal to (Length + Order). M * t: The new knot t to insert. M * Mult: The multiplicity that this knot should have in resulting M * knot vector. M * * * RETURN VALUE: M * CagdRType *: A new knot vector derived from KnotVector that has M * a mltiplicity of exacly Mult at the knot t. M * * * KEYWORDS: M * BspKnotInsertMult, knot vectors, knot insertion, refinement M *****************************************************************************/ CagdRType *BspKnotInsertMult(CagdRType *KnotVector, int Order, int *Len, CagdRType t, int Mult) { int i, CurrentMult, NewLen, FirstIndex; CagdRType *NewKnotVector; if (!BspKnotParamInDomain(KnotVector, *Len, Order, FALSE, t)) CAGD_FATAL_ERROR(CAGD_ERR_T_NOT_IN_CRV); CurrentMult = BspKnotFindMult(KnotVector, Order, *Len, t); NewLen = *Len + Mult - CurrentMult; NewKnotVector = (CagdRType *) IritMalloc(sizeof(CagdRType) * (NewLen + Order)); FirstIndex = BspKnotLastIndexL(KnotVector, *Len + Order, t) + 1; /* Copy all the knot before the knot t. */ CAGD_GEN_COPY(NewKnotVector, KnotVector, sizeof(CagdRType) * FirstIndex); /* Insert Mult knot of value t. */ for (i = FirstIndex; i < FirstIndex + Mult; i++) NewKnotVector[i] = t; /* And copy the second part. */ CAGD_GEN_COPY(&NewKnotVector[FirstIndex + Mult], &KnotVector[FirstIndex + CurrentMult], sizeof(CagdRType) * (*Len + Order - FirstIndex - Mult + 1)); *Len = NewLen; return NewKnotVector; } /***************************************************************************** * DESCRIPTION: M * Returns the multiplicity of knot t in knot vector KnotVector, zero if M * none. M * * * PARAMETERS: M * KnotVector: To test multiplicity of knot value t at. M * Order: Of geometry that exploits KnotVector. M * Len: Length of curve/surface using KnotVector. This is NOT the M * length of KnotVector which is equal to (Length + Order). M * t: The knot to verify the multiplicity of. M * * * RETURN VALUE: M * int: Multiplicity of t in KnotVector. M * * * KEYWORDS: M * BspKnotFindMult, knot vectors M *****************************************************************************/ int BspKnotFindMult(CagdRType *KnotVector, int Order, int Len, CagdRType t) { int LastIndex, FirstIndex; if (!BspKnotParamInDomain(KnotVector, Len, Order, FALSE, t)) CAGD_FATAL_ERROR(CAGD_ERR_T_NOT_IN_CRV); FirstIndex = BspKnotLastIndexL(KnotVector, Len, t) + 1; for (LastIndex = FirstIndex; LastIndex < Len && APX_EQ_EPS(KnotVector[LastIndex], t, IRIT_UEPS); LastIndex++); return LastIndex - FirstIndex; } /***************************************************************************** * DESCRIPTION: M * Scans the given knot vector to a potential C1 discontinuity. M * Returns TRUE if found one and set t to its parameter value. M * Assumes knot vector has open end condition. M * A knot vector with multiplicity of a knot of (Order - 1) can be C1 M * discontinuous at that knot. However, this is only a necessary condition M * for C1 discontinuity in the geometry. M * * * PARAMETERS: M * KnotVector: To test for potential C1 dicontinuities. M * Order: Of geometry that exploits KnotVector. M * Length: Length of curve/surface using KnotVector. This is NOT the M * length of KnotVector which is equal to (Length + Order). M * t: Where to put the parameter value (knot) that can be C2 M * discontinuous. M * * * RETURN VALUE: M * CagdBType: TRUE if found a potential C1 discontinuity, FALSE M * otherwise. M * * * KEYWORDS: M * BspKnotC1Discont, knot vectors, continuity, discontinuity M *****************************************************************************/ CagdBType BspKnotC1Discont(CagdRType *KnotVector, int Order, int Length, CagdRType *t) { int i, Count; CagdRType LastT = KnotVector[0] - 1.0; for (i = Order, Count = 0; i < Length; i++) { if (APX_EQ_EPS(LastT, KnotVector[i], IRIT_UEPS)) Count++; else { Count = 1; LastT = KnotVector[i]; } if (Count >= Order - 1) { *t = LastT; return TRUE; } } return FALSE; } /***************************************************************************** * DESCRIPTION: M * Scans the given knot vector for all potential C1 discontinuity. Returns M * a vector holding the parameter values of the potential C1 discontinuities, M * NULL of none found. M * Sets n to length of returned vector. M * Assumes knot vector has open end condition. M * A knot vector with multiplicity of a knot of (Order - 1) can be C1 M * discontinuous at that knot. However, this is only a necessary condition M * for C1 discontinuity in the geometry. M * * * PARAMETERS: M * KnotVector: To test for potential C1 dicontinuities. M * Order: Of geometry that exploits KnotVector. M * Length: Length of curve/surface using KnotVector. This is NOT the M * length of KnotVector which is equal to (Length + Order). M * n: Length of returned vector - number of potential C1 M * discontinuities found. M * * * RETURN VALUE: M * CagdRType *: Vector holding all parametr values with potential C1 M * discontinuities. M * * * KEYWORDS: M * BspKnotAllC1Discont, knot vectors, continuity, discontinuity M *****************************************************************************/ CagdRType *BspKnotAllC1Discont(CagdRType *KnotVector, int Order, int Length, int *n) { int i, Count, C1DiscontCount = 0; CagdRType LastT = KnotVector[0] - 1.0, *C1Disconts = (CagdRType *) IritMalloc(sizeof(CagdRType) * Length); for (i = Order, Count = 0; i < Length; i++) { if (APX_EQ_EPS(LastT, KnotVector[i], IRIT_UEPS)) Count++; else { Count = 1; LastT = KnotVector[i]; } if (Count >= Order - 1) { C1Disconts[C1DiscontCount++] = LastT; Count = 0; } } if ((*n = C1DiscontCount) > 0) { return C1Disconts; } else { IritFree((VoidPtr) C1Disconts); return NULL; } } /***************************************************************************** * DESCRIPTION: M * Routine to determine where to sample along the provided paramtric domain, M * given the C1 discontinuities along it. M * Returns a vector of length NumSamples. M * If C1Disconts != NULL (NumC1Disconts > 0), C1Discont is being freed. M * * * PARAMETERS: M * PMin: Minimum of parametric domain. M * PMax: Maximum of parametric domain. M * NumSamples: To allocate for the vector of samples. M * C1Disconts: A vector of potential C1 discontinuities in the M * (PMin, PMax) domain. This vector is freed by this M * routine, if it is not NULL. M * NumC1Disconts: Length of C1Discont. if zero then C1Discont == NULL. M * * * RETURN VALUE: M * CagdRType *: A vector of the suggested set of sampling locations. M * * * KEYWORDS: M * BspKnotParamValues, piecewise linear approximation, knot vectors M *****************************************************************************/ CagdRType *BspKnotParamValues(CagdRType PMin, CagdRType PMax, int NumSamples, CagdRType *C1Disconts, int NumC1Disconts) { int i, CrntIndex, NextIndex, CrntDiscont; CagdRType Step, *Samples = (CagdRType *) IritMalloc(sizeof(CagdRType) * NumSamples); Samples[0] = PMin; if (NumSamples <= 1) return Samples; Samples[NumSamples - 1] = PMax; if (NumSamples <= 2) return Samples; if (NumC1Disconts > NumSamples - 2) { /* More C1 discontinuities than we are sampling. Grab a subset of */ /* The discontinuities as the sampling set. */ Step = ((CagdRType) (NumC1Disconts - 1)) / (NumSamples - 2) - IRIT_UEPS; for (i = 0; i < NumSamples - 2; i++) Samples[i + 1] = C1Disconts[(int) (i * Step)]; } else { /* More samples than C1 discontinuites. Uniformly distribute the C1 */ /* discontinuites between the samples and linearly interpolate the */ /* sample values in between. */ Step = ((CagdRType) (NumC1Disconts + 1)) / (NumSamples - 2); CrntIndex = CrntDiscont = 0; while (CrntDiscont < NumC1Disconts) { NextIndex = (int) ((CrntDiscont + 1) / Step); Samples[NextIndex] = C1Disconts[CrntDiscont++]; for (i = CrntIndex + 1; i < NextIndex; i++) { CagdRType t = (i - CrntIndex) / ((CagdRType) (NextIndex - CrntIndex)), t1 = 1.0 - t; Samples[i] = Samples[CrntIndex] * t1 + Samples[NextIndex] * t; } CrntIndex = NextIndex; } /* Interpolate the last interval: */ for (i = CrntIndex + 1; i < NumSamples - 1; i++) { CagdRType t = (i - CrntIndex) / ((CagdRType) (NumSamples - 1 - CrntIndex)), t1 = 1.0 - t; Samples[i] = Samples[CrntIndex] * t1 + Samples[NumSamples - 1] * t; } } if (C1Disconts != NULL) IritFree((VoidPtr) C1Disconts); return Samples; } /***************************************************************************** * DESCRIPTION: M * Given a knot vector, make sure adjacent knots that are close "enough" are M * actually identical. Important for robustness of subdiv/refinement algs. M * M * M * * * PARAMETERS: M * KV: Knot vector to make robust, in place. M * Len: Length of knot vector KV. M * * * RETURN VALUE: M * void M * * * KEYWORDS: M * BspKnotMakeRobustKVm knot vectors M *****************************************************************************/ void BspKnotMakeRobustKV(CagdRType *KV, int Len) { int i; for (i = 1; i < Len; i++) if (KV[i] - KV[i - 1] < KNOT_IS_THE_SAME && KV[i] != KV[i - 1]) KV[i] = KV[i - 1]; }