/****************************************************************************** * CBsp-Aux.c - Bspline curve auxilary routines. * ******************************************************************************* * (C) Gershon Elber, Technion, Israel Institute of Technology * ******************************************************************************* * Written by Gershon Elber, Aug. 90. * ******************************************************************************/ #include #include #include #include "cagd_loc.h" #define EPS_ROUND_KNOT 1e-12 #define MOEBIUS_MEU(t) (1 + (t) * (1 - c) / c) #define MOEBIUS_REPARAM(t) ((t) / ((t) + c * (1 - (t)))) #define CAGD_SAME_PT_EPS 1e-6 STATIC_DATA CagdBType GlblDeriveScalar = FALSE; /***************************************************************************** * DESCRIPTION: M * Apply B-spline subdivision to the given curve Crv at parameter value t, M * and save the result in curves LPoints/RPoints. M * * * PARAMETERS: M * Crv: To subdivide at parametr value t. M * LPoints, RPoints: Where the results are kept. M * t: Parameter value to subdivide Crv at. M * * * RETURN VALUE: M * void M * * * SEE ALSO: M * BzrCrvSubdivCtlPoly, BspCrvSubdivAtParam M * * * KEYWORDS: M * BspCrvSubdivCtlPoly M *****************************************************************************/ void BspCrvSubdivCtlPoly(CagdCrvStruct *Crv, CagdRType **LPoints, CagdRType **RPoints, int LLength, int RLength, CagdRType t, BspKnotAlphaCoeffStruct *A, int Mult) { CagdBType IsNotRational; int i, j, Len, MaxCoord, k = Crv -> Order; CagdRType TMin, TMax, **Points; CagdCrvDomain(Crv, &TMin, &TMax); /* If it is a Bezier knot sequence, use the simpler Bezier subdivision. */ if (BspCrvHasBezierKV(Crv)) { BzrCrvSubdivCtlPoly(Crv -> Points, LPoints, RPoints, k, Crv -> PType, (t - TMin) / (TMax - TMin)); return; } Len = Crv -> Length; IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv); MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType); Points = Crv -> Points; if (Mult > 0) { CagdRType *NewKV = IritMalloc(sizeof(CagdRType) * Mult); for (i = 0; i < Mult; i++) NewKV[i] = (t == TMax ? t - CAGD_DOMAIN_IRIT_EPS : t); A = BspKnotEvalAlphaCoefMerge(k, Crv -> KnotVector, Len, NewKV, Mult, FALSE); IritFree(NewKV); } else { A = BspKnotEvalAlphaCoef(k, Crv -> KnotVector, Len, Crv -> KnotVector, Len, FALSE); } /* Note that Mult can be negative in cases where original */ /* multiplicity was order or more and we need to compensate */ /* here, since Alpha matrix will be just a unit matrix then. */ Mult = Mult >= 0 ? 0 : -Mult; /* Blend Crv into LCrv. */ for (j = IsNotRational; j <= MaxCoord; j++) { BspKnotAlphaLoopBlendNotPeriodic(A, 0, LLength, Points[j], LPoints[j]); } /* Blend Crv into RCrv. */ for (j = IsNotRational; j <= MaxCoord; j++) { BspKnotAlphaLoopBlendNotPeriodic(A, LLength - 1 + Mult, LLength + RLength - 1 + Mult, Points[j], RPoints[j]); } BspKnotFreeAlphaCoef(A); } /***************************************************************************** * DESCRIPTION: M * Given a Bspline curve - subdivides it into two sub-curves at the given M * parametric value. M * Returns pointer to first curve in a list of two subdivided curves. M * The subdivision is achieved by inserting (order-1) knot at the given M * parameter value t and spliting the control polygon and knot vector at that M * location. M * * * PARAMETERS: M * Crv: To subdivide at parametr value t. M * t: Parameter value to subdivide Crv at. M * * * RETURN VALUE: M * CagdCrvStruct *: A list of the two subdivided curves. M * * * KEYWORDS: M * BspCrvSubdivAtParam, subdivision, refinement M *****************************************************************************/ CagdCrvStruct *BspCrvSubdivAtParam(CagdCrvStruct *Crv, CagdRType t) { CagdBType NewCrv = FALSE; int i, j, Len, KVLen, Index1, Index2, Mult, k = Crv -> Order; CagdRType TMin, TMax; CagdCrvStruct *LCrv, *RCrv; if (CAGD_IS_PERIODIC_CRV(Crv)) { NewCrv = TRUE; Crv = CnvrtPeriodic2FloatCrv(Crv); } Len = Crv -> Length; KVLen = k + Len; i = BspKnotLastIndexLE(Crv -> KnotVector, KVLen, t); if (APX_EQ_EPS(t, Crv -> KnotVector[i], EPS_ROUND_KNOT)) t = Crv -> KnotVector[i]; else if (i < KVLen && APX_EQ_EPS(t, Crv -> KnotVector[i + 1], EPS_ROUND_KNOT)) t = Crv -> KnotVector[i + 1]; Index1 = BspKnotLastIndexL(Crv -> KnotVector, KVLen, t); if (Index1 + 1 < k) Index1 = k - 1; Index2 = BspKnotFirstIndexG(Crv -> KnotVector, KVLen, t); if (Index2 > Len) Index2 = Len; Mult = k - 1 - (Index2 - Index1 - 1); CAGD_DOMAIN_GET_AND_VERIFY_CRV(t, Crv, TMin, TMax); LCrv = BspCrvNew(Index1 + 1, k, Crv -> PType); RCrv = BspCrvNew(Len - Index2 + k, k, Crv -> PType); /* Update the new knot vectors. */ CAGD_GEN_COPY(LCrv -> KnotVector, Crv -> KnotVector, sizeof(CagdRType) * (Index1 + 1)); /* Close the knot vector with multiplicity Order: */ for (j = Index1 + 1; j <= Index1 + k; j++) LCrv -> KnotVector[j] = t; CAGD_GEN_COPY(&RCrv -> KnotVector[k], &Crv -> KnotVector[Index2], sizeof(CagdRType) * (Len + k - Index2)); /* Make sure knot vector starts with multiplicity Order: */ for (j = 0; j < k; j++) RCrv -> KnotVector[j] = t; BspKnotMakeRobustKV(RCrv -> KnotVector, RCrv -> Order + RCrv -> Length); BspKnotMakeRobustKV(LCrv -> KnotVector, LCrv -> Order + LCrv -> Length); /* Now handle the control polygon refinement. */ BspCrvSubdivCtlPoly(Crv, LCrv -> Points, RCrv -> Points, LCrv -> Length, RCrv -> Length, t, NULL, Mult); LCrv -> Pnext = RCrv; CAGD_PROPAGATE_ATTR(LCrv, Crv); CAGD_PROPAGATE_ATTR(RCrv, Crv); if (NewCrv) CagdCrvFree(Crv); return LCrv; } /***************************************************************************** * DESCRIPTION: M * Inserts n knot, all with the value t. In no case will the multiplicity of M * a knot be greater or equal to the curve order. M * * * PARAMETERS: M * Crv: To refine by insertion (upto) n knot of value t. M * t: Parameter value of new knot to insert. M * n: Maximum number of times t should be inserted. M * * * RETURN VALUE: M * CagdCrvStruct *: Refined Crv with n knots of value t. M * * * KEYWORDS: M * BspCrvKnotInsertNSame, refinement, subdivision M *****************************************************************************/ CagdCrvStruct *BspCrvKnotInsertNSame(CagdCrvStruct *Crv, CagdRType t, int n) { int i, CrntMult = BspKnotFindMult(Crv -> KnotVector, Crv -> Order, Crv -> Length, t), Mult = MIN(n, Crv -> Order - CrntMult - 1); CagdCrvStruct *RefinedCrv; if (Mult > 0) { CagdRType *NewKV = (CagdRType *) IritMalloc(sizeof(CagdRType) * Mult); for (i = 0; i < Mult; i++) NewKV[i] = t; RefinedCrv = BspCrvKnotInsertNDiff(Crv, FALSE, NewKV, Mult); IritFree(NewKV); } else { RefinedCrv = CagdCrvCopy(Crv); } return RefinedCrv; } /***************************************************************************** * DESCRIPTION: M * Inserts n knot with different values as defined by the vector t. If, M * however, Replace is TRUE, the knot are simply replacing the current knot M * vector. M * * * PARAMETERS: M * Crv: To refine by insertion (upto) n knot of value t. M * Replace: if TRUE, the n knots in t should replace the knot vector M * of size n of Crv. Sizes must match. If False, n new knots M * as defined by t will be introduced into Crv. M * t: New knots to introduce/replace knot vector of Crv. M * n: Size of t. M * * * RETURN VALUE: M * CagdCrvStruct *: Refined Crv with n new knots. M * * * KEYWORDS: M * BspCrvKnotInsertNDiff, refinement, subdivision M *****************************************************************************/ CagdCrvStruct *BspCrvKnotInsertNDiff(CagdCrvStruct *Crv, CagdBType Replace, CagdRType *t, int n) { CagdBType IsPeriodic = CAGD_IS_PERIODIC_CRV(Crv), IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv); CagdRType *KnotVector = Crv -> KnotVector; int Len = Crv -> Length, Length = CAGD_CRV_PT_LST_LEN(Crv), Order = Crv -> Order, MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType); CagdCrvStruct *RefinedCrv; if (Replace) { int i; if (Order + Length != n) CAGD_FATAL_ERROR(CAGD_ERR_NUM_KNOT_MISMATCH); for (i = 1; i < n; i++) if (t[i] < t[i - 1]) CAGD_FATAL_ERROR(CAGD_ERR_KNOT_NOT_ORDERED); RefinedCrv = CagdCrvCopy(Crv); for (i = 0; i < n; i++) RefinedCrv -> KnotVector[i] = *t++; } else if (n <= 0) { RefinedCrv = CagdCrvCopy(Crv); } else { int i, j, LengthKVt, NewLen; BspKnotAlphaCoeffStruct *A; CagdRType *MergedKVt, TMin, TMax; BspCrvDomain(Crv, &TMin, &TMax); for (i = 1; i < n; i++) if (t[i] < t[i - 1]) CAGD_FATAL_ERROR(CAGD_ERR_KNOT_NOT_ORDERED); for (i = 0; i < n; i++) if (t[i] >= TMax) t[i] = TMax - CAGD_DOMAIN_IRIT_EPS; /* Compute the Alpha refinement matrix. */ MergedKVt = BspKnotMergeTwo(KnotVector, Length + Order, t, n, 0, &LengthKVt); A = BspKnotEvalAlphaCoef(Order, KnotVector, Length, MergedKVt, LengthKVt - Order, IsPeriodic); NewLen = Len + n; RefinedCrv = BspPeriodicCrvNew(NewLen, Order, IsPeriodic, Crv -> PType); /* Update the knot vector. */ IritFree(RefinedCrv -> KnotVector); RefinedCrv -> KnotVector = MergedKVt; if (IsPeriodic) { /* Make sure the knot vector is indeed periodic. */ BspKnotVerifyPeriodicKV(RefinedCrv -> KnotVector, Order, CAGD_CRV_PT_LST_LEN(RefinedCrv)); } /* Update the control mesh */ for (j = IsNotRational; j <= MaxCoord; j++) { CagdRType *ROnePts = RefinedCrv -> Points[j], *OnePts = Crv -> Points[j]; if (CAGD_IS_PERIODIC_CRV(Crv)) { BspKnotAlphaLoopBlendPeriodic(A, 0, NewLen, OnePts, Len, ROnePts); } else { BspKnotAlphaLoopBlendNotPeriodic(A, 0, NewLen, OnePts, ROnePts); } } BspKnotFreeAlphaCoef(A); } BspKnotMakeRobustKV(RefinedCrv -> KnotVector, RefinedCrv -> Order + RefinedCrv -> Length); CAGD_PROPAGATE_ATTR(RefinedCrv, Crv); return RefinedCrv; } /***************************************************************************** * DESCRIPTION: M * Returns a new curve, identical to the original but with order N. M * Degree raise is computed by multiplying by a constant 1 curve of order M * * * PARAMETERS: M * Crv: To raise its degree to a NewOrder. M * NewOrder: NewOrder for Crv. M * * * RETURN VALUE: M * CagdCrvStruct *: A curve of order NewOrder representing the same M * geometry as Crv. M * * * KEYWORDS: M * BspCrvDegreeRaiseN, degree raising M *****************************************************************************/ CagdCrvStruct *BspCrvDegreeRaiseN(CagdCrvStruct *Crv, int NewOrder) { CagdBType NewCrv = FALSE; int i, j, RaisedOrder, Length, KvLen1, Order = Crv -> Order, MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType); CagdCrvStruct *RaisedCrv, *UnitCrv; CagdRType *Kv; if (CAGD_IS_PERIODIC_CRV(Crv)) { NewCrv = TRUE; Crv = CnvrtPeriodic2FloatCrv(Crv); } Length = Crv -> Length; KvLen1 = Order + Length - 1; Kv = Crv -> KnotVector; if (NewOrder < Order) { CAGD_FATAL_ERROR(CAGD_ERR_WRONG_ORDER); return NULL; } RaisedOrder = NewOrder - Order + 1; UnitCrv = BspCrvNew(RaisedOrder, RaisedOrder, CAGD_MAKE_PT_TYPE(FALSE, MaxCoord)); for (i = 0; i < RaisedOrder * 2; i++) UnitCrv -> KnotVector[i] = i >= RaisedOrder ? Kv[KvLen1] : Kv[0]; for (i = 1; i <= MaxCoord; i++) for (j = 0; j < RaisedOrder; j++) UnitCrv -> Points[i][j] = 1.0; RaisedCrv = BspCrvMult(Crv, UnitCrv); CagdCrvFree(UnitCrv); CAGD_PROPAGATE_ATTR(RaisedCrv, Crv); if (NewCrv) CagdCrvFree(Crv); return RaisedCrv; } /***************************************************************************** * DESCRIPTION: M * Returns a new curve, identical to the original but with one degree higher. M * * * PARAMETERS: M * Crv: To raise it degree by one. M * * * RETURN VALUE: M * CagdCrvStruct *: A curve with one degree higher representing the same M * geometry as Crv. M * * * KEYWORDS: M * BspCrvDegreeRaise, degree raising M *****************************************************************************/ CagdCrvStruct *BspCrvDegreeRaise(CagdCrvStruct *Crv) { CagdBType NewCrv = FALSE, IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv); int i, i2, j, RaisedLen, Length, Order = Crv -> Order, MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType); CagdCrvStruct *RaisedCrv; if (CAGD_IS_PERIODIC_CRV(Crv)) { NewCrv = TRUE; Crv = CnvrtPeriodic2FloatCrv(Crv); } Length = Crv -> Length; if (Order > 2) return BspCrvDegreeRaiseN(Crv, Order + 1); /* If curve is linear, degree raising means basically to increase the */ /* knot multiplicity of each segment by one and add a middle point for */ /* each such segment. */ RaisedLen = Length * 2 - 1; RaisedCrv = BspCrvNew(RaisedLen, Order + 1, Crv -> PType); /* Update the control polygon; */ for (j = IsNotRational; j <= MaxCoord; j++) /* First point. */ RaisedCrv -> Points[j][0] = Crv -> Points[j][0]; for (i = 1, i2 = 1; i < Length; i++, i2 += 2) for (j = IsNotRational; j <= MaxCoord; j++) { RaisedCrv -> Points[j][i2] = Crv -> Points[j][i-1] * 0.5 + Crv -> Points[j][i] * 0.5; RaisedCrv -> Points[j][i2 + 1] = Crv -> Points[j][i]; } /* Update the knot vector. */ for (i = 0; i < 3; i++) RaisedCrv -> KnotVector[i] = Crv -> KnotVector[0]; for (i = 2, j = 3; i < Length; i++, j += 2) RaisedCrv -> KnotVector[j] = RaisedCrv -> KnotVector[j + 1] = Crv -> KnotVector[i]; for (i = j; i < j + 3; i++) RaisedCrv -> KnotVector[i] = Crv -> KnotVector[Length]; CAGD_PROPAGATE_ATTR(RaisedCrv, Crv); if (NewCrv) CagdCrvFree(Crv); return RaisedCrv; } /***************************************************************************** * DESCRIPTION: M * Returns a (unit) vector, equal to the tangent to Crv at parameter value t. M * Algorithm: insert (order - 1) knots and return control polygon tangent. M * The unnormalized normal does not equal dC/dt in its magnitude, only in M * its direction. M * * * PARAMETERS: M * Crv: Crv for which to compute a (unit) tangent. M * t: The parameter at which to compute the unit tangent. M * Normalize: If TRUE, attempt is made to normalize the returned vector. M * If FALSE, returned is an unnormalized vector in the right M * direction of the tangent. M * * * RETURN VALUE: M * CagdVecStruct *: A pointer to a static vector holding the tangent M * information. M * * * KEYWORDS: M * BspCrvTangent, tangent M *****************************************************************************/ CagdVecStruct *BspCrvTangent(CagdCrvStruct *Crv, CagdRType t, CagdBType Normalize) { STATIC_DATA CagdVecStruct Tan; CagdBType IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv); CagdVecStruct *T; CagdRType TMin, TMax, DKnots, Pt0[4], Pt1[4]; CagdPType P0, P1; CagdCrvStruct *FCrv = CAGD_IS_PERIODIC_CRV(Crv) ? CnvrtPeriodic2FloatCrv(Crv) : Crv; int i, k = FCrv -> Order, FlipPts = FALSE, Len = FCrv -> Length, OpenEnd = BspCrvHasOpenEC(FCrv), MaxCoord = MIN(CAGD_NUM_OF_PT_COORD(Crv -> PType), 3), Index = BspKnotLastIndexL(FCrv -> KnotVector, k + Len, t); CagdPointType PType = FCrv -> PType; CagdCrvStruct *RefinedCrv; /* Can not compute for constant curves. */ if (k <= 1) return NULL; CAGD_DOMAIN_GET_AND_VERIFY_CRV(t, FCrv, TMin, TMax); if (APX_EQ(t, TMin) && OpenEnd) { /* Use Crv starting tangent direction. */ if (IsNotRational) { CagdCoerceToE3(&Pt0[1], FCrv -> Points, 0, PType); CagdCoerceToE3(&Pt1[1], FCrv -> Points, 1, PType); } else { CagdCoerceToP3(Pt0, FCrv -> Points, 0, PType); CagdCoerceToP3(Pt1, FCrv -> Points, 1, PType); } DKnots = FCrv -> KnotVector[k] - FCrv -> KnotVector[k - 1]; } else if (APX_EQ(t, TMax) && OpenEnd) { /* Use Crv ending tangent direction. */ if (IsNotRational) { CagdCoerceToE3(&Pt0[1], FCrv -> Points, Len - 2, PType); CagdCoerceToE3(&Pt1[1], FCrv -> Points, Len - 1, PType); } else { CagdCoerceToP3(Pt0, FCrv -> Points, Len - 2, PType); CagdCoerceToP3(Pt1, FCrv -> Points, Len - 1, PType); } DKnots = FCrv -> KnotVector[Len] - FCrv -> KnotVector[Len - 1]; FlipPts = TRUE; /* Pt1 is the point on Crv(t). */ } else { int j, Mult = BspKnotFindMult(FCrv -> KnotVector, k, Len, t); /* Note (k - Mult - 1) could be negative in which case no refinement */ /* will actually take place: RefinedCrv will be copied from FCrv. */ RefinedCrv = BspCrvKnotInsertNSame(FCrv, t, k - Mult - 1); for (j = IsNotRational; j <= MaxCoord; j++) { Pt0[j] = RefinedCrv -> Points[j][Index - 1]; Pt1[j] = RefinedCrv -> Points[j][Index]; } for (j = MaxCoord + 1; j <= 3; j++) Pt0[j] = Pt1[j] = 0.0; DKnots = RefinedCrv -> KnotVector[Index + 1] - RefinedCrv -> KnotVector[Index]; FlipPts = TRUE; /* Pt1 is the point on Crv(t). */ CagdCrvFree(RefinedCrv); } PT_COPY(P0, &Pt0[1]); PT_COPY(P1, &Pt1[1]); if (IsNotRational) { PT_SUB(Tan.Vec, P1, P0); } else if (Pt0[0] == 0.0 || Pt1[0] == 0.0) { VEC_RESET(Tan.Vec); } else { PT_SCALE(P0, 1.0 / Pt0[0]); PT_SCALE(P1, 1.0 / Pt1[0]); PT_SUB(Tan.Vec, P1, P0); } if (CAGD_SQR_LEN_VECTOR(Tan) < IRIT_UEPS) { if (AttrGetIntAttrib(Crv -> Attr, "_tan") != TRUE) { /* Try to move a little. This location has zero speed. However, */ /* do it only once since we can be here forever. The "_tan" */ /* attribute guarantee we will try to move IRIT_EPS only once. */ AttrSetIntAttrib(&Crv -> Attr, "_tan", TRUE); T = BspCrvTangent(Crv, t - TMin < TMax - t ? t + IRIT_EPS : t - IRIT_EPS, Normalize); AttrFreeOneAttribute(&Crv -> Attr, "_tan"); if (FCrv != Crv) CagdCrvFree(FCrv); return T; } else { /* A zero length vector signals failure to compute tangent. */ if (FCrv != Crv) CagdCrvFree(FCrv); return &Tan; } } else { if (IsNotRational) { if (Normalize) CAGD_NORMALIZE_VECTOR(Tan) /* Normalize the vector. */ else VEC_SCALE(Tan.Vec, (k - 1) / DKnots); } else { /* Make P1 hold the derivative. P0 already holds the position. */ if (FlipPts) { for (i = 0; i < 4; i++) { Pt0[i] = (Pt1[i] - Pt0[i]) * (k - 1) / DKnots; SWAP(CagdRType, Pt1[i], Pt0[i]); } } else { for (i = 0; i < 4; i++) Pt1[i] = (Pt1[i] - Pt0[i]) * (k - 1) / DKnots; } /* And use to quotient rule to get the exact tangent. */ for (i = 1; i <= 3; i++) Tan.Vec[i - 1] = (Pt1[i] * Pt0[0] - Pt1[0] * Pt0[i]) / SQR(Pt0[0]); if (Normalize) CAGD_NORMALIZE_VECTOR(Tan); /* Normalize the vector. */ } if (FCrv != Crv) CagdCrvFree(FCrv); return &Tan; } } /***************************************************************************** * DESCRIPTION: M * Returns a (unit) vector, equal to the binormal to Crv at parameter value t.M * Algorithm: insert (order - 1) knots and using 3 consecutive control M * points at the refined location (p1, p2, p3), compute to binormal to be the M * cross product of the two vectors (p1 - p2) and (p2 - p3). M * Since a curve may have not BiNormal at inflection points or if the 3 M * points are collinear, NULL will be returned at such cases. M * * * PARAMETERS: M * Crv: Crv for which to compute a (unit) binormal. M * t: The parameter at which to compute the unit binormal. M * Normalize: If TRUE, attempt is made to normalize the returned vector. M * If FALSE, length is a function of given parametrization. M * * * RETURN VALUE: M * CagdVecStruct *: A pointer to a static vector holding the binormal M * information. M * * * KEYWORDS: M * BspCrvBiNormal, binormal M *****************************************************************************/ CagdVecStruct *BspCrvBiNormal(CagdCrvStruct *Crv, CagdRType t, CagdBType Normalize) { STATIC_DATA CagdVecStruct P3; CagdVecStruct P1, P2, *B; CagdRType TMin, TMax; CagdCrvStruct *FCrv = CAGD_IS_PERIODIC_CRV(Crv) ? CnvrtPeriodic2FloatCrv(Crv) : Crv; int k = FCrv -> Order, Len = FCrv -> Length, OpenEnd = BspCrvHasOpenEC(FCrv), Index = BspKnotLastIndexL(FCrv -> KnotVector, k + Len, t); CagdPointType PType = FCrv -> PType; CagdCrvStruct *RefinedCrv; CAGD_DOMAIN_GET_AND_VERIFY_CRV(t, FCrv, TMin, TMax); /* Can not compute for linear curves. */ if (k <= 2) return NULL; /* For planar curves, B is trivially the Z axis. */ if (CAGD_NUM_OF_PT_COORD(FCrv -> PType) == 2) { P3.Vec[0] = P3.Vec[1] = 0.0; P3.Vec[2] = 1.0; return &P3; } if (APX_EQ(t, TMin) && OpenEnd) { /* Use Crv starting tangent direction. */ CagdCoerceToE3(P1.Vec, FCrv -> Points, 0, PType); CagdCoerceToE3(P2.Vec, FCrv -> Points, 1, PType); CagdCoerceToE3(P3.Vec, FCrv -> Points, 2, PType); } else if (APX_EQ(t, TMax) && OpenEnd) { /* Use Crv ending tangent direction. */ CagdCoerceToE3(P1.Vec, FCrv -> Points, Len - 3, PType); CagdCoerceToE3(P2.Vec, FCrv -> Points, Len - 2, PType); CagdCoerceToE3(P3.Vec, FCrv -> Points, Len - 1, PType); } else { int i, Mult = BspKnotFindMult(FCrv -> KnotVector, k, Len, t); /* Note (k - Mult - 1) could be negative in which case no refinement */ /* will actually take place: RefinedCrv will be copied from FCrv. */ RefinedCrv = BspCrvKnotInsertNSame(FCrv, t, k - Mult - 1); CagdCoerceToE3(P1.Vec, RefinedCrv -> Points, Index, PType); CagdCoerceToE3(P2.Vec, RefinedCrv -> Points, Index + 1, PType); CagdCoerceToE3(P3.Vec, RefinedCrv -> Points, Index + 2, PType); for (i = Index - 1; PT_APX_EQ(P1.Vec, P2.Vec) && i >= 0; i--) CagdCoerceToE3(P1.Vec, RefinedCrv -> Points, i, PType); for (i = Index + 3; PT_APX_EQ(P2.Vec, P3.Vec) && i < RefinedCrv -> Length; i++) CagdCoerceToE3(P3.Vec, RefinedCrv -> Points, i, PType); CagdCrvFree(RefinedCrv); } CAGD_SUB_VECTOR(P1, P2); CAGD_SUB_VECTOR(P2, P3); CROSS_PROD(P3.Vec, P1.Vec, P2.Vec); if (CAGD_SQR_LEN_VECTOR(P3) < IRIT_UEPS) { if (AttrGetIntAttrib(Crv -> Attr, "_bn") != TRUE) { /* Try to move a little. This location has zero speed. However, */ /* do it only once since we can be here forever. The "_bn" */ /* attribute guarantee we will try to move IRIT_EPS only once. */ AttrSetIntAttrib(&Crv -> Attr, "_bn", TRUE); B = BspCrvBiNormal(Crv, t - TMin < TMax - t ? t + IRIT_EPS : t - IRIT_EPS, Normalize); AttrFreeOneAttribute(&Crv -> Attr, "_bn"); if (FCrv != Crv) CagdCrvFree(FCrv); return B; } else { /* A zero length vector signals failure to compute tangent. */ if (FCrv != Crv) CagdCrvFree(FCrv); return &P3; } } else { if (Normalize) CAGD_NORMALIZE_VECTOR(P3); /* Normalize the vector. */ if (FCrv != Crv) CagdCrvFree(FCrv); return &P3; } } /***************************************************************************** * DESCRIPTION: M * Returns a (unit) vector, equal to the normal of Crv at parameter value t. M * Algorithm: returns the cross product of the curve tangent and binormal. M * * * PARAMETERS: M * Crv: Crv for which to compute a (unit) normal. M * t: The parameter at which to compute the unit normal. M * Normalize: If TRUE, attempt is made to normalize the returned vector. M * If FALSE, length is a function of given parametrization. M * * * RETURN VALUE: M * CagdVecStruct *: A pointer to a static vector holding the normal M * information. M * * * KEYWORDS: M * BspCrvNormal, normal M *****************************************************************************/ CagdVecStruct *BspCrvNormal(CagdCrvStruct *Crv, CagdRType t, CagdBType Normalize) { STATIC_DATA CagdVecStruct N, *T, *B; T = BspCrvTangent(Crv, t, FALSE); B = BspCrvBiNormal(Crv, t, FALSE); if (T == NULL || B == NULL) return NULL; CROSS_PROD(N.Vec, B -> Vec, T -> Vec); if (Normalize) CAGD_NORMALIZE_VECTOR(N); /* Normalize the vector. */ return &N; } /***************************************************************************** * DESCRIPTION: M * Returns a new curve, equal to the given curve, differentiated once. M * Let old control polygon be P(i), i = 0 to k-1, and Q(i) be new one then: M * Q(i) = Degree * (P(i+1) - P(i)) / (Kv(i + k) - Kv(i + 1)), i = 0 to k-2. V * * * PARAMETERS: M * Crv: To differentiate. M * * * RETURN VALUE: M * CagdCrvStruct *: Differentiated curve. M * * * SEE ALSO: M * BzrCrvDerive, CagdCrvDerive, BzrCrvDeriveRational, BspCrvDeriveRational, M * BspCrvDeriveScalar M * * * KEYWORDS: M * BspCrvDerive, derivatives M *****************************************************************************/ CagdCrvStruct *BspCrvDerive(CagdCrvStruct *Crv) { CagdBType NewCrv = FALSE, IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv); int i, j, Len, k = Crv -> Order, MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType); CagdRType *Kv; CagdCrvStruct *DerivedCrv; if (CAGD_IS_PERIODIC_CRV(Crv)) { NewCrv = TRUE; Crv = CnvrtPeriodic2FloatCrv(Crv); } if (!GlblDeriveScalar && !IsNotRational) { DerivedCrv = BspCrvDeriveRational(Crv); if (NewCrv) CagdCrvFree(Crv); return DerivedCrv; } Len = Crv -> Length; Kv = Crv -> KnotVector; DerivedCrv = BspCrvNew(MAX(Len - 1, 1), MAX(k - 1, 1), Crv -> PType); if (k >= 2) { for (i = 0; i < Len - 1; i++) { CagdRType Denom = Kv[i + k] - Kv[i + 1]; if (APX_EQ_EPS(Denom, 0.0, IRIT_UEPS)) Denom = IRIT_UEPS; for (j = IsNotRational; j <= MaxCoord; j++) DerivedCrv -> Points[j][i] = (k - 1) * (Crv -> Points[j][i + 1] - Crv -> Points[j][i]) / Denom; } } else { for (i = 0; i < MAX(Len - 1, 1); i++) for (j = IsNotRational; j <= MaxCoord; j++) DerivedCrv -> Points[j][i] = 0.0; } CAGD_GEN_COPY(DerivedCrv -> KnotVector, &Crv -> KnotVector[k < 2 ? 0 : 1], sizeof(CagdRType) * (MAX(k - 1, 1) + MAX(Len - 1, 1))); if (NewCrv) CagdCrvFree(Crv); return DerivedCrv; } /***************************************************************************** * DESCRIPTION: M * Returns a new curve, equal to the given curve, differentiated once. M * Let old control polygon be P(i), i = 0 to k-1, and Q(i) be new one then: M * Q(i) = (k - 1) * (P(i+1) - P(i)) / (Kv(i + k) - Kv(i + 1)), i = 0 to k-2. V * For a Euclidean curve this is the same as CagdCrvDerive but for a M * rational curve the returned curve is not the vector field but simply the M * derivatives of all the curve's coefficients, including the weights. M * * * PARAMETERS: M * Crv: To differentiate. M * * * RETURN VALUE: M * CagdCrvStruct *: Differentiated curve. M * * * SEE ALSO: M * BzrCrvDerive, CagdCrvDerive, BzrCrvDeriveRational, BspCrvDeriveRational M * BspCrvDerive, BzrCrvDeriveScalar, CagdCrvDeriveScalar M * * * KEYWORDS: M * BspCrvDeriveScalar, derivatives M *****************************************************************************/ CagdCrvStruct *BspCrvDeriveScalar(CagdCrvStruct *Crv) { GlblDeriveScalar = TRUE; Crv = BspCrvDerive(Crv); GlblDeriveScalar = FALSE; return Crv; } /***************************************************************************** * DESCRIPTION: M * Returns a new Bspline curve, equal to the integral of the given Bspline M * crv. M * The given Bspline curve should be nonrational. M * V * l l l l+1 V * / /- - / - P - V * | | \ n \ | n \ i \ n+1 V * | C(t) = | / P B (t) = / P | B (t) = / ----- / ( t - t ) B (t) = V * / / - i i - i / i - n + 1 - j+n j j V * i=0 i=0 i=0 j=i+1 V * V * l+1 j-1 V * - - V * 1 \ \ n+1 V * = ----- / / P ( t - t ) B (t) V * n + 1 - - i j+n j j V * j=1 i=0 V * M * * * PARAMETERS: M * Crv: Curve to integrate. M * * * RETURN VALUE: M * CagdCrvStruct *: Integrated curve. M * * * SEE ALSO: M * BzrCrvIntegrate, BspSrfIntegrate, CagdCrvIntegrate M * * * KEYWORDS: M * BspCrvIntegrate, integrals M *****************************************************************************/ CagdCrvStruct *BspCrvIntegrate(CagdCrvStruct *Crv) { CagdBType NewCrv = FALSE; int i, j, k, Len, n = Crv -> Order, MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType); CagdCrvStruct *IntCrv; CagdRType *Kv; if (CAGD_IS_PERIODIC_CRV(Crv)) { NewCrv = TRUE; Crv = CnvrtPeriodic2FloatCrv(Crv); } if (CAGD_IS_RATIONAL_CRV(Crv)) CAGD_FATAL_ERROR(CAGD_ERR_RATIONAL_NO_SUPPORT); Len = Crv -> Length; Kv = Crv -> KnotVector; IntCrv = BspCrvNew(Len + 1, n + 1, Crv -> PType); /* Copy the knot vector and duplicate the two end knots. */ CAGD_GEN_COPY(&IntCrv -> KnotVector[1], Kv, sizeof(CagdRType) * (Len + n)); IntCrv -> KnotVector[0] = Kv[0]; IntCrv -> KnotVector[Len + n + 1] = Kv[Len + n - 1]; Kv = IntCrv -> KnotVector; for (k = 1; k <= MaxCoord; k++) { CagdRType *Points = Crv -> Points[k], *IntPoints = IntCrv -> Points[k]; for (j = 0; j < Len + 1; j++) { IntPoints[j] = 0.0; for (i = 0; i < j; i++) IntPoints[j] += Points[i] * (Kv[i + n + 1] - Kv[i + 1]); IntPoints[j] /= n; } } if (NewCrv) CagdCrvFree(Crv); return IntCrv; } /***************************************************************************** * DESCRIPTION: M * Apply the Moebius transformation to a rational Bspline curve. M * See "Moebius reparametrization of rational Bsplines", by Lee & Lucian, M * CAGD 8 (1991) pp 213-215. M * * * PARAMETERS: M * Crv: Curve to apply the Moebius transformation to. M * c: The scaling coefficient - c^n is the ratio between the first M * and last weight of the curve. M * If c == 0, the first and last weights are made equal. M * * * RETURN VALUE: M * CagdCrvStruct *: The modified curve with the same shape but different M * speed. M * * * SEE ALSO: M * BzrCrvMoebiusTransform, BspSrfMoebiusTransform M * * * KEYWORDS: M * BspCrvMoebiusTransform M *****************************************************************************/ CagdCrvStruct *BspCrvMoebiusTransform(CagdCrvStruct *Crv, CagdRType c) { int i, j, Order = Crv -> Order, Length = Crv -> Length, MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType); CagdRType c0, TMin, TMax, **Points, *KV, MaxW = IRIT_UEPS; if (!CAGD_IS_BSPLINE_CRV(Crv)) { CAGD_FATAL_ERROR(CAGD_ERR_WRONG_CRV); return NULL; } if (!CAGD_IS_RATIONAL_CRV(Crv)) { if (c == 1.0) return CagdCrvCopy(Crv); else Crv = CagdCoerceCrvTo(Crv, CAGD_MAKE_PT_TYPE(TRUE, MaxCoord)); } else Crv = CagdCrvCopy(Crv); Points = Crv -> Points; /* Save original curve domain and temporary map to zero to one. */ CagdCrvDomain(Crv, &TMin, &TMax); KV = Crv -> KnotVector; BspKnotAffineTransOrder2(KV, Order, CAGD_CRV_PT_LST_LEN(Crv) + Order, 0.0, 1.0); if (Points[0][0] == 0 || Points[0][Length - 1] == 0) { CAGD_FATAL_ERROR(CAGD_ERR_W_ZERO); return NULL; } if (c == 0.0) { c = Points[0][0] / Points[0][Length - 1]; c = pow(c, 1.0 / (Order - 1.0)); } /* Update the control points of the curve. */ for (i = 1, c0 = 1.0; i < Order; i++) c0 *= MOEBIUS_MEU(KV[i]); for (i = 0; i < Length; i++) { for (j = 0; j <= MaxCoord; j++) Points[j][i] /= c0; c0 *= MOEBIUS_MEU(KV[Order + i]) / MOEBIUS_MEU(KV[i + 1]); } /* Normalize all weights so largest has magnitude of one. */ for (i = 0; i < Length; i++) { if (MaxW < FABS(Points[0][i])) MaxW = FABS(Points[0][i]); } for (i = 0; i < Length; i++) { for (j = 0; j <= MaxCoord; j++) Points[j][i] /= MaxW; } /* Update the knot sequence of the curve. */ for (i = 0; i < CAGD_CRV_PT_LST_LEN(Crv) + Order; i++) KV[i] = MOEBIUS_REPARAM(KV[i]); BspKnotAffineTransOrder2(KV, Order, CAGD_CRV_PT_LST_LEN(Crv) + Order, TMin, TMax); return Crv; } /***************************************************************************** * DESCRIPTION: M * Returns a new linear Bspline curve constructed from the given polyline. M * * * PARAMETERS: M * Poly: To convert to a linear bspline curve. M * * * RETURN VALUE: M * CagdCrvStruct *: A linear Bspline curve representing Poly. M * * * SEE ALSO: M * UserPolylines2LinBsplineCrvs, CnvrtPolyline2LinBsplineCrv, M * CnvrtLinBsplineCrv2Polyline M * * * KEYWORDS: M * CnvrtPolyline2LinBsplineCrv, linear curves, conversion M *****************************************************************************/ CagdCrvStruct *CnvrtPolyline2LinBsplineCrv(CagdPolylineStruct *Poly) { int i, j, Length = Poly -> Length; CagdCrvStruct *NCrv, *Crv = BspCrvNew(Length, 2, CAGD_PT_E3_TYPE); CagdRType **Points = Crv -> Points; CagdPolylnStruct *Pts = Poly -> Polyline; for (i = j = 0; i < Length; i++, Pts++) { if (j > 0 && APX_EQ_EPS(Points[1][j - 1], Pts -> Pt[0], CAGD_SAME_PT_EPS) && APX_EQ_EPS(Points[2][j - 1], Pts -> Pt[1], CAGD_SAME_PT_EPS) && APX_EQ_EPS(Points[3][j - 1], Pts -> Pt[2], CAGD_SAME_PT_EPS)) { /* Must skip identical points. */ continue; } Points[1][j] = Pts -> Pt[0]; Points[2][j] = Pts -> Pt[1]; Points[3][j++] = Pts -> Pt[2]; } if (Crv -> Length != j) { int Size = sizeof(CagdRType) * j; NCrv = BspCrvNew(j, j > 1 ? 2 : 1, CAGD_PT_E3_TYPE); for (i = 1; i <= 3; i++) CAGD_GEN_COPY(NCrv -> Points[i], Crv -> Points[i], Size); CagdCrvFree(Crv); Crv = NCrv; } if (j == 1) { Crv -> Order = 1; BspKnotUniformOpen(1, 1, Crv -> KnotVector); } else { BspKnotUniformOpen(j, 2, Crv -> KnotVector); } return Crv; } /***************************************************************************** * DESCRIPTION: M * Returns a new polyline representing same geometry as the given linear M * Bspline curve. M * * * PARAMETERS: M * Crv: A linear Bspline curve to convert to a polyline. M * * * RETURN VALUE: M * CagdPolylineStruct *: A polyline same as linear curve Crv. M * * * SEE ALSO: M * UserPolylines2LinBsplineCrvs, UserPolyline2LinBsplineCrv, M * CnvrtPolyline2LinBsplineCrv M * * * KEYWORDS: M * CnvrtLinBsplineCrv2Polyline, linear curves, conversion M *****************************************************************************/ CagdPolylineStruct *CnvrtLinBsplineCrv2Polyline(CagdCrvStruct *Crv) { int i, j, Length = Crv -> Length; CagdPolylineStruct *Poly = CagdPolylineNew(Length); CagdRType **CrvPts = Crv -> Points; CagdPolylnStruct *PlPts = Poly -> Polyline; for (i = j = 0; i < Length; i++) { CagdPType Pt; CagdCoerceToE3(Pt, CrvPts, i, Crv -> PType); if (j > 0 && APX_EQ_EPS(Pt[0], PlPts[-1].Pt[0], CAGD_SAME_PT_EPS) && APX_EQ_EPS(Pt[1], PlPts[-1].Pt[1], CAGD_SAME_PT_EPS) && APX_EQ_EPS(Pt[2], PlPts[-1].Pt[2], CAGD_SAME_PT_EPS)) { /* Must skip identical points. */ continue; } PlPts -> Pt[0] = Pt[0]; PlPts -> Pt[1] = Pt[1]; PlPts -> Pt[2] = Pt[2]; PlPts++; j++; } Poly -> Length = j; return Poly; } /***************************************************************************** * DESCRIPTION: M * Converts a Bspline curve to a Bspline curve with floating end conditions. M * * * PARAMETERS: M * Crv: Bspline curve to convert to floating end conditions. Assume M * Crv is either periodic or has floating end condition. M * * * RETURN VALUE: M * CagdCrvStruct *: A Bspline curve with floating end conditions, M * representing the same geometry as Crv. M * * * KEYWORDS: M * CnvrtPeriodic2FloatCrv, conversion M *****************************************************************************/ CagdCrvStruct *CnvrtPeriodic2FloatCrv(CagdCrvStruct *Crv) { int i, Order = Crv -> Order, Length = Crv -> Length, MaxAxis = CAGD_NUM_OF_PT_COORD(Crv -> PType); CagdCrvStruct *NewCrv; if (!CAGD_IS_BSPLINE_CRV(Crv)) { CAGD_FATAL_ERROR(CAGD_ERR_BSP_CRV_EXPECT); return NULL; } if (!CAGD_IS_PERIODIC_CRV(Crv)) { CAGD_FATAL_ERROR(CAGD_ERR_PERIODIC_EXPECTED); return NULL; } NewCrv = BspCrvNew(Length + Order - 1, Order, Crv -> PType); CAGD_GEN_COPY(NewCrv -> KnotVector, Crv -> KnotVector, sizeof(CagdRType) * (Length + Order + Order - 1)); for (i = !CAGD_IS_RATIONAL_PT(Crv -> PType); i <= MaxAxis; i++) { CAGD_GEN_COPY(NewCrv -> Points[i], Crv -> Points[i], sizeof(CagdRType) * Length); CAGD_GEN_COPY(&NewCrv -> Points[i][Length], Crv -> Points[i], sizeof(CagdRType) * (Order - 1)); } for (i = MaxAxis + 1; i <= CAGD_MAX_PT_COORD; i++) NewCrv -> Points[i] = NULL; CAGD_PROPAGATE_ATTR(NewCrv, Crv); return NewCrv; } /***************************************************************************** * DESCRIPTION: M * Converts a float Bspline curve to a Bspline curve with open end M * conditions. M * * * PARAMETERS: M * Crv: Bspline curve to convert to open end conditions. M * * * RETURN VALUE: M * CagdCrvStruct *: A Bspline curve with open end conditions, M * representing the same geometry as Crv. M * * * KEYWORDS: M * CnvrtFloat2OpenCrv, conversion M *****************************************************************************/ CagdCrvStruct *CnvrtFloat2OpenCrv(CagdCrvStruct *Crv) { CagdRType TMin, TMax; CagdCrvStruct *NewCrv; if (!CAGD_IS_BSPLINE_CRV(Crv)) { CAGD_FATAL_ERROR(CAGD_ERR_BSP_CRV_EXPECT); return NULL; } CagdCrvDomain(Crv, &TMin, &TMax); NewCrv = CagdCrvRegionFromCrv(Crv, TMin, TMax); CAGD_PROPAGATE_ATTR(NewCrv, Crv); return NewCrv; } /***************************************************************************** * DESCRIPTION: M * Converts a Bspline curve to a Bspline curve with open end conditions. M * * * PARAMETERS: M * Crv: Bspline curve to convert to open end conditions. M * * * RETURN VALUE: M * CagdCrvStruct *: A Bspline curve with open end conditions, M * representing the same geometry as Crv. M * * * KEYWORDS: M * CnvrtBsp2OpenCrv, conversion M *****************************************************************************/ CagdCrvStruct *CnvrtBsp2OpenCrv(CagdCrvStruct *Crv) { CagdRType TMin, TMax; CagdCrvStruct *NewCrv; if (!CAGD_IS_BSPLINE_CRV(Crv)) { CAGD_FATAL_ERROR(CAGD_ERR_BSP_CRV_EXPECT); return NULL; } if (CAGD_IS_PERIODIC_CRV(Crv)) { CagdCrvStruct *TCrv = CnvrtPeriodic2FloatCrv(Crv); CagdCrvDomain(TCrv, &TMin, &TMax); NewCrv = CagdCrvRegionFromCrv(TCrv, TMin, TMax); CagdCrvFree(TCrv); } else if (!BspCrvHasOpenEC(Crv)) { CagdCrvDomain(Crv, &TMin, &TMax); NewCrv = CagdCrvRegionFromCrv(Crv, TMin, TMax); } else NewCrv = CagdCrvCopy(Crv); CAGD_PROPAGATE_ATTR(NewCrv, Crv); return NewCrv; } /***************************************************************************** * DESCRIPTION: M * Convert a curve into a list of control points. M * * * PARAMETERS: M * Crv: To the curve to convert to list of control points. M * * * RETURN VALUE: M * CagdCtlPtStruct *: List of control points of curve Crv. M * * * KEYWORDS: M * CnvrtCrvToCtlPts M *****************************************************************************/ CagdCtlPtStruct *CnvrtCrvToCtlPts(CagdCrvStruct *Crv) { CagdPointType PtType = Crv -> PType; int i, j, NumCoords = CAGD_NUM_OF_PT_COORD(PtType); CagdBType IsRational = CAGD_IS_RATIONAL_PT(PtType); CagdRType **Points = Crv -> Points; CagdCtlPtStruct *Pt, *CtlPts = NULL; for (i = 0; i < Crv -> Length; i++) { Pt = CagdCtlPtNew(PtType); for (j = !IsRational; j <= NumCoords; j++) { Pt -> Coords[j] = Points[j][i]; } LIST_PUSH(Pt, CtlPts); } return CagdListReverse(CtlPts); } /***************************************************************************** * DESCRIPTION: M * Reparameterize a curve to follow a desired parametrization. M * * * PARAMETERS: M * Crv: The curve to update its paraametrization. M * ParamType: The desired parametrization type: uniform, chord len., etc. M * * * RETURN VALUE: M * void M * * * SEE ALSO: M * BspCrvInterpBuildKVs, BspReparameterizeSrf M * * * KEYWORDS: M * BspReparameterizeCrv M *****************************************************************************/ void BspReparameterizeCrv(CagdCrvStruct *Crv, CagdParametrizationType ParamType) { CagdRType *NewPtKnots, *NewKV; CagdCtlPtStruct *CtlPts; if (!CAGD_IS_BSPLINE_CRV(Crv)) { CAGD_FATAL_ERROR(CAGD_ERR_BSP_CRV_EXPECT); return; } CtlPts = CnvrtCrvToCtlPts(Crv); BspCrvInterpBuildKVs(CtlPts, Crv -> Order, Crv -> Length, ParamType, Crv -> Periodic, &NewPtKnots, &NewKV); IritFree(NewPtKnots); CagdCtlPtFreeList(CtlPts); IritFree(Crv -> KnotVector); Crv -> KnotVector = NewKV; }