/****************************************************************************** * Mvar_Der.c - Compute derivatives of multi-variates. * ******************************************************************************* * (C) Gershon Elber, Technion, Israel Institute of Technology * ******************************************************************************* * Written by Gershon Elber, May. 97. * ******************************************************************************/ #include #include #include #include "mvar_loc.h" /***************************************************************************** * DESCRIPTION: M * Given a multi-variate, computes its partial derivative multi-variate in M * direction Dir. M * * * PARAMETERS: M * MV: Multi-Variate to differentiate. M * Dir: Direction of differentiation. M * * * RETURN VALUE: M * MvarMVStruct *: Differentiated multi-variate in direction Dir. M * * * KEYWORDS: M * MvarMVDerive, multi-variates M *****************************************************************************/ MvarMVStruct *MvarMVDerive(MvarMVStruct *MV, MvarMVDirType Dir) { switch (MV -> GType) { case MVAR_BEZIER_TYPE: return MvarBzrMVDerive(MV, Dir); case MVAR_BSPLINE_TYPE: return MvarBspMVDerive(MV, Dir); default: MVAR_FATAL_ERROR(MVAR_ERR_UNDEF_GEOM); return NULL; } } /***************************************************************************** * DESCRIPTION: M * Given a Bezier multi-variate, computes its partial derivative M * multi-variate in direction Dir. M * Let old control polygon be P(i), i = 0 to k-1, and Q(i) be new one. M * Then: M * Q(i) = (k - 1) * (P(i+1) - P(i)), i = 0 to k-2. V * This function computes the derivative of a rational function component- M * wise with out taking into consideration the quotient rule. M * * * PARAMETERS: M * MV: Multi-Variate to differentiate. M * Dir: Direction of differentiation. M * * * RETURN VALUE: M * MvarMVStruct *: Differentiated multi-variate in direction Dir. A M * Bezier multi-variate. M * * * KEYWORDS: M * MvarBzrMVDerive, multi-variates M *****************************************************************************/ MvarMVStruct *MvarBzrMVDerive(MvarMVStruct *MV, MvarMVDirType Dir) { CagdBType IsNotRational = !MVAR_IS_RATIONAL_MV(MV); int l, DIndex, Length = MV -> Lengths[Dir], *Indices = (int *) IritMalloc(sizeof(int) * MV -> Dim), MaxCoord = CAGD_NUM_OF_PT_COORD(MV -> PType); CagdRType **DPoints; MvarMVStruct *DerivedMV = NULL; ZAP_MEM(Indices, sizeof(int) * MV -> Dim); if (Length > 1) MV -> Lengths[Dir]--; DerivedMV = MvarBzrMVNew(MV -> Dim, MV -> Lengths, MV -> PType); if (Length > 1) MV -> Lengths[Dir]++; DPoints = DerivedMV -> Points; DIndex = 0; do { int Index = MvarGetPointsMeshIndices(MV, Indices), NextIndex = Index + MVAR_NEXT_DIM(MV, Dir); for (l = IsNotRational; l <= MaxCoord; l++) DPoints[l][DIndex] = Length < 2 ? 0.0 : (Length - 1) * (MV -> Points[l][NextIndex] - MV -> Points[l][Index]); } while (MvarIncrementMeshIndices2(DerivedMV, Indices, &DIndex)); IritFree(Indices); return DerivedMV; } /***************************************************************************** * DESCRIPTION: M * Given a Bspline multi-variate, computes its partial derivative M * multi-variate in direction Dir. M * Let old control polygon be P(i), i = 0 to k-1, and Q(i) be new one. M * Then: M * Q(i) = (k - 1) * (P(i+1) - P(i)) / (Kv(i + k) - Kv(i + 1)), i = 0 to k-2. V * This function computes the derivative of a rational function component- M * wise with out taking into consideration the quotient rule. M * * * PARAMETERS: M * MV: Multi-Variate to differentiate. M * Dir: Direction of differentiation. M * * * RETURN VALUE: M * MvarMVStruct *: Differentiated multi-variate in direction Dir. A M * Bspline multi-variate. M * * * KEYWORDS: M * MvarBspMVDerive, multi-variates M *****************************************************************************/ MvarMVStruct *MvarBspMVDerive(MvarMVStruct *MV, MvarMVDirType Dir) { CagdBType IsNotRational = !MVAR_IS_RATIONAL_MV(MV); int i, DIndex, NewLength, NewOrder, Length = MV -> Lengths[Dir], Order = MV -> Orders[Dir], *Indices = (int *) IritMalloc(sizeof(int) * MV -> Dim), MaxCoord = CAGD_NUM_OF_PT_COORD(MV -> PType); CagdRType **DPoints, *KV = MV -> KnotVectors[Dir], **Points = MV -> Points; MvarMVStruct *DerivedMV = NULL; ZAP_MEM(Indices, sizeof(int) * MV -> Dim); NewLength = Order < 2 ? Length : Length - 1; NewOrder = MAX(Order - 1, 1); MV -> Lengths[Dir] = NewLength; MV -> Orders[Dir] = NewOrder; DerivedMV = MvarBspMVNew(MV -> Dim, MV -> Lengths, MV -> Orders, MV -> PType); MV -> Lengths[Dir] = Length; MV -> Orders[Dir] = Order; for (i = 0; i < MV -> Dim; i++) { if (i == Dir) CAGD_GEN_COPY(DerivedMV -> KnotVectors[i], &MV -> KnotVectors[i][Order < 2 ? 0 : 1], sizeof(CagdRType) * (NewLength + NewOrder)); else CAGD_GEN_COPY(DerivedMV -> KnotVectors[i], MV -> KnotVectors[i], sizeof(CagdRType) * (MV -> Lengths[i] + MV -> Orders[i])); } DPoints = DerivedMV -> Points; DIndex = 0; do { int l, Index = MvarGetPointsMeshIndices(MV, Indices), NextIndex = Index + MVAR_NEXT_DIM(MV, Dir); CagdRType Denom = KV[Indices[Dir] + Order] - KV[Indices[Dir] + 1]; if (APX_EQ(Denom, 0.0)) Denom = IRIT_UEPS; for (l = IsNotRational; l <= MaxCoord; l++) { DPoints[l][DIndex] = Order < 2 ? 0.0 : (Order - 1) * (Points[l][NextIndex] - Points[l][Index]) / Denom; } } while (MvarIncrementMeshIndices2(DerivedMV, Indices, &DIndex)); IritFree(Indices); return DerivedMV; } /***************************************************************************** * DESCRIPTION: M * Builds a gradient for the given scalar multivariate. M * If the input is rational, returned is a dynamically allocated vector of M * scalar mulitvariate functions each representing Dm/Dui, i from 1 to Dim. M * The returned partial derivative are differentiated directly without the M * quotient rule which must be applied manually. M * Otherwise, if the input is polynomial, the gradient is returned as one M * vector function. M * * * PARAMETERS: M * MV: Input scalar field to compute its gradient function. M * * * RETURN VALUE: M * MvarMVGradientStruct *: The gradient function of the input scalar field. M * * * SEE ALSO: M * MvarMVDerive, MvarMVFreeGradient, MvarMVEvalGradient M * * * KEYWORDS: M * MvarMVPrepGradient M *****************************************************************************/ MvarMVGradientStruct *MvarMVPrepGradient(MvarMVStruct *MV) { CagdBType IsRational = MVAR_IS_RATIONAL_MV(MV); int i, MaxCoord = CAGD_NUM_OF_PT_COORD(MV -> PType); MvarMVGradientStruct *MVGrad; if (MaxCoord != 1) { MVAR_FATAL_ERROR(MVAR_ERR_SCALAR_EXPECTED); return NULL; } if (MV -> Dim >= CAGD_MAX_PT_COORD) { MVAR_FATAL_ERROR(MVAR_ERR_DIM_TOO_HIGH); return NULL; } MVGrad = (MvarMVGradientStruct *) IritMalloc(sizeof(MvarMVGradientStruct)); ZAP_MEM(MVGrad, sizeof(MvarMVGradientStruct)); MVGrad -> Dim = MV -> Dim; MVGrad -> IsRational = IsRational; MVGrad -> MV = MvarMVCopy(MV); for (i = 0; i < MV -> Dim; i++) MVGrad -> MVRGrad[i + 1] = MvarMVDerive(MV, i); if (!IsRational) { MVGrad -> MVGrad = MvarMVMergeScalar(MVGrad -> MVRGrad); for (i = 0; i < MV -> Dim; i++) MvarMVFree(MVGrad -> MVRGrad[i + 1]); } return MVGrad; } /***************************************************************************** * DESCRIPTION: M * Free an gradient function. M * * * PARAMETERS: M * MVGrad: Gradient function to free. M * * * RETURN VALUE: M * void M * * * SEE ALSO: M * MvarMVDerive, MvarMVPrepGradient, MvarMVEvalGradient M * * * KEYWORDS: M * MvarMVFreeGradient M *****************************************************************************/ void MvarMVFreeGradient(MvarMVGradientStruct *MVGrad) { if (MVGrad -> IsRational) { int i; for (i = 0; i < MVGrad -> Dim; i++) MvarMVFree(MVGrad -> MVRGrad[i + 1]); } else { MvarMVFree(MVGrad -> MVGrad); } MvarMVFree(MVGrad -> MV); IritFree(MVGrad); } /***************************************************************************** * DESCRIPTION: M * Evaluates the gradient function at the given parameteric location. M * * * PARAMETERS: M * MVGrad: Input gradient function to evaluate at. M * Params: Parametric location to evaluate gradient at. M * * * RETURN VALUE: M * CagdRType *: The gradient at Params. M * * * SEE ALSO: M * MvarMVDerive, MvarMVFreeGradient, MvarMVPrepGradient, M * MvarMVEvalGradient2 M * * * KEYWORDS: M * MvarMVEvalGradient M *****************************************************************************/ CagdRType *MvarMVEvalGradient(MvarMVGradientStruct *MVGrad, CagdRType *Params) { STATIC_DATA CagdRType Grad[CAGD_MAX_PT_COORD]; CagdRType *R; if (MVGrad -> IsRational) { int i; CagdRType Pos[2]; /* Evaluate the multivariate itself - must be a scalar field. */ R = MvarMVEval(MVGrad -> MV, Params); CAGD_GEN_COPY(Pos, R, sizeof(CagdRType) * 2); for (i = 0; i < MVGrad -> Dim; i++) { R = MvarMVEval(MVGrad -> MVRGrad[i + 1], Params); /* Apply the quotient rule: */ Grad[i] = (R[1] * Pos[0] - R[0] * Pos[1]) / SQR(Pos[0]); } } else { R = MvarMVEval(MVGrad -> MVGrad, Params); CAGD_GEN_COPY(Grad, &R[1], sizeof(CagdRType) * MVGrad -> Dim); } return Grad; }