/****************************************************************************** * SBsp_Sym.c - Bspline surface symbolic computation. * ******************************************************************************* * (C) Gershon Elber, Technion, Israel Institute of Technology * ******************************************************************************* * Written by Gershon Elber, Sep. 92. * ******************************************************************************/ #include "symb_loc.h" #define NODE_EQUAL_SHIFT 0.8 STATIC_DATA int BspMultUsingInterpolation = TRUE; STATIC_DATA CagdCrvStruct *GlblInnerProdCrv1 = NULL, *GlblInnerProdCrv2 = NULL; static CagdCrvStruct *BspCrvMultAux(CagdCrvStruct *Crv1, CagdCrvStruct *Crv2); static CagdSrfStruct *BspSrfMultAux(CagdSrfStruct *Srf1, CagdSrfStruct *Srf2); /***************************************************************************** * DESCRIPTION: M * Sets method of Bspline product computation. M * * * PARAMETERS: M * BspMultUsingInter: If TRUE, Bspline product is computed by setting an M * interpolation problem. Otherwise, by decomposing the M * Bspline geometry to Bezier geometry. M * * * RETURN VALUE: M * int: Previous setting. M * * * KEYWORDS: M * BspMultInterpFlag, product M *****************************************************************************/ int BspMultInterpFlag(int BspMultUsingInter) { int i = BspMultUsingInterpolation; BspMultUsingInterpolation = BspMultUsingInter; return i; } /***************************************************************************** * DESCRIPTION: M * Given two Bspline curves - multiply them coordinatewise. M * The two curves are promoted to same point type before the multiplication M * can take place. M * See also BspMultInterpFlag. M * * * PARAMETERS: M * Crv1, Crv2: The two curves to multiply. M * * * RETURN VALUE: M * CagdCrvStruct *: The product Crv1 * Crv2 coordinatewise. M * * * KEYWORDS: M * BspCrvMult, product M *****************************************************************************/ CagdCrvStruct *BspCrvMult(CagdCrvStruct *Crv1, CagdCrvStruct *Crv2) { CagdCrvStruct *ProdCrv, *TCrv; Crv1 = CagdCrvCopy(Crv1); Crv2 = CagdCrvCopy(Crv2); if (!CagdMakeCrvsCompatible(&Crv1, &Crv2, FALSE, FALSE)) { SYMB_FATAL_ERROR(SYMB_ERR_CRV_FAIL_CMPT); return NULL; } /* Use interpolation to compute the product only if so requested and at */ /* least one of the two curves is non constant. */ if (BspMultUsingInterpolation && (Crv1 -> Order > 1 || Crv2 -> Order > 1)) { CagdPointType PType = Crv1 -> PType; CagdBType IsRational = CAGD_IS_RATIONAL_PT(PType); int i, j, KVLen, ResLength, NumCoords = CAGD_NUM_OF_PT_COORD(PType), ResOrder = Crv1 -> Order + Crv2 -> Order - 1; CagdRType *R, *KV = BspKnotContinuityMergeTwo(Crv1 -> KnotVector, Crv1 -> Length + Crv1 -> Order, Crv1 -> Order, Crv2 -> KnotVector, Crv2 -> Length + Crv2 -> Order, Crv2 -> Order, ResOrder, &KVLen), *KVNodes = BspKnotNodes(KV, KVLen, ResOrder); CagdCtlPtStruct *CtlPt = NULL, *CtlPtList = NULL; ResLength = KVLen - ResOrder; /* Verify that all nodes are distinct. */ for (i = 0, R = KVNodes; i < ResLength - 1; i++, R++) { if (APX_EQ_EPS(R[0], R[1], IRIT_UEPS)) { if (i > 0) R[0] = R[0] * NODE_EQUAL_SHIFT + R[-1] * (1 - NODE_EQUAL_SHIFT); } } /* Evaluate the multiplication at the node values. */ for (i = 0, R = KVNodes; i < ResLength; i++, R++) { CagdRType *Evl; if (CtlPt == NULL) CtlPt = CtlPtList = CagdCtlPtNew(PType); else { CtlPt -> Pnext = CagdCtlPtNew(PType); CtlPt = CtlPt -> Pnext; } Evl = CagdCrvEval(Crv1, *R); SYMB_GEN_COPY(CtlPt -> Coords, Evl, CAGD_MAX_PT_SIZE * sizeof(CagdRType)); Evl = CagdCrvEval(Crv2, *R); for (j = !IsRational; j <= NumCoords; j++) CtlPt -> Coords[j] *= Evl[j]; } if ((ProdCrv = BspCrvInterpolate(CtlPtList, ResLength, KVNodes, KV, ResLength, ResOrder, FALSE)) == NULL) { SYMB_FATAL_ERROR(SYMB_ERR_SPL_PROD_FAILED); } IritFree(KVNodes); IritFree(KV); CagdCtlPtFreeList(CtlPtList); } else { TCrv = BspCrvMultAux(Crv1, Crv2); if (CAGD_IS_BEZIER_CRV(TCrv)) { ProdCrv = CnvrtBezier2BsplineCrv(TCrv); CagdCrvFree(TCrv); } else ProdCrv = TCrv; } CagdCrvFree(Crv1); CagdCrvFree(Crv2); return ProdCrv; } /***************************************************************************** * DESCRIPTION: M * Computes a matrix of size (Len x (Len - (Order1 - Order2)) of inner M * products, SymbBspBasisInnerProd style. matrix is allocated dynamically. M * * * PARAMETERS: M * KV: Knot vector of the basis functions (of Order). M * Len: Length of knot vector KV. M * Order1: Order of first basis function. M * Order2: Order of second basis function, <= Order1. M * * * RETURN VALUE: M * CagdRType **: The allocated matrix and values of inner products. M * * * SEE ALSO: M * SymbBspBasisInnerProd, SymbBspBasisInnerProdPrep M * * * KEYWORDS: M * SymbBspBasisInnerProdMat M *****************************************************************************/ CagdRType **SymbBspBasisInnerProdMat(CagdRType *KV, int Len, int Order1, int Order2) { int i, j, PtsLen = Len - Order1, DOrder = Order1 - Order2; CagdRType **M; M = (CagdRType **) IritMalloc(sizeof(CagdRType *) * PtsLen); for (i = 0; i < PtsLen; i++) M[i] = (CagdRType *) IritMalloc(sizeof(CagdRType) * (PtsLen - DOrder)); SymbBspBasisInnerProdPrep(KV, Len, Order1, Order2); for (i = 0; i < PtsLen; i++) for (j = 0; j < PtsLen - DOrder; j++) M[i][j] = SymbBspBasisInnerProd(i, j); #ifdef DEBUG { IRIT_SET_IF_DEBUG_ON_PARAMETER(_DebugPrintInnerProdMat, FALSE) { fprintf(stderr, "Inner product matrix of size %d x %d\n", PtsLen, PtsLen - DOrder); for (i = 0; i < PtsLen; i++) { for (j = 0; j < PtsLen - DOrder; j++) fprintf(stderr, "%5.3f ", M[i][j]); fprintf(stderr, "\n"); } } } #endif /* DEBUG */ return M; } /***************************************************************************** * DESCRIPTION: M * Prepares for the computation of the inner product of pair of basis M * functions over a similar function space. M * The inner product is defined as "int( B1(t) * B2(t) )" where "int ( . )" M * denotes the integral of the function over all the domain. M * * * PARAMETERS: M * KV: Knot vector of the basis functions (of Order). M * Len: Length of knot vector KV. M * Order1: Order of first basis function. M * Order2: Order of second basis function, <= Order1. M * * * RETURN VALUE: M * void M * * * SEE ALSO: M * SymbBspBasisInnerProd, SymbBspBasisInnerProdPrep2 M * * * KEYWORDS: M * SymbBspBasisInnerProdPrep M *****************************************************************************/ void SymbBspBasisInnerProdPrep(CagdRType *KV, int Len, int Order1, int Order2) { if (GlblInnerProdCrv1 != NULL) CagdCrvFree(GlblInnerProdCrv1); if (GlblInnerProdCrv2 != NULL) CagdCrvFree(GlblInnerProdCrv2); /* Construct two zero curves with the proper orders and knot vectors. */ GlblInnerProdCrv1 = BspCrvNew(Len - Order1, Order1, CAGD_PT_E1_TYPE); ZAP_MEM(GlblInnerProdCrv1 -> Points[1], sizeof(CagdRType) * GlblInnerProdCrv1 -> Length); CAGD_GEN_COPY(GlblInnerProdCrv1 -> KnotVector, KV, sizeof(CagdRType) * Len); GlblInnerProdCrv2 = BspCrvNew(Len - Order2 - (Order1 - Order2) * 2, Order2, CAGD_PT_E1_TYPE); ZAP_MEM(GlblInnerProdCrv2 -> Points[1], sizeof(CagdRType) * GlblInnerProdCrv2 -> Length); CAGD_GEN_COPY(GlblInnerProdCrv2 -> KnotVector, &KV[Order1 - Order2], sizeof(CagdRType) * (Len - (Order1 - Order2) * 2)); } /***************************************************************************** * DESCRIPTION: M * Prepares for the computation of the inner product of pair of basis M * functions over a similar function space. M * The inner product is defined as "int( B1(t) * B2(t) )" where "int ( . )" M * denotes the integral of the function over all the domain. M * * * PARAMETERS: M * KV1: Knot vector of the first basis functions (of Order1). M * KV2: Knot vector of the second basis functions (of Order2). M * Len1: Length of knot vector KV1. M * Len2: Length of knot vector KV2. M * Order1: Order of first basis function. M * Order2: Order of second basis function, <= Order1. M * * * RETURN VALUE: M * void M * * * SEE ALSO: M * SymbBspBasisInnerProd2, SymbBspBasisInnerProdPrep M * * * KEYWORDS: M * SymbBspBasisInnerProdPrep2 M *****************************************************************************/ void SymbBspBasisInnerProdPrep2(CagdRType *KV1, CagdRType *KV2, int Len1, int Len2, int Order1, int Order2) { if (GlblInnerProdCrv1 != NULL) CagdCrvFree(GlblInnerProdCrv1); if (GlblInnerProdCrv2 != NULL) CagdCrvFree(GlblInnerProdCrv2); /* Construct two zero curves with the proper orders and knot vectors. */ GlblInnerProdCrv1 = BspCrvNew(Len1 - Order1, Order1, CAGD_PT_E1_TYPE); ZAP_MEM(GlblInnerProdCrv1 -> Points[1], sizeof(CagdRType) * GlblInnerProdCrv1 -> Length); CAGD_GEN_COPY(GlblInnerProdCrv1 -> KnotVector, KV1, sizeof(CagdRType) * Len1); GlblInnerProdCrv2 = BspCrvNew(Len2 - Order2, Order2, CAGD_PT_E1_TYPE); ZAP_MEM(GlblInnerProdCrv2 -> Points[1], sizeof(CagdRType) * GlblInnerProdCrv2 -> Length); CAGD_GEN_COPY(GlblInnerProdCrv2 -> KnotVector, KV2, sizeof(CagdRType) * Len2); } /***************************************************************************** * DESCRIPTION: M * Computes the inner product of two basis functions over a similar M * function space, as created via SymbBspBasisInnerProdPrep. M * The inner product is defined as "int( B1(t) * B2(t) )" where "int ( . )" M * denotes the integral of the function over all the domain. M * * * PARAMETERS: M * Index1: Index of first basis function. M * Index2: Index of second basis function. M * * * RETURN VALUE: M * CagdRType: The value of the inner product. M * * * SEE ALSO: M * SymbBspBasisInnerProdPrep, SymbBspBasisInnerProdMat M * SymbBspBasisInnerProd2 M * * * KEYWORDS: M * SymbBspBasisInnerProd M *****************************************************************************/ CagdRType SymbBspBasisInnerProd(int Index1, int Index2) { int Order1, Order2; CagdRType *R, *KV1, *KV2, FirstVal, LastVal, MinVal, MaxVal; CagdCrvStruct *Crv1, *Crv2, *TCrv, *ICrv; if (GlblInnerProdCrv1 == NULL || GlblInnerProdCrv2 == NULL) return 0.0; KV1 = GlblInnerProdCrv1 -> KnotVector; KV2 = GlblInnerProdCrv2 -> KnotVector; Order1 = GlblInnerProdCrv1 -> Order; Order2 = GlblInnerProdCrv2 -> Order; CagdCrvDomain(GlblInnerProdCrv1, &MinVal, &MaxVal); if (Index1 < 0 || Index1 >= GlblInnerProdCrv1 -> Length || Index2 < 0 || Index2 >= GlblInnerProdCrv2 -> Length || KV1[Index1 + Order1] <= KV2[Index2] || KV2[Index2 + Order2] <= KV1[Index1]) return 0.0; /* No common domain to basis functions. */ GlblInnerProdCrv1 -> Points[1][Index1] = 1.0; GlblInnerProdCrv2 -> Points[1][Index2] = 1.0; /* Compute indices. */ FirstVal = MAX(KV1[Index1], KV2[Index2]); FirstVal = MAX(FirstVal, MinVal); LastVal = MIN(KV1[Index1 + Order1], KV2[Index2 + Order2]); LastVal = MIN(LastVal, MaxVal); /* Construct the two curves with minimal supporting domain. */ Crv1 = CagdCrvRegionFromCrv(GlblInnerProdCrv1, FirstVal, LastVal); Crv2 = CagdCrvRegionFromCrv(GlblInnerProdCrv2, FirstVal, LastVal); TCrv = BspCrvMultAux(Crv1, Crv2); CagdCrvFree(Crv1); CagdCrvFree(Crv2); ICrv = BspCrvIntegrate(TCrv); CagdCrvFree(TCrv); R = CagdCrvEval(ICrv, LastVal); CagdCrvFree(ICrv); GlblInnerProdCrv1 -> Points[1][Index1] = 0.0; GlblInnerProdCrv2 -> Points[1][Index2] = 0.0; return R[1]; } /***************************************************************************** * DESCRIPTION: * * Auxiliary routine. Subdivides the curves into Bezier curves, multiply * * the Bezier curves and merge them back. All is done simultaneously. * * * * PARAMETERS: * * Crv1, Crv2: The two curves to multiply. * * * * RETURN VALUE: * * CagdCrvStruct *: The product Crv1 * Crv2 coordinatewise. * *****************************************************************************/ static CagdCrvStruct *BspCrvMultAux(CagdCrvStruct *Crv1, CagdCrvStruct *Crv2) { CagdCrvStruct *Crv1a, *Crv1b, *Crv2a, *Crv2b, *CrvA, *CrvB, *ProdCrv; if (Crv1 -> Order != Crv1 -> Length || Crv2 -> Order != Crv2 -> Length) { CagdRType SubdivVal = Crv1 -> Order != Crv1 -> Length ? Crv1 -> KnotVector[(Crv1 -> Length + Crv1 -> Order) >> 1] : Crv2 -> KnotVector[(Crv2 -> Length + Crv2 -> Order) >> 1]; /* Subdivide. */ Crv1a = BspCrvSubdivAtParam(Crv1, SubdivVal); Crv1b = Crv1a -> Pnext; Crv1a -> Pnext = NULL; Crv2a = BspCrvSubdivAtParam(Crv2, SubdivVal); Crv2b = Crv2a -> Pnext; Crv2a -> Pnext = NULL; CrvA = BspCrvMultAux(Crv1a, Crv2a); CrvB = BspCrvMultAux(Crv1b, Crv2b); CagdCrvFree(Crv1a); CagdCrvFree(Crv1b); CagdCrvFree(Crv2a); CagdCrvFree(Crv2b); ProdCrv = CagdMergeCrvCrv(CrvA, CrvB, FALSE); CagdCrvFree(CrvA); CagdCrvFree(CrvB); } else { int i; CagdRType TMin, TMax; CagdCrvStruct *Crv1Bzr = CnvrtBspline2BezierCrv(Crv1), *Crv2Bzr = CnvrtBspline2BezierCrv(Crv2), *ProdCrvAux = BzrCrvMult(Crv1Bzr, Crv2Bzr); CagdCrvDomain(Crv1, &TMin, &TMax); ProdCrv = CnvrtBezier2BsplineCrv(ProdCrvAux); for (i = 0; i < ProdCrv -> Order; i++) { ProdCrv -> KnotVector[i] = TMin; ProdCrv -> KnotVector[i + ProdCrv -> Order] = TMax; } CagdCrvFree(Crv1Bzr); CagdCrvFree(Crv2Bzr); CagdCrvFree(ProdCrvAux); } return ProdCrv; } /***************************************************************************** * DESCRIPTION: M * Given two Bspline surfaces - multiply them coordinatewise. M * The two surfaces are promoted to same point type before multiplication M * can take place. M * See also BspMultInterpFlag. M * * * PARAMETERS: M * Srf1, Srf2: The two surfaces to multiply. M * * * RETURN VALUE: M * CagdSrfStruct *: The product Srf1 * Srf2 coordinatewise. M * * * KEYWORDS: M * BspSrfMult, product M *****************************************************************************/ CagdSrfStruct *BspSrfMult(CagdSrfStruct *Srf1, CagdSrfStruct *Srf2) { CagdSrfStruct *ProdSrf, *TSrf; Srf1 = CagdSrfCopy(Srf1); Srf2 = CagdSrfCopy(Srf2); if (!CagdMakeSrfsCompatible(&Srf1, &Srf2, FALSE, FALSE, FALSE, FALSE)) { SYMB_FATAL_ERROR(SYMB_ERR_SRF_FAIL_CMPT); return NULL; } TSrf = BspSrfMultAux(Srf1, Srf2); if (CAGD_IS_BEZIER_SRF(TSrf)) { ProdSrf = CnvrtBezier2BsplineSrf(TSrf); CagdSrfFree(TSrf); } else ProdSrf = TSrf; CagdSrfFree(Srf1); CagdSrfFree(Srf2); return ProdSrf; } /***************************************************************************** * DESCRIPTION: * * Auxiliary routine. Subdivides the surfaces into Bezier surfaces, multiply * * the Bezier surfaces and merge them back. All is done simultaneously. * * * * PARAMETERS: * * Srf1, Srf2: The two surfaces to multiply. * * * * RETURN VALUE: * * CagdSrfStruct *: The product Srf1 * Srf2 coordinatewise. * *****************************************************************************/ static CagdSrfStruct *BspSrfMultAux(CagdSrfStruct *Srf1, CagdSrfStruct *Srf2) { CagdSrfStruct *Srf1a, *Srf1b, *Srf2a, *Srf2b, *SrfA, *SrfB, *ProdSrf; if (Srf1 -> UOrder != Srf1 -> ULength || Srf2 -> UOrder != Srf2 -> ULength) { CagdRType SubdivVal = Srf1 -> UOrder != Srf1 -> ULength ? Srf1 -> UKnotVector[(Srf1 -> ULength + Srf1 -> UOrder) >> 1] : Srf2 -> UKnotVector[(Srf2 -> ULength + Srf2 -> UOrder) >> 1]; int Mult1 = BspKnotFindMult(Srf1 -> UKnotVector, Srf1 -> UOrder, Srf1 -> ULength, SubdivVal), Mult2 = BspKnotFindMult(Srf2 -> UKnotVector, Srf2 -> UOrder, Srf2 -> ULength, SubdivVal), C0Discont = Mult1 >= Srf1 -> UOrder || Mult2 >= Srf2 -> UOrder; /* Subdivide along U. */ Srf1a = BspSrfSubdivAtParam(Srf1, SubdivVal, CAGD_CONST_U_DIR); Srf1b = Srf1a -> Pnext; Srf1a -> Pnext = NULL; Srf2a = BspSrfSubdivAtParam(Srf2, SubdivVal, CAGD_CONST_U_DIR); Srf2b = Srf2a -> Pnext; Srf2a -> Pnext = NULL; SrfA = BspSrfMultAux(Srf1a, Srf2a); SrfB = BspSrfMultAux(Srf1b, Srf2b); CagdSrfFree(Srf1a); CagdSrfFree(Srf1b); CagdSrfFree(Srf2a); CagdSrfFree(Srf2b); ProdSrf = CagdMergeSrfSrf(SrfA, SrfB, CAGD_CONST_U_DIR, !C0Discont, FALSE); CagdSrfFree(SrfA); CagdSrfFree(SrfB); } else if (Srf1 -> VOrder != Srf1 -> VLength || Srf2 -> VOrder != Srf2 -> VLength) { CagdRType SubdivVal = Srf1 -> VOrder != Srf1 -> VLength ? Srf1 -> VKnotVector[(Srf1 -> VLength + Srf1 -> VOrder) >> 1] : Srf2 -> VKnotVector[(Srf2 -> VLength + Srf2 -> VOrder) >> 1]; int Mult1 = BspKnotFindMult(Srf1 -> VKnotVector, Srf1 -> VOrder, Srf1 -> VLength, SubdivVal), Mult2 = BspKnotFindMult(Srf2 -> VKnotVector, Srf2 -> VOrder, Srf2 -> VLength, SubdivVal), C0Discont = Mult1 >= Srf1 -> VOrder || Mult2 >= Srf2 -> VOrder; /* Subdivide along V. */ Srf1a = BspSrfSubdivAtParam(Srf1, SubdivVal, CAGD_CONST_V_DIR); Srf1b = Srf1a -> Pnext; Srf1a -> Pnext = NULL; Srf2a = BspSrfSubdivAtParam(Srf2, SubdivVal, CAGD_CONST_V_DIR); Srf2b = Srf2a -> Pnext; Srf2a -> Pnext = NULL; SrfA = BspSrfMultAux(Srf1a, Srf2a); SrfB = BspSrfMultAux(Srf1b, Srf2b); CagdSrfFree(Srf1a); CagdSrfFree(Srf1b); CagdSrfFree(Srf2a); CagdSrfFree(Srf2b); ProdSrf = CagdMergeSrfSrf(SrfA, SrfB, CAGD_CONST_V_DIR, !C0Discont, FALSE); CagdSrfFree(SrfA); CagdSrfFree(SrfB); } else { int i; CagdRType UMin, UMax, VMin, VMax; CagdSrfStruct *Srf1Bzr = CnvrtBspline2BezierSrf(Srf1), *Srf2Bzr = CnvrtBspline2BezierSrf(Srf2), *ProdSrfAux = BzrSrfMult(Srf1Bzr, Srf2Bzr); CagdSrfDomain(Srf1, &UMin, &UMax, &VMin, &VMax); ProdSrf = CnvrtBezier2BsplineSrf(ProdSrfAux); for (i = 0; i < ProdSrf -> UOrder; i++) { ProdSrf -> UKnotVector[i] = UMin; ProdSrf -> UKnotVector[i + ProdSrf -> UOrder] = UMax; } for (i = 0; i < ProdSrf -> VOrder; i++) { ProdSrf -> VKnotVector[i] = VMin; ProdSrf -> VKnotVector[i + ProdSrf -> VOrder] = VMax; } CagdSrfFree(Srf1Bzr); CagdSrfFree(Srf2Bzr); CagdSrfFree(ProdSrfAux); } return ProdSrf; } /***************************************************************************** * DESCRIPTION: M * Factors out a (u - v) term from a surface, assuming it has one. M * * * PARAMETERS: M * Srf: To factor out a (u - v) term from. M * * * RETURN VALUE: M * CagdSrfStruct *: Srf / (u - v). M * * * SEE ALSO: M * BzrSrfFactorUMinusV, BzrSrfFactorBilinear M * * * KEYWORDS: M * BspSrfFactorUMinusV M *****************************************************************************/ CagdSrfStruct *BspSrfFactorUMinusV(CagdSrfStruct *Srf) { CagdSrfStruct *Srf1 , *Srf2, *Srf1f, *Srf2f, *FactorSrf; if (Srf -> UOrder != Srf -> ULength) { CagdRType SubdivVal = Srf -> UKnotVector[(Srf -> ULength + Srf -> UOrder) >> 1]; int Mult = BspKnotFindMult(Srf -> UKnotVector, Srf -> UOrder, Srf -> ULength, SubdivVal), C0Discont = Mult >= Srf -> UOrder; /* Subdivide along U. */ Srf1 = BspSrfSubdivAtParam(Srf, SubdivVal, CAGD_CONST_U_DIR); Srf2 = Srf1 -> Pnext; Srf1 -> Pnext = NULL; Srf1f = BspSrfFactorUMinusV(Srf1); Srf2f = BspSrfFactorUMinusV(Srf2); CagdSrfFree(Srf1); CagdSrfFree(Srf2); FactorSrf = CagdMergeSrfSrf(Srf1f, Srf2f, CAGD_CONST_U_DIR, !C0Discont, FALSE); CagdSrfFree(Srf1f); CagdSrfFree(Srf2f); } else if (Srf -> VOrder != Srf -> VLength) { CagdRType SubdivVal = Srf -> VKnotVector[(Srf -> VLength + Srf -> VOrder) >> 1]; int Mult = BspKnotFindMult(Srf -> VKnotVector, Srf -> VOrder, Srf -> VLength, SubdivVal), C0Discont = Mult >= Srf -> VOrder; /* Subdivide along V. */ Srf1 = BspSrfSubdivAtParam(Srf, SubdivVal, CAGD_CONST_V_DIR); Srf2 = Srf1 -> Pnext; Srf1 -> Pnext = NULL; Srf1f = BspSrfFactorUMinusV(Srf1); Srf2f = BspSrfFactorUMinusV(Srf2); CagdSrfFree(Srf1); CagdSrfFree(Srf2); FactorSrf = CagdMergeSrfSrf(Srf1f, Srf2f, CAGD_CONST_V_DIR, !C0Discont, FALSE); CagdSrfFree(Srf1f); CagdSrfFree(Srf2f); } else { int i; CagdRType UMin, UMax, VMin, VMax, A[4]; CagdSrfStruct *FactorSrfBzr, *SrfBzr = CnvrtBspline2BezierSrf(Srf); /* Bring to the appropriate domain. */ CagdSrfDomain(Srf, &UMin, &UMax, &VMin, &VMax); /* * Srf[u,v]/(u,v) = * SrfBzr[s,t] / {(u1-u0) s - (v1-v0) t + (u0-v0) * (u1-u0) s - (v1-v0) t + (u0-v0) * = (u0-v0)(1-s)(1-t) + (u1-v0)s(1-t) * (u0-v1)(1-s)t + (u1-v1)st */ A[0] = UMin - VMin; A[1] = UMax - VMin; A[2] = UMin - VMax; A[3] = UMax - VMax; FactorSrfBzr = BzrSrfFactorBilinear(SrfBzr, A); FactorSrf = CnvrtBezier2BsplineSrf(FactorSrfBzr); for (i = 0; i < FactorSrf -> UOrder; i++) { FactorSrf -> UKnotVector[i] = UMin; FactorSrf -> UKnotVector[i + FactorSrf -> UOrder] = UMax; } for (i = 0; i < FactorSrf -> VOrder; i++) { FactorSrf -> VKnotVector[i] = VMin; FactorSrf -> VKnotVector[i + FactorSrf -> VOrder] = VMax; } CagdSrfFree(SrfBzr); CagdSrfFree(FactorSrfBzr); } return FactorSrf; } /***************************************************************************** * DESCRIPTION: M * Given a rational Bspline curve - computes its derivative curve (Hodograph) M * using the quotient rule for differentiation. M * * * PARAMETERS: M * Crv: Bspline curve to differentiate. M * * * RETURN VALUE: M * CagdCrvStruct *: Differentiated rational Bspline curve. M * * * SEE ALSO: M * BzrCrvDerive, BspCrvDerive, BzrCrvDeriveRational, CagdCrvDerive M * * * KEYWORDS: M * BspCrvDeriveRational, derivatives M *****************************************************************************/ CagdCrvStruct *BspCrvDeriveRational(CagdCrvStruct *Crv) { CagdCrvStruct *CrvW, *CrvX, *CrvY, *CrvZ, *TCrv1, *TCrv2, *DeriveCrv, *DCrvW = NULL, *DCrvX = NULL, *DCrvY = NULL, *DCrvZ = NULL; SymbCrvSplitScalar(Crv, &CrvW, &CrvX, &CrvY, &CrvZ); if (CrvW) DCrvW = BspCrvDerive(CrvW); else { SYMB_FATAL_ERROR(SYMB_ERR_RATIONAL_EXPECTED); return NULL; } if (CrvX) { DCrvX = BspCrvDerive(CrvX); TCrv1 = BspCrvMult(DCrvX, CrvW); TCrv2 = BspCrvMult(CrvX, DCrvW); CagdCrvFree(CrvX); CrvX = SymbCrvSub(TCrv1, TCrv2); CagdCrvFree(TCrv1); CagdCrvFree(TCrv2); } if (CrvY) { DCrvY = BspCrvDerive(CrvY); TCrv1 = BspCrvMult(DCrvY, CrvW); TCrv2 = BspCrvMult(CrvY, DCrvW); CagdCrvFree(CrvY); CrvY = SymbCrvSub(TCrv1, TCrv2); CagdCrvFree(TCrv1); CagdCrvFree(TCrv2); } if (CrvZ) { DCrvZ = BspCrvDerive(CrvZ); TCrv1 = BspCrvMult(DCrvZ, CrvW); TCrv2 = BspCrvMult(CrvZ, DCrvW); CagdCrvFree(CrvZ); CrvZ = SymbCrvSub(TCrv1, TCrv2); CagdCrvFree(TCrv1); CagdCrvFree(TCrv2); } TCrv1 = BspCrvMult(CrvW, CrvW); CagdCrvFree(CrvW); CrvW = TCrv1; if (!CagdMakeCrvsCompatible(&CrvW, &CrvX, TRUE, TRUE) || !CagdMakeCrvsCompatible(&CrvW, &CrvY, TRUE, TRUE) || !CagdMakeCrvsCompatible(&CrvW, &CrvZ, TRUE, TRUE)) { SYMB_FATAL_ERROR(SYMB_ERR_CRV_FAIL_CMPT); return NULL; } /* Note CrvX/Y/Z might be different due to possible C0 discontinuities. */ if (!CagdMakeCrvsCompatible(&CrvX, &CrvY, TRUE, TRUE) || !CagdMakeCrvsCompatible(&CrvX, &CrvZ, TRUE, TRUE) || !CagdMakeCrvsCompatible(&CrvY, &CrvZ, TRUE, TRUE) || !CagdMakeCrvsCompatible(&CrvW, &CrvX, TRUE, TRUE) || !CagdMakeCrvsCompatible(&CrvW, &CrvY, TRUE, TRUE) || !CagdMakeCrvsCompatible(&CrvW, &CrvZ, TRUE, TRUE)) { SYMB_FATAL_ERROR(SYMB_ERR_CRV_FAIL_CMPT); return NULL; } DeriveCrv = SymbCrvMergeScalar(CrvW, CrvX, CrvY, CrvZ); if (CrvX) { CagdCrvFree(CrvX); CagdCrvFree(DCrvX); } if (CrvY) { CagdCrvFree(CrvY); CagdCrvFree(DCrvY); } if (CrvZ) { CagdCrvFree(CrvZ); CagdCrvFree(DCrvZ); } if (CrvW) { CagdCrvFree(CrvW); CagdCrvFree(DCrvW); } return DeriveCrv; } /***************************************************************************** * DESCRIPTION: M * Given a rational Bspline surface - computes its derivative surface in M * direction Dir, using the quotient rule for differentiation. M * * * PARAMETERS: M * Srf: Bspline surface to differentiate. M * Dir: Direction of Differentiation. Either U or V. M * * * RETURN VALUE: M * CagdSrfStruct *: Differentiated rational Bspline surface. M * * * SEE ALSO: M * CagdSrfDerive, BzrSrfDerive, BspSrfDerive, BzrSrfDeriveRational M * * * KEYWORDS: M * BspSrfDeriveRational, derivatives M *****************************************************************************/ CagdSrfStruct *BspSrfDeriveRational(CagdSrfStruct *Srf, CagdSrfDirType Dir) { CagdSrfStruct *SrfW, *SrfX, *SrfY, *SrfZ, *TSrf1, *TSrf2, *DeriveSrf, *DSrfW = NULL, *DSrfX = NULL, *DSrfY = NULL, *DSrfZ = NULL; SymbSrfSplitScalar(Srf, &SrfW, &SrfX, &SrfY, &SrfZ); if (SrfW) DSrfW = BspSrfDerive(SrfW, Dir); else { SYMB_FATAL_ERROR(SYMB_ERR_RATIONAL_EXPECTED); return NULL; } if (SrfX) { DSrfX = BspSrfDerive(SrfX, Dir); TSrf1 = BspSrfMult(DSrfX, SrfW); TSrf2 = BspSrfMult(SrfX, DSrfW); CagdSrfFree(SrfX); SrfX = SymbSrfSub(TSrf1, TSrf2); CagdSrfFree(TSrf1); CagdSrfFree(TSrf2); } if (SrfY) { DSrfY = BspSrfDerive(SrfY, Dir); TSrf1 = BspSrfMult(DSrfY, SrfW); TSrf2 = BspSrfMult(SrfY, DSrfW); CagdSrfFree(SrfY); SrfY = SymbSrfSub(TSrf1, TSrf2); CagdSrfFree(TSrf1); CagdSrfFree(TSrf2); } if (SrfZ) { DSrfZ = BspSrfDerive(SrfZ, Dir); TSrf1 = BspSrfMult(DSrfZ, SrfW); TSrf2 = BspSrfMult(SrfZ, DSrfW); CagdSrfFree(SrfZ); SrfZ = SymbSrfSub(TSrf1, TSrf2); CagdSrfFree(TSrf1); CagdSrfFree(TSrf2); } TSrf1 = BspSrfMult(SrfW, SrfW); CagdSrfFree(SrfW); SrfW = TSrf1; if (!CagdMakeSrfsCompatible(&SrfW, &SrfX, TRUE, TRUE, TRUE, TRUE) || !CagdMakeSrfsCompatible(&SrfW, &SrfY, TRUE, TRUE, TRUE, TRUE) || !CagdMakeSrfsCompatible(&SrfW, &SrfZ, TRUE, TRUE, TRUE, TRUE)) { SYMB_FATAL_ERROR(SYMB_ERR_SRF_FAIL_CMPT); return NULL; } /* Note SrfX/Y/Z might be different due to possible C0 discontinuities. */ if (!CagdMakeSrfsCompatible(&SrfX, &SrfY, TRUE, TRUE, TRUE, TRUE) || !CagdMakeSrfsCompatible(&SrfX, &SrfZ, TRUE, TRUE, TRUE, TRUE) || !CagdMakeSrfsCompatible(&SrfY, &SrfZ, TRUE, TRUE, TRUE, TRUE) || !CagdMakeSrfsCompatible(&SrfW, &SrfX, TRUE, TRUE, TRUE, TRUE) || !CagdMakeSrfsCompatible(&SrfW, &SrfY, TRUE, TRUE, TRUE, TRUE) || !CagdMakeSrfsCompatible(&SrfW, &SrfZ, TRUE, TRUE, TRUE, TRUE)) { SYMB_FATAL_ERROR(SYMB_ERR_SRF_FAIL_CMPT); return NULL; } DeriveSrf = SymbSrfMergeScalar(SrfW, SrfX, SrfY, SrfZ); if (SrfX) { CagdSrfFree(SrfX); CagdSrfFree(DSrfX); } if (SrfY) { CagdSrfFree(SrfY); CagdSrfFree(DSrfY); } if (SrfZ) { CagdSrfFree(SrfZ); CagdSrfFree(DSrfZ); } if (SrfW) { CagdSrfFree(SrfW); CagdSrfFree(DSrfW); } return DeriveSrf; }