File indexing completed on 2026-09-28 09:20:13
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
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
0032
0033
0034 namespace LProp_SurfaceUtils
0035 {
0036
0037
0038
0039
0040 template <typename T>
0041 T& Deref(T& theObj)
0042 {
0043 return theObj;
0044 }
0045
0046
0047 template <typename T>
0048 T& Deref(occ::handle<T>& theHandle)
0049 {
0050 return *theHandle;
0051 }
0052
0053
0054 template <typename T>
0055 const T& Deref(const occ::handle<T>& theHandle)
0056 {
0057 return *theHandle;
0058 }
0059
0060
0061
0062
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
0073
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
0087
0088
0089
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
0131
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
0174
0175
0176
0177
0178
0179
0180
0181
0182
0183
0184
0185
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
0213
0214
0215
0216
0217
0218
0219
0220
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
0254
0255
0256
0257
0258
0259
0260
0261
0262
0263
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
0345
0346
0347
0348
0349
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
0361
0362
0363
0364
0365
0366
0367
0368
0369
0370
0371
0372
0373
0374
0375
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
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
0483
0484
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
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
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
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
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
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
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
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
0639
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
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
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
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 }
0732
0733 #endif