/****************************************************************************** * CBzr-Aux.c - Bezier curve auxilary routines. * ******************************************************************************* * (C) Gershon Elber, Technion, Israel Institute of Technology * ******************************************************************************* * Written by Gershon Elber, Mar. 90. * ******************************************************************************/ #include #include #include #include "cagd_loc.h" /***************************************************************************** * DESCRIPTION: M * Given a Bezier 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 * * * 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 * BzrCrvSubdivAtParam, subdivision, refinement M *****************************************************************************/ CagdCrvStruct *BzrCrvSubdivAtParam(CagdCrvStruct *Crv, CagdRType t) { CagdBType IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv); int i, j, l, k = Crv -> Length, MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType); CagdRType t1 = 1.0 - t; CagdCrvStruct *LCrv, *RCrv; LCrv = BzrCrvNew(k, Crv -> PType); RCrv = BzrCrvNew(k, Crv -> PType); /* Copy Curve into RCrv, so we can apply the recursive algo. to it. */ for (i = 0; i < k; i++) for (j = IsNotRational; j <= MaxCoord; j++) RCrv -> Points[j][i] = Crv -> Points[j][i]; for (j = IsNotRational; j <= MaxCoord; j++) LCrv -> Points[j][0] = Crv -> Points[j][0]; /* Apply the recursive algorithm to RCrv, and update LCrv with the */ /* temporary results. Note we updated the first point of LCrv above. */ for (i = 1; i < k; i++) { for (l = 0; l < k - i; l++) for (j = IsNotRational; j <= MaxCoord; j++) RCrv -> Points[j][l] = RCrv -> Points[j][l] * t1 + RCrv -> Points[j][l+1] * t; /* Copy temporary result to LCrv: */ for (j = IsNotRational; j <= MaxCoord; j++) LCrv -> Points[j][i] = RCrv -> Points[j][0]; } LCrv -> Pnext = RCrv; return LCrv; } /***************************************************************************** * 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 * BzrCrvDegreeRaiseN, degree raising M *****************************************************************************/ CagdCrvStruct *BzrCrvDegreeRaiseN(CagdCrvStruct *Crv, int NewOrder) { int i, j, RaisedOrder, Order = Crv -> Order, MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType); CagdCrvStruct *RaisedCrv, *UnitCrv; if (NewOrder < Order) { CAGD_FATAL_ERROR(CAGD_ERR_WRONG_ORDER); return NULL; } RaisedOrder = NewOrder - Order + 1; UnitCrv = BzrCrvNew(RaisedOrder, CAGD_MAKE_PT_TYPE(FALSE, MaxCoord)); for (i = 1; i <= MaxCoord; i++) for (j = 0; j < RaisedOrder; j++) UnitCrv -> Points[i][j] = 1.0; RaisedCrv = BzrCrvMult(Crv, UnitCrv); CagdCrvFree(UnitCrv); return RaisedCrv; } /***************************************************************************** * DESCRIPTION: M * Returns a new curve, identical to the original but with one degree higher. M * Let old control polygon be P(i), i = 0 to k-1, and Q(i) be new one then: M * i k-i V * Q(0) = P(0), Q(i) = --- P(i-1) + (---) P(i), Q(k) = P(k-1). V * k k V * * * PARAMETERS: M * Crv: To raise it degree by one. M * * * RETURN VALUE: M * CagdCrvStruct *: A curve of one order higher representing the same M * geometry as Crv. M * * * KEYWORDS: M * BzrCrvDegreeRaise, degree raising M *****************************************************************************/ CagdCrvStruct *BzrCrvDegreeRaise(CagdCrvStruct *Crv) { CagdBType IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv); int i, j, k = Crv -> Length, MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType); CagdCrvStruct *RaisedCrv = BzrCrvNew(k + 1, Crv -> PType); for (j = IsNotRational; j <= MaxCoord; j++) /* Q(0). */ RaisedCrv -> Points[j][0] = Crv -> Points[j][0]; for (i = 1; i < k; i++) /* Q(i). */ for (j = IsNotRational; j <= MaxCoord; j++) RaisedCrv -> Points[j][i] = Crv -> Points[j][i-1] * (i / ((CagdRType) k)) + Crv -> Points[j][i] * ((k - i) / ((CagdRType) k)); for (j = IsNotRational; j <= MaxCoord; j++) /* Q(k). */ RaisedCrv -> Points[j][k] = Crv -> Points[j][k-1]; return RaisedCrv; } /***************************************************************************** * DESCRIPTION: M * Returns a (unit) vector, equal to the tangent to Crv at parameter value t. M * Algorithm: pseudo subdivide Crv at t and using control point of M * subdivided curve find the tangent as the difference of the 2 end points. 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, length is a function of given parametrization. M * * * RETURN VALUE: M * CagdVecStruct *: A pointer to a static vector holding the tangent M * information. M * * * KEYWORDS: M * BzrCrvTangent, tangent M *****************************************************************************/ CagdVecStruct *BzrCrvTangent(CagdCrvStruct *Crv, CagdRType t, CagdBType Normalize) { static CagdVecStruct P2; CagdVecStruct P1, *T; CagdBType IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv); int i, j, l, k = Crv -> Length, MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType); CagdRType t1 = 1.0 - t; CagdCrvStruct *RCrv; if (APX_EQ(t, 0.0)) { /* Use Crv starting tangent direction. */ CagdCoerceToE3(P1.Vec, Crv -> Points, 0, Crv -> PType); CagdCoerceToE3(P2.Vec, Crv -> Points, 1, Crv -> PType); } else if (APX_EQ(t, 1.0)) { /* Use Crv ending tangent direction. */ CagdCoerceToE3(P1.Vec, Crv -> Points, k - 2, Crv -> PType); CagdCoerceToE3(P2.Vec, Crv -> Points, k - 1, Crv -> PType); } else { if (t < 0.0 || t > 1.0) CAGD_FATAL_ERROR(CAGD_ERR_T_NOT_IN_CRV); RCrv = BzrCrvNew(k, Crv -> PType); /* Copy Crv into RCrv, so we can apply the recursive algo. to it. */ for (i = 0; i < k; i++) for (j = IsNotRational; j <= MaxCoord; j++) RCrv -> Points[j][i] = Crv -> Points[j][i]; /* Apply the recursive algorithm to RCrv. */ for (i = 1; i < k; i++) for (l = 0; l < k - i; l++) for (j = IsNotRational; j <= MaxCoord; j++) RCrv -> Points[j][l] = RCrv -> Points[j][l] * t1 + RCrv -> Points[j][l+1] * t; CagdCoerceToE3(P1.Vec, RCrv -> Points, 0, RCrv -> PType); CagdCoerceToE3(P2.Vec, RCrv -> Points, 1, RCrv -> PType); CagdCrvFree(RCrv); } CAGD_SUB_VECTOR(P2, P1); if (!Normalize) return &P2; if (CAGD_LEN_VECTOR(P2) < 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 = BzrCrvTangent(Crv, t < 0.5 ? t + IRIT_EPS : t - IRIT_EPS, Normalize); AttrFreeOneAttribute(&Crv -> Attr, "_tan"); return T; } else { /* A zero length vector signals failure to compute tangent. */ return &P2; } } else { CAGD_NORMALIZE_VECTOR(P2); /* Normalize the vector. */ return &P2; } } /***************************************************************************** * 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 colinear, 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 * BzrCrvBiNormal, binormal M *****************************************************************************/ CagdVecStruct *BzrCrvBiNormal(CagdCrvStruct *Crv, CagdRType t, CagdBType Normalize) { static CagdVecStruct P3; CagdVecStruct P1, P2; CagdBType IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv); int i, j, l, k = Crv -> Length, MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType); CagdRType t1 = 1.0 - t; CagdCrvStruct *RCrv; /* 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(Crv -> PType) == 2) { P3.Vec[0] = P3.Vec[1] = 0.0; P3.Vec[2] = 1.0; return &P3; } if (APX_EQ(t, 0.0)) { /* Use Crv starting tangent direction. */ CagdCoerceToE3(P1.Vec, Crv -> Points, 0, Crv -> PType); CagdCoerceToE3(P2.Vec, Crv -> Points, 1, Crv -> PType); CagdCoerceToE3(P3.Vec, Crv -> Points, 2, Crv -> PType); } else if (APX_EQ(t, 1.0)) { /* Use Crv ending tangent direction. */ CagdCoerceToE3(P1.Vec, Crv -> Points, k - 3, Crv -> PType); CagdCoerceToE3(P2.Vec, Crv -> Points, k - 2, Crv -> PType); CagdCoerceToE3(P3.Vec, Crv -> Points, k - 1, Crv -> PType); } else { if (t < 0.0 || t > 1.0) CAGD_FATAL_ERROR(CAGD_ERR_T_NOT_IN_CRV); RCrv = BzrCrvNew(k, Crv -> PType); /* Copy Crv into RCrv, so we can apply the recursive algo. to it. */ for (i = 0; i < k; i++) for (j = IsNotRational; j <= MaxCoord; j++) RCrv -> Points[j][i] = Crv -> Points[j][i]; /* Apply the recursive algorithm to RCrv. */ for (i = 1; i < k; i++) for (l = 0; l < k - i; l++) for (j = IsNotRational; j <= MaxCoord; j++) RCrv -> Points[j][l] = RCrv -> Points[j][l] * t1 + RCrv -> Points[j][l+1] * t; CagdCoerceToE3(P1.Vec, RCrv -> Points, 0, RCrv -> PType); CagdCoerceToE3(P2.Vec, RCrv -> Points, 1, RCrv -> PType); CagdCoerceToE3(P3.Vec, RCrv -> Points, 2, RCrv -> PType); CagdCrvFree(RCrv); } CAGD_SUB_VECTOR(P1, P2); CAGD_SUB_VECTOR(P2, P3); CROSS_PROD(P3.Vec, P1.Vec, P2.Vec); if (Normalize) { if ((t = CAGD_LEN_VECTOR(P3)) < IRIT_UEPS) return NULL; else CAGD_DIV_VECTOR(P3, t); /* Normalize the vector. */ } 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 * BzrCrvNoraml, normal M *****************************************************************************/ CagdVecStruct *BzrCrvNormal(CagdCrvStruct *Crv, CagdRType t, CagdBType Normalize) { static CagdVecStruct N; CagdVecStruct *T, *B; T = BzrCrvTangent(Crv, t, FALSE); B = BzrCrvBiNormal(Crv, t, FALSE); if (T == NULL || B == NULL) return NULL; CROSS_PROD(N.Vec, T -> Vec, B -> 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) = (k - 1) * (P(i+1) - P(i)), i = 0 to k-2. M * * * PARAMETERS: M * Crv: To differentiate. M * * * RETURN VALUE: M * CagdCrvStruct *: Differentiated curve. M * * * SEE ALSO: M * CagdCrvDerive, BspCrvDerive, BzrCrvDeriveRational, BspCrvDeriveRational M * * * KEYWORDS: M * BzrCrvDerive, derivatives M *****************************************************************************/ CagdCrvStruct *BzrCrvDerive(CagdCrvStruct *Crv) { CagdBType IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv); int i, j, k = Crv -> Length, MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType); CagdCrvStruct *DerivedCrv; if (!IsNotRational) return BzrCrvDeriveRational(Crv); DerivedCrv = BzrCrvNew(MAX(1, k - 1), Crv -> PType); if (k >= 2) { for (i = 0; i < k - 1; i++) for (j = IsNotRational; j <= MaxCoord; j++) DerivedCrv -> Points[j][i] = (k - 1) * (Crv -> Points[j][i+1] - Crv -> Points[j][i]); } else { for (j = IsNotRational; j <= MaxCoord; j++) DerivedCrv -> Points[j][0] = 0.0; } return DerivedCrv; } /***************************************************************************** * DESCRIPTION: M * Returns a new Bezier curve, equal to the integral of the given Bezier M * crv. M * The given Bezier curve should be nonrational. M * V * n n n n+1 V * / /- - / - P - V * | | \ n \ | n \ i \ n+1 V * | C(t) = | / P B (t) = / P | B (t) = / ----- / B (t) = V * / / - i i - i / i - n + 1 - j V * i=0 i=0 i=0 j=i+1 V * V * n+1 j-1 V * - - V * 1 \ \ n+1 V * = ----- / / P B (t) V * n + 1 - - i j V * j=1 i=0 V * V * M * * * PARAMETERS: M * Crv: Curve to integrate. M * * * RETURN VALUE: M * CagdCrvStruct *: Integrated curve. M * * * SEE ALSO: M * BspCrvIntegrate, CagdCrvIntegrate M * * * KEYWORDS: M * BzrCrvIntegrate, integrals M *****************************************************************************/ CagdCrvStruct *BzrCrvIntegrate(CagdCrvStruct *Crv) { int i, j, k, n = Crv -> Length, MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType); CagdCrvStruct *IntCrv; if (CAGD_IS_RATIONAL_CRV(Crv)) CAGD_FATAL_ERROR(CAGD_ERR_RATIONAL_NO_SUPPORT); IntCrv = BzrCrvNew(n + 1, Crv -> PType); for (k = 1; k <= MaxCoord; k++) { CagdRType *Points = Crv -> Points[k], *IntPoints = IntCrv -> Points[k]; for (j = 0; j < n + 1; j++) { IntPoints[j] = 0.0; for (i = 0; i < j; i++) IntPoints[j] += Points[i]; IntPoints[j] /= n; } } return IntCrv; } /***************************************************************************** * DESCRIPTION: M * Converts a Bezier curve into Bspline curve by adding an open knot vector. M * * * PARAMETERS: M * Crv: A Bezier curve to convert to a Bspline curve. M * * * RETURN VALUE: M * CagdCrvStruct *: A Bspline curve representing Bezier curve Crv. M * * * SEE ALSO: M * CnvrtBspline2BezierCrv, CnvrtBezier2PowerCrv, CnvrtPower2BezierCrv M * * * KEYWORDS: M * CnvrtBezier2BsplineCrv, conversion M *****************************************************************************/ CagdCrvStruct *CnvrtBezier2BsplineCrv(CagdCrvStruct *Crv) { CagdCrvStruct *BspCrv; if (Crv -> GType != CAGD_CBEZIER_TYPE) { CAGD_FATAL_ERROR(CAGD_ERR_WRONG_CRV); return NULL; } BspCrv = CagdCrvCopy(Crv); BspCrv -> Order = BspCrv -> Length; BspCrv -> KnotVector = BspKnotUniformOpen(BspCrv -> Length, BspCrv -> Order, NULL); BspCrv -> GType = CAGD_CBSPLINE_TYPE; return BspCrv; } /***************************************************************************** * DESCRIPTION: M * Converts a Bspline curve into a set of Bezier curves by subdividing the M * Bspline curve at all its internal knots. M * Returned is a list of Bezier curves. M * * * PARAMETERS: M * Crv: A Bspline curve to convert to a Bezier curve. M * * * RETURN VALUE: M * CagdCrvStruct *: A list of Bezier curves representing the Bspline M * curve Crv. M * * * SEE ALSO: M * CnvrtBezier2BsplineCrv, CnvrtBezier2PowerCrv, CnvrtPower2BezierCrv M * * * KEYWORDS: M * CnvrtBezier2BsplineCrv, conversion M *****************************************************************************/ CagdCrvStruct *CnvrtBspline2BezierCrv(CagdCrvStruct *Crv) { CagdBType NewCrv = FALSE; int i, Order, Length; CagdRType LastT, *KnotVector; CagdCrvStruct *OrigCrv, *BezierCrvs = NULL; if (Crv -> GType != CAGD_CBSPLINE_TYPE) { CAGD_FATAL_ERROR(CAGD_ERR_WRONG_CRV); return NULL; } if (CAGD_IS_PERIODIC_CRV(Crv)) { NewCrv = TRUE; Crv = CnvrtPeriodic2FloatCrv(Crv); } if (CAGD_IS_BSPLINE_CRV(Crv) && !BspCrvHasOpenEC(Crv)) { CagdCrvStruct *TCrv = BspCrvOpenEnd(Crv); if (NewCrv) CagdCrvFree(Crv); Crv = TCrv; NewCrv = TRUE; } Order = Crv -> Order, Length = Crv -> Length; KnotVector = Crv -> KnotVector; OrigCrv = Crv; for (i = Length - 1, LastT = KnotVector[Length]; i >= Order; i--) { CagdRType t = KnotVector[i]; if (!APX_EQ(LastT, t)) { CagdCrvStruct *Crvs = BspCrvSubdivAtParam(Crv, t); if (Crv != OrigCrv) CagdCrvFree(Crv); Crvs -> Pnext -> Pnext = BezierCrvs; BezierCrvs = Crvs -> Pnext; Crv = Crvs; Crv -> Pnext = NULL; LastT = t; } } if (Crv == OrigCrv) { /* No interior knots in this curve - just copy it: */ BezierCrvs = CagdCrvCopy(Crv); } else { Crv -> Pnext = BezierCrvs; BezierCrvs = Crv; } for (Crv = BezierCrvs; Crv != NULL; Crv = Crv -> Pnext) { Crv -> GType = CAGD_CBEZIER_TYPE; Crv -> Length = Crv -> Order; IritFree((VoidPtr) Crv -> KnotVector); Crv -> KnotVector = NULL; } if (NewCrv) CagdCrvFree(Crv); return BezierCrvs; }