/****************************************************************************** * MvCones.c - Tools to construct and intersect MV (anti-)cones & vectors. * ******************************************************************************* * (C) Gershon Elber, Technion, Israel Institute of Technology * ******************************************************************************* * Written by Iddo Haniel and Gershon Elber, May 2005. * ******************************************************************************/ #include "mvar_loc.h" #include #include /***************************************************************************** * DESCRIPTION: M * Add two multivariate vectors. M * * * PARAMETERS: M * VRes: Result. Can be one of V1 or V2. M * V1,V2: Two input vectors. M * * * RETURN VALUE: M * void M * * * SEE ALSO: M * MvarVecDotProd, MvarVecSqrLength, MvarVecLength, MvarVecScale M * MvarVecNormalize M * * * KEYWORDS: M * MvarVecAdd M *****************************************************************************/ void MvarVecAdd(MvarVecStruct *VRes, MvarVecStruct *V1, MvarVecStruct *V2) { int i; assert(V1 -> Dim == V2 -> Dim && VRes -> Dim == V2 -> Dim); for (i = 0; i < V1 -> Dim; i++) VRes -> Vec[i] = V1 -> Vec[i] + V2 -> Vec[i]; } /***************************************************************************** * DESCRIPTION: M * Compute the dot product of two multivariate vectors. M * * * PARAMETERS: M * V1,V2: Two input vectors. M * * * RETURN VALUE: M * CagdRType: The dot product. M * * * SEE ALSO: M * MvarVecAdd, MvarVecSqrLength, MvarVecLength, MvarVecScale M * MvarVecNormalize M * * * KEYWORDS: M * MvarVecDotProd M *****************************************************************************/ CagdRType MvarVecDotProd(MvarVecStruct *V1, MvarVecStruct *V2) { int i; CagdRType DotProd = 0.0; assert(V1 -> Dim == V2 -> Dim); for (i = 0; i < V1 -> Dim; i++) DotProd += V1 -> Vec[i] * V2 -> Vec[i]; return DotProd; } /***************************************************************************** * DESCRIPTION: M * Compute the length squared of a multivariate vector. M * * * PARAMETERS: M * V: Vector to compute its length. M * * * RETURN VALUE: M * CagdRType: Computed length squared. M * * * SEE ALSO: M * MvarVecAdd, MvarVecDotProd, MvarVecLength, MvarVecScale M * MvarVecNormalize M * * * KEYWORDS: M * MvarVecSqrLength M *****************************************************************************/ CagdRType MvarVecSqrLength(MvarVecStruct *V) { return MvarVecDotProd(V, V); } /***************************************************************************** * DESCRIPTION: M * Compute the length of a multivariate vector. M * * * PARAMETERS: M * V: Vector to compute its length. M * * * RETURN VALUE: M * CagdRType: Computed length. M * * * SEE ALSO: M * MvarVecAdd, MvarVecDotProd, MvarVecSqrLength, MvarVecScale M * MvarVecNormalize M * * * KEYWORDS: M * MvarVecSqrLength M *****************************************************************************/ CagdRType MvarVecLength(MvarVecStruct *V) { return sqrt(MvarVecSqrLength(V)); } /***************************************************************************** * DESCRIPTION: M * Scale a given multivariate vector V by a scaling factor ScaleFactor. M * * * PARAMETERS: M * V: Vector to scale, in place. M * ScaleFactor: Scaling factor to use. M * * * RETURN VALUE: M * void M * * * SEE ALSO: M * MvarVecAdd, MvarVecDotProd, MvarVecSqrLength, MvarVecLength, M * MvarVecNormalize M * * * KEYWORDS: M * MvarVecScale M *****************************************************************************/ void MvarVecScale(MvarVecStruct *V, CagdRType ScaleFactor) { int i; for (i = 0; i < V -> Dim; i++) V -> Vec[i] *= ScaleFactor; } /***************************************************************************** * DESCRIPTION: M * Normalize a given multivariate vector to a unit length, in place. M * * * PARAMETERS: M * V: Vector to normalize. M * * * RETURN VALUE: M * int: TRUE if successful, FALSE if the input is the ZERO vector. M * * * SEE ALSO: M * MvarVecAdd, MvarVecDotProd, MvarVecSqrLength, MvarVecLength, M * MvarVecScale M * * * KEYWORDS: M * MvarVecNormalize M *****************************************************************************/ int MvarVecNormalize(MvarVecStruct *V) { CagdRType VecLengthFactor = MvarVecLength(V); if (VecLengthFactor != 0.0) { VecLengthFactor = 1.0 / MvarVecLength(V); MvarVecScale(V, VecLengthFactor); return TRUE; } else return FALSE; } /***************************************************************************** * DESCRIPTION: M * Constructs a multivariate normal cone structure. M * * * PARAMETERS: M * Dim: Dimension of the cone. M * * * RETURN VALUE: M * MvarNormalConeStruct *: Constructed cone. M * * * SEE ALSO: M * MvarNormalConeFree, MvarNormalConeFreeList M * * * KEYWORDS: M * MvarNormalConeNew M *****************************************************************************/ MvarNormalConeStruct *MvarNormalConeNew(int Dim) { MvarNormalConeStruct *NewCone = (MvarNormalConeStruct *) IritMalloc(sizeof(MvarNormalConeStruct)); NewCone -> ConeAxis = MvarVecNew(Dim); ZAP_MEM((NewCone -> ConeAxis -> Vec), sizeof(CagdRType) * Dim); NewCone -> Attr = NULL; NewCone -> Pnext = NULL; return NewCone; } /***************************************************************************** * DESCRIPTION: M * Free a multivariate normal cone structure. M * * * PARAMETERS: M * NormalCone: Normal cone to free. M * * * RETURN VALUE: M * void M * * * SEE ALSO: M * MvarNormalConeNew, MvarNormalConeFreeList M * * * KEYWORDS: M * MvarNormalConeFree M *****************************************************************************/ void MvarNormalConeFree(MvarNormalConeStruct *NormalCone) { IP_ATTR_FREE_ATTRS(NormalCone -> Attr); MvarVecFree(NormalCone -> ConeAxis); IritFree(NormalCone); } /***************************************************************************** * DESCRIPTION: M * Deallocates and frees a list of multi-variate cone structures. M * * * PARAMETERS: M * NormalConeList: Multi-Variate cone list to free. M * * * RETURN VALUE: M * void M * * * SEE ALSO: M * MvarNormalConeFree, MvarNormalConeNew M * * * KEYWORDS: M * MvarNormalConeFreeList M *****************************************************************************/ void MvarNormalConeFreeList(MvarNormalConeStruct *NormalConeList) { MvarNormalConeStruct *NormalConeTemp; while (NormalConeList) { NormalConeTemp = NormalConeList -> Pnext; MvarNormalConeFree(NormalConeList); NormalConeList = NormalConeTemp; } } /***************************************************************************** * DESCRIPTION: M * Constructs a normal cone to a scalar multivariate MV. M * Note the normal space of the trivariate is assumed of dimension one, and M * the gradient of the multivariate is assumed to span the normal space. M * If the input multivariate is not scalar it is assumed to be of point M * type E(1+Dim), where Dim is the dimension of the MV. This scalar holds M * the gradient already in the Dim locations of the Points array in MV, in M * slots Points[2] to Points[Dim + 1], with Points[1] still holding the M * scalar field. M * * * PARAMETERS: M * MV: Multivariate to derive a cone bounding its normal spaace. M * * * RETURN VALUE: M * MvarNormalConeStruct *: A cone bounding the normal space ov MV, or NULL M * if failed (i.e. angular span of normal space too large). M * * * SEE ALSO: M * MvarMVConesOverlap M * * * KEYWORDS: M * MVarMVNormalCone M *****************************************************************************/ MvarNormalConeStruct *MVarMVNormalCone(MvarMVStruct *MV) { CagdBType IsRational = MVAR_IS_RATIONAL_MV(MV); int i, TotalLength = MVAR_CTL_MESH_LENGTH(MV), Dim = MV -> Dim, Index = 0; CagdRType **GradPoints, ConeAngleCosine, VecLength; MvarVecStruct *UnitVec = MvarVecNew(Dim); MvarNormalConeStruct *NormalCone = MvarNormalConeNew(Dim); MvarMVGradientStruct *Grad; if (IsRational) { MVAR_FATAL_ERROR(MVAR_ERR_RATIONAL_NO_SUPPORT); return NULL; } if (MVAR_NUM_OF_PT_COORD(MV) == 1) { /* No precomputed gradient. */ Grad = MvarMVPrepGradient(MV); GradPoints = &Grad -> MVGrad -> Points[1]; } else if (MV -> Dim == MVAR_NUM_OF_PT_COORD(MV) - 1) { /* Gradient is embedded in Points[2] to Points[Dim + 1]. */ Grad = NULL; GradPoints = &MV -> Points[2]; } else { MVAR_FATAL_ERROR(MVAR_ERR_DIM_TOO_HIGH); return NULL; } /* Construct the central axis of the cone, ConeAxis, as an average of */ /* the control points of the gradient. */ for (Index = 0; Index < TotalLength; Index++) { if (IsRational) { MVAR_FATAL_ERROR(MVAR_ERR_RATIONAL_NO_SUPPORT); } else { for (i = 0; i < Dim; i++) UnitVec -> Vec[i] = GradPoints[i][Index]; } MvarVecNormalize(UnitVec); MvarVecAdd(NormalCone -> ConeAxis, NormalCone -> ConeAxis, UnitVec); } MvarVecNormalize(NormalCone -> ConeAxis); /* Find the maximal angle, ConeAngleCosine, between ConeAxis and some */ /* vector in mesh. */ ConeAngleCosine = 1.0; for (Index = 0; Index < TotalLength; Index++) { CagdRType InnerProd; if (IsRational) { MVAR_FATAL_ERROR(MVAR_ERR_RATIONAL_NO_SUPPORT); } else { for (i = 0; i < Dim; ++i) UnitVec -> Vec[i] = GradPoints[i][Index]; } VecLength = MvarVecLength(UnitVec); if (VecLength > IRIT_UEPS) { MvarVecScale(UnitVec, 1.0 / VecLength); InnerProd = MvarVecDotProd(NormalCone -> ConeAxis, UnitVec); if (InnerProd < 0.0) { /* More than 90 degrees from axis of cone - abort. */ MvarVecFree(UnitVec); if (Grad != NULL) MvarMVFreeGradient(Grad); MvarNormalConeFree(NormalCone); return NULL; } if (ConeAngleCosine > InnerProd) ConeAngleCosine = InnerProd; } } if (Grad != NULL) MvarMVFreeGradient(Grad); MvarVecFree(UnitVec); NormalCone -> ConeAngleCosine = ConeAngleCosine; return NormalCone; } /***************************************************************************** * DESCRIPTION: M * Computes the tangency anti-cones of the set of multivariate constraints, M * and returns whether they overlap or not. M * * * PARAMETERS: M * MVs: Multivariates to derive their tangentcy anti-cones. M * NumOfZeroMVs: Size of the vector MVs. M * * * RETURN VALUE: M * CagdBType: M * * * SEE ALSO: M * MVarMVNormalCone M * * * KEYWORDS: M * MvarMVConesOverlap M *****************************************************************************/ CagdBType MvarMVConesOverlap(MvarMVStruct **MVs, int NumOfZeroMVs) { static int AllocDim = 0; static CagdRType *A = NULL, *x = NULL, *b = NULL, *bCopy = NULL; int i = 0, j = 0, k = 0, Dim = MVs[0] -> Dim, PowerTwoDim = (int) pow(2, Dim - 1); /* In order to save alloc/free run-time we use static variables. */ if (AllocDim < Dim) { if (AllocDim > 0) { IritFree(A); IritFree(x); IritFree(b); IritFree(bCopy); } AllocDim = Dim * 2; A = (CagdRType *) IritMalloc(sizeof(CagdRType) * AllocDim * AllocDim); x = (CagdRType *) IritMalloc(sizeof(CagdRType) * AllocDim); b = (CagdRType *) IritMalloc(sizeof(CagdRType) * AllocDim); bCopy = (CagdRType *) IritMalloc(sizeof(CagdRType) * AllocDim); } for (i = 0; i < NumOfZeroMVs; i++) { MvarNormalConeStruct *Cone; if ((Cone = MVarMVNormalCone(MVs[i])) == NULL) { /* Angular span is too large - return cones overlap. */ return TRUE; } /* Cone is valid: add plane orthogonal to cone axis to matrix A. */ CAGD_GEN_COPY(&A[i * Dim], Cone -> ConeAxis -> Vec, sizeof(CagdRType) * Dim); /* Take the anti-cone = cos(90 - angle) = sqrt(1-cos^2). */ b[i] = sqrt(1.0 - SQR(Cone -> ConeAngleCosine)); MvarNormalConeFree(Cone); } /* We now have the matrix A containing planes orthogonal to cone's axes */ /* and the vector b of the expected solutions. */ /* Compute QR decomposition of matrix A. */ if (IritQRUnderdetermined(A, NULL, NULL, Dim, Dim)) { return TRUE; /* Something went wrong - return cones are overlapping. */ } /* Loop over 2^(d-1) combinations of b vector (000 -> ---, 111 -> +++). */ /* If Qx=b returns a point that is out of unit hyper-sphere return TRUE */ /* - cones overlap. */ for (i = 0; i < PowerTwoDim; i++) { /* Construct relevant copy of b (+/- of b[j] defined by binary rep. */ CagdRType VecSqrLength; k = i; for (j = 0; j < Dim; ++j) { bCopy[j] = k & 1 ? b[j] : -b[j]; k >>= 1; } IritQRUnderdetermined(NULL, x, bCopy, Dim, Dim); for (VecSqrLength = 0.0, j = 0; j < Dim; ++j) { VecSqrLength += SQR(x[j]); } if (VecSqrLength >= 1.0) return TRUE; } return FALSE; }