/****************************************************************************** * SBsp-Aux.c - Bspline surface auxilary routines. * ******************************************************************************* * (C) Gershon Elber, Technion, Israel Institute of Technology * ******************************************************************************* * Written by Gershon Elber, July. 90. * ******************************************************************************/ #include #include #include #include "cagd_loc.h" /* Define some marcos to make some of the routines below look better. They */ /* calculate the index of the U, V point of the control mesh in Points. */ #define DERIVED_SRF(U, V) CAGD_MESH_UV(DerivedSrf, U, V) #define RAISED_SRF(U, V) CAGD_MESH_UV(RaisedSrf, U, V) #define SRF(U, V) CAGD_MESH_UV(Srf, U, V) #define NORMAL_IRIT_EPS 1e-4 static CagdBType EvalTangentVector(CagdVType Vec, CagdCrvStruct *Crv, CagdRType t, CagdRType TMin, CagdRType TMax); /***************************************************************************** * DESCRIPTION: M * Given a Bspline surface - subdivides it into two sub-surfaces at the given M * parametric value. M * Returns pointer to first surface in a list of two subdivided surfaces. M * * * PARAMETERS: M * Srf: To subdivide at parameter value t. M * t: Parameter value to subdivide Srf at. M * Dir: Direction of subdivision. Either U or V. M * * * RETURN VALUE: M * CagdSrfStruct *: A list of the two subdivided surfaces. M * * * SEE ALSO: M * CagdSrfSubdivAtParam, BzrSrfSubdivAtParam, TrimSrfSubdivAtParam M * * * KEYWORDS: M * BspSrfSubdivAtParam, subdivision, refinement M *****************************************************************************/ CagdSrfStruct *BspSrfSubdivAtParam(CagdSrfStruct *Srf, CagdRType t, CagdSrfDirType Dir) { CagdBType NewSrf = FALSE, IsNotRational = !CAGD_IS_RATIONAL_CRV(Srf); int i, j, Row, Col, KVLen, Index1, Index2, Mult, LULength, RULength, LVLength, RVLength, ULength, VLength, UOrder = Srf -> UOrder, VOrder = Srf -> VOrder, MaxCoord = CAGD_NUM_OF_PT_COORD(Srf -> PType); CagdRType *RefKV, **Pts, **LPts, **RPts, *LOnePts, *ROnePts, *OnePts, UMin, UMax, VMin, VMax; CagdSrfStruct *RSrf, *LSrf; BspKnotAlphaCoeffType *A; if (CAGD_IS_PERIODIC_SRF(Srf)) { NewSrf = TRUE; Srf = CnvrtPeriodic2FloatSrf(Srf); } ULength = Srf -> ULength; VLength = Srf -> VLength; switch (Dir) { case CAGD_CONST_U_DIR: RefKV = Srf -> UKnotVector; KVLen = UOrder + ULength; Index1 = BspKnotLastIndexL(RefKV, KVLen, t); if (Index1 + 1 < UOrder) Index1 = UOrder - 1; Index2 = BspKnotFirstIndexG(RefKV, KVLen, t); if (Index2 > ULength) Index2 = ULength; LSrf = BspSrfNew(Index1 + 1, VLength, UOrder, VOrder, Srf -> PType); RSrf = BspSrfNew(ULength - Index2 + UOrder, VLength, UOrder, VOrder, Srf -> PType); Mult = UOrder - 1 - (Index2 - Index1 - 1); /* Update the new knot vectors. */ CAGD_GEN_COPY(LSrf -> UKnotVector, Srf -> UKnotVector, sizeof(CagdRType) * (Index1 + 1)); /* Close the knot vector with multiplicity Order: */ for (j = Index1 + 1; j <= Index1 + UOrder; j++) LSrf -> UKnotVector[j] = t; CAGD_GEN_COPY(&RSrf -> UKnotVector[UOrder], &Srf -> UKnotVector[Index2], sizeof(CagdRType) * (ULength + UOrder - Index2)); /* Make sure knot vector starts with multiplicity Order: */ for (j = 0; j < UOrder; j++) RSrf -> UKnotVector[j] = t; /* And copy the other direction knot vectors. */ CAGD_GEN_COPY(LSrf -> VKnotVector, Srf -> VKnotVector, sizeof(CagdRType) * (VOrder + VLength)); CAGD_GEN_COPY(RSrf -> VKnotVector, Srf -> VKnotVector, sizeof(CagdRType) * (VOrder + VLength)); break; case CAGD_CONST_V_DIR: RefKV = Srf -> VKnotVector; KVLen = VOrder + VLength; Index1 = BspKnotLastIndexL(RefKV, KVLen, t); if (Index1 + 1 < VOrder) Index1 = VOrder - 1; Index2 = BspKnotFirstIndexG(RefKV, KVLen, t); if (Index2 > VLength) Index2 = VLength; LSrf = BspSrfNew(ULength, Index1 + 1, UOrder, VOrder, Srf -> PType); RSrf = BspSrfNew(ULength, VLength - Index2 + VOrder, UOrder, VOrder, Srf -> PType); Mult = VOrder - 1 - (Index2 - Index1 - 1); /* Update the new knot vectors. */ CAGD_GEN_COPY(LSrf -> VKnotVector, Srf -> VKnotVector, sizeof(CagdRType) * (Index1 + 1)); /* Close the knot vector with multiplicity Order: */ for (j = Index1 + 1; j <= Index1 + VOrder; j++) LSrf -> VKnotVector[j] = t; CAGD_GEN_COPY(&RSrf -> VKnotVector[VOrder], &Srf -> VKnotVector[Index2], sizeof(CagdRType) * (VLength + VOrder - Index2)); /* Make sure knot vector starts with multiplicity Order: */ for (j = 0; j < VOrder; j++) RSrf -> VKnotVector[j] = t; /* And copy the other direction knot vectors. */ CAGD_GEN_COPY(LSrf -> UKnotVector, Srf -> UKnotVector, sizeof(CagdRType) * (UOrder + ULength)); CAGD_GEN_COPY(RSrf -> UKnotVector, Srf -> UKnotVector, sizeof(CagdRType) * (UOrder + ULength)); break; default: Mult = 1; RefKV = NULL; LSrf = RSrf = NULL; CAGD_FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV); break; } Pts = Srf -> Points; LPts = LSrf -> Points; RPts = RSrf -> Points; LULength = LSrf -> ULength, RULength = RSrf -> ULength; LVLength = LSrf -> VLength, RVLength = RSrf -> VLength; CagdSrfDomain(Srf, &UMin, &UMax, &VMin, &VMax); switch (Dir) { case CAGD_CONST_U_DIR: /* Compute the Alpha refinement matrix. */ if (Mult > 0) { CagdRType *NewKV = (CagdRType *) IritMalloc(sizeof(CagdRType) * Mult); CAGD_DOMAIN_T_VERIFY(t, UMin, UMax); if (t == UMax) t -= CAGD_DOMAIN_IRIT_EPS; for (i = 0; i < Mult; i++) NewKV[i] = t; A = BspKnotEvalAlphaCoefMerge(UOrder, RefKV, ULength, NewKV, Mult); IritFree((VoidPtr) NewKV); } else A = BspKnotEvalAlphaCoefMerge(UOrder, RefKV, ULength, NULL, 0); /* Now work on the two new surfaces meshes. */ /* 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; for (Row = 0; Row < VLength; Row++) { /* Blend Srf into LSrf. */ for (j = IsNotRational; j <= MaxCoord; j++) { LOnePts = &LPts[j][Row * LULength]; OnePts = &Pts[j][Row * ULength]; for (i = 0; i < LULength; i++, LOnePts++) CAGD_ALPHA_BLEND(A, i, OnePts, Srf -> ULength, LOnePts); } /* Blend Srf into LSrf. */ for (j = IsNotRational; j <= MaxCoord; j++) { ROnePts = &RPts[j][Row * RULength]; OnePts = &Pts[j][Row * ULength]; for (i = LSrf -> ULength - 1 + Mult; i < LSrf -> ULength + RSrf -> ULength - 1 + Mult; i++, ROnePts++) CAGD_ALPHA_BLEND(A, i, OnePts, Srf -> ULength, ROnePts); } } BspKnotFreeAlphaCoef(A); break; case CAGD_CONST_V_DIR: /* Compute the Alpha refinement matrix. */ if (Mult > 0) { CagdRType *NewKV = (CagdRType *) IritMalloc(sizeof(CagdRType) * Mult); CAGD_DOMAIN_T_VERIFY(t, VMin, VMax); if (t == VMax) t -= CAGD_DOMAIN_IRIT_EPS; for (i = 0; i < Mult; i++) NewKV[i] = t; A = BspKnotEvalAlphaCoefMerge(VOrder, RefKV, VLength, NewKV, Mult); IritFree((VoidPtr) NewKV); } else A = BspKnotEvalAlphaCoefMerge(VOrder, RefKV, VLength, NULL, 0); /* Now work on the two new surfaces meshes. */ /* 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; for (Col = 0; Col < ULength; Col++) { LULength = LSrf -> ULength; RULength = RSrf -> ULength; /* Blend Srf into LSrf. */ for (j = IsNotRational; j <= MaxCoord; j++) { LOnePts = &LPts[j][Col]; OnePts = &Pts[j][Col]; for (i = 0; i < LVLength; i++, LOnePts += LULength) CAGD_ALPHA_BLEND_STEP(A, i, OnePts, Srf -> VLength, LOnePts, ULength); } /* Blend Srf into RSrf. */ for (j = IsNotRational; j <= MaxCoord; j++) { ROnePts = &RPts[j][Col]; OnePts = &Pts[j][Col]; for (i = LVLength - 1 + Mult; i < LVLength + RVLength - 1 + Mult; i++, ROnePts += RULength) CAGD_ALPHA_BLEND_STEP(A, i, OnePts, Srf -> VLength, ROnePts, ULength); } } BspKnotFreeAlphaCoef(A); break; default: CAGD_FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV); break; } BspKnotMakeRobustKV(RSrf -> UKnotVector, RSrf -> UOrder + RSrf -> ULength); BspKnotMakeRobustKV(RSrf -> VKnotVector, RSrf -> VOrder + RSrf -> VLength); BspKnotMakeRobustKV(LSrf -> UKnotVector, LSrf -> UOrder + LSrf -> ULength); BspKnotMakeRobustKV(LSrf -> VKnotVector, LSrf -> VOrder + LSrf -> VLength); LSrf -> Pnext = RSrf; if (NewSrf) CagdSrfFree(Srf); return LSrf; } /***************************************************************************** * DESCRIPTION: M * Inserts n knot, all with the value t in direction Dir. In no case will M * the multiplicity of a knot be greater or equal to the curve order. M * * * PARAMETERS: M * Srf: To refine by insertion (upto) n knot of value t. M * Dir: Direction of refinement. Either U or V. M * t: Parameter value of new knot to insert. M * n: Maximum number of times t should be inserted. M * * * RETURN VALUE: M * CagdSrfStruct *: Refined Srf with n knots of value t in direction Dir. M * * * KEYWORDS: M * BspSrfKnotInsertNSame, refinement, subdivision M *****************************************************************************/ CagdSrfStruct *BspSrfKnotInsertNSame(CagdSrfStruct *Srf, CagdSrfDirType Dir, CagdRType t, int n) { CagdBType NewSrf = FALSE; int i, CrntMult, Mult; CagdSrfStruct *RefinedSrf; if (CAGD_IS_PERIODIC_SRF(Srf)) { NewSrf = TRUE; Srf = CnvrtPeriodic2FloatSrf(Srf); } switch (Dir) { case CAGD_CONST_U_DIR: CrntMult = BspKnotFindMult(Srf -> UKnotVector, Srf -> UOrder, Srf -> ULength, t), Mult = MIN(n, Srf -> UOrder - CrntMult - 1); break; case CAGD_CONST_V_DIR: CrntMult = BspKnotFindMult(Srf -> VKnotVector, Srf -> VOrder, Srf -> VLength, t), Mult = MIN(n, Srf -> VOrder - CrntMult - 1); break; default: Mult = 0; CAGD_FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV); break; } if (Mult > 0) { CagdRType *NewKV = (CagdRType *) IritMalloc(sizeof(CagdRType) * Mult); for (i = 0; i < Mult; i++) NewKV[i] = t; RefinedSrf = BspSrfKnotInsertNDiff(Srf, Dir, FALSE, NewKV, Mult); IritFree((VoidPtr) NewKV); } else { RefinedSrf = CagdSrfCopy(Srf); } if (NewSrf) CagdSrfFree(Srf); return RefinedSrf; } /***************************************************************************** * 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 * Srf: To refine by insertion (upto) n knot of value t. M * Dir: Direction of refinement. Either U or V. M * Replace: if TRUE, the n knots in t should replace the knot vector M * of size n of Srf. Sizes must match. If False, n new knots M * as defined by t will be introduced into Srf. M * t: New knots to introduce/replace knot vector of Srf. M * n: Size of t. M * * * RETURN VALUE: M * CagdSrfStruct *: Refined Srf with n new knots in direction Dir. M * * * KEYWORDS: M * BspSrfKnotInsertNDiff, refinement, subdivision M *****************************************************************************/ CagdSrfStruct *BspSrfKnotInsertNDiff(CagdSrfStruct *Srf, CagdSrfDirType Dir, int Replace, CagdRType *t, int n) { CagdBType NewSrf = FALSE, IsNotRational = !CAGD_IS_RATIONAL_SRF(Srf); int i, Row, Col, ULength, VLength, UOrder = Srf -> UOrder, VOrder = Srf -> VOrder, MaxCoord = CAGD_NUM_OF_PT_COORD(Srf -> PType); CagdSrfStruct *RefSrf = NULL; if (CAGD_IS_PERIODIC_SRF(Srf)) { NewSrf = TRUE; Srf = CnvrtPeriodic2FloatSrf(Srf); } ULength = Srf -> ULength; VLength = Srf -> VLength; if (Replace) { for (i = 1; i < n; i++) if (t[i] < t[i - 1]) CAGD_FATAL_ERROR(CAGD_ERR_KNOT_NOT_ORDERED); switch (Dir) { case CAGD_CONST_U_DIR: if (Srf -> UOrder + Srf -> ULength != n) CAGD_FATAL_ERROR(CAGD_ERR_NUM_KNOT_MISMATCH); RefSrf = CagdSrfCopy(Srf); for (i = 0; i < n; i++) RefSrf -> UKnotVector[i] = *t++; break; case CAGD_CONST_V_DIR: if (Srf -> VOrder + Srf -> VLength != n) CAGD_FATAL_ERROR(CAGD_ERR_NUM_KNOT_MISMATCH); RefSrf = CagdSrfCopy(Srf); for (i = 0; i < n; i++) RefSrf -> VKnotVector[i] = *t++; break; default: CAGD_FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV); break; } } else if (n == 0) { RefSrf = CagdSrfCopy(Srf); } else { int j, LengthKVt, RULength, RVLength; BspKnotAlphaCoeffType *A; CagdRType *MergedKVt, UMin, UMax, VMin, VMax, *UKnotVector = Srf -> UKnotVector, *VKnotVector = Srf -> VKnotVector; CagdSrfDomain(Srf, &UMin, &UMax, &VMin, &VMax); for (i = 1; i < n; i++) if (t[i] < t[i - 1]) CAGD_FATAL_ERROR(CAGD_ERR_KNOT_NOT_ORDERED); switch (Dir) { case CAGD_CONST_U_DIR: for (i = 0; i < n; i++) { CAGD_DOMAIN_T_VERIFY(t[i], UMin, UMax); if (t[i] == UMax) t[i] -= CAGD_DOMAIN_IRIT_EPS; } /* Compute the Alpha refinement matrix. */ MergedKVt = BspKnotMergeTwo(UKnotVector, ULength + UOrder, t, n, 0, &LengthKVt); A = BspKnotEvalAlphaCoef(UOrder, UKnotVector, ULength, MergedKVt, LengthKVt - UOrder); RefSrf = BspSrfNew(ULength + n, VLength, UOrder, VOrder, Srf -> PType); IritFree((VoidPtr) RefSrf -> UKnotVector); IritFree((VoidPtr) RefSrf -> VKnotVector); RefSrf -> UKnotVector = MergedKVt; RefSrf -> VKnotVector = BspKnotCopy(Srf -> VKnotVector, Srf -> VLength + Srf -> VOrder); RULength = RefSrf -> ULength; /* Update the control mesh */ for (Row = 0; Row < VLength; Row++) { for (j = IsNotRational; j <= MaxCoord; j++) { CagdRType *ROnePts = &RefSrf -> Points[j][Row * RULength], *OnePts = &Srf -> Points[j][Row * ULength]; for (i = 0; i < RULength; i++, ROnePts++) CAGD_ALPHA_BLEND(A, i, OnePts, Srf -> ULength, ROnePts); } } BspKnotFreeAlphaCoef(A); break; case CAGD_CONST_V_DIR: for (i = 0; i < n; i++) { CAGD_DOMAIN_T_VERIFY(t[i], VMin, VMax); if (t[i] == VMax) t[i] -= CAGD_DOMAIN_IRIT_EPS; } /* Compute the Alpha refinement matrix. */ MergedKVt = BspKnotMergeTwo(VKnotVector, VLength + VOrder, t, n, 0, &LengthKVt); A = BspKnotEvalAlphaCoef(VOrder, VKnotVector, VLength, MergedKVt, LengthKVt - VOrder); RefSrf = BspSrfNew(ULength, VLength + n, UOrder, VOrder, Srf -> PType); IritFree((VoidPtr) RefSrf -> UKnotVector); IritFree((VoidPtr) RefSrf -> VKnotVector); RefSrf -> UKnotVector = BspKnotCopy(Srf -> UKnotVector, Srf -> ULength + Srf -> UOrder); RefSrf -> VKnotVector = MergedKVt; RULength = RefSrf -> ULength; RVLength = RefSrf -> VLength; /* Update the control mesh */ for (Col = 0; Col < ULength; Col++) { for (j = IsNotRational; j <= MaxCoord; j++) { CagdRType *ROnePts = &RefSrf -> Points[j][Col], *OnePts = &Srf -> Points[j][Col]; for (i = 0; i < RVLength; i++, ROnePts += RULength) CAGD_ALPHA_BLEND_STEP(A, i, OnePts, Srf -> VLength, ROnePts, ULength); } } BspKnotFreeAlphaCoef(A); break; default: CAGD_FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV); break; } } BspKnotMakeRobustKV(RefSrf -> UKnotVector, RefSrf -> UOrder + RefSrf -> ULength); BspKnotMakeRobustKV(RefSrf -> VKnotVector, RefSrf -> VOrder + RefSrf -> VLength); if (NewSrf) CagdSrfFree(Srf); return RefSrf; } /***************************************************************************** * DESCRIPTION: M * Returns a new Bspline surface, identical to the original but with one M * degree higher, in the requested direction Dir. M * * * PARAMETERS: M * Srf: To raise it degree by one. M * Dir: Direction of degree raising. Either U or V. M * * * RETURN VALUE: M * CagdSrfStruct *: A surface with one degree higher in direction Dir, M * representing the same geometry as Srf. M * * * SEE ALSO: M * CagdSrfDegreeRaise, BzrSrfDegreeRaise, TrimSrfDegreeRaise M * * * KEYWORDS: M * BspSrfDegreeRaise, degree raising M *****************************************************************************/ CagdSrfStruct *BspSrfDegreeRaise(CagdSrfStruct *Srf, CagdSrfDirType Dir) { CagdBType NewSrf = FALSE, IsNotRational = !CAGD_IS_RATIONAL_SRF(Srf); int i, i2, j, RaisedLen, Row, Col, Length, Order = Dir == CAGD_CONST_V_DIR ? Srf -> UOrder : Srf -> VOrder, MaxCoord = CAGD_NUM_OF_PT_COORD(Srf -> PType); CagdSrfStruct *RaisedSrf = NULL; if (CAGD_IS_PERIODIC_SRF(Srf)) { NewSrf = TRUE; Srf = CnvrtPeriodic2FloatSrf(Srf); } Length = Dir == CAGD_CONST_V_DIR ? Srf -> ULength : Srf -> VLength; if (Order > 2) { CagdSrfStruct *UnitSrf; int UKvLen1 = Srf -> UOrder + Srf -> ULength - 1, VKvLen1 = Srf -> VOrder + Srf -> VLength - 1; CagdRType *UKv = Srf -> UKnotVector, *VKv = Srf -> VKnotVector; /* Degree raise by multiplying by a constant 1 linear surface in the */ /* raised direction and constant 1 constant surface in the other. */ switch (Dir) { case CAGD_CONST_U_DIR: UnitSrf = BspSrfNew(1, 2, 1, 2, CAGD_MAKE_PT_TYPE(FALSE, MaxCoord)); for (i = 0; i < 2; i++) UnitSrf -> UKnotVector[i] = i > 0 ? UKv[UKvLen1] : UKv[0]; for (i = 0; i < 4; i++) UnitSrf -> VKnotVector[i] = i > 1 ? VKv[VKvLen1] : VKv[0]; break; case CAGD_CONST_V_DIR: UnitSrf = BspSrfNew(2, 1, 2, 1, CAGD_MAKE_PT_TYPE(FALSE, MaxCoord)); for (i = 0; i < 4; i++) UnitSrf -> UKnotVector[i] = i > 1 ? UKv[UKvLen1] : UKv[0]; for (i = 0; i < 2; i++) UnitSrf -> VKnotVector[i] = i > 0 ? VKv[VKvLen1] : VKv[0]; break; default: UnitSrf = NULL; CAGD_FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV); break; } for (i = 1; i <= MaxCoord; i++) UnitSrf -> Points[i][0] = UnitSrf -> Points[i][1] = 1.0; RaisedSrf = BspSrfMult(Srf, UnitSrf); CagdSrfFree(UnitSrf); if (NewSrf) CagdSrfFree(Srf); return RaisedSrf; } /* If surface 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; switch (Dir) { case CAGD_CONST_U_DIR: RaisedSrf = BspSrfNew(Srf -> ULength, RaisedLen, Srf -> UOrder, Order + 1, Srf -> PType); /* Update the knot vectors. */ CAGD_GEN_COPY(RaisedSrf -> UKnotVector, Srf -> UKnotVector, sizeof(CagdRType) * (Srf -> ULength + Srf -> UOrder)); for (i = 0; i < 3; i++) RaisedSrf -> VKnotVector[i] = Srf -> VKnotVector[0]; for (i = 2, j = 3; i < Length; i++, j += 2) RaisedSrf -> VKnotVector[j] = RaisedSrf -> VKnotVector[j + 1] = Srf -> VKnotVector[i]; for (i = j; i < j + 3; i++) RaisedSrf -> VKnotVector[i] = Srf -> VKnotVector[Length]; /* Update the mesh. */ for (Col = 0; Col < Srf -> ULength; Col++) { for (j = IsNotRational; j <= MaxCoord; j++) /* First point. */ RaisedSrf -> Points[j][RAISED_SRF(Col, 0)] = Srf -> Points[j][SRF(Col, 0)]; for (i = 1, i2 = 1; i < Length; i++, i2 += 2) for (j = IsNotRational; j <= MaxCoord; j++) { RaisedSrf -> Points[j][RAISED_SRF(Col, i2)] = Srf -> Points[j][SRF(Col, i - 1)] * 0.5 + Srf -> Points[j][SRF(Col, i)] * 0.5; RaisedSrf -> Points[j][RAISED_SRF(Col, i2 + 1)] = Srf -> Points[j][SRF(Col, i)]; } } break; case CAGD_CONST_V_DIR: RaisedSrf = BspSrfNew(RaisedLen, Srf -> VLength, Order + 1, Srf -> VOrder, Srf -> PType); /* Update the knot vectors. */ CAGD_GEN_COPY(RaisedSrf -> VKnotVector, Srf -> VKnotVector, sizeof(CagdRType) * (Srf -> VLength + Srf -> VOrder)); for (i = 0; i < 3; i++) RaisedSrf -> UKnotVector[i] = Srf -> UKnotVector[0]; for (i = 2, j = 3; i < Length; i++, j += 2) RaisedSrf -> UKnotVector[j] = RaisedSrf -> UKnotVector[j + 1] = Srf -> UKnotVector[i]; for (i = j; i < j + 3; i++) RaisedSrf -> UKnotVector[i] = Srf -> UKnotVector[Length]; /* Update the mesh. */ for (Row = 0; Row < Srf -> VLength; Row++) { for (j = IsNotRational; j <= MaxCoord; j++) /* First point. */ RaisedSrf -> Points[j][RAISED_SRF(0, Row)] = Srf -> Points[j][SRF(0, Row)]; for (i = 1, i2 = 1; i < Length; i++, i2 += 2) for (j = IsNotRational; j <= MaxCoord; j++) { RaisedSrf -> Points[j][RAISED_SRF(i2, Row)] = Srf -> Points[j][SRF(i - 1, Row)] * 0.5 + Srf -> Points[j][SRF(i, Row)] * 0.5; RaisedSrf -> Points[j][RAISED_SRF(i2 + 1, Row)] = Srf -> Points[j][SRF(i, Row)]; } } break; default: CAGD_FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV); break; } if (NewSrf) CagdSrfFree(Srf); return RaisedSrf; } /***************************************************************************** * DESCRIPTION: M * Returns a new surface equal to the given surface, differentiated once in M * the direction Dir. 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 * This is applied to all rows/cols of the surface. M * * * PARAMETERS: M * Srf: To differentiate. M * Dir: Direction of differentiation. Either U or V. M * * * RETURN VALUE: M * CagdSrfStruct *: Differentiated surface. M * * * SEE ALSO: M * CagdSrfDerive, BzrSrfDerive, BzrSrfDeriveRational, BspSrfDeriveRational M * * * KEYWORDS: M * BspSrfDerive, derivatives, partial derivatives M *****************************************************************************/ CagdSrfStruct *BspSrfDerive(CagdSrfStruct *Srf, CagdSrfDirType Dir) { CagdBType NewSrf = FALSE, IsNotRational = !CAGD_IS_RATIONAL_SRF(Srf); int i, j, Row, Col, NewUOrder, NewVOrder, NewULength, NewVLength, ULength, VLength, UOrder = Srf -> UOrder, VOrder = Srf -> VOrder, MaxCoord = CAGD_NUM_OF_PT_COORD(Srf -> PType); CagdRType **DPoints, *UKv, *VKv, **Points; CagdSrfStruct *DerivedSrf = NULL; if (CAGD_IS_PERIODIC_SRF(Srf)) { NewSrf = TRUE; Srf = CnvrtPeriodic2FloatSrf(Srf); } if (!IsNotRational) return BspSrfDeriveRational(Srf, Dir); ULength = Srf -> ULength; VLength = Srf -> VLength; UKv = Srf -> UKnotVector; VKv = Srf -> VKnotVector; Points = Srf -> Points; switch (Dir) { case CAGD_CONST_U_DIR: NewULength = UOrder < 2 ? ULength : ULength - 1; NewUOrder = MAX(UOrder - 1, 1); DerivedSrf = BspSrfNew(NewULength, VLength, NewUOrder, VOrder, Srf -> PType); CAGD_GEN_COPY(DerivedSrf -> UKnotVector, &UKv[UOrder < 2 ? 0 : 1], sizeof(CagdRType) * (NewULength + NewUOrder)); CAGD_GEN_COPY(DerivedSrf -> VKnotVector, VKv, sizeof(CagdRType) * (VLength + VOrder)); DPoints = DerivedSrf -> Points; for (Row = 0; Row < VLength; Row++) for (i = 0; i < NewULength; i++) { CagdRType Denom = UKv[i + UOrder] - UKv[i + 1]; for (j = IsNotRational; j <= MaxCoord; j++) { if (APX_EQ(Denom, 0.0)) Denom = IRIT_INFNTY; DPoints[j][DERIVED_SRF(i, Row)] = UOrder < 2 ? 0.0 : (UOrder - 1) * (Points[j][SRF(i + 1, Row)] - Points[j][SRF(i, Row)]) / Denom; } } break; case CAGD_CONST_V_DIR: NewVLength = VOrder < 2 ? VLength : VLength - 1; NewVOrder = MAX(VOrder - 1, 1); DerivedSrf = BspSrfNew(ULength, NewVLength, UOrder, NewVOrder, Srf -> PType); CAGD_GEN_COPY(DerivedSrf -> UKnotVector, UKv, sizeof(CagdRType) * (ULength + UOrder)); CAGD_GEN_COPY(DerivedSrf -> VKnotVector, &VKv[VOrder < 2 ? 0 : 1], sizeof(CagdRType) * (NewVLength + NewVOrder)); DPoints = DerivedSrf -> Points; for (Col = 0; Col < ULength; Col++) for (i = 0; i < NewVLength; i++) { CagdRType Denom = VKv[i + VOrder] - VKv[i + 1]; for (j = IsNotRational; j <= MaxCoord; j++) { if (APX_EQ(Denom, 0.0)) Denom = IRIT_INFNTY; DPoints[j][DERIVED_SRF(Col, i)] = VOrder < 2 ? 0.0 : (VOrder - 1) * (Points[j][SRF(Col, i + 1)] - Points[j][SRF(Col, i)]) / Denom; } } break; default: CAGD_FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV); break; } if (NewSrf) CagdSrfFree(Srf); return DerivedSrf; } /***************************************************************************** * DESCRIPTION: M * Evaluates the (unit) tangent to a surface at a given parametric location M * (u, v) and given direction Dir. M * * * PARAMETERS: M * Srf: Bspline surface to evaluate (unit) tangent vector for. M * u, v: Parametric location of required (unit) tangent. M * Dir: Direction of tangent vector. Either U or V. 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 (unit) M * tangent information. M * * * SEE ALSO: M * CagdSrfTangent, BzrSrfTangent M * * * KEYWORDS: M * BspSrfTangent, tangent M *****************************************************************************/ CagdVecStruct *BspSrfTangent(CagdSrfStruct *Srf, CagdRType u, CagdRType v, CagdSrfDirType Dir, CagdBType Normalize) { CagdVecStruct *Tangent = NULL; CagdCrvStruct *Crv; switch (Dir) { case CAGD_CONST_V_DIR: Crv = BspSrfCrvFromSrf(Srf, v, Dir); Tangent = BspCrvTangent(Crv, u, Normalize); CagdCrvFree(Crv); break; case CAGD_CONST_U_DIR: Crv = BspSrfCrvFromSrf(Srf, u, Dir); Tangent = BspCrvTangent(Crv, v, Normalize); CagdCrvFree(Crv); break; default: CAGD_FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV); break; } return Tangent; } /***************************************************************************** * DESCRIPTION: M * Evaluate the (unit) normal of a surface at a given parametric location. M * If we fail to compute the normal at given location we retry by moving a M * tad. M * * * PARAMETERS: M * Srf: Bspline surface to evaluate (unit) normal vector for. M * u, v: Parametric location of required (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 (unit) normal M * information. M * * * SEE ALSO: M * CagdSrfNormal, BzrSrfNormal, BspSrfMeshNormals, SymbSrfNormalSrf M * * * KEYWORDS: M * BspSrfNormal, normal M *****************************************************************************/ CagdVecStruct *BspSrfNormal(CagdSrfStruct *Srf, CagdRType u, CagdRType v, CagdBType Normalize) { static CagdVecStruct Normal; CagdVecStruct *V, T1, T2; CagdRType UMin, UMax, VMin, VMax; CAGD_DOMAIN_GET_AND_VERIFY_SRF(u, v, Srf, UMin, UMax, VMin, VMax); V = BspSrfTangent(Srf, u, v, CAGD_CONST_U_DIR, FALSE); if (CAGD_LEN_VECTOR(*V) < IRIT_UEPS) V = BspSrfTangent(Srf, u > UMin + IRIT_EPS ? u - IRIT_EPS : u + IRIT_EPS, v > VMin + IRIT_EPS ? v - IRIT_EPS : v + IRIT_EPS, CAGD_CONST_U_DIR, FALSE); CAGD_COPY_VECTOR(T1, *V); V = BspSrfTangent(Srf, u, v, CAGD_CONST_V_DIR, FALSE); if (CAGD_LEN_VECTOR(*V) < IRIT_UEPS) V = BspSrfTangent(Srf, u > UMin + IRIT_EPS ? u - IRIT_EPS : u + IRIT_EPS, v > VMin + IRIT_EPS ? v - IRIT_EPS : v + IRIT_EPS, CAGD_CONST_V_DIR, FALSE); CAGD_COPY_VECTOR(T2, *V); /* The normal is the cross product of T1 and T2: */ Normal.Vec[0] = T1.Vec[1] * T2.Vec[2] - T1.Vec[2] * T2.Vec[1]; Normal.Vec[1] = T1.Vec[2] * T2.Vec[0] - T1.Vec[0] * T2.Vec[2]; Normal.Vec[2] = T1.Vec[0] * T2.Vec[1] - T1.Vec[1] * T2.Vec[0]; if (Normalize) CAGD_NORMALIZE_VECTOR(Normal); /* Normalize the vector. */ return &Normal; } /***************************************************************************** * DESCRIPTION: M * Evaluates the unit normals of a surface at a mesh defined by subdividing M * the parametric space into a grid of size UFineNess by VFineNess. M * The normals are saved in a linear CagdVecStruct vector which is M * allocated dynamically. Data is saved u inc. first. M * This routine is much faster than evaluating normal for each point, M * individually. M * * * PARAMETERS: M * Srf: To compute normals on a grid of its parametric domain. M * UFineNess: U Fineness of imposed grid on Srf's parametric domain. M * VFineNess: V Fineness of imposed grid on Srf's parametric domain. M * * * RETURN VALUE: M * CagdVecStruct *: An vector of unit normals (u increments first). M * * * SEE ALSO: M * CagdSrfNormal, BspSrfNormal, SymbSrfNormalSrf M * * * KEYWORDS: M * BspSrfMeshNormals, normal M *****************************************************************************/ CagdVecStruct *BspSrfMeshNormals(CagdSrfStruct *Srf, int UFineNess, int VFineNess) { int i, j; CagdRType UMin, UMax, VMin, VMax; CagdVecStruct *Normals, *NPtr; CagdSrfStruct *DuSrf = CagdSrfDerive(Srf, CAGD_CONST_U_DIR), *DvSrf = CagdSrfDerive(Srf, CAGD_CONST_V_DIR); UFineNess = MAX(2, UFineNess); VFineNess = MAX(2, VFineNess); Normals = CagdVecArrayNew(UFineNess * VFineNess); CagdSrfDomain(Srf, &UMin, &UMax, &VMin, &VMax); NPtr = Normals; for (i = 0; i < UFineNess; i++) { CagdRType U = UMin + (UMax - UMin) * i / (UFineNess - 1); CagdCrvStruct *DuCrv, *DvCrv, *DuCrv2 = NULL, *DvCrv2 = NULL; if (U > UMax) /* Due to floating point round off. */ U = UMax; DuCrv = CagdCrvFromSrf(DuSrf, U, CAGD_CONST_U_DIR), DvCrv = CagdCrvFromSrf(DvSrf, U, CAGD_CONST_U_DIR); for (j = 0; j < VFineNess; j++) { CagdVType Du, Dv; CagdRType V = VMin + (VMax - VMin) * j / (VFineNess - 1); if (V > VMax) /* Due to floating point round off. */ V = VMax; if (!EvalTangentVector(Du, DuCrv, V, VMin, VMax)) { if (DuCrv2 == NULL) { DuCrv2 = CagdCrvFromSrf(DuSrf, U > UMin + NORMAL_IRIT_EPS ? U - NORMAL_IRIT_EPS : U + NORMAL_IRIT_EPS, CAGD_CONST_U_DIR); } EvalTangentVector(Du, DuCrv2, V, VMin, VMax); } if (!EvalTangentVector(Dv, DvCrv, V, VMin, VMax)) { if (DvCrv2 == NULL) { DvCrv2 = CagdCrvFromSrf(DvSrf, U > UMin + NORMAL_IRIT_EPS ? U - NORMAL_IRIT_EPS : U + NORMAL_IRIT_EPS, CAGD_CONST_U_DIR); } EvalTangentVector(Dv, DvCrv2, V, VMin, VMax); } CROSS_PROD(NPtr -> Vec, Dv, Du); NPtr++; } if (DuCrv2 != NULL) CagdCrvFree(DuCrv2); if (DvCrv2 != NULL) CagdCrvFree(DvCrv2); CagdCrvFree(DuCrv); CagdCrvFree(DvCrv); } /* Normalize the results. */ NPtr = Normals; for (i = 0; i < UFineNess; i++) { for (j = 0; j < VFineNess; j++) { CagdRType Len = PT_LENGTH(NPtr -> Vec); if (Len < PT_NORMALIZE_ZERO) { CagdVecStruct *NPtrTmp; /* Try its neighbors. */ if (i > 0) { NPtrTmp = NPtr - VFineNess; PT_COPY(NPtr -> Vec, NPtrTmp -> Vec); Len = PT_LENGTH(NPtr -> Vec); } if (Len < PT_NORMALIZE_ZERO && i < UFineNess - 1) { NPtrTmp = NPtr + VFineNess; PT_COPY(NPtr -> Vec, NPtrTmp -> Vec); Len = PT_LENGTH(NPtr -> Vec); } if (Len < PT_NORMALIZE_ZERO && j > 0) { NPtrTmp = NPtr - 1; PT_COPY(NPtr -> Vec, NPtrTmp -> Vec); Len = PT_LENGTH(NPtr -> Vec); } if (Len < PT_NORMALIZE_ZERO && j < VFineNess - 1) { NPtrTmp = NPtr + 1; PT_COPY(NPtr -> Vec, NPtrTmp -> Vec); Len = PT_LENGTH(NPtr -> Vec); } if (Len > PT_NORMALIZE_ZERO) { Len = 1 / Len; PT_SCALE(NPtr -> Vec, Len); } else { /* Do something. */ NPtr -> Vec[0] = NPtr -> Vec[1] = 0.0; NPtr -> Vec[2] = 1.0; } } else { Len = 1 / Len; PT_SCALE(NPtr -> Vec, Len); } NPtr++; } } CagdSrfFree(DuSrf); CagdSrfFree(DvSrf); return Normals; } /***************************************************************************** * DESCRIPTION: * * Evaluates a vector field Crv at param t. If zero length, tries to move a * * tad. * * * * PARAMETERS: * * Vec: Where to place the result. * * Crv: Vector field to evaluate. * * t: Parameter value where to evaluate the vector field. * * TMin: Minimum of vector field domain. * * TMax: Maximum of vector field domain. * * * * RETURN VALUE: * * CagdBType: TRUE, if successful, FALSE otherwise. * *****************************************************************************/ static CagdBType EvalTangentVector(CagdVType Vec, CagdCrvStruct *Crv, CagdRType t, CagdRType TMin, CagdRType TMax) { CagdRType Len, *R = CagdCrvEval(Crv, t); CagdCoerceToE3(Vec, &R, -1, Crv -> PType); Len = PT_LENGTH(Vec); if (Len < PT_NORMALIZE_ZERO && t - NORMAL_IRIT_EPS > TMin) { CagdCrvEval(Crv, t - NORMAL_IRIT_EPS); CagdCoerceToE3(Vec, &R, -1, Crv -> PType); Len = PT_LENGTH(Vec); } if (Len < PT_NORMALIZE_ZERO && t + NORMAL_IRIT_EPS < TMax) { CagdCrvEval(Crv, t + NORMAL_IRIT_EPS); CagdCoerceToE3(Vec, &R, -1, Crv -> PType); Len = PT_LENGTH(Vec); } if (Len > PT_NORMALIZE_ZERO) { Len = 1.0 / Len; PT_SCALE(Vec, Len); return TRUE; } else return FALSE; } /***************************************************************************** * DESCRIPTION: M * Evaluates the unit normals of a surface at a mesh defined by subdividing M * the parametric space into a grid of size UFineNess by VFineNess. M * The normals are saved in a linear CagdVecStruct vector which is M * allocated dynamically. Data is saved u inc. first. M * This routine is much faster than evaluating normal for each point, M * individually. M * * * PARAMETERS: M * Srf: To compute normals on a grid of its parametric domain. M * UFineNess: U Fineness of imposed grid on Srf's parametric domain. M * VFineNess: V Fineness of imposed grid on Srf's parametric domain. M * * * RETURN VALUE: M * CagdVecStruct *: An vector of unit normals (u increments first). M * * * KEYWORDS: M * BspSrfMeshNormalsSymb, normal M *****************************************************************************/ CagdVecStruct *BspSrfMeshNormalsSymb(CagdSrfStruct *Srf, int UFineNess, int VFineNess) { int i, j; CagdRType U, V, UMin, UMax, VMin, VMax, **Points; CagdVecStruct *Normals, *NPtr; CagdSrfStruct *TstNormalSrf, *NormalSrf = SymbSrfNormalSrf(Srf); /* Verify we have a valid regular surface. */ TstNormalSrf = CagdCoerceSrfTo(NormalSrf, CAGD_PT_E3_TYPE); for (Points = TstNormalSrf -> Points, i = 0; i < TstNormalSrf -> ULength * TstNormalSrf -> VLength; i++) if (!APX_EQ(Points[1][i], 0.0 ) || !APX_EQ(Points[2][i], 0.0 ) || !APX_EQ(Points[3][i], 0.0 )) break; CagdSrfFree(TstNormalSrf); if (i >= TstNormalSrf -> ULength * TstNormalSrf -> VLength) { /* Not a regular surface - ignore it, */ return NULL; } UFineNess = MAX(2, UFineNess); VFineNess = MAX(2, VFineNess); Normals = CagdVecArrayNew(UFineNess * VFineNess); CagdSrfDomain(Srf, &UMin, &UMax, &VMin, &VMax); NPtr = Normals; for (i = 0; i < UFineNess; i++) { for (j = 0; j < VFineNess; j++) { U = UMin + (UMax - UMin) * i / (UFineNess - 1); V = VMin + (VMax - VMin) * j / (VFineNess - 1); CagdEvaluateSurfaceVecField(NPtr -> Vec, NormalSrf, U, V); NPtr++; } } CagdSrfFree(NormalSrf); return Normals; } /***************************************************************************** * DESCRIPTION: M * Converts a Bspline surface into a Bspline surface with floating end M * conditions. M * * * PARAMETERS: M * Srf: Bspline surface to convert to floating end conditions. Assume M * Crv is either periodic or has floating end condition. M * * * RETURN VALUE: M * CagdSrfStruct *: A Bspline surface with floating end conditions, M * representing the same geometry as Srf. M * * * SEE ALSO: M * CnvrtFloat2OpenSrf M * * * KEYWORDS: M * CnvrtPeriodic2FloatSrf, conversion M *****************************************************************************/ CagdSrfStruct *CnvrtPeriodic2FloatSrf(CagdSrfStruct *Srf) { int i, j, UOrder = Srf -> UOrder, VOrder = Srf -> VOrder, ULength = Srf -> ULength, VLength = Srf -> VLength, MaxAxis = CAGD_NUM_OF_PT_COORD(Srf -> PType); CagdSrfStruct *NewSrf; if (!CAGD_IS_UPERIODIC_SRF(Srf) && !CAGD_IS_VPERIODIC_SRF(Srf)) { CAGD_FATAL_ERROR(CAGD_ERR_PERIODIC_EXPECTED); return NULL; } NewSrf = BspSrfNew(CAGD_SRF_UPT_LST_LEN(Srf), CAGD_SRF_VPT_LST_LEN(Srf), UOrder, VOrder, Srf -> PType); NewSrf -> UKnotVector = BspKnotCopy(Srf -> UKnotVector, CAGD_SRF_UPT_LST_LEN(Srf) + UOrder); NewSrf -> VKnotVector = BspKnotCopy(Srf -> VKnotVector, CAGD_SRF_VPT_LST_LEN(Srf) + UOrder); for (i = !CAGD_IS_RATIONAL_PT(Srf -> PType); i <= MaxAxis; i++) { CagdRType *NewPts, *Pts = Srf -> Points[i]; NewSrf -> Points[i] = (CagdRType *) IritMalloc(sizeof(CagdRType) * CAGD_SRF_UPT_LST_LEN(Srf) * CAGD_SRF_VPT_LST_LEN(Srf)); NewPts = NewSrf -> Points[i]; for (j = 0; j < VLength; j++, Pts += ULength) { CAGD_GEN_COPY(NewPts, Pts, sizeof(CagdRType) * ULength); NewPts += ULength; if (Srf -> UPeriodic) { CAGD_GEN_COPY(NewPts, Pts, sizeof(CagdRType) * (UOrder - 1)); NewPts += UOrder - 1; } } if (Srf -> VPeriodic) { CAGD_GEN_COPY(NewPts, NewSrf -> Points[i], sizeof(CagdRType) * (VOrder - 1) * CAGD_SRF_UPT_LST_LEN(Srf)); } } for (i = MaxAxis + 1; i <= CAGD_MAX_PT_COORD; i++) NewSrf -> Points[i] = NULL; return NewSrf; } /***************************************************************************** * DESCRIPTION: M * Converts a Bspline surface to a Bspline surface with open end conditions. M * * * PARAMETERS: M * Srf: Bspline surface to convert to open end conditions. M * * * RETURN VALUE: M * CagdSrfStruct *: A Bspline surface with open end conditions, M * representing the same geometry as Srf. M * * * SEE ALSO: M * CnvrtPeriodic2FloatSrf M * * * KEYWORDS: M * CnvrtFloat2OpenSrf, conversion M *****************************************************************************/ CagdSrfStruct *CnvrtFloat2OpenSrf(CagdSrfStruct *Srf) { CagdRType UMin, UMax, VMin, VMax; CagdSrfStruct *TSrf1, *TSrf2; if (!CAGD_IS_BSPLINE_SRF(Srf)) { CAGD_FATAL_ERROR(CAGD_ERR_BSP_SRF_EXPECT); return NULL; } CagdSrfDomain(Srf, &UMin, &UMax, &VMin, &VMax); TSrf1 = CagdSrfRegionFromSrf(Srf, UMin, UMax, CAGD_CONST_U_DIR); TSrf2 = CagdSrfRegionFromSrf(TSrf1, VMin, VMax, CAGD_CONST_V_DIR); CagdSrfFree(TSrf1); return TSrf2; }