Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-09-28 09:20:13

0001 // Copyright (c) 2025 OPEN CASCADE SAS
0002 //
0003 // This file is part of Open CASCADE Technology software library.
0004 //
0005 // This library is free software; you can redistribute it and/or modify it under
0006 // the terms of the GNU Lesser General Public License version 2.1 as published
0007 // by the Free Software Foundation, with special exception defined in the file
0008 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
0009 // distribution for complete text of the license and disclaimer of any warranty.
0010 //
0011 // Alternatively, this file may be used under the terms of Open CASCADE
0012 // commercial license or contractual agreement.
0013 
0014 #ifndef _GeomLProp_SurfaceUtils_HeaderFile
0015 #define _GeomLProp_SurfaceUtils_HeaderFile
0016 
0017 #include <Adaptor3d_Surface.hxx>
0018 #include <CSLib.hxx>
0019 #include <CSLib_DerivativeStatus.hxx>
0020 #include <Geom_Surface.hxx>
0021 #include <LProp_NotDefined.hxx>
0022 #include <LProp_Status.hxx>
0023 #include <Standard_Handle.hxx>
0024 #include <gp_Dir.hxx>
0025 #include <gp_Pnt.hxx>
0026 #include <gp_Vec.hxx>
0027 #include <math_DirectPolynomialRoots.hxx>
0028 
0029 #include <cmath>
0030 
0031 //! Template utility functions for SLProps surface local property computation.
0032 //! Provides access policies (DirectAccess, ToolAccess), surface bounds helpers,
0033 //! and shared algorithms for evaluating derivatives, tangents, normals, and curvatures.
0034 namespace LProp_SurfaceUtils
0035 {
0036 
0037 // ==================== Deref helpers ====================
0038 
0039 //! Dereference by-value or reference types (identity).
0040 template <typename T>
0041 T& Deref(T& theObj)
0042 {
0043   return theObj;
0044 }
0045 
0046 //! Dereference occ::handle types.
0047 template <typename T>
0048 T& Deref(occ::handle<T>& theHandle)
0049 {
0050   return *theHandle;
0051 }
0052 
0053 //! Dereference const occ::handle types.
0054 template <typename T>
0055 const T& Deref(const occ::handle<T>& theHandle)
0056 {
0057   return *theHandle;
0058 }
0059 
0060 // ==================== Surface Bounds Helpers ====================
0061 
0062 //! Get bounds from Geom_Surface (uses Bounds method with U1, U2, V1, V2 order).
0063 inline void GetSurfBounds(const Geom_Surface& theSurf,
0064                           double&             theU1,
0065                           double&             theV1,
0066                           double&             theU2,
0067                           double&             theV2)
0068 {
0069   theSurf.Bounds(theU1, theU2, theV1, theV2);
0070 }
0071 
0072 //! Get bounds from Adaptor3d_Surface (uses individual parameter methods).
0073 //! Also works for BRepAdaptor_Surface which inherits from Adaptor3d_Surface.
0074 inline void GetSurfBounds(const Adaptor3d_Surface& theSurf,
0075                           double&                  theU1,
0076                           double&                  theV1,
0077                           double&                  theU2,
0078                           double&                  theV2)
0079 {
0080   theU1 = theSurf.FirstUParameter();
0081   theV1 = theSurf.FirstVParameter();
0082   theU2 = theSurf.LastUParameter();
0083   theV2 = theSurf.LastVParameter();
0084 }
0085 
0086 // ==================== Access Policies ====================
0087 
0088 //! Direct access policy: calls D0/D1/D2 methods on the surface object.
0089 //! Works with occ::handle<T> and by-value surface types.
0090 struct DirectAccess
0091 {
0092   template <typename S>
0093   static void D0(S& theSurf, double theU, double theV, gp_Pnt& thePnt)
0094   {
0095     Deref(theSurf).D0(theU, theV, thePnt);
0096   }
0097 
0098   template <typename S>
0099   static void D1(S&      theSurf,
0100                  double  theU,
0101                  double  theV,
0102                  gp_Pnt& thePnt,
0103                  gp_Vec& theD1u,
0104                  gp_Vec& theD1v)
0105   {
0106     Deref(theSurf).D1(theU, theV, thePnt, theD1u, theD1v);
0107   }
0108 
0109   template <typename S>
0110   static void D2(S&      theSurf,
0111                  double  theU,
0112                  double  theV,
0113                  gp_Pnt& thePnt,
0114                  gp_Vec& theD1u,
0115                  gp_Vec& theD1v,
0116                  gp_Vec& theD2u,
0117                  gp_Vec& theD2v,
0118                  gp_Vec& theDuv)
0119   {
0120     Deref(theSurf).D2(theU, theV, thePnt, theD1u, theD1v, theD2u, theD2v, theDuv);
0121   }
0122 
0123   template <typename S>
0124   static void Bounds(S& theSurf, double& theU1, double& theV1, double& theU2, double& theV2)
0125   {
0126     GetSurfBounds(Deref(theSurf), theU1, theV1, theU2, theV2);
0127   }
0128 };
0129 
0130 //! Tool-based access policy: delegates to static Tool methods.
0131 //! Used for HLRBRep types where Tool class provides the interface.
0132 template <typename Tool>
0133 struct ToolAccess
0134 {
0135   template <typename S>
0136   static void D0(S& theSurf, double theU, double theV, gp_Pnt& thePnt)
0137   {
0138     Tool::Value(theSurf, theU, theV, thePnt);
0139   }
0140 
0141   template <typename S>
0142   static void D1(S&      theSurf,
0143                  double  theU,
0144                  double  theV,
0145                  gp_Pnt& thePnt,
0146                  gp_Vec& theD1u,
0147                  gp_Vec& theD1v)
0148   {
0149     Tool::D1(theSurf, theU, theV, thePnt, theD1u, theD1v);
0150   }
0151 
0152   template <typename S>
0153   static void D2(S&      theSurf,
0154                  double  theU,
0155                  double  theV,
0156                  gp_Pnt& thePnt,
0157                  gp_Vec& theD1u,
0158                  gp_Vec& theD1v,
0159                  gp_Vec& theD2u,
0160                  gp_Vec& theD2v,
0161                  gp_Vec& theDuv)
0162   {
0163     Tool::D2(theSurf, theU, theV, thePnt, theD1u, theD1v, theD2u, theD2v, theDuv);
0164   }
0165 
0166   template <typename S>
0167   static void Bounds(S& theSurf, double& theU1, double& theV1, double& theU2, double& theV2)
0168   {
0169     Tool::Bounds(theSurf, theU1, theV1, theU2, theV2);
0170   }
0171 };
0172 
0173 // ==================== Algorithm Utilities ====================
0174 
0175 //! Evaluate surface derivatives at parameters (theU, theV) up to specified order.
0176 //! @param[in,out] theSurf   surface object
0177 //! @param[in]     theU      U parameter value
0178 //! @param[in]     theV      V parameter value
0179 //! @param[in]     theOrder  maximum derivative order (0-2)
0180 //! @param[out]    thePnt    evaluated point
0181 //! @param[out]    theD1u    first U derivative
0182 //! @param[out]    theD1v    first V derivative
0183 //! @param[out]    theD2u    second U derivative
0184 //! @param[out]    theD2v    second V derivative
0185 //! @param[out]    theDuv    mixed UV derivative
0186 template <typename Access, typename Surface>
0187 void EvalSurfDerivatives(Surface& theSurf,
0188                          double   theU,
0189                          double   theV,
0190                          int      theOrder,
0191                          gp_Pnt&  thePnt,
0192                          gp_Vec&  theD1u,
0193                          gp_Vec&  theD1v,
0194                          gp_Vec&  theD2u,
0195                          gp_Vec&  theD2v,
0196                          gp_Vec&  theDuv)
0197 {
0198   switch (theOrder)
0199   {
0200     case 0:
0201       Access::D0(theSurf, theU, theV, thePnt);
0202       break;
0203     case 1:
0204       Access::D1(theSurf, theU, theV, thePnt, theD1u, theD1v);
0205       break;
0206     case 2:
0207       Access::D2(theSurf, theU, theV, thePnt, theD1u, theD1v, theD2u, theD2v, theDuv);
0208       break;
0209   }
0210 }
0211 
0212 //! Check if the surface tangent is defined for a given derivative direction.
0213 //! Searches for the first non-null derivative in either U or V direction.
0214 //! @param[in]  theD1        first derivative (D1U for U direction, D1V for V direction)
0215 //! @param[in]  theD2        second derivative (D2U for U direction, D2V for V direction)
0216 //! @param[in]  theCN        continuity order
0217 //! @param[in]  theTolSq     squared linear tolerance
0218 //! @param[out] theOrder     order of first significant derivative
0219 //! @param[out] theStatus    resulting tangent status
0220 //! @return true if tangent is defined
0221 inline bool FindSurfTangentOrder(const gp_Vec& theD1,
0222                                  const gp_Vec& theD2,
0223                                  int           theCN,
0224                                  double        theTolSq,
0225                                  int&          theOrder,
0226                                  LProp_Status& theStatus)
0227 {
0228   const gp_Vec* aDerivs[2] = {&theD1, &theD2};
0229   theOrder                 = 0;
0230 
0231   while (theOrder < 3)
0232   {
0233     theOrder++;
0234     if (theCN >= theOrder)
0235     {
0236       if (theOrder <= 2 && aDerivs[theOrder - 1]->SquareMagnitude() > theTolSq)
0237       {
0238         theStatus = LProp_Defined;
0239         return true;
0240       }
0241     }
0242     else
0243     {
0244       theStatus = LProp_Undefined;
0245       return false;
0246     }
0247   }
0248 
0249   theStatus = LProp_Undefined;
0250   return false;
0251 }
0252 
0253 //! Compute surface tangent direction with sign correction for higher-order derivatives.
0254 //! When the first significant derivative has order > 1, the sign is determined
0255 //! by comparing with the chord direction near the point.
0256 //! @param[in,out] theSurf       surface object
0257 //! @param[in]     theU          current U parameter
0258 //! @param[in]     theV          current V parameter
0259 //! @param[in]     theFirstDeriv first derivative in this direction (D1U or D1V)
0260 //! @param[in]     theSecDeriv   second derivative in this direction (D2U or D2V)
0261 //! @param[in]     theSigOrder   order of first significant derivative
0262 //! @param[in]     theIsU        true for U direction, false for V direction
0263 //! @param[out]    theDir        computed tangent direction
0264 template <typename Access, typename Surface>
0265 void ComputeSurfTangent(Surface&      theSurf,
0266                         double        theU,
0267                         double        theV,
0268                         const gp_Vec& theFirstDeriv,
0269                         const gp_Vec& theSecDeriv,
0270                         int           theSigOrder,
0271                         bool          theIsU,
0272                         gp_Dir&       theDir)
0273 {
0274   if (theSigOrder == 1)
0275   {
0276     theDir = gp_Dir(theFirstDeriv);
0277     return;
0278   }
0279 
0280   constexpr double THE_DIVISION_FACTOR = 1.0e-3;
0281   constexpr double THE_MIN_STEP        = 1.0e-7;
0282 
0283   double anUinfimum, anVinfimum, anUsupremum, anVsupremum;
0284   Access::Bounds(theSurf, anUinfimum, anVinfimum, anUsupremum, anVsupremum);
0285 
0286   if (theIsU)
0287   {
0288     double aDu;
0289     if ((anUsupremum >= RealLast()) || (anUinfimum <= RealFirst()))
0290       aDu = 0.0;
0291     else
0292       aDu = anUsupremum - anUinfimum;
0293 
0294     const double aDelta = std::max(aDu * THE_DIVISION_FACTOR, THE_MIN_STEP);
0295 
0296     gp_Vec aV = theSecDeriv;
0297 
0298     double anOtherU;
0299     if (theU - anUinfimum < aDelta)
0300       anOtherU = theU + aDelta;
0301     else
0302       anOtherU = theU - aDelta;
0303 
0304     gp_Pnt aP1, aP2;
0305     Access::D0(theSurf, std::min(theU, anOtherU), theV, aP1);
0306     Access::D0(theSurf, std::max(theU, anOtherU), theV, aP2);
0307 
0308     gp_Vec aChord(aP1, aP2);
0309     if (aV.Dot(aChord) < 0.0)
0310       aV = -aV;
0311 
0312     theDir = gp_Dir(aV);
0313   }
0314   else
0315   {
0316     double aDv;
0317     if ((anVsupremum >= RealLast()) || (anVinfimum <= RealFirst()))
0318       aDv = 0.0;
0319     else
0320       aDv = anVsupremum - anVinfimum;
0321 
0322     const double aDelta = std::max(aDv * THE_DIVISION_FACTOR, THE_MIN_STEP);
0323 
0324     gp_Vec aV = theSecDeriv;
0325 
0326     double anOtherV;
0327     if (theV - anVinfimum < aDelta)
0328       anOtherV = theV + aDelta;
0329     else
0330       anOtherV = theV - aDelta;
0331 
0332     gp_Pnt aP1, aP2;
0333     Access::D0(theSurf, theU, std::min(theV, anOtherV), aP1);
0334     Access::D0(theSurf, theU, std::max(theV, anOtherV), aP2);
0335 
0336     gp_Vec aChord(aP1, aP2);
0337     if (aV.Dot(aChord) < 0.0)
0338       aV = -aV;
0339 
0340     theDir = gp_Dir(aV);
0341   }
0342 }
0343 
0344 //! Check if surface normal is defined, and compute it via CSLib::Normal.
0345 //! @param[in]  theD1u     first U derivative
0346 //! @param[in]  theD1v     first V derivative
0347 //! @param[in]  theLinTol  linear tolerance
0348 //! @param[out] theNormal  computed normal direction (if defined)
0349 //! @return true if normal is defined
0350 inline bool ComputeSurfNormal(const gp_Vec& theD1u,
0351                               const gp_Vec& theD1v,
0352                               double        theLinTol,
0353                               gp_Dir&       theNormal)
0354 {
0355   CSLib_DerivativeStatus aStatus = CSLib_Done;
0356   CSLib::Normal(theD1u, theD1v, theLinTol, aStatus, theNormal);
0357   return aStatus == CSLib_Done;
0358 }
0359 
0360 //! Compute principal curvatures and directions via fundamental forms.
0361 //! Solves the eigenvalue problem for the shape operator using
0362 //! first and second fundamental form coefficients.
0363 //! @param[in]  theD1u      first U derivative
0364 //! @param[in]  theD1v      first V derivative
0365 //! @param[in]  theD2u      second U derivative
0366 //! @param[in]  theD2v      second V derivative
0367 //! @param[in]  theDuv      mixed UV derivative
0368 //! @param[in]  theNormal   surface normal direction
0369 //! @param[out] theMinCurv  minimum principal curvature
0370 //! @param[out] theMaxCurv  maximum principal curvature
0371 //! @param[out] theDirMin   direction of minimum curvature
0372 //! @param[out] theDirMax   direction of maximum curvature
0373 //! @param[out] theMeanCurv mean curvature
0374 //! @param[out] theGausCurv Gaussian curvature
0375 //! @return true if curvatures are successfully computed
0376 inline bool ComputeSurfCurvatures(const gp_Vec& theD1u,
0377                                   const gp_Vec& theD1v,
0378                                   const gp_Vec& theD2u,
0379                                   const gp_Vec& theD2v,
0380                                   const gp_Vec& theDuv,
0381                                   const gp_Dir& theNormal,
0382                                   double&       theMinCurv,
0383                                   double&       theMaxCurv,
0384                                   gp_Dir&       theDirMin,
0385                                   gp_Dir&       theDirMax,
0386                                   double&       theMeanCurv,
0387                                   double&       theGausCurv)
0388 {
0389   const gp_Vec aNorm(theNormal);
0390 
0391   const double anE = theD1u.SquareMagnitude();
0392   const double anF = theD1u.Dot(theD1v);
0393   const double aG  = theD1v.SquareMagnitude();
0394 
0395   const double aL = aNorm.Dot(theD2u);
0396   const double aM = aNorm.Dot(theDuv);
0397   const double aN = aNorm.Dot(theD2v);
0398 
0399   const double anA0 = anE * aM - anF * aL;
0400   const double aB0  = anE * aN - aG * aL;
0401   const double aC0  = anF * aN - aG * aM;
0402 
0403   const double aMaxABC = std::max(std::max(std::abs(anA0), std::abs(aB0)), std::abs(aC0));
0404   if (aMaxABC < RealEpsilon())
0405   {
0406     // Umbilic point
0407     if (aG < RealEpsilon())
0408     {
0409       return false;
0410     }
0411     theMinCurv  = aN / aG;
0412     theMaxCurv  = theMinCurv;
0413     theDirMin   = gp_Dir(theD1u);
0414     theDirMax   = gp_Dir(theD1u.Crossed(aNorm));
0415     theMeanCurv = theMinCurv;
0416     theGausCurv = theMinCurv * theMinCurv;
0417     return true;
0418   }
0419 
0420   const double anA = anA0 / aMaxABC;
0421   const double aB  = aB0 / aMaxABC;
0422   const double aC  = aC0 / aMaxABC;
0423 
0424   double aCurv1, aCurv2;
0425   gp_Vec aVectCurv1, aVectCurv2;
0426 
0427   if (std::abs(anA) > RealEpsilon())
0428   {
0429     math_DirectPolynomialRoots aRoot(anA, aB, aC);
0430     if (aRoot.NbSolutions() != 2)
0431       return false;
0432 
0433     const double aRoot1 = aRoot.Value(1);
0434     const double aRoot2 = aRoot.Value(2);
0435     aCurv1 = ((aL * aRoot1 + 2.0 * aM) * aRoot1 + aN) / ((anE * aRoot1 + 2.0 * anF) * aRoot1 + aG);
0436     aCurv2 = ((aL * aRoot2 + 2.0 * aM) * aRoot2 + aN) / ((anE * aRoot2 + 2.0 * anF) * aRoot2 + aG);
0437     aVectCurv1 = aRoot1 * theD1u + theD1v;
0438     aVectCurv2 = aRoot2 * theD1u + theD1v;
0439   }
0440   else if (std::abs(aC) > RealEpsilon())
0441   {
0442     math_DirectPolynomialRoots aRoot(aC, aB, anA);
0443     if (aRoot.NbSolutions() != 2)
0444       return false;
0445 
0446     const double aRoot1 = aRoot.Value(1);
0447     const double aRoot2 = aRoot.Value(2);
0448     aCurv1 = ((aN * aRoot1 + 2.0 * aM) * aRoot1 + aL) / ((aG * aRoot1 + 2.0 * anF) * aRoot1 + anE);
0449     aCurv2 = ((aN * aRoot2 + 2.0 * aM) * aRoot2 + aL) / ((aG * aRoot2 + 2.0 * anF) * aRoot2 + anE);
0450     aVectCurv1 = theD1u + aRoot1 * theD1v;
0451     aVectCurv2 = theD1u + aRoot2 * theD1v;
0452   }
0453   else
0454   {
0455     aCurv1     = aL / anE;
0456     aCurv2     = aN / aG;
0457     aVectCurv1 = theD1u;
0458     aVectCurv2 = theD1v;
0459   }
0460 
0461   if (aCurv1 < aCurv2)
0462   {
0463     theMinCurv = aCurv1;
0464     theMaxCurv = aCurv2;
0465     theDirMin  = gp_Dir(aVectCurv1);
0466     theDirMax  = gp_Dir(aVectCurv2);
0467   }
0468   else
0469   {
0470     theMinCurv = aCurv2;
0471     theMaxCurv = aCurv1;
0472     theDirMin  = gp_Dir(aVectCurv2);
0473     theDirMax  = gp_Dir(aVectCurv1);
0474   }
0475 
0476   const double anEG_FF = (anE * aG) - (anF * anF);
0477   theMeanCurv          = ((aN * anE) - (2.0 * aM * anF) + (aL * aG)) / (2.0 * anEG_FF);
0478   theGausCurv          = ((aL * aN) - (aM * aM)) / anEG_FF;
0479   return true;
0480 }
0481 
0482 // ==================== Higher-Level Method Wrappers ====================
0483 
0484 //! SetParameters: sets parameters, evaluates derivatives, resets all statuses.
0485 template <typename Access, typename Surface>
0486 void SetParameters(Surface&      theSurf,
0487                    double        theU,
0488                    double        theV,
0489                    double&       theStoredU,
0490                    double&       theStoredV,
0491                    int           theDerOrder,
0492                    gp_Pnt&       thePnt,
0493                    gp_Vec&       theD1u,
0494                    gp_Vec&       theD1v,
0495                    gp_Vec&       theD2u,
0496                    gp_Vec&       theD2v,
0497                    gp_Vec&       theDuv,
0498                    LProp_Status& theUTanSt,
0499                    LProp_Status& theVTanSt,
0500                    LProp_Status& theNormSt,
0501                    LProp_Status& theCurvSt)
0502 {
0503   theStoredU = theU;
0504   theStoredV = theV;
0505   EvalSurfDerivatives<
0506     Access>(theSurf, theU, theV, theDerOrder, thePnt, theD1u, theD1v, theD2u, theD2v, theDuv);
0507   theUTanSt = LProp_Undecided;
0508   theVTanSt = LProp_Undecided;
0509   theNormSt = LProp_Undecided;
0510   theCurvSt = LProp_Undecided;
0511 }
0512 
0513 //! Ensure surface derivatives up to the required order. Returns the specified result field.
0514 template <typename Access, typename Surface>
0515 const gp_Vec& EnsureSurfDeriv(Surface&      theSurf,
0516                               double        theU,
0517                               double        theV,
0518                               int&          theDerOrder,
0519                               int           theRequired,
0520                               gp_Pnt&       thePnt,
0521                               gp_Vec&       theD1u,
0522                               gp_Vec&       theD1v,
0523                               gp_Vec&       theD2u,
0524                               gp_Vec&       theD2v,
0525                               gp_Vec&       theDuv,
0526                               const gp_Vec& theResult)
0527 {
0528   if (theDerOrder < theRequired)
0529   {
0530     theDerOrder = theRequired;
0531     EvalSurfDerivatives<
0532       Access>(theSurf, theU, theV, theDerOrder, thePnt, theD1u, theD1v, theD2u, theD2v, theDuv);
0533   }
0534   return theResult;
0535 }
0536 
0537 //! IsTangentUDefined: checks U tangent status, then searches via D1U/D2U.
0538 template <typename Props>
0539 bool IsTangentUDefined(Props&        theProps,
0540                        int           theCN,
0541                        double        theLinTol,
0542                        int&          theSigOrder,
0543                        LProp_Status& theTanStatus)
0544 {
0545   if (theTanStatus == LProp_Undefined)
0546     return false;
0547   if (theTanStatus >= LProp_Defined)
0548     return true;
0549   return FindSurfTangentOrder(theProps.D1U(),
0550                               theProps.D2U(),
0551                               theCN,
0552                               theLinTol * theLinTol,
0553                               theSigOrder,
0554                               theTanStatus);
0555 }
0556 
0557 //! IsTangentVDefined: checks V tangent status, then searches via D1V/D2V.
0558 template <typename Props>
0559 bool IsTangentVDefined(Props&        theProps,
0560                        int           theCN,
0561                        double        theLinTol,
0562                        int&          theSigOrder,
0563                        LProp_Status& theTanStatus)
0564 {
0565   if (theTanStatus == LProp_Undefined)
0566     return false;
0567   if (theTanStatus >= LProp_Defined)
0568     return true;
0569   return FindSurfTangentOrder(theProps.D1V(),
0570                               theProps.D2V(),
0571                               theCN,
0572                               theLinTol * theLinTol,
0573                               theSigOrder,
0574                               theTanStatus);
0575 }
0576 
0577 //! TangentU: checks IsTangentUDefined, then computes U tangent direction.
0578 template <typename Access, typename Props, typename Surface>
0579 void TangentU(Props&        theProps,
0580               Surface&      theSurf,
0581               double        theU,
0582               double        theV,
0583               const gp_Vec& theD1u,
0584               const gp_Vec& theD2u,
0585               int           theSigOrder,
0586               gp_Dir&       theDir)
0587 {
0588   if (!theProps.IsTangentUDefined())
0589     throw LProp_NotDefined();
0590   ComputeSurfTangent<Access>(theSurf, theU, theV, theD1u, theD2u, theSigOrder, true, theDir);
0591 }
0592 
0593 //! TangentV: checks IsTangentVDefined, then computes V tangent direction.
0594 template <typename Access, typename Props, typename Surface>
0595 void TangentV(Props&        theProps,
0596               Surface&      theSurf,
0597               double        theU,
0598               double        theV,
0599               const gp_Vec& theD1v,
0600               const gp_Vec& theD2v,
0601               int           theSigOrder,
0602               gp_Dir&       theDir)
0603 {
0604   if (!theProps.IsTangentVDefined())
0605     throw LProp_NotDefined();
0606   ComputeSurfTangent<Access>(theSurf, theU, theV, theD1v, theD2v, theSigOrder, false, theDir);
0607 }
0608 
0609 //! IsNormalDefined: checks normal status, then computes via CSLib.
0610 inline bool IsNormalDefined(const gp_Vec& theD1u,
0611                             const gp_Vec& theD1v,
0612                             double        theLinTol,
0613                             gp_Dir&       theNormal,
0614                             LProp_Status& theNormStatus)
0615 {
0616   if (theNormStatus == LProp_Undefined)
0617     return false;
0618   if (theNormStatus >= LProp_Defined)
0619     return true;
0620   if (ComputeSurfNormal(theD1u, theD1v, theLinTol, theNormal))
0621   {
0622     theNormStatus = LProp_Computed;
0623     return true;
0624   }
0625   theNormStatus = LProp_Undefined;
0626   return false;
0627 }
0628 
0629 //! Normal: checks IsNormalDefined (via theProps), then returns normal.
0630 template <typename Props>
0631 const gp_Dir& Normal(Props& theProps, const gp_Dir& theNormal)
0632 {
0633   if (!theProps.IsNormalDefined())
0634     throw LProp_NotDefined();
0635   return theNormal;
0636 }
0637 
0638 //! IsCurvatureDefined: full check including continuity, normal, tangents, and curvature
0639 //! computation.
0640 template <typename Props>
0641 bool IsCurvatureDefined(Props&        theProps,
0642                         int           theCN,
0643                         int&          theDerOrder,
0644                         const gp_Vec& theD1u,
0645                         const gp_Vec& theD1v,
0646                         const gp_Vec& theD2u,
0647                         const gp_Vec& theD2v,
0648                         const gp_Vec& theDuv,
0649                         const gp_Dir& theNormal,
0650                         double&       theMinCurv,
0651                         double&       theMaxCurv,
0652                         gp_Dir&       theDirMin,
0653                         gp_Dir&       theDirMax,
0654                         double&       theMeanCurv,
0655                         double&       theGausCurv,
0656                         LProp_Status& theCurvStatus)
0657 {
0658   if (theCurvStatus == LProp_Undefined)
0659     return false;
0660   if (theCurvStatus >= LProp_Defined)
0661     return true;
0662   if (theCN < 2)
0663   {
0664     theCurvStatus = LProp_Undefined;
0665     return false;
0666   }
0667   if (!theProps.IsNormalDefined())
0668   {
0669     theCurvStatus = LProp_Undefined;
0670     return false;
0671   }
0672   if (!theProps.IsTangentUDefined() || !theProps.IsTangentVDefined())
0673   {
0674     theCurvStatus = LProp_Undefined;
0675     return false;
0676   }
0677   if (theDerOrder < 2)
0678     theProps.D2U();
0679   if (ComputeSurfCurvatures(theD1u,
0680                             theD1v,
0681                             theD2u,
0682                             theD2v,
0683                             theDuv,
0684                             theNormal,
0685                             theMinCurv,
0686                             theMaxCurv,
0687                             theDirMin,
0688                             theDirMax,
0689                             theMeanCurv,
0690                             theGausCurv))
0691   {
0692     theCurvStatus = LProp_Computed;
0693     return true;
0694   }
0695   theCurvStatus = LProp_Undefined;
0696   return false;
0697 }
0698 
0699 //! Require curvature to be defined, throw if not. Returns the given value.
0700 template <typename Props>
0701 double RequireCurvature(Props& theProps, double theValue)
0702 {
0703   if (!theProps.IsCurvatureDefined())
0704     throw LProp_NotDefined();
0705   return theValue;
0706 }
0707 
0708 //! IsUmbilic: checks curvature, then compares min/max.
0709 template <typename Props>
0710 bool IsUmbilic(Props& theProps, double theMaxCurv, double theMinCurv)
0711 {
0712   if (!theProps.IsCurvatureDefined())
0713     throw LProp_NotDefined();
0714   return std::abs(theMaxCurv - theMinCurv) < std::abs(Epsilon(theMaxCurv));
0715 }
0716 
0717 //! CurvatureDirections: checks curvature, then returns max/min directions.
0718 template <typename Props>
0719 void CurvatureDirections(Props&        theProps,
0720                          const gp_Dir& theDirMax,
0721                          const gp_Dir& theDirMin,
0722                          gp_Dir&       theMax,
0723                          gp_Dir&       theMin)
0724 {
0725   if (!theProps.IsCurvatureDefined())
0726     throw LProp_NotDefined();
0727   theMax = theDirMax;
0728   theMin = theDirMin;
0729 }
0730 
0731 } // namespace LProp_SurfaceUtils
0732 
0733 #endif // _GeomLProp_SurfaceUtils_HeaderFile