/****************************************************************************** * Symbolic.c - Generic symbolic computation. * ******************************************************************************* * (C) Gershon Elber, Technion, Israel Institute of Technology * ******************************************************************************* * Written by Gershon Elber, Nov. 92. * ******************************************************************************/ #include #include #include #include "symb_loc.h" static CagdCrvStruct *SymbCrvAddSubAux(CagdCrvStruct *Crv1, CagdCrvStruct *Crv2, CagdBType OperationAdd); static CagdSrfStruct *SymbSrfAddSubAux(CagdSrfStruct *Srf1, CagdSrfStruct *Srf2, CagdBType OperationAdd); /***************************************************************************** * DESCRIPTION: M * Given two curves - add them coordinatewise. M * The two curves are promoted to same point type before the multiplication M * can take place. Furthermore, order and continuity are matched as well. M * * * PARAMETERS: M * Crv1, Crv2: Two curve to add up coordinatewise. M * * * RETURN VALUE: M * CagdCrvStruct *: The summation of Crv1 + Crv2 coordinatewise. M * * * SEE ALSO: M * SymbCrvSub, SymbCrvMult M * * * KEYWORDS: M * SymbCrvAdd, addition, symbolic computation M *****************************************************************************/ CagdCrvStruct *SymbCrvAdd(CagdCrvStruct *Crv1, CagdCrvStruct *Crv2) { return SymbCrvAddSubAux(Crv1, Crv2, TRUE); } /***************************************************************************** * DESCRIPTION: M * Given two curves - subtract them coordinatewise. M * The two curves are promoted to same point type before the multiplication M * can take place. Furthermore, order and continuity are matched as well. M * * * PARAMETERS: M * Crv1, Crv2: Two curve to subtract coordinatewise. M * * * RETURN VALUE: M * CagdCrvStruct *: The difference of Crv1 - Crv2 coordinatewise. M * * * SEE ALSO: M * SymbCrvAdd, SymbCrvMult M * * * KEYWORDS: M * SymbCrvSub, subtraction, symbolic computation M *****************************************************************************/ CagdCrvStruct *SymbCrvSub(CagdCrvStruct *Crv1, CagdCrvStruct *Crv2) { return SymbCrvAddSubAux(Crv1, Crv2, FALSE); } /***************************************************************************** * DESCRIPTION: M * Given two curves - multiply them coordinatewise. M * The two curves are promoted to same point type before the multiplication M * can take place. M * * * PARAMETERS: M * Crv1, Crv2: Two curve to multiply coordinatewise. M * * * RETURN VALUE: M * CagdCrvStruct *: The product of Crv1 * Crv2 coordinatewise. M * * * SEE ALSO: M * SymbCrvDotProd, SymbCrvVecDotProd, SymbCrvInvert, SymbCrvMultScalar M * * * KEYWORDS: M * SymbCrvMult, product, symbolic computation M *****************************************************************************/ CagdCrvStruct *SymbCrvMult(CagdCrvStruct *Crv1, CagdCrvStruct *Crv2) { CagdCrvStruct *ProdCrv = NULL; if (Crv1 -> GType == CAGD_CBEZIER_TYPE && Crv2 -> GType == CAGD_CBEZIER_TYPE) ProdCrv = BzrCrvMult(Crv1, Crv2); else if ((Crv1 -> GType == CAGD_CBEZIER_TYPE || Crv1 -> GType == CAGD_CBSPLINE_TYPE) && (Crv2 -> GType == CAGD_CBEZIER_TYPE || Crv2 -> GType == CAGD_CBSPLINE_TYPE)) ProdCrv = BspCrvMult(Crv1, Crv2); else SYMB_FATAL_ERROR(SYMB_ERR_UNDEF_CRV); return ProdCrv; } /***************************************************************************** * DESCRIPTION: M * Given a scalar curve, returns a scalar curve representing the reciprocal M * values, by making it rational (if was not one) and flipping the numerator M * and the denominator. M * * * PARAMETERS: M * Crv: A scalar curve to compute a reciprocal value for. M * * * RETURN VALUE: M * CagdCrvStruct *: A rational scalar curve that is equalto the M * reciprocal value of Crv. M * * * SEE ALSO: M * SymbCrvDotProd, SymbCrvVecDotProd, SymbCrvMult, SymbCrvMultScalar M * * * KEYWORDS: M * SymbCrvInvert, division, symbolic computation, reciprocal value M *****************************************************************************/ CagdCrvStruct *SymbCrvInvert(CagdCrvStruct *Crv) { int i; CagdRType *R, **Points; Crv = CagdCrvCopy(Crv); Points = Crv -> Points; switch (Crv -> PType) { case CAGD_PT_E1_TYPE: Points[0] = Points[1]; Points[1] = R = (CagdRType *) IritMalloc(sizeof(CagdRType) * Crv -> Length); for (i = 0; i < Crv -> Length; i++) *R++ = 1.0; Crv -> PType = CAGD_PT_P1_TYPE; break; case CAGD_PT_P1_TYPE: R = Points[0]; Points[0] = Points[1]; Points[1] = R; break; default: SYMB_FATAL_ERROR(SYMB_ERR_UNSUPPORT_PT); break; } return Crv; } /***************************************************************************** * DESCRIPTION: M * Given a curve, scale it by Scale. M * * * PARAMETERS: M * Crv: A curve to scale by magnitude Scale. M * Scale: Scaling factor. M * * * RETURN VALUE: M * CagdCrvStruct *: A curves scaled by Scale compared to Crv. M * * * SEE ALSO: M * SymbCrvDotProd, SymbCrvVecDotProd, SymbCrvMult, SymbCrvMultScalar M * * * KEYWORDS: M * SymbCrvScalarScale, scaling, symbolic computation M *****************************************************************************/ CagdCrvStruct *SymbCrvScalarScale(CagdCrvStruct *Crv, CagdRType Scale) { int i, j; CagdRType *R; Crv = CagdCrvCopy(Crv); for (j = 1; j <= CAGD_NUM_OF_PT_COORD(Crv -> PType); j++) for (i = 0, R = Crv -> Points[j]; i < Crv -> Length; i++) *R++ *= Scale; return Crv; } /***************************************************************************** * DESCRIPTION: M * Given two curves - computes their dot product. M * Returned curve is a scalar curve representing the dot product of the M * two given curves. M * * * PARAMETERS: M * Crv1, Crv2: Two curve to multiply and compute a dot product for. M * * * RETURN VALUE: M * CagdCrvStruct *: A scalar curve representing the dot product of M * Crv1 . Crv2. M * * * SEE ALSO: M * SymbCrvScalarScale, SymbCrvVecDotProd, SymbCrvMult, SymbCrvMultScalar M * * * KEYWORDS: M * SymbCrvDotProd, product, dot product, symbolic computation M *****************************************************************************/ CagdCrvStruct *SymbCrvDotProd(CagdCrvStruct *Crv1, CagdCrvStruct *Crv2) { CagdCrvStruct *PCrvW, *PCrvX, *PCrvY, *PCrvZ, *TCrv1, *TCrv2, *DotProdCrv, *ProdCrv = SymbCrvMult(Crv1, Crv2); SymbCrvSplitScalar(ProdCrv, &PCrvW, &PCrvX, &PCrvY, &PCrvZ); CagdCrvFree(ProdCrv); if (PCrvY != NULL) { TCrv1 = SymbCrvAdd(PCrvX, PCrvY); CagdCrvFree(PCrvX); CagdCrvFree(PCrvY); } else TCrv1 = PCrvX; if (PCrvZ != NULL) { TCrv2 = SymbCrvAdd(TCrv1, PCrvZ); CagdCrvFree(TCrv1); CagdCrvFree(PCrvZ); TCrv1 = TCrv2; } DotProdCrv = SymbCrvMergeScalar(PCrvW, TCrv1, NULL, NULL); CagdCrvFree(PCrvW); CagdCrvFree(TCrv1); return DotProdCrv; } /***************************************************************************** * DESCRIPTION: M * Given a curve and a vector - computes their dot product. M * Returned curve is a scalar curve representing the dot product. M * * * PARAMETERS: M * Crv: Curve to multiply and compute a dot product for. M * Vec: Vector to project Crv onto. M * * * RETURN VALUE: M * CagdCrvStruct *: A scalar curve representing the dot product of M * Crv . Vec. M * * * SEE ALSO: M * SymbCrvDotProd, SymbCrvMult, SymbCrvMultScalar, SymbCrvCrossProd M * * * KEYWORDS: M * SymbCrvVecDotProd, product, dot product, symbolic computation M *****************************************************************************/ CagdCrvStruct *SymbCrvVecDotProd(CagdCrvStruct *Crv, CagdVType Vec) { CagdCrvStruct *PCrvW, *PCrvX, *PCrvY, *PCrvZ, *TCrv, *DotProdCrv; SymbCrvSplitScalar(Crv, &PCrvW, &PCrvX, &PCrvY, &PCrvZ); TCrv = SymbCrvScalarScale(PCrvX, Vec[0]); CagdCrvFree(PCrvX); PCrvX = TCrv; if (PCrvY != NULL) { TCrv = SymbCrvScalarScale(PCrvY, Vec[1]); CagdCrvFree(PCrvY); PCrvY = TCrv; TCrv = SymbCrvAdd(PCrvX, PCrvY); CagdCrvFree(PCrvX); CagdCrvFree(PCrvY); PCrvX = TCrv; } if (PCrvZ != NULL) { TCrv = SymbCrvScalarScale(PCrvZ, Vec[2]); CagdCrvFree(PCrvZ); PCrvZ = TCrv; TCrv = SymbCrvAdd(PCrvX, PCrvZ); CagdCrvFree(PCrvX); CagdCrvFree(PCrvZ); PCrvX = TCrv; } DotProdCrv = SymbCrvMergeScalar(PCrvW, PCrvX, NULL, NULL); CagdCrvFree(PCrvW); CagdCrvFree(PCrvX); return DotProdCrv; } /***************************************************************************** * DESCRIPTION: M * Given two curves - a vector curve Crv1 and a scalar curve Crv2, multiply M * all Crv1's coordinates by the scalar curve Crv2. M * Returned curve is a curve representing the product of the two given M * curves. M * * * PARAMETERS: M * Crv1, Crv2: Two curve to multiply. M * * * RETURN VALUE: M * CagdCrvStruct *: A curve representing the product of Crv1 and Crv2. M * * * SEE ALSO: M * SymbCrvDotProd, SymbCrvVecDotProd, SymbCrvMult, SymbCrvCrossProd, M * SymbSrfMultScalar M * * * KEYWORDS: M * SymbCrvMultScalar, product, symbolic computation M *****************************************************************************/ CagdCrvStruct *SymbCrvMultScalar(CagdCrvStruct *Crv1, CagdCrvStruct *Crv2) { CagdCrvStruct *PCrv1W, *PCrv1X, *PCrv1Y, *PCrv1Z, *PCrv2W, *PCrv2X, *PCrv2Y, *PCrv2Z, *TCrv, *ProdCrv; SymbCrvSplitScalar(Crv1, &PCrv1W, &PCrv1X, &PCrv1Y, &PCrv1Z); SymbCrvSplitScalar(Crv2, &PCrv2W, &PCrv2X, &PCrv2Y, &PCrv2Z); TCrv = SymbCrvMult(PCrv1X, PCrv2X); CagdCrvFree(PCrv1X); PCrv1X = TCrv; if (PCrv1Y != NULL) { TCrv = SymbCrvMult(PCrv1Y, PCrv2X); CagdCrvFree(PCrv1Y); PCrv1Y = TCrv; } if (PCrv1Z != NULL) { TCrv = SymbCrvMult(PCrv1Z, PCrv2X); CagdCrvFree(PCrv1Z); PCrv1Z = TCrv; } if (PCrv1W != NULL && PCrv2W != NULL) { TCrv = SymbCrvMult(PCrv1W, PCrv2W); CagdCrvFree(PCrv1W); PCrv1W = TCrv; } else if (PCrv2W != NULL) { PCrv1W = PCrv2W; PCrv2W = NULL; } ProdCrv = SymbCrvMergeScalar(PCrv1W, PCrv1X, PCrv1Y, PCrv1Z); CagdCrvFree(PCrv1W); CagdCrvFree(PCrv1X); CagdCrvFree(PCrv1Y); CagdCrvFree(PCrv1Z); CagdCrvFree(PCrv2W); CagdCrvFree(PCrv2X); CagdCrvFree(PCrv2Y); CagdCrvFree(PCrv2Z); return ProdCrv; } /***************************************************************************** * DESCRIPTION: M * Given two curves - computes their cross product. M * Returned curve is a scalar curve representing the cross product of the M * two given curves. M * * * PARAMETERS: M * Crv1, Crv2: Two curve to multiply and compute a cross product for. M * * * RETURN VALUE: M * CagdCrvStruct *: A scalar curve representing the cross product of M * Crv1 x Crv2. M * * * SEE ALSO: M * SymbCrvDotProd, SymbCrvVecDotProd, SymbCrvMult, SymbCrvMultScalar M * * * KEYWORDS: M * SymbCrvCrossProd, product, cross product, symbolic computation M *****************************************************************************/ CagdCrvStruct *SymbCrvCrossProd(CagdCrvStruct *Crv1, CagdCrvStruct *Crv2) { CagdCrvStruct *Crv1W, *Crv1X, *Crv1Y, *Crv1Z, *Crv2W, *Crv2X, *Crv2Y, *Crv2Z, *TCrv1, *TCrv2, *CrossProdCrv, *PCrvW, *PCrvX, *PCrvY, *PCrvZ; CagdBType Crv1New = FALSE, Crv2New = FALSE; if (Crv1 -> GType != CAGD_PT_E3_TYPE && Crv1 -> GType != CAGD_PT_P3_TYPE) { Crv1 = CagdCoerceCrvTo(Crv1, CAGD_IS_RATIONAL_CRV(Crv1) ? CAGD_PT_P3_TYPE : CAGD_PT_E3_TYPE); Crv1New = TRUE; } if (Crv2 -> GType != CAGD_PT_E3_TYPE && Crv2 -> GType != CAGD_PT_P3_TYPE) { Crv2 = CagdCoerceCrvTo(Crv2, CAGD_IS_RATIONAL_CRV(Crv2) ? CAGD_PT_P3_TYPE : CAGD_PT_E3_TYPE); Crv2New = TRUE; } SymbCrvSplitScalar(Crv1, &Crv1W, &Crv1X, &Crv1Y, &Crv1Z); SymbCrvSplitScalar(Crv2, &Crv2W, &Crv2X, &Crv2Y, &Crv2Z); if (Crv1New) CagdCrvFree(Crv1); if (Crv2New) CagdCrvFree(Crv2); /* Cross product X axis. */ TCrv1 = SymbCrvMult(Crv1Y, Crv2Z); TCrv2 = SymbCrvMult(Crv2Y, Crv1Z); PCrvX = SymbCrvSub(TCrv1, TCrv2); CagdCrvFree(TCrv2); CagdCrvFree(TCrv1); /* Cross product Y axis. */ TCrv1 = SymbCrvMult(Crv1Z, Crv2X); TCrv2 = SymbCrvMult(Crv2Z, Crv1X); PCrvY = SymbCrvSub(TCrv1, TCrv2); CagdCrvFree(TCrv2); CagdCrvFree(TCrv1); /* Cross product Z axis. */ TCrv1 = SymbCrvMult(Crv1X, Crv2Y); TCrv2 = SymbCrvMult(Crv2X, Crv1Y); PCrvZ = SymbCrvSub(TCrv1, TCrv2); CagdCrvFree(TCrv1); CagdCrvFree(TCrv2); /* Cross product W axis. */ if (Crv1W || Crv2W) { if (Crv1W == NULL) PCrvW = CagdCrvCopy(Crv2W); else if (Crv2W == NULL) PCrvW = CagdCrvCopy(Crv1W); else PCrvW = SymbCrvMult(Crv1W, Crv2W); } else PCrvW = NULL; CagdCrvFree(Crv1X); CagdCrvFree(Crv1Y); CagdCrvFree(Crv1Z); CagdCrvFree(Crv1W); CagdCrvFree(Crv2X); CagdCrvFree(Crv2Y); CagdCrvFree(Crv2Z); CagdCrvFree(Crv2W); if (!CagdMakeCrvsCompatible(&PCrvX, &PCrvY, TRUE, TRUE) || !CagdMakeCrvsCompatible(&PCrvX, &PCrvZ, TRUE, TRUE) || !CagdMakeCrvsCompatible(&PCrvY, &PCrvZ, TRUE, TRUE) || !CagdMakeCrvsCompatible(&PCrvW, &PCrvX, TRUE, TRUE) || !CagdMakeCrvsCompatible(&PCrvW, &PCrvY, TRUE, TRUE) || !CagdMakeCrvsCompatible(&PCrvW, &PCrvZ, TRUE, TRUE)) SYMB_FATAL_ERROR(SYMB_ERR_CRV_FAIL_CMPT); CrossProdCrv = SymbCrvMergeScalar(PCrvW, PCrvX, PCrvY, PCrvZ); CagdCrvFree(PCrvX); CagdCrvFree(PCrvY); CagdCrvFree(PCrvZ); CagdCrvFree(PCrvW); return CrossProdCrv; } /***************************************************************************** * DESCRIPTION: M * Given two curves - multiply them using the quotient product rule: M * X = X1 W2 +/- X2 W1 V * All provided curves are assumed to be non rational scalar curves. M * Returned is a non rational scalar curve (CAGD_PT_E1_TYPE). M * * * PARAMETERS: M * Crv1X: Numerator of first curve. M * Crv1W: Denominator of first curve. M * Crv2X: Numerator of second curve. M * Crv2W: Denominator of second curve. M * OperationAdd: TRUE for addition, FALSE for subtraction. M * * * RETURN VALUE: M * CagdCrvStruct *: The result of Crv1X Crv2W +/- Crv2X Crv1W. M * * * SEE ALSO: M * SymbCrvDotProd, SymbCrvVecDotProd, SymbCrvMult, SymbCrvMultScalar M * * * KEYWORDS: M * SymbCrvRtnlMult, product, symbolic computation M *****************************************************************************/ CagdCrvStruct *SymbCrvRtnlMult(CagdCrvStruct *Crv1X, CagdCrvStruct *Crv1W, CagdCrvStruct *Crv2X, CagdCrvStruct *Crv2W, CagdBType OperationAdd) { CagdCrvStruct *CTmp1, *CTmp2, *CTmp3; if (Crv1X == NULL || Crv2X == NULL) return NULL; /* Make the two curves - same order and point type. */ Crv1X = CagdCrvCopy(Crv1X); Crv1W = CagdCrvCopy(Crv1W); Crv2X = CagdCrvCopy(Crv2X); Crv2W = CagdCrvCopy(Crv2W); if (!CagdMakeCrvsCompatible(&Crv1X, &Crv2X, FALSE, FALSE) || !CagdMakeCrvsCompatible(&Crv1W, &Crv2W, FALSE, FALSE) || !CagdMakeCrvsCompatible(&Crv1X, &Crv2W, FALSE, FALSE) || !CagdMakeCrvsCompatible(&Crv1W, &Crv2X, FALSE, FALSE)) SYMB_FATAL_ERROR(SYMB_ERR_CRV_FAIL_CMPT); CTmp1 = SymbCrvMult(Crv1X, Crv2W); CTmp2 = SymbCrvMult(Crv2X, Crv1W); CTmp3 = SymbCrvAddSubAux(CTmp1, CTmp2, OperationAdd); CagdCrvFree(CTmp1); CagdCrvFree(CTmp2); CagdCrvFree(Crv1X); CagdCrvFree(Crv1W); CagdCrvFree(Crv2X); CagdCrvFree(Crv2W); return CTmp3; } /***************************************************************************** * DESCRIPTION: * * Given two curve - add or subtract them as prescribed by OperationAdd, * * coordinatewise. * * The two curves are promoted to same type, point type, and order before * * addition can take place. * * Returned is a curve representing their sum or difference. * * * * PARAMETERS: * * Crv1, Crv2: Two curve to subtract coordinatewise. * * OperationAdd: TRUE of addition, FALSE for subtraction. * * * * RETURN VALUE: * * CagdCrvStruct *: The summation or difference of Crv1 and Crv2. * *****************************************************************************/ static CagdCrvStruct *SymbCrvAddSubAux(CagdCrvStruct *Crv1, CagdCrvStruct *Crv2, CagdBType OperationAdd) { CagdBType IsNotRational; int i, k; CagdCrvStruct *SumCrv; CagdRType **Points1, **Points2; /* Make the two curves - same order and point type. */ Crv1 = CagdCrvCopy(Crv1); Crv2 = CagdCrvCopy(Crv2); if (!CagdMakeCrvsCompatible(&Crv1, &Crv2, TRUE, TRUE)) SYMB_FATAL_ERROR(SYMB_ERR_CRV_FAIL_CMPT); SumCrv = CagdCrvNew(Crv1 -> GType, Crv1 -> PType, k = Crv1 -> Length); SumCrv -> Order = Crv1 -> Order; if (CAGD_IS_BSPLINE_CRV(SumCrv)) { SumCrv -> KnotVector = BspKnotCopy(Crv1 -> KnotVector, Crv1 -> Length + Crv1 -> Order); } IsNotRational = !CAGD_IS_RATIONAL_CRV(SumCrv); Points1 = Crv1 -> Points; Points2 = Crv2 -> Points; if (IsNotRational) { /* Simply add the respective control polygons. */ SymbMeshAddSub(SumCrv -> Points, Crv1 -> Points, Crv2 -> Points, SumCrv -> PType, SumCrv -> Length, OperationAdd); } else { /* Maybe the weights are identical, in which we can still add the */ /* the respective control polygons. */ for (i = 0; i < k; i++) if (!APX_EQ(Points1[0][i], Points2[0][i])) break; if (i < k) { CagdCrvStruct *Crv1W, *Crv1X, *Crv1Y, *Crv1Z, *Crv2W, *Crv2X, *Crv2Y, *Crv2Z, *SumCrvW, *SumCrvX, *SumCrvY, *SumCrvZ; /* Weights are different. Must use the addition of rationals */ /* rule ( we invoke SymbCrvMult here): */ /* */ /* x1 x2 x1 w2 +/- x2 w1 */ /* -- +/- -- = --------------- */ /* w1 w2 w1 w2 */ /* */ SymbCrvSplitScalar(Crv1, &Crv1W, &Crv1X, &Crv1Y, &Crv1Z); SymbCrvSplitScalar(Crv2, &Crv2W, &Crv2X, &Crv2Y, &Crv2Z); SumCrvW = SymbCrvMult(Crv1W, Crv2W); SumCrvX = SymbCrvRtnlMult(Crv1X, Crv1W, Crv2X, Crv2W, OperationAdd); SumCrvY = SymbCrvRtnlMult(Crv1Y, Crv1W, Crv2Y, Crv2W, OperationAdd); SumCrvZ = SymbCrvRtnlMult(Crv1Z, Crv1W, Crv2Z, Crv2W, OperationAdd); CagdCrvFree(Crv1W); CagdCrvFree(Crv1X); CagdCrvFree(Crv1Y); CagdCrvFree(Crv1Z); CagdCrvFree(Crv2W); CagdCrvFree(Crv2X); CagdCrvFree(Crv2Y); CagdCrvFree(Crv2Z); CagdCrvFree(SumCrv); if (!CagdMakeCrvsCompatible(&SumCrvW, &SumCrvX, TRUE, TRUE) || (SumCrvY != NULL && !CagdMakeCrvsCompatible(&SumCrvW, &SumCrvY, TRUE, TRUE)) || (SumCrvZ != NULL && !CagdMakeCrvsCompatible(&SumCrvW, &SumCrvZ, TRUE, TRUE))) SYMB_FATAL_ERROR(SYMB_ERR_CRV_FAIL_CMPT); SumCrv = SymbCrvMergeScalar(SumCrvW, SumCrvX, SumCrvY, SumCrvZ); CagdCrvFree(SumCrvW); CagdCrvFree(SumCrvX); CagdCrvFree(SumCrvY); CagdCrvFree(SumCrvZ); } else { /* Simply add respective control polygons (w is just copied). */ SymbMeshAddSub(SumCrv -> Points, Crv1 -> Points, Crv2 -> Points, SumCrv -> PType, SumCrv -> Length, OperationAdd); } } CagdCrvFree(Crv1); CagdCrvFree(Crv2); return SumCrv; } /***************************************************************************** * DESCRIPTION: M * Given a planar curve, compute its enclosed area field curve. M * This has little meaning unless Crv is closed, in which by evaluation M * the resulting area filed curve at the end points, the area enclosed by M * Crv can be computed. M * * * PARAMETERS: M * Crv: A curve to compute area filed curve for. M * * * RETURN VALUE: M * CagdCrvStruct *: The area field curve. M * * * KEYWORDS: M * SymbCrvEnclosedArea, area, symbolic computation M *****************************************************************************/ CagdCrvStruct *SymbCrvEnclosedArea(CagdCrvStruct *Crv) { CagdVType Trans; CagdCrvStruct *CrvX, *CrvY, *CrvZ, *CrvW, *DCrvX, *DCrvY, *CTmp1, *CTmp2, *CTmp3; switch (Crv -> GType) { case CAGD_CBEZIER_TYPE: case CAGD_CBSPLINE_TYPE: break; case CAGD_CPOWER_TYPE: SYMB_FATAL_ERROR(SYMB_ERR_POWER_NO_SUPPORT); return NULL; default: SYMB_FATAL_ERROR(SYMB_ERR_UNDEF_CRV); return NULL; } /* Using Green's theorem, if C(t1) = C(t2) the area enclosed is equal to */ /* */ /* t2 */ /* / */ /* 1 | */ /* - | -x'(t) y(t) + x(t) y'(t) dt */ /* 2 | */ /* / */ /* t1 */ /* */ SymbCrvSplitScalar(Crv, &CrvW, &CrvX, &CrvY, &CrvZ); if (CrvZ) CagdCrvFree(CrvZ); if (CrvW) { SYMB_FATAL_ERROR(SYMB_ERR_RATIONAL_NO_SUPPORT); CagdCrvFree(CrvW); } DCrvX = CagdCrvDerive(CrvX); DCrvY = CagdCrvDerive(CrvY); CTmp1 = SymbCrvMult(CrvX, DCrvY); CTmp2 = SymbCrvMult(DCrvX, CrvY); CagdCrvFree(CrvX); CagdCrvFree(CrvY); CagdCrvFree(DCrvX); CagdCrvFree(DCrvY); CTmp3 = SymbCrvSub(CTmp1, CTmp2); CagdCrvFree(CTmp1); CagdCrvFree(CTmp2); CTmp1 = CagdCrvIntegrate(CTmp3); CagdCrvFree(CTmp3); /* Scale by a half. */ Trans[0] = Trans[1] = Trans[2] = 0.0; CagdCrvTransform(CTmp1, Trans, 0.5); return CTmp1; } /***************************************************************************** * DESCRIPTION: M * Given two surfaces - add them coordinatewise. M * The two surfaces are promoted to same point type before the M * multiplication can take place. Furthermore, order and continuity are M * matched as well. M * * * PARAMETERS: M * Srf1, Srf2: Two surface to add up coordinatewise. M * * * RETURN VALUE: M * CagdSrfStruct *: The summation of Srf1 + Srf2 coordinatewise. M * * * SEE ALSO: M * SymbSrfSub, SymbMeshAddSub, SymbSrfMult M * * * KEYWORDS: M * SymbSrfAdd, addition, symbolic computation M *****************************************************************************/ CagdSrfStruct *SymbSrfAdd(CagdSrfStruct *Srf1, CagdSrfStruct *Srf2) { return SymbSrfAddSubAux(Srf1, Srf2, TRUE); } /***************************************************************************** * DESCRIPTION: M * Given two surfaces - subtract them coordinatewise. M * The two surfaces are promoted to same point type before the M * multiplication can take place. Furthermore, order and continuity are M * matched as well. M * * * PARAMETERS: M * Srf1, Srf2: Two surface to subtract coordinatewise. M * * * RETURN VALUE: M * CagdSrfStruct *: The difference of Srf1 - Srf2 coordinatewise. M * * * SEE ALSO: M * SymbSrfAdd, SymbMeshAddSub, SymbSrfMult M * * * KEYWORDS: M * SymbSrfSub, subtraction, symbolic computation M *****************************************************************************/ CagdSrfStruct *SymbSrfSub(CagdSrfStruct *Srf1, CagdSrfStruct *Srf2) { return SymbSrfAddSubAux(Srf1, Srf2, FALSE); } /***************************************************************************** * DESCRIPTION: M * Given two surfaces - multiply them coordinatewise. M * The two surfaces are promoted to same point type before the M * multiplication can take place. M * * * PARAMETERS: M * Srf1, Srf2: Two surface to multiply coordinatewise. M * * * RETURN VALUE: M * CagdSrfStruct *: The product of Srf1 * Srf2 coordinatewise. M * * * SEE ALSO: M * SymbSrfDotProd, SymbSrfVecDotProd, SymbSrfScalarScale, SymbSrfMultScalar,M * SymbSrfInvert M * * * KEYWORDS: M * SymbSrfMult, product, symbolic computation M *****************************************************************************/ CagdSrfStruct *SymbSrfMult(CagdSrfStruct *Srf1, CagdSrfStruct *Srf2) { CagdSrfStruct *ProdSrf = NULL; if (Srf1 -> GType == CAGD_SBEZIER_TYPE && Srf2 -> GType == CAGD_SBEZIER_TYPE) ProdSrf = BzrSrfMult(Srf1, Srf2); else if ((Srf1 -> GType == CAGD_SBEZIER_TYPE || Srf1 -> GType == CAGD_SBSPLINE_TYPE) && (Srf2 -> GType == CAGD_SBEZIER_TYPE || Srf2 -> GType == CAGD_SBSPLINE_TYPE)) ProdSrf = BspSrfMult(Srf1, Srf2); else SYMB_FATAL_ERROR(SYMB_ERR_UNDEF_SRF); return ProdSrf; } /***************************************************************************** * DESCRIPTION: M * Given a scalar surface, returns a scalar surface representing the M * reciprocal values, by making it rational (if was not one) and flipping M * the numerator and the denominator. M * * * PARAMETERS: M * Srf: A scalar surface to compute a reciprocal value for. M * * * RETURN VALUE: M * CagdSrfStruct *: A rational scalar surface that is equalto the M * reciprocal value of Srf. M * * * SEE ALSO: M * SymbSrfDotProd, SymbSrfVecDotProd, SymbSrfScalarScale, SymbSrfMultScalar,M * SymbSrfMult, SymbSrfCrossProd M * * * KEYWORDS: M * SymbSrfInvert, division, symbolic computation, reciprocal value M *****************************************************************************/ CagdSrfStruct *SymbSrfInvert(CagdSrfStruct *Srf) { int i; CagdRType *R, **Points; Srf = CagdSrfCopy(Srf); Points = Srf -> Points; switch (Srf -> PType) { case CAGD_PT_E1_TYPE: Points[0] = Points[1]; Points[1] = R = (CagdRType *) IritMalloc(sizeof(CagdRType) * Srf -> ULength * Srf -> VLength); for (i = 0; i < Srf -> ULength * Srf -> VLength; i++) *R++ = 1.0; Srf -> PType = CAGD_PT_P1_TYPE; break; case CAGD_PT_P1_TYPE: R = Points[0]; Points[0] = Points[1]; Points[1] = R; break; default: SYMB_FATAL_ERROR(SYMB_ERR_UNSUPPORT_PT); break; } return Srf; } /***************************************************************************** * DESCRIPTION: M * Given a surface, scale it by Scale. M * * * PARAMETERS: M * Srf: A surface to scale by magnitude Scale. M * Scale: Scaling factor. M * * * RETURN VALUE: M * CagdSrfStruct *: A surfaces scaled by Scale compared to Srf. M * * * SEE ALSO: M * SymbSrfDotProd, SymbSrfVecDotProd, SymbSrfMult, SymbSrfMultScalar, M * SymbSrfInvert, SymbSrfCrossProd M * * * KEYWORDS: M * SymbSrfScalarScale, scaling, symbolic computation M *****************************************************************************/ CagdSrfStruct *SymbSrfScalarScale(CagdSrfStruct *Srf, CagdRType Scale) { int i; CagdRType *R; Srf = CagdSrfCopy(Srf); switch (Srf -> PType) { case CAGD_PT_P1_TYPE: case CAGD_PT_E1_TYPE: for (i = 0, R = Srf -> Points[1]; i < Srf -> ULength * Srf -> VLength; i++) *R++ *= Scale; break; default: SYMB_FATAL_ERROR(SYMB_ERR_UNSUPPORT_PT); break; } return Srf; } /***************************************************************************** * DESCRIPTION: M * Given two surface - a vector curve Srf1 and a scalar curve Srf2, multiply M * all Srf1's coordinates by the scalar curve Srf2. M * Returned surface is a surface representing the product of the two given M * surfaces. M * * * PARAMETERS: M * Srf1, Srf2: Two surfaces to multiply. M * * * RETURN VALUE: M * CagdSrfStruct *: A surface representing the product of Srf1 and Srf2. M * * * SEE ALSO: M * SymbSrfDotProd, SymbSrfVecDotProd, SymbSrfMult, SymbSrfCrossProd, M * SymbCrvMultScalar M * * * KEYWORDS: M * SymbSrfMultScalar, product, symbolic computation M *****************************************************************************/ CagdSrfStruct *SymbSrfMultScalar(CagdSrfStruct *Srf1, CagdSrfStruct *Srf2) { CagdSrfStruct *PSrf1W, *PSrf1X, *PSrf1Y, *PSrf1Z, *PSrf2W, *PSrf2X, *PSrf2Y, *PSrf2Z, *TSrf, *ProdSrf; SymbSrfSplitScalar(Srf1, &PSrf1W, &PSrf1X, &PSrf1Y, &PSrf1Z); SymbSrfSplitScalar(Srf2, &PSrf2W, &PSrf2X, &PSrf2Y, &PSrf2Z); TSrf = SymbSrfMult(PSrf1X, PSrf2X); CagdSrfFree(PSrf1X); PSrf1X = TSrf; if (PSrf1Y != NULL) { TSrf = SymbSrfMult(PSrf1Y, PSrf2X); CagdSrfFree(PSrf1Y); PSrf1Y = TSrf; } if (PSrf1Z != NULL) { TSrf = SymbSrfMult(PSrf1Z, PSrf2X); CagdSrfFree(PSrf1Z); PSrf1Z = TSrf; } if (PSrf1W != NULL && PSrf2W != NULL) { TSrf = SymbSrfMult(PSrf1W, PSrf2W); CagdSrfFree(PSrf1W); PSrf1W = TSrf; } else if (PSrf2W != NULL) { PSrf1W = PSrf2W; PSrf2W = NULL; } ProdSrf = SymbSrfMergeScalar(PSrf1W, PSrf1X, PSrf1Y, PSrf1Z); CagdSrfFree(PSrf1W); CagdSrfFree(PSrf1X); CagdSrfFree(PSrf1Y); CagdSrfFree(PSrf1Z); CagdSrfFree(PSrf2W); CagdSrfFree(PSrf2X); CagdSrfFree(PSrf2Y); CagdSrfFree(PSrf2Z); return ProdSrf; } /***************************************************************************** * DESCRIPTION: M * Given two surfaces - computes their dot product. M * Returned surface is a scalar surface representing the dot product of the M * two given surfaces. M * * * PARAMETERS: M * Srf1, Srf2: Two surface to multiply and compute a dot product for. M * * * RETURN VALUE: M * CagdSrfStruct *: A scalar surface representing the dot product of M * Srf1 . Srf2. M * * * SEE ALSO: M * SymbSrfMult, SymbSrfVecDotProd, SymbSrfScalarScale, SymbSrfMultScalar, M * SymbSrfInvert, SymbSrfCrossProd M * * * KEYWORDS: M * SymbSrfDotProd, product, dot product, symbolic computation M *****************************************************************************/ CagdSrfStruct *SymbSrfDotProd(CagdSrfStruct *Srf1, CagdSrfStruct *Srf2) { CagdSrfStruct *PSrfW, *PSrfX, *PSrfY, *PSrfZ, *TSrf1, *TSrf2, *DotProdSrf, *ProdSrf = SymbSrfMult(Srf1, Srf2); SymbSrfSplitScalar(ProdSrf, &PSrfW, &PSrfX, &PSrfY, &PSrfZ); CagdSrfFree(ProdSrf); if (PSrfY != NULL) { TSrf1 = SymbSrfAdd(PSrfX, PSrfY); CagdSrfFree(PSrfX); CagdSrfFree(PSrfY); } else TSrf1 = PSrfX; if (PSrfZ != NULL) { TSrf2 = SymbSrfAdd(TSrf1, PSrfZ); CagdSrfFree(TSrf1); CagdSrfFree(PSrfZ); TSrf1 = TSrf2; } DotProdSrf = SymbSrfMergeScalar(PSrfW, TSrf1, NULL, NULL); CagdSrfFree(PSrfW); CagdSrfFree(TSrf1); return DotProdSrf; } /***************************************************************************** * DESCRIPTION: M * Given a surface and a vector - computes their dot product. M * Returned surface is a scalar surface representing the dot product. M * * * PARAMETERS: M * Srf: Curve to multiply and compute a dot product for. M * Vec: Vector to project Srf onto. M * * * RETURN VALUE: M * CagdSrfStruct *: A scalar surface representing the dot product of M * Srf . Vec. M * * * SEE ALSO: M * SymbSrfDotProd, SymbSrfMult, SymbSrfScalarScale, SymbSrfMultScalar, M * SymbSrfInvert, SymbSrfCrossProd M * * * KEYWORDS: M * SymbSrfVecDotProd, product, dot product, symbolic computation M *****************************************************************************/ CagdSrfStruct *SymbSrfVecDotProd(CagdSrfStruct *Srf, CagdVType Vec) { CagdSrfStruct *PSrfW, *PSrfX, *PSrfY, *PSrfZ, *TSrf, *DotProdSrf; SymbSrfSplitScalar(Srf, &PSrfW, &PSrfX, &PSrfY, &PSrfZ); TSrf = SymbSrfScalarScale(PSrfX, Vec[0]); CagdSrfFree(PSrfX); PSrfX = TSrf; if (PSrfY != NULL) { TSrf = SymbSrfScalarScale(PSrfY, Vec[1]); CagdSrfFree(PSrfY); PSrfY = TSrf; TSrf = SymbSrfAdd(PSrfX, PSrfY); CagdSrfFree(PSrfX); CagdSrfFree(PSrfY); PSrfX = TSrf; } if (PSrfZ != NULL) { TSrf = SymbSrfScalarScale(PSrfZ, Vec[2]); CagdSrfFree(PSrfZ); PSrfZ = TSrf; TSrf = SymbSrfAdd(PSrfX, PSrfZ); CagdSrfFree(PSrfX); CagdSrfFree(PSrfZ); PSrfX = TSrf; } DotProdSrf = SymbSrfMergeScalar(PSrfW, PSrfX, NULL, NULL); CagdSrfFree(PSrfW); CagdSrfFree(PSrfX); return DotProdSrf; } /***************************************************************************** * DESCRIPTION: M * Given two surfaces - computes their cross product. M * Returned surface is a scalar surface representing the cross product of M * the two given surfaces. M * * * PARAMETERS: M * Srf1, Srf2: Two surface to multiply and compute a cross product for. M * * * RETURN VALUE: M * CagdSrfStruct *: A scalar surface representing the cross product of M * Srf1 x Srf2. M * * * SEE ALSO: M * SymbSrfDotProd, SymbSrfVecDotProd, SymbSrfScalarScale, SymbSrfMultScalar,M * SymbSrfInvert M * * * KEYWORDS: M * SymbSrfCrossProd, product, cross product, symbolic computation M *****************************************************************************/ CagdSrfStruct *SymbSrfCrossProd(CagdSrfStruct *Srf1, CagdSrfStruct *Srf2) { CagdSrfStruct *Srf1W, *Srf1X, *Srf1Y, *Srf1Z, *Srf2W, *Srf2X, *Srf2Y, *Srf2Z, *TSrf1, *TSrf2, *CrossProdSrf, *PSrfW = NULL, *PSrfX = NULL, *PSrfY = NULL, *PSrfZ = NULL; SymbSrfSplitScalar(Srf1, &Srf1W, &Srf1X, &Srf1Y, &Srf1Z); SymbSrfSplitScalar(Srf2, &Srf2W, &Srf2X, &Srf2Y, &Srf2Z); if (Srf1X == NULL || Srf1Y == NULL || Srf2X == NULL || Srf2Y == NULL) SYMB_FATAL_ERROR(SYMB_ERR_NO_CROSS_PROD); /* Cross product X axis. */ TSrf1 = Srf2Z ? SymbSrfMult(Srf1Y, Srf2Z) : NULL; TSrf2 = Srf1Z ? SymbSrfMult(Srf2Y, Srf1Z) : NULL; if (TSrf1) { if (TSrf2) { PSrfX = SymbSrfSub(TSrf1, TSrf2); CagdSrfFree(TSrf2); } CagdSrfFree(TSrf1); } /* Cross product Y axis. */ TSrf1 = Srf1Z ? SymbSrfMult(Srf1Z, Srf2X) : NULL; TSrf2 = Srf2Z ? SymbSrfMult(Srf2Z, Srf1X) : NULL; if (TSrf1) { if (TSrf2) { PSrfY = SymbSrfSub(TSrf1, TSrf2); CagdSrfFree(TSrf2); } CagdSrfFree(TSrf1); } /* Cross product Z axis. */ TSrf1 = SymbSrfMult(Srf1X, Srf2Y); TSrf2 = SymbSrfMult(Srf2X, Srf1Y); PSrfZ = SymbSrfSub(TSrf1, TSrf2); CagdSrfFree(TSrf1); CagdSrfFree(TSrf2); /* Cross product W axis. */ if (Srf1W || Srf2W) { if (Srf1W == NULL) PSrfW = CagdSrfCopy(Srf2W); else if (Srf2W == NULL) PSrfW = CagdSrfCopy(Srf1W); else PSrfW = SymbSrfMult(Srf1W, Srf2W); } CagdSrfFree(Srf1X); CagdSrfFree(Srf1Y); CagdSrfFree(Srf1Z); CagdSrfFree(Srf1W); CagdSrfFree(Srf2X); CagdSrfFree(Srf2Y); CagdSrfFree(Srf2Z); CagdSrfFree(Srf2W); if (!CagdMakeSrfsCompatible(&PSrfW, &PSrfX, TRUE, TRUE, TRUE, TRUE) || !CagdMakeSrfsCompatible(&PSrfW, &PSrfY, TRUE, TRUE, TRUE, TRUE) || !CagdMakeSrfsCompatible(&PSrfW, &PSrfZ, TRUE, TRUE, TRUE, TRUE)) SYMB_FATAL_ERROR(SYMB_ERR_SRF_FAIL_CMPT); CrossProdSrf = SymbSrfMergeScalar(PSrfW, PSrfX, PSrfY, PSrfZ); CagdSrfFree(PSrfX); CagdSrfFree(PSrfY); CagdSrfFree(PSrfZ); CagdSrfFree(PSrfW); return CrossProdSrf; } /***************************************************************************** * DESCRIPTION: M * Given two surfaces - multiply them using the quotient product rule: M * X = X1 W2 +/- X2 W1 V * All provided surfaces are assumed to be non rational scalar surfaces. M * Returned is a non rational scalar surface (CAGD_PT_E1_TYPE). M * * * PARAMETERS: M * Srf1X: Numerator of first surface. M * Srf1W: Denominator of first surface. M * Srf2X: Numerator of second surface. M * Srf2W: Denominator of second surface. M * OperationAdd: TRUE for addition, FALSE for subtraction. M * * * RETURN VALUE: M * CagdSrfStruct *: The result of Srf1X Srf2W +/- Srf2X Srf1W. M * * * SEE ALSO: M * SymbSrfDotProd, SymbSrfVecDotProd, SymbSrfScalarScale, SymbSrfMultScalar,M * SymbSrfInvert M * * * KEYWORDS: M * SymbSrfRtnlMult, product, symbolic computation M *****************************************************************************/ CagdSrfStruct *SymbSrfRtnlMult(CagdSrfStruct *Srf1X, CagdSrfStruct *Srf1W, CagdSrfStruct *Srf2X, CagdSrfStruct *Srf2W, CagdBType OperationAdd) { CagdSrfStruct *CTmp1, *CTmp2, *CTmp3; /* Make the two surfaces - same order and point type. */ Srf1X = CagdSrfCopy(Srf1X); Srf1W = CagdSrfCopy(Srf1W); Srf2X = CagdSrfCopy(Srf2X); Srf2W = CagdSrfCopy(Srf2W); if (!CagdMakeSrfsCompatible(&Srf1X, &Srf2X, FALSE, FALSE, FALSE, FALSE) || !CagdMakeSrfsCompatible(&Srf1W, &Srf2W, FALSE, FALSE, FALSE, FALSE) || !CagdMakeSrfsCompatible(&Srf1X, &Srf2W, FALSE, FALSE, FALSE, FALSE) || !CagdMakeSrfsCompatible(&Srf1W, &Srf2X, FALSE, FALSE, FALSE, FALSE)) SYMB_FATAL_ERROR(SYMB_ERR_SRF_FAIL_CMPT); CTmp1 = SymbSrfMult(Srf1X, Srf2W); CTmp2 = SymbSrfMult(Srf2X, Srf1W); CTmp3 = SymbSrfAddSubAux(CTmp1, CTmp2, OperationAdd); CagdSrfFree(CTmp1); CagdSrfFree(CTmp2); CagdSrfFree(Srf1X); CagdSrfFree(Srf1W); CagdSrfFree(Srf2X); CagdSrfFree(Srf2W); return CTmp3; } /***************************************************************************** * DESCRIPTION: M * Given a surface - compute its unnormalized norrmal vectorfield surface, as M * the cross product if this partial derivatives. M * * * PARAMETERS: M * Srf: To compute an unnormalized normal vector field for. M * * * RETURN VALUE: M * CagdSrfStruct *: A vector field representing the unnormalized normal M * vector field of Srf. M * * * KEYWORDS: M * SymbSrfNormalSrf, normal, symbolic computation M *****************************************************************************/ CagdSrfStruct *SymbSrfNormalSrf(CagdSrfStruct *Srf) { CagdSrfStruct *SrfDU = CagdSrfDerive(Srf, CAGD_CONST_U_DIR), *SrfDV = CagdSrfDerive(Srf, CAGD_CONST_V_DIR), *NormalSrf = SymbSrfCrossProd(SrfDV, SrfDU); CagdSrfFree(SrfDU); CagdSrfFree(SrfDV); return NormalSrf; } /***************************************************************************** * DESCRIPTION: * * Given two surface - add or subtract them as prescribed by OperationAdd, * * coordinatewise. * * The two surfaces are promoted to same type, point type, and order before * * addition can take place. * * Returned is a surface representing their sum or difference. * * * * PARAMETERS: * * Srf1, Srf2: Two surface to subtract coordinatewise. * * OperationAdd: TRUE of addition, FALSE for subtraction. * * * * RETURN VALUE: * * CagdSrfStruct *: The summation or difference of Srf1 and Srf2. * *****************************************************************************/ static CagdSrfStruct *SymbSrfAddSubAux(CagdSrfStruct *Srf1, CagdSrfStruct *Srf2, CagdBType OperationAdd) { CagdBType IsNotRational; int i, Len; CagdSrfStruct *SumSrf; CagdRType **Points1, **Points2; /* Make the two surfaces - same order and point type. */ Srf1 = CagdSrfCopy(Srf1); Srf2 = CagdSrfCopy(Srf2); if (!CagdMakeSrfsCompatible(&Srf1, &Srf2, TRUE, TRUE, TRUE, TRUE)) SYMB_FATAL_ERROR(SYMB_ERR_SRF_FAIL_CMPT); Len = Srf1 -> ULength * Srf1 -> VLength; SumSrf = CagdSrfNew(Srf1 -> GType, Srf1 -> PType, Srf1 -> ULength, Srf1 ->VLength); SumSrf -> UOrder = Srf1 -> UOrder; SumSrf -> VOrder = Srf1 -> VOrder; if (CAGD_IS_BSPLINE_SRF(SumSrf)) { SumSrf -> UKnotVector = BspKnotCopy(Srf1 -> UKnotVector, Srf1 -> ULength + Srf1 -> UOrder); SumSrf -> VKnotVector = BspKnotCopy(Srf1 -> VKnotVector, Srf1 -> VLength + Srf1 -> VOrder); } IsNotRational = !CAGD_IS_RATIONAL_SRF(SumSrf); Points1 = Srf1 -> Points; Points2 = Srf2 -> Points; if (IsNotRational) { /* Simply add the respective control polygons. */ SymbMeshAddSub(SumSrf -> Points, Srf1 -> Points, Srf2 -> Points, SumSrf -> PType, Len, OperationAdd); } else { /* Maybe the weights are identical, in which we can still add the */ /* the respective control polygons. */ for (i = 0; i < Len; i++) if (!APX_EQ(Points1[0][i], Points2[0][i])) break; if (i < Len) { CagdSrfStruct *Srf1W, *Srf1X, *Srf1Y, *Srf1Z, *Srf2W, *Srf2X, *Srf2Y, *Srf2Z, *SumSrfW, *SumSrfX, *SumSrfY, *SumSrfZ; /* Weights are different. Must use the addition of rationals */ /* rule ( we invoke SymbSrfMult here): */ /* */ /* x1 x2 x1 w2 +/- x2 w1 */ /* -- +/- -- = --------------- */ /* w1 w2 w1 w2 */ /* */ SymbSrfSplitScalar(Srf1, &Srf1W, &Srf1X, &Srf1Y, &Srf1Z); SymbSrfSplitScalar(Srf2, &Srf2W, &Srf2X, &Srf2Y, &Srf2Z); SumSrfW = SymbSrfMult(Srf1W, Srf2W); SumSrfX = SymbSrfRtnlMult(Srf1X, Srf1W, Srf2X, Srf2W, OperationAdd); SumSrfY = SymbSrfRtnlMult(Srf1Y, Srf1W, Srf2Y, Srf2W, OperationAdd); SumSrfZ = SymbSrfRtnlMult(Srf1Z, Srf1W, Srf2Z, Srf2W, OperationAdd); CagdSrfFree(Srf1W); CagdSrfFree(Srf1X); CagdSrfFree(Srf1Y); CagdSrfFree(Srf1Z); CagdSrfFree(Srf2W); CagdSrfFree(Srf2X); CagdSrfFree(Srf2Y); CagdSrfFree(Srf2Z); CagdSrfFree(SumSrf); if (!CagdMakeSrfsCompatible(&SumSrfW, &SumSrfX, TRUE, TRUE, TRUE, TRUE) || (SumSrfY != NULL && !CagdMakeSrfsCompatible(&SumSrfW, &SumSrfY, TRUE, TRUE, TRUE, TRUE)) || (SumSrfZ != NULL && !CagdMakeSrfsCompatible(&SumSrfW, &SumSrfZ, TRUE, TRUE, TRUE, TRUE))) SYMB_FATAL_ERROR(SYMB_ERR_SRF_FAIL_CMPT); SumSrf = SymbSrfMergeScalar(SumSrfW, SumSrfX, SumSrfY, SumSrfZ); CagdSrfFree(SumSrfW); CagdSrfFree(SumSrfX); CagdSrfFree(SumSrfY); CagdSrfFree(SumSrfZ); } else { /* Simply add respective control polygons (w is just copied). */ SymbMeshAddSub(SumSrf -> Points, Srf1 -> Points, Srf2 -> Points, SumSrf -> PType, Len, OperationAdd); } } CagdSrfFree(Srf1); CagdSrfFree(Srf2); return SumSrf; } /****************************************************************************** * DESCRIPTION: M * Given two control polygons/meshes - add them coordinatewise. M * If mesh is rational, weights are assumed identical and are just copied. M * * * PARAMETERS: M * DestPoints: Where addition or difference result should go to. M * Points1: First control polygon/mesh. M * Points2: Second control polygon/mesh. M * PType: Type of points we are dealing with. M * Size: Length of each vector in Points1/2. M * OperationAdd: TRUE of addition, FALSE for subtraction. M * * * RETURN VALUE: M * void M * * * SEE ALSO: M * SymbSrfSub, SymbSrfAdd M * * * KEYWORDS: M * SymbMeshAddSub, addition, subtraction, symbolic computation M *****************************************************************************/ void SymbMeshAddSub(CagdRType **DestPoints, CagdRType **Points1, CagdRType **Points2, CagdPointType PType, int Size, CagdBType OperationAdd) { CagdBType IsRational = CAGD_IS_RATIONAL_PT(PType); int i, j, NumCoords = CAGD_NUM_OF_PT_COORD(PType); for (i = 1; i <= NumCoords; i++) { CagdRType *DPts = DestPoints[i], *Pts1 = Points1[i], *Pts2 = Points2[i]; for (j = 0; j < Size; j++) *DPts++ = OperationAdd ? *Pts1++ + *Pts2++ : *Pts1++ - *Pts2++; } if (IsRational) { /* Copy the weights (should be identical in both). */ CagdRType *DPts = DestPoints[0], *Pts1 = Points1[0], *Pts2 = Points2[0]; for (j = 0; j < Size; j++) { if (!APX_EQ(*Pts1, *Pts2)) SYMB_FATAL_ERROR(SYMB_ERR_W_NOT_SAME); *DPts++ = *Pts1++; Pts2++; } } } /***************************************************************************** * DESCRIPTION: M * Given a curve splits it to its scalar component curves. Ignores all M * dimensions beyond the third, Z, dimension. M * * * PARAMETERS: M * Crv: Curve to split. M * CrvW: The weight component of Crv, if have any. M * CrvX: The X component of Crv. M * CrvY: The Y component of Crv, if have any. M * CrvZ: The Z component of Crv, if have any. M * * * RETURN VALUE: M * void M * * * SEE ALSO: M * SymbSrfSplitScalar, SymbCrvMergeScalar M * * * KEYWORDS: M * SymbCrvSplitScalar, split, symbolic computation M *****************************************************************************/ void SymbCrvSplitScalar(CagdCrvStruct *Crv, CagdCrvStruct **CrvW, CagdCrvStruct **CrvX, CagdCrvStruct **CrvY, CagdCrvStruct **CrvZ) { CagdBType IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv); int i, Length = Crv -> Length, NumCoords = CAGD_NUM_OF_PT_COORD(Crv -> PType); CagdCrvStruct *Crvs[CAGD_MAX_PT_SIZE]; for (i = 0; i < CAGD_MAX_PT_SIZE; i++) Crvs[i] = NULL; for (i = IsNotRational; i <= NumCoords; i++) { Crvs[i] = CagdCrvNew(Crv -> GType, CAGD_PT_E1_TYPE, Length); Crvs[i] -> Order = Crv -> Order; if (Crv -> KnotVector != NULL) Crvs[i] -> KnotVector = BspKnotCopy(Crv -> KnotVector, Crv -> Length + Crv -> Order); SYMB_GEN_COPY(Crvs[i] -> Points[1], Crv -> Points[i], sizeof(CagdRType) * Length); } *CrvW = Crvs[0]; *CrvX = Crvs[1]; *CrvY = Crvs[2]; *CrvZ = Crvs[3]; } /***************************************************************************** * DESCRIPTION: M * Given a set of scalar curves, treat them as coordinates into a new curve. M * Assumes at least CrvX is not NULL in which a scalar curve is returned. M * Assumes CrvX/Y/Z/W are either E1 or P1 in which the weights are assumed M * to be identical and can be ignored if CrvW exists or copied otheriwse. M * * * PARAMETERS: M * CrvW: The weight component of new constructed curve, if have any. M * CrvX: The X component of new constructed curve. M * CrvY: The Y component of new constructed curve, if have any. M * CrvZ: The Z component of new constructed curve, if have any. M * * * RETURN VALUE: M * CagdCrvStruct *: A new curve constructed from given scalar curves. M * * * SEE ALSO: M * SymbSrfMergeScalar, SymbCrvSplitScalar M * * * KEYWORDS: M * SymbCrvMergeScalar, merge, symbolic computation M *****************************************************************************/ CagdCrvStruct *SymbCrvMergeScalar(CagdCrvStruct *CrvW, CagdCrvStruct *CrvX, CagdCrvStruct *CrvY, CagdCrvStruct *CrvZ) { CagdBType WeightCopied = FALSE, IsRational = CrvW != NULL; int i, Length, NumCoords = (CrvX != NULL) + (CrvY != NULL) + (CrvZ != NULL); CagdPointType PType = CAGD_MAKE_PT_TYPE(IsRational, NumCoords); CagdCrvStruct *Crvs[CAGD_MAX_PT_SIZE], *Crv; Crvs[0] = CrvW == NULL ? NULL : CagdCrvCopy(CrvW); Crvs[1] = CrvX == NULL ? NULL : CagdCrvCopy(CrvX); Crvs[2] = CrvY == NULL ? NULL : CagdCrvCopy(CrvY); Crvs[3] = CrvZ == NULL ? NULL : CagdCrvCopy(CrvZ); for (i = 0; i < 4; i++) { int j; for (j = i + 1; j < 4; j++) if (Crvs[i] != NULL && Crvs[j] != NULL) CagdMakeCrvsCompatible(&Crvs[i], &Crvs[j], TRUE, TRUE); } Length = CrvX -> Length; Crv = CagdCrvNew(CrvX -> GType, PType, Length); Crv -> Order = CrvX -> Order; if (CrvX -> KnotVector != NULL) Crv -> KnotVector = BspKnotCopy(CrvX -> KnotVector, Length + CrvX -> Order); for (i = !IsRational; i <= NumCoords; i++) { if (Crvs[i] != NULL) { if (Crvs[i] -> PType != CAGD_PT_E1_TYPE) { if (Crvs[i] -> PType != CAGD_PT_P1_TYPE) SYMB_FATAL_ERROR(SYMB_ERR_SCALAR_EXPECTED); else if (CrvW == NULL && WeightCopied == FALSE) { SYMB_GEN_COPY(Crv -> Points[0], Crvs[i] -> Points[0], sizeof(CagdRType) * Length); WeightCopied = TRUE; } } SYMB_GEN_COPY(Crv -> Points[i], Crvs[i] -> Points[1], sizeof(CagdRType) * Length); } } for (i = 0; i < 4; i++) CagdCrvFree(Crvs[i]); return Crv; } /***************************************************************************** * DESCRIPTION: M * Given a surface splits it to its scalar component surfaces. Ignores all M * dimensions beyond the third, Z, dimension. M * * * PARAMETERS: M * Srf: Surface to split. M * SrfW: The weight component of Srf, if have any. M * SrfX: The X component of Srf. M * SrfY: The Y component of Srf, if have any. M * SrfZ: The Z component of Srf, if have any. M * * * RETURN VALUE: M * void M * * * SEE ALSO: M * SymbSrfMergeScalar, SymbCrvSplitScalar M * * * KEYWORDS: M * SymbSrfSplitScalar, split, symbolic computation M *****************************************************************************/ void SymbSrfSplitScalar(CagdSrfStruct *Srf, CagdSrfStruct **SrfW, CagdSrfStruct **SrfX, CagdSrfStruct **SrfY, CagdSrfStruct **SrfZ) { CagdBType IsNotRational = !CAGD_IS_RATIONAL_SRF(Srf); int i, ULength = Srf -> ULength, VLength = Srf -> VLength, NumCoords = CAGD_NUM_OF_PT_COORD(Srf -> PType); CagdSrfStruct *Srfs[CAGD_MAX_PT_SIZE]; for (i = 0; i < CAGD_MAX_PT_SIZE; i++) Srfs[i] = NULL; for (i = IsNotRational; i <= NumCoords; i++) { Srfs[i] = CagdSrfNew(Srf -> GType, CAGD_PT_E1_TYPE, ULength, VLength); Srfs[i] -> UOrder = Srf -> UOrder; Srfs[i] -> VOrder = Srf -> VOrder; if (Srf -> UKnotVector != NULL) Srfs[i] -> UKnotVector = BspKnotCopy(Srf -> UKnotVector, ULength + Srf -> UOrder); if (Srf -> VKnotVector != NULL) Srfs[i] -> VKnotVector = BspKnotCopy(Srf -> VKnotVector, VLength + Srf -> VOrder); SYMB_GEN_COPY(Srfs[i] -> Points[1], Srf -> Points[i], sizeof(CagdRType) * Srf -> ULength * Srf -> VLength); } *SrfW = Srfs[0]; *SrfX = Srfs[1]; *SrfY = Srfs[2]; *SrfZ = Srfs[3]; } /***************************************************************************** * DESCRIPTION: M * Given a set of scalar surfaces, treat them as coordinates into a new M * surface. M * Assumes at least SrfX is not NULL in which a scalar surface is returned. M * Assumes SrfX/Y/Z/W are either E1 or P1 in which the weights are assumed M * to be identical and can be ignored if SrfW exists or copied otheriwse. M * * * PARAMETERS: M * SrfW: The weight component of new constructed surface, if have any. M * SrfX: The X component of new constructed surface. M * SrfY: The Y component of new constructed surface, if have any. M * SrfZ: The Z component of new constructed surface, if have any. M * * * RETURN VALUE: M * CagdSrfStruct *: A new surface constructed from given scalar surfaces. M * * * SEE ALSO: M * SymbSrfSplitScalar, SymbCrvMergeScalar M * * * KEYWORDS: M * SymbSrfMergeScalar, merge, symbolic computation M *****************************************************************************/ CagdSrfStruct *SymbSrfMergeScalar(CagdSrfStruct *SrfW, CagdSrfStruct *SrfX, CagdSrfStruct *SrfY, CagdSrfStruct *SrfZ) { CagdBType WeightCopied = FALSE, IsRational = SrfW != NULL; int i, ULength, VLength, NumCoords = (SrfX != NULL) + (SrfY != NULL) + (SrfZ != NULL); CagdPointType PType = CAGD_MAKE_PT_TYPE(IsRational, NumCoords); CagdSrfStruct *Srfs[CAGD_MAX_PT_SIZE], *Srf; Srfs[0] = SrfW == NULL ? NULL : CagdSrfCopy(SrfW); Srfs[1] = SrfX == NULL ? NULL : CagdSrfCopy(SrfX); Srfs[2] = SrfY == NULL ? NULL : CagdSrfCopy(SrfY); Srfs[3] = SrfZ == NULL ? NULL : CagdSrfCopy(SrfZ); for (i = 0; i < 4; i++) { int j; for (j = i + 1; j < 4; j++) if (Srfs[i] != NULL && Srfs[j] != NULL) CagdMakeSrfsCompatible(&Srfs[i], &Srfs[j], TRUE, TRUE, TRUE, TRUE); } ULength = SrfX -> ULength; VLength = SrfX -> VLength; Srf = CagdSrfNew(SrfX -> GType, PType, ULength, VLength); Srf -> UOrder = SrfX -> UOrder; Srf -> VOrder = SrfX -> VOrder; if (SrfX -> UKnotVector != NULL) Srf -> UKnotVector = BspKnotCopy(SrfX -> UKnotVector, ULength + SrfX -> UOrder); if (SrfX -> VKnotVector != NULL) Srf -> VKnotVector = BspKnotCopy(SrfX -> VKnotVector, VLength + SrfX -> VOrder); for (i = !IsRational; i <= NumCoords; i++) { if (Srfs[i] != NULL) { if (Srfs[i] -> PType != CAGD_PT_E1_TYPE) { if (Srfs[i] -> PType != CAGD_PT_P1_TYPE) SYMB_FATAL_ERROR(SYMB_ERR_SCALAR_EXPECTED); else if (SrfW == NULL && WeightCopied == FALSE) { SYMB_GEN_COPY(Srf -> Points[0], Srfs[i] -> Points[0], sizeof(CagdRType) * ULength * VLength); WeightCopied = TRUE; } } SYMB_GEN_COPY(Srf -> Points[i], Srfs[i] -> Points[1], sizeof(CagdRType) * ULength * VLength); } } for (i = 0; i < 4; i++) CagdSrfFree(Srfs[i]); return Srf; } /***************************************************************************** * DESCRIPTION: M * Promote a scalar curve to two dimensions by moving the scalar axis to be M * the Y axis and adding monotone X axis. M * * * PARAMETERS: M * Crv: Scalar curve to promose to a two dimensional one. M * Min: Minimum of new monotone X axis. M * Max: Maximum of new monotone X axis. M * * * RETURN VALUE: M * CagdCrvStruct *: A two dimensional curve. M * * * SEE ALSO: M * SymbPrmtSclrSrfTo3D M * * * KEYWORDS: M * SymbPrmtSclrCrvTo2D, promotion, conversion, symbolic computation M *****************************************************************************/ CagdCrvStruct *SymbPrmtSclrCrvTo2D(CagdCrvStruct *Crv, CagdRType Min, CagdRType Max) { int i, Length = Crv -> Length; CagdBType IsRational = CAGD_IS_RATIONAL_CRV(Crv); CagdRType *R, *Points, *WPoints, Step = (Max - Min) / (Length - 1); CagdCrvStruct *PrmtCrv = CagdCoerceCrvTo(Crv, IsRational ? CAGD_PT_P2_TYPE : CAGD_PT_E2_TYPE); R = PrmtCrv -> Points[1]; PrmtCrv -> Points[1] = PrmtCrv -> Points[2]; PrmtCrv -> Points[2] = R; Points = PrmtCrv -> Points[1]; WPoints = IsRational ? PrmtCrv -> Points[0] : NULL; for (i = 0; i < Length; i++) *Points++ = (Min + Step * i) * (IsRational ? *WPoints++ : 1.0); return PrmtCrv; } /***************************************************************************** * DESCRIPTION: M * Promote a scalar surface to three dimensions by moving the scalar axis to M * be the Z axis and adding monotone X and Y axes. M * * * PARAMETERS: M * Srf: M * UMin: Minimum of new monotone X axis. M * UMax: Maximum of new monotone X axis. M * VMin: Minimum of new monotone Y axis. M * VMax: Maximum of new monotone X axis. M * * * RETURN VALUE: M * CagdSrfStruct *: A three dimensional surface. M * * * SEE ALSO: M * SymbPrmtSclrCrvTo2D M * * * KEYWORDS: M * SymbPrmtSclrSrfTo3D, promotion, conversion, symbolic computation M *****************************************************************************/ CagdSrfStruct *SymbPrmtSclrSrfTo3D(CagdSrfStruct *Srf, CagdRType UMin, CagdRType UMax, CagdRType VMin, CagdRType VMax) { int i, j, ULength = Srf -> ULength, VLength = Srf -> VLength; CagdBType IsRational = CAGD_IS_RATIONAL_SRF(Srf); CagdRType *R, *Points, *WPoints, UStep = (UMax - UMin) / (ULength - 1), VStep = (VMax - VMin) / (VLength - 1); CagdSrfStruct *PrmtSrf = CagdCoerceSrfTo(Srf, IsRational ? CAGD_PT_P3_TYPE : CAGD_PT_E3_TYPE); R = PrmtSrf -> Points[1]; PrmtSrf -> Points[1] = PrmtSrf -> Points[3]; PrmtSrf -> Points[3] = R; Points = PrmtSrf -> Points[1]; WPoints = IsRational ? PrmtSrf -> Points[0] : NULL; for (j = 0; j < VLength; j++) for (i = 0; i < ULength; i++) *Points++ = (UMin + UStep * i) * (IsRational ? *WPoints++ : 1.0); Points = PrmtSrf -> Points[2]; WPoints = IsRational ? PrmtSrf -> Points[0] : NULL; for (j = 0; j < VLength; j++) for (i = 0; i < ULength; i++) *Points++ = (VMin + VStep * j) * (IsRational ? *WPoints++ : 1.0); return PrmtSrf; } /***************************************************************************** * DESCRIPTION: M * Given a control polygon/mesh, computes the extremum values of them all. M * * * PARAMETERS: M * Points: To scan for extremum values. M * Length: Length of each vector in Points. M * FindMinimum: TRUE for minimum, FALSE for maximum. M * * * RETURN VALUE: M * CagdRType *: A vector holding PType point with the extremum values of M * each axis independently. M * * * KEYWORDS: M * SymbExtremumCntPtVals, extremum values, symbolic computation M *****************************************************************************/ CagdRType *SymbExtremumCntPtVals(CagdRType **Points, int Length, CagdBType FindMinimum) { static CagdRType Extremum[CAGD_MAX_PT_SIZE]; int i, j; for (i = 1; Points[i] != NULL && i <= CAGD_MAX_PT_COORD; i++) { CagdRType *R = Points[0], *P = Points[i]; Extremum[i] = FindMinimum ? IRIT_INFNTY : -IRIT_INFNTY; for (j = 0; j < Length; j++) { CagdRType Val = R == NULL ? *P++ : *P++ / *R++; if (FindMinimum) { if (Extremum[i] > Val) Extremum[i] = Val; } else { if (Extremum[i] < Val) Extremum[i] = Val; } } } return Extremum; }