/****************************************************************************** * MvBisect.c - Compute bisectors of multivariates. * ******************************************************************************* * (C) Gershon Elber, Technion, Israel Institute of Technology * ******************************************************************************* * Written by Gershon Elber, May. 97. * ******************************************************************************/ #include #include #include #include "geom_lib.h" #include "mvar_loc.h" /***************************************************************************** * DESCRIPTION: M * Compute bisector to two given multivariates. M * * * PARAMETERS: M * MV1, MV2: The two multivariates to compute the bisector for. M * * * RETURN VALUE: M * MvarMVStruct *: The result bisector. M * * * SEE ALSO: M * SymbSrfPtBisectorSrf3D, SymbCrvPtBisectorSrf3D, SymbCrvCrvBisectorSrf3D M * * * KEYWORDS: M * MvarMVsBisector, bisectors M *****************************************************************************/ MvarMVStruct *MvarMVsBisector(MvarMVStruct *MV1, MvarMVStruct *MV2) { MvarMVStruct *MVBisect = NULL; /* Find the case in hand and dispatch to the proper routine. */ if (MV1 -> Dim == 1 && MV2 -> Dim == 1) { /* Crv-Crv bisector. */ if (CAGD_NUM_OF_PT_COORD(MV1 -> PType) != 3 && CAGD_NUM_OF_PT_COORD(MV2 -> PType) != 3) { MVAR_FATAL_ERROR(MVAR_ERR_PT_OR_LEN_MISMATCH); return NULL; } else { CagdCrvStruct *Crv1 = MvarMVToCrv(MV1), *Crv2 = MvarMVToCrv(MV2); CagdSrfStruct *Srf = SymbCrvCrvBisectorSrf3D(Crv1, Crv2, 0.5); CagdCrvFree(Crv1); CagdCrvFree(Crv2); MVBisect = MvarSrfToMV(Srf); CagdSrfFree(Srf); } } else if ((MV1 -> Dim == 1 && MV2 -> Dim == 2) || (MV1 -> Dim == 2 && MV2 -> Dim == 1)) { /* Crv-Srf bisector. */ if (MV1 -> Dim == 2 && MV2 -> Dim == 1) { MvarMVStruct *MVTmp1 = MvarCrvSrfBisector(MV1, MV2), *MVTmp2 = MvarMVReverse(MVTmp1, 0, 2); MVBisect = MvarMVReverse(MVTmp2, 1, 2); MvarMVFree(MVTmp1); MvarMVFree(MVTmp2); } else MVBisect = MvarCrvSrfBisector(MV1, MV2); } else if (MV1 -> Dim == 2 && MV2 -> Dim == 2) { /* Srf-Srf bisector. */ MVBisect = MvarSrfSrfBisector(MV1, MV2); } else { MVAR_FATAL_ERROR(MVAR_ERR_GEOM_NO_SUPPORT); return NULL; } return MVBisect; } /***************************************************************************** * DESCRIPTION: M * Computes the bisectors of a curve and a surface in R^4. M * * * PARAMETERS: M * MV1: The univariate (curve) in R^4. M * MV2: The bivariate (surface) in R^4. M * * * RETURN VALUE: M * MvarMVStruct *: The resulting bisector. M * * * SEE ALSO: M * MvarMVsBisector, MvarSrfSrfBisector, MvarCrvSrfBisectorApprox M * * * KEYWORDS: M * MvarCrvSrfBisector, bisectors M *****************************************************************************/ MvarMVStruct *MvarCrvSrfBisector(MvarMVStruct *MV1, MvarMVStruct *MV2) { int i, j; STATIC_DATA CagdVType Translate = { 0.0, 0.0, 0.0 }; MvarMVStruct *MVTmp, *MVTmp2, *A[4][4], *B[4], *MVs[CAGD_MAX_PT_SIZE], **MVSplit, *MVBisect; if (CAGD_NUM_OF_PT_COORD(MV1 -> PType) != 4 && CAGD_NUM_OF_PT_COORD(MV2 -> PType) != 4) { MVAR_FATAL_ERROR(MVAR_ERR_PT_OR_LEN_MISMATCH); return NULL; } if (MV1 -> GType != MV2 -> GType) { MVAR_FATAL_ERROR(MVAR_ERR_SAME_GTYPE_EXPECTED); return NULL; } /* Bring both the curve and the surface into a trivariate form. */ if (MV1 -> Dim == 1 && MV2 -> Dim == 2) { /* Crv-Srf bisector. */ CagdRType Min, Max; MV1 = MvarPromoteMVToMV2(MV1, 3, 0); /* Trivariate at axis 0. */ MV2 = MvarPromoteMVToMV2(MV2, 3, 1); /* Trivariate at axes 1,2. */ /* Make sure domain are the same. */ if (MV1 -> GType == MVAR_BSPLINE_TYPE) { MvarMVDomain(MV1, &Min, &Max, 0); BspKnotAffineTrans2(MV2 -> KnotVectors[0], MV2 -> Lengths[0] + MV2 -> Orders[0], Min, Max); MvarMVDomain(MV2, &Min, &Max, 1); BspKnotAffineTrans2(MV1 -> KnotVectors[1], MV1 -> Lengths[1] + MV1 -> Orders[1], Min, Max); MvarMVDomain(MV2, &Min, &Max, 2); BspKnotAffineTrans2(MV1 -> KnotVectors[2], MV1 -> Lengths[2] + MV1 -> Orders[2], Min, Max); } } else { MVAR_FATAL_ERROR(MVAR_ERR_GEOM_NO_SUPPORT); return NULL; } /* Compute the derivative of the curve. */ MVTmp = MvarMVDerive(MV1, 0); MVSplit = MvarMVSplitScalar(MVTmp); for (i = 0; i < 4; i++) A[0][i] = MVSplit[i + 1]; B[0] = MvarMVDotProd(MVTmp, MV1); MvarMVFree(MVTmp); /* Compute the partial derivatives of the surface. */ MVTmp = MvarMVDerive(MV2, 1); MVSplit = MvarMVSplitScalar(MVTmp); for (i = 0; i < 4; i++) A[1][i] = MVSplit[i + 1]; B[1] = MvarMVDotProd(MVTmp, MV2); MvarMVFree(MVTmp); MVTmp = MvarMVDerive(MV2, 2); MVSplit = MvarMVSplitScalar(MVTmp); for (i = 0; i < 4; i++) A[2][i] = MVSplit[i + 1]; B[2] = MvarMVDotProd(MVTmp, MV2); MvarMVFree(MVTmp); /* Compute the distance constraint. */ MVTmp = MvarMVSub(MV1, MV2); MVSplit = MvarMVSplitScalar(MVTmp); for (i = 0; i < 4; i++) A[3][i] = MVSplit[i + 1]; MVTmp2 = MvarMVAdd(MV1, MV2); MvarMVTransform(MVTmp2, Translate, 0.5); B[3] = MvarMVDotProd(MVTmp, MVTmp2); MvarMVFree(MVTmp); MvarMVFree(MVTmp2); /* Done with preparations - compute solution of the 4x4 linear system. */ MVAR_CLEAR_SCALARS(MVs); MVs[0] = MvarMVDeterminant4(A[0][0], A[0][1], A[0][2], A[0][3], A[1][0], A[1][1], A[1][2], A[1][3], A[2][0], A[2][1], A[2][2], A[2][3], A[3][0], A[3][1], A[3][2], A[3][3]); MVs[1] = MvarMVDeterminant4(B[0], A[0][1], A[0][2], A[0][3], B[1], A[1][1], A[1][2], A[1][3], B[2], A[2][1], A[2][2], A[2][3], B[3], A[3][1], A[3][2], A[3][3]); MVs[2] = MvarMVDeterminant4(A[0][0], B[0], A[0][2], A[0][3], A[1][0], B[1], A[1][2], A[1][3], A[2][0], B[2], A[2][2], A[2][3], A[3][0], B[3], A[3][2], A[3][3]); MVs[3] = MvarMVDeterminant4(A[0][0], A[0][1], B[0], A[0][3], A[1][0], A[1][1], B[1], A[1][3], A[2][0], A[2][1], B[2], A[2][3], A[3][0], A[3][1], B[3], A[3][3]); MVs[4] = MvarMVDeterminant4(A[0][0], A[0][1], A[0][2], B[0], A[1][0], A[1][1], A[1][2], B[1], A[2][0], A[2][1], A[2][2], B[2], A[3][0], A[3][1], A[3][2], B[3]); for (i = 0; i < 4; i++) { MvarMVFree(B[i]); for (j = 0; j < 4; j++) MvarMVFree(A[i][j]); } MVBisect = MvarMVMergeScalar(MVs); MVAR_FREE_SCALARS(MVs); MvarMVFree(MV1); MvarMVFree(MV2); return MVBisect; } /***************************************************************************** * DESCRIPTION: M * Computes the bisectors of two surfaces in R^5. M * * * PARAMETERS: M * MV1, MV2: The two bivariates (surfaces) in R^5. M * * * RETURN VALUE: M * MvarMVStruct *: The resulting bisector. M * * * SEE ALSO: M * MvarMVsBisector, MvarCrvSrfBisector, MvarSrfSrfBisectorApprox M * * * KEYWORDS: M * MvarSrfSrfBisector, bisectors M *****************************************************************************/ MvarMVStruct *MvarSrfSrfBisector(MvarMVStruct *MV1, MvarMVStruct *MV2) { int i, j; STATIC_DATA CagdVType Translate = { 0.0, 0.0, 0.0 }; MvarMVStruct *MVTmp, *MVTmp2, *A[5][5], *B[5], *MVs[CAGD_MAX_PT_SIZE], **MVSplit, *MVBisect; if (CAGD_NUM_OF_PT_COORD(MV1 -> PType) != 5 && CAGD_NUM_OF_PT_COORD(MV2 -> PType) != 5) { MVAR_FATAL_ERROR(MVAR_ERR_PT_OR_LEN_MISMATCH); return NULL; } if (MV1 -> GType != MV2 -> GType) { MVAR_FATAL_ERROR(MVAR_ERR_SAME_GTYPE_EXPECTED); return NULL; } /* Bring both surfaces into a four-variate form. */ if (MV1 -> Dim == 2 && MV2 -> Dim == 2) { CagdRType Min, Max; MV1 = MvarPromoteMVToMV2(MV1, 4, 0); /* Four variate at axes 0,1. */ MV2 = MvarPromoteMVToMV2(MV2, 4, 2); /* Four variate at axes 2,3. */ /* Make sure domain are the same. */ if (MV1 -> GType == MVAR_BSPLINE_TYPE) { MvarMVDomain(MV1, &Min, &Max, 0); BspKnotAffineTrans2(MV2 -> KnotVectors[0], MV2 -> Lengths[0] + MV2 -> Orders[0], Min, Max); MvarMVDomain(MV1, &Min, &Max, 1); BspKnotAffineTrans2(MV2 -> KnotVectors[1], MV2 -> Lengths[1] + MV2 -> Orders[1], Min, Max); MvarMVDomain(MV2, &Min, &Max, 2); BspKnotAffineTrans2(MV1 -> KnotVectors[2], MV1 -> Lengths[2] + MV1 -> Orders[2], Min, Max); MvarMVDomain(MV2, &Min, &Max, 3); BspKnotAffineTrans2(MV1 -> KnotVectors[3], MV1 -> Lengths[3] + MV1 -> Orders[3], Min, Max); } } else { MVAR_FATAL_ERROR(MVAR_ERR_GEOM_NO_SUPPORT); return NULL; } /* Compute the partial derivatives of the first surface. */ MVTmp = MvarMVDerive(MV1, 0); MVSplit = MvarMVSplitScalar(MVTmp); for (i = 0; i < 5; i++) A[0][i] = MVSplit[i + 1]; B[0] = MvarMVDotProd(MVTmp, MV1); MvarMVFree(MVTmp); MVTmp = MvarMVDerive(MV1, 1); MVSplit = MvarMVSplitScalar(MVTmp); for (i = 0; i < 5; i++) A[1][i] = MVSplit[i + 1]; B[1] = MvarMVDotProd(MVTmp, MV1); MvarMVFree(MVTmp); /* Compute the partial derivatives of the second surface. */ MVTmp = MvarMVDerive(MV2, 2); MVSplit = MvarMVSplitScalar(MVTmp); for (i = 0; i < 5; i++) A[2][i] = MVSplit[i + 1]; B[2] = MvarMVDotProd(MVTmp, MV2); MvarMVFree(MVTmp); MVTmp = MvarMVDerive(MV2, 3); MVSplit = MvarMVSplitScalar(MVTmp); for (i = 0; i < 5; i++) A[3][i] = MVSplit[i + 1]; B[3] = MvarMVDotProd(MVTmp, MV2); MvarMVFree(MVTmp); /* Compute the distance constraint. */ MVTmp = MvarMVSub(MV1, MV2); MVSplit = MvarMVSplitScalar(MVTmp); for (i = 0; i < 5; i++) A[4][i] = MVSplit[i + 1]; MVTmp2 = MvarMVAdd(MV1, MV2); MvarMVTransform(MVTmp2, Translate, 0.5); B[4] = MvarMVDotProd(MVTmp, MVTmp2); MvarMVFree(MVTmp); MvarMVFree(MVTmp2); /* Done with preparations - compute solution of the 4x4 linear system. */ MVAR_CLEAR_SCALARS(MVs); MVs[0] = MvarMVDeterminant5(A[0][0], A[0][1], A[0][2], A[0][3], A[0][4], A[1][0], A[1][1], A[1][2], A[1][3], A[1][4], A[2][0], A[2][1], A[2][2], A[2][3], A[2][4], A[3][0], A[3][1], A[3][2], A[3][3], A[3][4], A[4][0], A[4][1], A[4][2], A[4][3], A[4][4]); MVs[1] = MvarMVDeterminant5(B[0], A[0][1], A[0][2], A[0][3], A[0][4], B[1], A[1][1], A[1][2], A[1][3], A[1][4], B[2], A[2][1], A[2][2], A[2][3], A[2][4], B[3], A[3][1], A[3][2], A[3][3], A[3][4], B[4], A[4][1], A[4][2], A[4][3], A[4][4]); MVs[2] = MvarMVDeterminant5(A[0][0], B[0], A[0][2], A[0][3], A[0][4], A[1][0], B[1], A[1][2], A[1][3], A[1][4], A[2][0], B[2], A[2][2], A[2][3], A[2][4], A[3][0], B[3], A[3][2], A[3][3], A[3][4], A[4][0], B[4], A[4][2], A[4][3], A[4][4]); MVs[3] = MvarMVDeterminant5(A[0][0], A[0][1], B[0], A[0][3], A[0][4], A[1][0], A[1][1], B[1], A[1][3], A[1][4], A[2][0], A[2][1], B[2], A[2][3], A[2][4], A[3][0], A[3][1], B[3], A[3][3], A[3][4], A[4][0], A[4][1], B[4], A[4][3], A[4][4]); MVs[4] = MvarMVDeterminant5(A[0][0], A[0][1], A[0][2], B[0], A[0][4], A[1][0], A[1][1], A[1][2], B[1], A[1][4], A[2][0], A[2][1], A[2][2], B[2], A[2][4], A[3][0], A[3][1], A[3][2], B[3], A[3][4], A[4][0], A[4][1], A[4][2], B[4], A[4][4]); MVs[0] = MvarMVDeterminant5(A[0][0], A[0][1], A[0][2], A[0][3], B[0], A[1][0], A[1][1], A[1][2], A[1][3], B[1], A[2][0], A[2][1], A[2][2], A[2][3], B[2], A[3][0], A[3][1], A[3][2], A[3][3], B[3], A[4][0], A[4][1], A[4][2], A[4][3], B[4]); for (i = 0; i < 5; i++) { MvarMVFree(B[i]); for (j = 0; j < 5; j++) MvarMVFree(A[i][j]); } MVBisect = MvarMVMergeScalar(MVs); MVAR_FREE_SCALARS(MVs); MvarMVFree(MV1); MvarMVFree(MV2); return MVBisect; } /***************************************************************************** * DESCRIPTION: M * Computes an approximation to the bisector of a curve and a surface. M * Let C(t) be the parametric curve and T(t) its unnormalized tangent field. M * Let S(u,v) be the parametric surface and n(u, v) its unnormalized normal. M * Then: M * M * < P - C(t), T(t) > = 0 defines the normal plane of T(t), V * M * and the solution of M * M * < S(u, v) + n(u, v) Alpha - C(t), T(t) > = 0 V * M * finds the intersection point of of the normal of the surface with the M * normal plane of the curve. Then, M * M * < C(t) - S(u, v), T(t) > V * Alpha(u, v, t) = ------------------------ V * < n(u, v), T(t) > V * M * We now can define this intersection point, P, as M * M * P(u, v, t) = S(u, v) + Alpha(u, v, t) n(u, v) M * M * and end up with a single function we must extract its zero set M * M * < C(t) - P(u, v, t), C(t) - P(u, v, t) > - V * < S(u, v) - P(u, v, t), S(u, v) - P(u, v, t) > = 0 V * M * or M * M * < C(t) - P(u, v, t), C(t) - P(u, v, t) > - V * < Alpha(u, v, t) n(u, v), Alpha(u, v, t) n(u, v) > = 0 V * M * Finding the zero set of the last equation provides the correspondance M * between the (u, v) location and the surface and (t) locations on the curve M * that serve as mutual foot point for some bisector point. M * * * PARAMETERS: M * MV1: The univariate (curve) in R^3. M * MV2: The bivariate (surface) in R^3. M * OutputType: Expected output type: M * 1. For the computed multivariate constraints. M * 2. For the computed point cloud on the bisector. M * 3. Points in a form of (u, v, x, y, z) where (u, v) are M * the parameter space of the surface. M * SubdivTol: Tolerance of the first zero set finding subdivision stage. M * NumericTol: Tolerance of the second zero set finding numeric stage. M * * * RETURN VALUE: M * VoidPtr: Following OutputType, either a set of multivariates (as a M * linked list of MvarMVStruct), or a cloud of points on the M * bisector (as a linked list of MvarPtStruct). M * * * SEE ALSO: M * MvarMVsBisector, MvarSrfSrfBisector, MvarCrvSrfBisector M * MvarSrfSrfBisectorApprox M * * * KEYWORDS: M * MvarCrvSrfBisectorApprox, bisectors M *****************************************************************************/ VoidPtr MvarCrvSrfBisectorApprox(MvarMVStruct *MV1, MvarMVStruct *MV2, int OutputType, CagdRType SubdivTol, CagdRType NumericTol) { STATIC_DATA MvarConstraintType Constraint = MVAR_CNSTRNT_ZERO; MvarMVStruct *MVTmp1, *MVTmp2, *MVTmp3, *Alpha, *P, *DMV1, *NrmlMV2, *Final; MvarPtStruct *Pts, *Pt; if (CAGD_NUM_OF_PT_COORD(MV1 -> PType) != 3 && CAGD_NUM_OF_PT_COORD(MV2 -> PType) != 3) { MVAR_FATAL_ERROR(MVAR_ERR_PT_OR_LEN_MISMATCH); return NULL; } if (MV1 -> GType != MV2 -> GType) { MVAR_FATAL_ERROR(MVAR_ERR_SAME_GTYPE_EXPECTED); return NULL; } /* Bring both the curve and the surface into a trivariate form. */ if (MV1 -> Dim == 1 && MV2 -> Dim == 2) { /* Crv-Srf bisector. */ CagdRType Min, Max; MV1 = MvarPromoteMVToMV2(MV1, 3, 0); /* Trivariate at axis 0. */ MV2 = MvarPromoteMVToMV2(MV2, 3, 1); /* Trivariate at axes 1,2. */ /* Make sure domains are the same. */ if (MV1 -> GType == MVAR_BSPLINE_TYPE) { MvarMVDomain(MV1, &Min, &Max, 0); BspKnotAffineTrans2(MV2 -> KnotVectors[0], MV2 -> Lengths[0] + MV2 -> Orders[0], Min, Max); MvarMVDomain(MV2, &Min, &Max, 1); BspKnotAffineTrans2(MV1 -> KnotVectors[1], MV1 -> Lengths[1] + MV1 -> Orders[1], Min, Max); MvarMVDomain(MV2, &Min, &Max, 2); BspKnotAffineTrans2(MV1 -> KnotVectors[2], MV1 -> Lengths[2] + MV1 -> Orders[2], Min, Max); } } else { MVAR_FATAL_ERROR(MVAR_ERR_GEOM_NO_SUPPORT); return NULL; } DMV1 = MvarMVDerive(MV1, 0); MVTmp1 = MvarMVDerive(MV2, 1); MVTmp2 = MvarMVDerive(MV2, 2); NrmlMV2 = MvarMVCrossProd(MVTmp1, MVTmp2); MvarMVFree(MVTmp1); MvarMVFree(MVTmp2); /* Compute Alpha(u, v, t). */ MVTmp1 = MvarMVSub(MV1, MV2); MVTmp2 = MvarMVDotProd(MVTmp1, DMV1); /* We have the numerator! */ MvarMVFree(MVTmp1); MVTmp1 = MvarMVDotProd(NrmlMV2, DMV1); /* We have the denominator! */ if (MVAR_IS_RATIONAL_MV(MVTmp1) || MVAR_IS_RATIONAL_MV(MVTmp2)) { MvarMVStruct **ScalarMVs, *MVs[CAGD_MAX_PT_SIZE], *MVTmp1W, *MVTmp1X, *MVTmp2W, *MVTmp2X; ScalarMVs = MvarMVSplitScalar(MVTmp1); MVTmp1W = ScalarMVs[0]; MVTmp1X = ScalarMVs[1]; ScalarMVs = MvarMVSplitScalar(MVTmp2); MVTmp2W = ScalarMVs[0]; MVTmp2X = ScalarMVs[1]; if (MVTmp1W != NULL) { MvarMVStruct *MVTmp = MvarMVMult(MVTmp1W, MVTmp2X); MvarMVFree(MVTmp2X); MVTmp2X = MVTmp; } if (MVTmp2W != NULL) { MvarMVStruct *MVTmp = MvarMVMult(MVTmp2W, MVTmp1X); MvarMVFree(MVTmp1X); MVTmp1X = MVTmp; } MvarMVFree(MVTmp1W); MvarMVFree(MVTmp2W); MVAR_CLEAR_SCALARS(MVs); MVs[0] = MVTmp1X; MVs[1] = MVTmp2X; Alpha = MvarMVMergeScalar(MVs); MvarMVFree(MVTmp1X); MvarMVFree(MVTmp2X); } else { MvarMVStruct *MVs[CAGD_MAX_PT_SIZE]; MVAR_CLEAR_SCALARS(MVs); MVs[0] = MVTmp1; MVs[1] = MVTmp2; Alpha = MvarMVMergeScalar(MVs); } MvarMVFree(MVTmp1); MvarMVFree(MVTmp2); /* Compute S(u, v) + Alpha(u, v, t) n(u, v). */ MVTmp1 = MvarMVMultScalar(NrmlMV2, Alpha); P = MvarMVAdd(MV2, MVTmp1); MvarMVFree(MVTmp1); /* Compute the final function, we seek its zero set. */ MVTmp1 = MvarMVSub(MV1, P); MVTmp2 = MvarMVDotProd(MVTmp1, MVTmp1); MvarMVFree(MVTmp1); MVTmp1 = MvarMVSub(MV2, P); MVTmp3 = MvarMVDotProd(MVTmp1, MVTmp1); MvarMVFree(MVTmp1); Final = MvarMVSub(MVTmp2, MVTmp3); MvarMVFree(MVTmp2); MvarMVFree(MVTmp3); MvarMVFree(DMV1); MvarMVFree(NrmlMV2); MvarMVFree(Alpha); MvarMVFree(MV1); MvarMVFree(MV2); if (OutputType == 1) { /* Place "S(u, v) + Alpha(u, v, t) n(u, v)" as "MVEuclid" attrib. */ AttrSetPtrAttrib(&Final -> Attr, "MVEuclid", P); return Final; } /* Compute the Zero set for Final. */ Pts = MvarMVsZeros(&Final, &Constraint, 1, SubdivTol, NumericTol); MvarMVFree(Final); /* Evaluate these points at P for the points on the bisector surface. */ for (Pt = Pts; Pt != NULL; Pt = Pt -> Pnext) { RealType *R; R = MvarMVEval(P, Pt -> Pt); if (OutputType == 2) { CagdCoerceToE3(Pt -> Pt, &R, -1, P -> PType); Pt -> Dim = 3; } else { /* OutputType == 3 */ /* Allocate a longer vector to hold (u, v, x, y, z): */ MvarPtRealloc(Pt, 5); Pt -> Pt[0] = Pt -> Pt[1]; Pt -> Pt[1] = Pt -> Pt[2]; CagdCoerceToE3(&Pt -> Pt[2], &R, -1, P -> PType); } } MvarMVFree(P); return Pts; } /***************************************************************************** * DESCRIPTION: M * Computes an approximation to the bisector of two surfaces. M * Let S1(u, v) and S2(r, s) be two parametric surfaces and let n1(u, v) and M * n2(r, s) be their unnormalized normal fields. M * Becuase the two normals of the two surfaces must be coplanar we M * introduce the following constraint, forcing the three vectors n1, n2, and M * S1 - S2 to all be in the same plane. M * M * < ( S1(u, v) - S2(r, s) ) x n1(u, v), n2(r, s) > = 0. V * M * To make sure the distance to the intersection point of the normals, from M * both surface's foot points we also coerces these three vectors to form a M * isosceles triangle: M * M * || n2(r, s) ||^2 < S1(u, v) - S2(r, s), n1(u, v) > ^ 2 - V * || n1(r, s) ||^2 < S1(u, v) - S2(r, s), n2(u, v) > ^ 2 V * M * Finding the zero set of the last equation provides the correspondance M * between the (u, v) location and the first surface and (r, s) locations on M * the second surface that serve as mutual foot point for some bisector M * point. M * * * PARAMETERS: M * MV1, MV2: The two bivariates (surfaces) in R^3. M * OutputType: Expected output type: M * 1. For the computed multivariate constraints. M * 2. For the computed point cloud on the bisector. M * 3. Points in a form of (u1, v2, x, y, z) where (u1, v1) are M * the parameter space of the first surface. M * SubdivTol: Tolerance of the first zero set finding subdivision stage. M * NumericTol: Tolerance of the second zero set finding numeric stage. M * * * RETURN VALUE: M * VoidPtr: Following OutputType, either a set of multivariates (as a M * linked list of MvarMVStruct), or a cloud of points on the M * bisector (as a linked list of MvarPtStruct). M * * * SEE ALSO: M * MvarMVsBisector, MvarCrvSrfBisector, MvarSrfSrfBisector M * MvarCrvSrfBisectorApprox M * * * KEYWORDS: M * MvarSrfSrfBisectorApprox, bisectors M *****************************************************************************/ VoidPtr MvarSrfSrfBisectorApprox(MvarMVStruct *MV1, MvarMVStruct *MV2, int OutputType, CagdRType SubdivTol, CagdRType NumericTol) { STATIC_DATA MvarConstraintType Constraints[2] = { MVAR_CNSTRNT_ZERO, MVAR_CNSTRNT_ZERO }; MvarMVStruct *MVTmp1, *MVTmp2, *MVTmp3, *MVTmp4, *MVTmp5, *NrmlMV1, *NrmlMV2, *AlphaNumer, *AlphaDenom, *Finals[2]; MvarPtStruct *Pts, *Pt; if (CAGD_NUM_OF_PT_COORD(MV1 -> PType) != 3 && CAGD_NUM_OF_PT_COORD(MV2 -> PType) != 3) { MVAR_FATAL_ERROR(MVAR_ERR_PT_OR_LEN_MISMATCH); return NULL; } if (MV1 -> GType != MV2 -> GType) { MVAR_FATAL_ERROR(MVAR_ERR_SAME_GTYPE_EXPECTED); return NULL; } /* Bring both surfaces into a four-variate form. */ if (MV1 -> Dim == 2 && MV2 -> Dim == 2) { CagdRType Min, Max; MV1 = MvarPromoteMVToMV2(MV1, 4, 0); /* Four variate at axes 0,1. */ MV2 = MvarPromoteMVToMV2(MV2, 4, 2); /* Four variate at axes 2,3. */ /* Make sure domains are the same. */ if (MV1 -> GType == MVAR_BSPLINE_TYPE) { MvarMVDomain(MV1, &Min, &Max, 0); BspKnotAffineTrans2(MV2 -> KnotVectors[0], MV2 -> Lengths[0] + MV2 -> Orders[0], Min, Max); MvarMVDomain(MV1, &Min, &Max, 1); BspKnotAffineTrans2(MV2 -> KnotVectors[1], MV2 -> Lengths[1] + MV2 -> Orders[1], Min, Max); MvarMVDomain(MV2, &Min, &Max, 2); BspKnotAffineTrans2(MV1 -> KnotVectors[2], MV1 -> Lengths[2] + MV1 -> Orders[2], Min, Max); MvarMVDomain(MV2, &Min, &Max, 3); BspKnotAffineTrans2(MV1 -> KnotVectors[3], MV1 -> Lengths[3] + MV1 -> Orders[3], Min, Max); } } else { MVAR_FATAL_ERROR(MVAR_ERR_GEOM_NO_SUPPORT); return NULL; } /* Compute normal fields. */ MVTmp1 = MvarMVDerive(MV1, 0); MVTmp2 = MvarMVDerive(MV1, 1); NrmlMV1 = MvarMVCrossProd(MVTmp1, MVTmp2); MvarMVFree(MVTmp1); MvarMVFree(MVTmp2); MVTmp4 = MvarMVDerive(MV2, 2); MVTmp5 = MvarMVDerive(MV2, 3); NrmlMV2 = MvarMVCrossProd(MVTmp4, MVTmp5); /* Compute Alpha(u, v, r, s) = / -2 .*/ MVTmp1 = MvarMVSub(MV1, MV2); AlphaNumer = MvarMVDotProd(MVTmp1, MVTmp1); MVTmp2 = MvarMVDotProd(NrmlMV1, MVTmp1); AlphaDenom = MvarMVScalarScale(MVTmp2, -2); MvarMVFree(MVTmp2); /* Compute Delta = -2 (S1 - S2) + n1 . */ MVTmp2 = MvarMVMultScalar(MVTmp1, AlphaDenom); MVTmp3 = MvarMVMultScalar(NrmlMV1, AlphaNumer); MvarMVFree(MVTmp1); MVTmp1 = MvarMVAdd(MVTmp2, MVTmp3); MvarMVFree(MVTmp2); MvarMVFree(MVTmp3); /* Compute , and . */ Finals[0] = MvarMVDotProd(MVTmp1, MVTmp4); Finals[1] = MvarMVDotProd(MVTmp1, MVTmp5); MvarMVFree(MVTmp1); MvarMVFree(MVTmp4); MvarMVFree(MVTmp5); if (OutputType == 1) { MvarMVFree(AlphaDenom); MvarMVFree(AlphaNumer); MvarMVFree(NrmlMV1); MvarMVFree(NrmlMV2); MvarMVFree(MV1); MvarMVFree(MV2); Finals[0] -> Pnext = Finals[1]; return Finals[0]; } /* Compute the Zero set for Finals. */ Pts = MvarMVsZeros(Finals, Constraints, 2, SubdivTol, NumericTol); MvarMVFree(Finals[0]); MvarMVFree(Finals[1]); /* Evaluate these points for the points on the bisector surface. */ for (Pt = Pts; Pt != NULL; Pt = Pt -> Pnext) { RealType *R, t1; PointType Pt1E3; VectorType Nrml1E3; R = MvarMVEval(MV1, Pt -> Pt); CagdCoerceToE3(Pt1E3, &R, -1, MV1 -> PType); R = MvarMVEval(AlphaNumer, Pt -> Pt); t1 = R[1]; R = MvarMVEval(AlphaDenom, Pt -> Pt); t1 /= R[1]; R = MvarMVEval(NrmlMV1, Pt -> Pt); CagdCoerceToE3(Nrml1E3, &R, -1, NrmlMV1 -> PType); PT_SCALE(Nrml1E3, t1); PT_ADD(Pt1E3, Pt1E3, Nrml1E3); if (OutputType == 2) { PT_COPY(Pt -> Pt, Pt1E3); Pt -> Dim = 3; } else { /* OutputType == 3 */ /* Allocate a longer vector to hold (u, v, x, y, z): */ MvarPtRealloc(Pt, 5); PT_COPY(&Pt -> Pt[2], Pt1E3); } } MvarMVFree(AlphaDenom); MvarMVFree(AlphaNumer); MvarMVFree(NrmlMV1); MvarMVFree(NrmlMV2); MvarMVFree(MV1); MvarMVFree(MV2); return Pts; }