/****************************************************************************** * EvalCurv.c - Evaluate curvature properties of curves and surfaces. * ******************************************************************************* * (C) Gershon Elber, Technion, Israel Institute of Technology * ******************************************************************************* * Written by Gershon Elber, Nov. 96. * ******************************************************************************/ #include #include "symb_loc.h" #include "geom_lib.h" #include "miscattr.h" /***************************************************************************** * DESCRIPTION: M * Preprocess a given surface so we can evaluate curvature properties from M * it efficiently, at every point. See SymbEvalCurvature for actual M * curvature at surface point evaluations M * * * PARAMETERS: M * Srf: Surface to preprocess. M * Init: TRUE for initializing, FALSE for clearing out. M * * * RETURN VALUE: M * void M * * * SEE ALSO: M * SymbEvalCurvature M * * * KEYWORDS: M * SymbEvalSrfCurvPrep, curvature M *****************************************************************************/ void SymbEvalSrfCurvPrep(CagdSrfStruct *Srf, CagdBType Init) { CagdSrfStruct **DSrfs; if (Init) { DSrfs = (CagdSrfStruct **) IritMalloc(6 * sizeof(CagdSrfStruct *)); DSrfs[0] = SymbSrfNormalSrf(Srf); DSrfs[1] = CagdSrfDerive(Srf, CAGD_CONST_U_DIR); DSrfs[2] = CagdSrfDerive(Srf, CAGD_CONST_V_DIR); DSrfs[3] = CagdSrfDerive(DSrfs[1], CAGD_CONST_U_DIR); DSrfs[4] = CagdSrfDerive(DSrfs[2], CAGD_CONST_V_DIR); DSrfs[5] = CagdSrfDerive(DSrfs[2], CAGD_CONST_U_DIR); AttrSetPtrAttrib(&Srf -> Attr, "_EvalCurv", DSrfs); } else { int i; if ((DSrfs = (CagdSrfStruct **) AttrGetPtrAttrib(Srf -> Attr, "_EvalCurv")) == NULL) return; /* Release the three preprocessed derivative surfaces. */ for (i = 0; i < 6; i++) CagdSrfFree(DSrfs[i]); IritFree(DSrfs); } } /***************************************************************************** * DESCRIPTION: M * Evaluate a given surface's curvature properties. Returns the principal M * curvatures and directions for the given surface location. M * This function must be invoked after SymbEvalSrfCurvPrep was called to M * initialize the proper data structures, for fast curvature evaluations at M * many points. SymbEvalSrfCurvPrep should be called at the end to release M * these data structures. M * * * PARAMETERS: M * Srf: Surface to evaluate its curvature properties. M * U, V: Location of evaluation. M * DirInUV: If TRUE principal directions are given in UV, otherwise in M * Euclidean 3-space. M * K1, K2: Principal curvatures. M * D1, D2: Principal directions. M * * * RETURN VALUE: M * int: TRUE if succesful, FALSE otherwise. M * * * SEE ALSO: M * SymbEvalSrfCurvPrep M * * * KEYWORDS: M * SymbEvalSrfCurvature, curvature M *****************************************************************************/ int SymbEvalSrfCurvature(CagdSrfStruct *Srf, CagdRType U, CagdRType V, CagdBType DirInUV, CagdRType *K1, CagdRType *K2, CagdVType D1, CagdVType D2) { int RetVal = FALSE; CagdRType *R, E, F, G, L, M, N, A, B, C, Disc, v1, v2; CagdVType Nrml, DSrfU, DSrfV, DSrfUU, DSrfVV, DSrfUV, Vec; CagdSrfStruct **DSrfs; if ((DSrfs = (CagdSrfStruct **) AttrGetPtrAttrib(Srf -> Attr, "_EvalCurv")) == NULL) return RetVal; R = CagdSrfEval(DSrfs[0], U, V); CagdCoerceToE3(Nrml, &R, -1, DSrfs[0] -> PType); PT_NORMALIZE(Nrml); R = CagdSrfEval(DSrfs[1], U, V); CagdCoerceToE3(DSrfU, &R, -1, DSrfs[1] -> PType); R = CagdSrfEval(DSrfs[2], U, V); CagdCoerceToE3(DSrfV, &R, -1, DSrfs[2] -> PType); R = CagdSrfEval(DSrfs[3], U, V); CagdCoerceToE3(DSrfUU, &R, -1, DSrfs[3] -> PType); R = CagdSrfEval(DSrfs[4], U, V); CagdCoerceToE3(DSrfVV, &R, -1, DSrfs[4] -> PType); R = CagdSrfEval(DSrfs[5], U, V); CagdCoerceToE3(DSrfUV, &R, -1, DSrfs[5] -> PType); /* Compute the coefficients of the first and second fundamental forms. */ E = DOT_PROD(DSrfU, DSrfU); F = DOT_PROD(DSrfU, DSrfV); G = DOT_PROD(DSrfV, DSrfV); L = DOT_PROD(DSrfUU, Nrml); M = DOT_PROD(DSrfUV, Nrml); N = DOT_PROD(DSrfVV, Nrml); /* Solve the quadratic equation for the principal curvature: */ /* (EG - F^2) k^2 - (GL + EN - 2FM) k + (LN - M^2) = 0. */ A = E * G - F * F; B = -(G * L + E * N - 2 * F * M); C = L * N - M * M; if ((Disc = B * B - 4 * A * C) < 0.0) { if (Disc < -IRIT_EPS) RetVal = FALSE; Disc = 0.0; } else RetVal = TRUE; Disc = sqrt(Disc); *K1 = (-B - Disc) / (2 * A); *K2 = (-B + Disc) / (2 * A); /* Compute the principal directions - K1. */ A = L - *K1 * E; B = M - *K1 * F; C = N - *K1 * G; if (FABS(A) > FABS(C)) { v1 = B; v2 = -A; } else { v1 = C; v2 = -B; } if (DirInUV) { D1[0] = v1; D1[1] = v2; D1[2] = 0.0; } else { PT_COPY(D1, DSrfU); PT_SCALE(D1, v1); PT_COPY(Vec, DSrfV); PT_SCALE(Vec, v2); PT_ADD(D1, D1, Vec); PT_NORMALIZE(D1); } /* Compute the principal directions - K2. */ A = L - *K2 * E; B = M - *K2 * F; C = N - *K2 * G; if (FABS(A) > FABS(C)) { v1 = B; v2 = -A; } else { v1 = C; v2 = -B; } if (DirInUV) { D2[0] = v1; D2[1] = v2; D2[2] = 0.0; } else { PT_COPY(D2, DSrfU); PT_SCALE(D2, v1); PT_COPY(Vec, DSrfV); PT_SCALE(Vec, v2); PT_ADD(D2, D2, Vec); PT_NORMALIZE(D2); } return RetVal; } /***************************************************************************** * DESCRIPTION: M * Computes the asymptotic directions of the surface in the given U, V M * location, if exists. If DirInUV, the returned asymptotic direction is M * in UV space, otherwise, in Euclidean surface tangent plane space. M * This function must be invoked after SymbEvalSrfCurvPrep was called to M * initialize the proper data structures, for fast curvature evaluations at M * many points. SymbEvalSrfCurvPrep should be called at the end to release M * these data structures. M * The asymptotic direction(s) is the direction for which the normal M * curvature vanish and hence can exist for hyperbolic regions. We solve: M * M * [t (1-t)] [ L M ] [t ] V * [ ] [ ] V * [ M N ] [1-t] V * M * and look for solution of t between zero and one as: M * M * t = (N-M +/- sqrt(M^2-LN)) / (L-2M+N). V * * * PARAMETERS: M * Srf: To compute asymptotic directions for. M * U, V: Location of evaluation. M * DirInUV: If TRUE asymptotic direction is given in UV, otherwise in M * Euclidean 3-space. M * AsympDir1, AsympDir2: Returned values. M * * * RETURN VALUE: M * int: Number of asymptotic directions found, zero for none. M * * * SEE ALSO: M * SymbEvalSrfCurvature, SymbEvalSrfCurvPrep M * * * KEYWORDS: M * SymbEvalSrfAsympDir, curvature M *****************************************************************************/ int SymbEvalSrfAsympDir(CagdSrfStruct *Srf, CagdRType U, CagdRType V, CagdBType DirInUV, CagdVType AsympDir1, CagdVType AsympDir2) { int NumSol = 0; CagdRType *R, L, M, N, Desc, t1, t2; CagdVType Nrml, DSrfU, DSrfV, DSrfUU, DSrfVV, DSrfUV; CagdSrfStruct **DSrfs; if ((DSrfs = (CagdSrfStruct **) AttrGetPtrAttrib(Srf -> Attr, "_EvalCurv")) == NULL) return NumSol; R = CagdSrfEval(DSrfs[0], U, V); CagdCoerceToE3(Nrml, &R, -1, DSrfs[0] -> PType); PT_NORMALIZE(Nrml); R = CagdSrfEval(DSrfs[1], U, V); CagdCoerceToE3(DSrfU, &R, -1, DSrfs[1] -> PType); R = CagdSrfEval(DSrfs[2], U, V); CagdCoerceToE3(DSrfV, &R, -1, DSrfs[2] -> PType); R = CagdSrfEval(DSrfs[3], U, V); CagdCoerceToE3(DSrfUU, &R, -1, DSrfs[3] -> PType); R = CagdSrfEval(DSrfs[4], U, V); CagdCoerceToE3(DSrfVV, &R, -1, DSrfs[4] -> PType); R = CagdSrfEval(DSrfs[5], U, V); CagdCoerceToE3(DSrfUV, &R, -1, DSrfs[5] -> PType); /* Compute the coefficients of the second fundamental forms. */ L = DOT_PROD(DSrfUU, Nrml); M = DOT_PROD(DSrfUV, Nrml); N = DOT_PROD(DSrfVV, Nrml); Desc = SQR(M) - L * N; if (Desc < 0) return NumSol; /* No asymptotic direction. */ t1 = (N - M + sqrt(Desc)) / (L - 2 * M + N); t2 = (N - M - sqrt(Desc)) / (L - 2 * M + N); /* Valid solution. */ AsympDir1[0] = t1; AsympDir1[1] = 1 - t1; AsympDir1[2] = 0; NumSol++; if (!APX_EQ(t1, t2)) { AsympDir2[0] = t2; AsympDir2[1] = 1 - t2; AsympDir2[2] = 0; NumSol++; } if (DirInUV) return NumSol; if (NumSol >= 1) { VEC_BLEND(AsympDir1, DSrfU, DSrfV, t1); VEC_NORMALIZE(AsympDir1); } if (NumSol >= 2) { VEC_BLEND(AsympDir2, DSrfU, DSrfV, t2); VEC_NORMALIZE(AsympDir2); } return NumSol; }