/****************************************************************************** * SBsp-Aux.c - Bspline surface auxilary routines. * ******************************************************************************* * 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) /****************************************************************************** * Given a bspline surface - subdivide it into two at given parametric value. * * Returns pointer to first surface in a list of two srfs (subdivided ones). * * The subdivision is exact result of evaluating the surface int a curve at t * * using the recursive algorithm - the left resulting points is left surface, * * and the right resulting points is right surface (left is below t). * ******************************************************************************/ CagdSrfStruct *BspSrfSubdivAtParam(CagdSrfStruct *Srf, CagdRType t, CagdSrfDirType Dir) { CagdBType IsNotRational = !CAGD_IS_RATIONAL_CRV(Srf); int i, j, Row, Col, KVLen, Index1, Index2, Mult, LULength, RULength, LVLength, RVLength, ULength = Srf -> ULength, VLength = Srf -> VLength, UOrder = Srf -> UOrder, VOrder = Srf -> VOrder, MaxCoord = CAGD_NUM_OF_PT_COORD(Srf -> PType); CagdRType *RefKV, **Pts, **LPts, **RPts, *LOnePts, *ROnePts, *OnePts; CagdSrfStruct *RSrf, *LSrf; BspKnotAlphaCoeffType *A; switch (Dir) { case CAGD_CONST_U_DIR: RefKV = Srf -> UKnotVector; KVLen = UOrder + ULength; Index1 = BspKnotLastIndexL(RefKV, KVLen, t); Index2 = BspKnotFirstIndexG(RefKV, KVLen, t); 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); Index2 = BspKnotFirstIndexG(RefKV, KVLen, t); 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; 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; switch (Dir) { case CAGD_CONST_U_DIR: /* Compute the Alpha refinement matrix. */ if (Mult > 0) { CagdRType *NewKV = (CagdRType *) IritMalloc(sizeof(CagdRType) * Mult); 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, 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, ROnePts); } } BspKnotFreeAlphaCoef(A); break; case CAGD_CONST_V_DIR: /* Compute the Alpha refinement matrix. */ if (Mult > 0) { CagdRType *NewKV = (CagdRType *) IritMalloc(sizeof(CagdRType) * Mult); 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, 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, ROnePts, ULength); } } BspKnotFreeAlphaCoef(A); break; default: 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; return LSrf; } /****************************************************************************** * Insert n knot all with the value t in direction Dir. In no case will the * * multiplicity of knot be greater or equal to the curve order. * ******************************************************************************/ CagdSrfStruct *BspSrfKnotInsertNSame(CagdSrfStruct *BspSrf, CagdSrfDirType Dir, CagdRType t, int n) { int i, CrntMult, Mult; CagdSrfStruct *RefinedSrf; switch (Dir) { case CAGD_CONST_U_DIR: CrntMult = BspKnotFindMult(BspSrf -> UKnotVector, BspSrf -> UOrder, BspSrf -> ULength, t), Mult = MIN(n, BspSrf -> UOrder - CrntMult - 1); break; case CAGD_CONST_V_DIR: CrntMult = BspKnotFindMult(BspSrf -> VKnotVector, BspSrf -> VOrder, BspSrf -> VLength, t), Mult = MIN(n, BspSrf -> VOrder - CrntMult - 1); break; default: Mult = 0; 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(BspSrf, Dir, FALSE, NewKV, Mult); IritFree((VoidPtr) NewKV); } else { RefinedSrf = CagdSrfCopy(BspSrf); } return RefinedSrf; } /****************************************************************************** * Insert n knot with different values as defined by t. If however Replace is * * TRUE, the knot are simply replacing the current ones. * ******************************************************************************/ CagdSrfStruct *BspSrfKnotInsertNDiff(CagdSrfStruct *Srf, CagdSrfDirType Dir, int Replace, CagdRType *t, int n) { CagdBType IsNotRational = !CAGD_IS_RATIONAL_CRV(Srf); int i, Row, Col, ULength = Srf -> ULength, VLength = Srf -> VLength, UOrder = Srf -> UOrder, VOrder = Srf -> VOrder, MaxCoord = CAGD_NUM_OF_PT_COORD(Srf -> PType); CagdSrfStruct *RefSrf = NULL; if (Replace) { for (i = 1; i < n; i++) if (t[i] < t[i - 1]) FATAL_ERROR(CAGD_ERR_KNOT_NOT_ORDERED); switch (Dir) { case CAGD_CONST_U_DIR: if (Srf -> UOrder + Srf -> ULength != n) 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) FATAL_ERROR(CAGD_ERR_NUM_KNOT_MISMATCH); RefSrf = CagdSrfCopy(Srf); for (i = 0; i < n; i++) RefSrf -> VKnotVector[i] = *t++; break; default: 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, *UKnotVector = Srf -> UKnotVector, *VKnotVector = Srf -> VKnotVector; switch (Dir) { case CAGD_CONST_U_DIR: /* 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, ROnePts ); } } BspKnotFreeAlphaCoef(A); break; case CAGD_CONST_V_DIR: /* 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, ROnePts, ULength); } } BspKnotFreeAlphaCoef(A); break; default: FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV); break; } } BspKnotMakeRobustKV(RefSrf -> UKnotVector, RefSrf -> UOrder + RefSrf -> ULength); BspKnotMakeRobustKV(RefSrf -> VKnotVector, RefSrf -> VOrder + RefSrf -> VLength); return RefSrf; } /****************************************************************************** * Return a new surface, identical to the original but with one degree higher * * in the given direction. * ******************************************************************************/ CagdSrfStruct *BspSrfDegreeRaise(CagdSrfStruct *Srf, CagdSrfDirType Dir) { CagdBType IsNotRational = !CAGD_IS_RATIONAL_SRF(Srf); int i, i2, j, RaisedLen, Row, Col, Order = Dir == CAGD_CONST_V_DIR ? Srf -> UOrder : Srf -> VOrder, Length = Dir == CAGD_CONST_V_DIR ? Srf -> ULength : Srf -> VLength, MaxCoord = CAGD_NUM_OF_PT_COORD(Srf -> PType); CagdSrfStruct *RaisedSrf = NULL; 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; 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); 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: FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV); break; } return RaisedSrf; } /****************************************************************************** * Return a new surface equal to the derived surface in the direction Dir. * * Let old control polygon be P(i), i = 0 to k-1, and Q(i) be new one then: * * Q(i) = (k - 1) * P(i+1) - P(i), i = 0 to k-2. * * This is applied to all rows/cols of the surface. * ******************************************************************************/ CagdSrfStruct *BspSrfDerive(CagdSrfStruct *Srf, CagdSrfDirType Dir) { CagdBType IsNotRational = !CAGD_IS_RATIONAL_SRF(Srf); int i, j, Row, Col, ULength = Srf -> ULength, VLength = Srf -> VLength, UOrder = Srf -> UOrder, VOrder = Srf -> VOrder, MaxCoord = CAGD_NUM_OF_PT_COORD(Srf -> PType); CagdRType **DPoints, *UKv = Srf -> UKnotVector, *VKv = Srf -> VKnotVector, **Points = Srf -> Points; CagdSrfStruct *DerivedSrf = NULL; if (!IsNotRational) return BspSrfDeriveRational(Srf, Dir); switch (Dir) { case CAGD_CONST_V_DIR: if (UOrder < 2) FATAL_ERROR(CAGD_ERR_LIN_NO_SUPPORT); DerivedSrf = BspSrfNew(ULength - 1, VLength, UOrder - 1, VOrder, Srf -> PType); CAGD_GEN_COPY(DerivedSrf -> UKnotVector, &UKv[1], sizeof(CagdRType) * (ULength + UOrder - 2)); CAGD_GEN_COPY(DerivedSrf -> VKnotVector, VKv, sizeof(CagdRType) * (VLength + VOrder)); DPoints = DerivedSrf -> Points; for (Row = 0; Row < VLength; Row++) for (i = 0; i < ULength - 1; i++) { CagdRType Denom = UKv[i + UOrder] - UKv[i + 1]; for (j = IsNotRational; j <= MaxCoord; j++) { if (APX_EQ(Denom, 0.0)) Denom = INFINITY; DPoints[j][DERIVED_SRF(i, Row)] = (UOrder - 1) * (Points[j][SRF(i + 1, Row)] - Points[j][SRF(i, Row)]) / Denom; } } break; case CAGD_CONST_U_DIR: if (VOrder < 2) FATAL_ERROR(CAGD_ERR_LIN_NO_SUPPORT); DerivedSrf = BspSrfNew(ULength, VLength - 1, UOrder, VOrder - 1, Srf -> PType); CAGD_GEN_COPY(DerivedSrf -> UKnotVector, UKv, sizeof(CagdRType) * (ULength + UOrder)); CAGD_GEN_COPY(DerivedSrf -> VKnotVector, &VKv[1], sizeof(CagdRType) * (VLength + VOrder - 2)); DPoints = DerivedSrf -> Points; for (Col = 0; Col < ULength; Col++) for (i = 0; i < VLength - 1; i++) { CagdRType Denom = VKv[i + VOrder] - VKv[i + 1]; for (j = IsNotRational; j <= MaxCoord; j++) { if (APX_EQ(Denom, 0.0)) Denom = INFINITY; DPoints[j][DERIVED_SRF(Col, i)] = (VOrder - 1) * (Points[j][SRF(Col, i + 1)] - Points[j][SRF(Col, i)]) / Denom; } } break; default: FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV); break; } return DerivedSrf; } /****************************************************************************** * Evaluate the tangent to a surface at a given point and given direction. * ******************************************************************************/ CagdVecStruct *BspSrfTangent(CagdSrfStruct *Srf, CagdRType u, CagdRType v, CagdSrfDirType Dir) { CagdVecStruct *Tangent = NULL; CagdCrvStruct *Crv; switch (Dir) { case CAGD_CONST_V_DIR: Crv = BspSrfCrvFromSrf(Srf, v, Dir); Tangent = BspCrvTangent(Crv, u); CagdCrvFree(Crv); break; case CAGD_CONST_U_DIR: Crv = BspSrfCrvFromSrf(Srf, u, Dir); Tangent = BspCrvTangent(Crv, v); CagdCrvFree(Crv); break; default: FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV); break; } return Tangent; } /****************************************************************************** * Evaluate the normal of a surface at a given point. * * If we fail to compute the normal at given location we try by moving a tad. * ******************************************************************************/ CagdVecStruct *BspSrfNormal(CagdSrfStruct *Srf, CagdRType u, CagdRType v) { static CagdVecStruct Normal; CagdVecStruct *V, T1, T2; CagdRType UMin, UMax, VMin, VMax; CagdSrfDomain(Srf, &UMin, &UMax, &VMin, &VMax); V = BspSrfTangent(Srf, u, v, CAGD_CONST_U_DIR); if (CAGD_LEN_VECTOR(*V) < EPSILON) V = BspSrfTangent(Srf, u > UMin + EPSILON ? u - EPSILON : u + EPSILON, v > VMin + EPSILON ? v - EPSILON : v + EPSILON, CAGD_CONST_U_DIR); CAGD_COPY_VECTOR(T1, *V); V = BspSrfTangent(Srf, u, v, CAGD_CONST_V_DIR); if (CAGD_LEN_VECTOR(*V) < EPSILON) V = BspSrfTangent(Srf, u > UMin + EPSILON ? u - EPSILON : u + EPSILON, v > VMin + EPSILON ? v - EPSILON : v + EPSILON, CAGD_CONST_V_DIR); 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]; CAGD_NORMALIZE_VECTOR(Normal); /* Normalize the vector. */ return &Normal; } /****************************************************************************** * Evaluate the normals of a surface at a mesh defined by subdividing the * * parametric space into a grid of size UFineNess by VFineNess. * * The normals are saved in a linear CagdVecStruct vector which is allocated * * dynamically. Data is saved u inc. first. * * This routine is much faster than evaluating normal per each point. * ******************************************************************************/ CagdVecStruct *BspSrfMeshNormals(CagdSrfStruct *Srf, int UFineNess, int VFineNess) { int i, j; CagdRType u, v, UMin, UMax, VMin, VMax; CagdVecStruct *Normals, *NPtr, *T, T1, T2; CagdCrvStruct **UCurves, **VCurves, *UCrv, *VCrv; BspSrfDomain(Srf, &UMin, &UMax, &VMin, &VMax); Normals = (CagdVecStruct *) IritMalloc(sizeof(CagdVecStruct) * UFineNess * VFineNess); UCurves = (CagdCrvStruct **) IritMalloc(sizeof(CagdCrvStruct *) * UFineNess); VCurves = (CagdCrvStruct **) IritMalloc(sizeof(CagdCrvStruct *) * VFineNess); UFineNess--; VFineNess--; for (i = 0; i <= UFineNess; i++) /* Prepare Iso U curves. */ UCurves[i] = BspSrfCrvFromSrf(Srf, UMin + (UMax - UMin) * ((CagdRType) i) / UFineNess, CAGD_CONST_U_DIR); for (j = 0; j <= VFineNess; j++) /* Prepare Iso V curves. */ VCurves[j] = BspSrfCrvFromSrf(Srf, VMin + (VMax - VMin) * ((CagdRType) j) / VFineNess, CAGD_CONST_V_DIR); NPtr = Normals; for (i = 0; i <= UFineNess; i++) { UCrv = UCurves[i]; u = UMin + (UMax - UMin) * ((CagdRType) i) / UFineNess; for (j = 0; j <= VFineNess; j++) { VCrv = VCurves[j]; v = VMin + (VMax - VMin) * ((CagdRType) j) / VFineNess; /* We need to copy the tangents as BspCrvTangent save it in as */ /* static so the second call will overwrite first call value. */ /* If we fail to compute the tangent, we try adjacent isoline. */ T = BspCrvTangent(UCrv, v); if (CAGD_LEN_VECTOR(*T) < EPSILON) T = BspCrvTangent(UCurves[i == 0 ? i + 1 : i - 1], v); CAGD_COPY_VECTOR(T1, *T); T = BspCrvTangent(VCrv, u); if (CAGD_LEN_VECTOR(*T) < EPSILON) T = BspCrvTangent(VCurves[j == 0 ? j + 1 : j - 1], u); CAGD_COPY_VECTOR(T2, *T); /* The normal is the cross product of T1 and T2: */ NPtr -> Vec[0] = T1.Vec[1] * T2.Vec[2] - T1.Vec[2] * T2.Vec[1]; NPtr -> Vec[1] = T1.Vec[2] * T2.Vec[0] - T1.Vec[0] * T2.Vec[2]; NPtr -> Vec[2] = T1.Vec[0] * T2.Vec[1] - T1.Vec[1] * T2.Vec[0]; CAGD_NORMALIZE_VECTOR(*NPtr); /* Normalize the vector. */ NPtr++; } } for (i = 0; i <= UFineNess; i++) CagdCrvFree(UCurves[i]); IritFree(UCurves); for (j = 0; j <= VFineNess; j++) CagdCrvFree(VCurves[j]); IritFree(VCurves); return Normals; }