Back to home page

EIC code displayed by LXR
 
 

    

Warning, file /include/opencascade/Extrema_GFuncExtPC.hxx was not indexed or was modified since last indexation (in which case cross-reference links may be missing, inaccurate or erroneous).

0001 // Created on: 1991-02-26
0002 // Created by: Isabelle GRIGNON
0003 // Copyright (c) 1991-1999 Matra Datavision
0004 // Copyright (c) 1999-2014 OPEN CASCADE SAS
0005 //
0006 // This file is part of Open CASCADE Technology software library.
0007 //
0008 // This library is free software; you can redistribute it and/or modify it under
0009 // the terms of the GNU Lesser General Public License version 2.1 as published
0010 // by the Free Software Foundation, with special exception defined in the file
0011 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
0012 // distribution for complete text of the license and disclaimer of any warranty.
0013 //
0014 // Alternatively, this file may be used under the terms of Open CASCADE
0015 // commercial license or contractual agreement.
0016 
0017 #ifndef _Extrema_GFuncExtPC_HeaderFile
0018 #define _Extrema_GFuncExtPC_HeaderFile
0019 
0020 #include <Standard.hxx>
0021 #include <Standard_DefineAlloc.hxx>
0022 
0023 #include <GeomAbs_CurveType.hxx>
0024 #include <math_FunctionWithDerivative.hxx>
0025 #include <Precision.hxx>
0026 #include <Standard_TypeMismatch.hxx>
0027 #include <Standard_Integer.hxx>
0028 #include <NCollection_Sequence.hxx>
0029 
0030 #include <cmath>
0031 
0032 //! Template class for computing extremal distance function between a point and a curve.
0033 //! Searches for a parameter value u such that dist(P, C(u)) passes through an extremum.
0034 //! Inherits from math_FunctionWithDerivative and is used by math_FunctionRoot and
0035 //! math_FunctionRoots algorithms.
0036 //!
0037 //! If D1c and D2c are the first and second derivatives:
0038 //! F(u) = (C(u)-P).D1c(u) / ||D1c||
0039 //! DF(u) = ||D1c|| + (C(u)-P).D2c(u)/||D1c|| - F(u)*D2c.D1c/||D1c||^2
0040 //!
0041 //! @tparam TheCurve    Curve type (e.g., Adaptor3d_Curve, Adaptor2d_Curve2d)
0042 //! @tparam TheCurveTool Tool for curve operations
0043 //! @tparam ThePOnC     Point on curve type
0044 //! @tparam ThePoint    Point type (e.g., gp_Pnt, gp_Pnt2d)
0045 //! @tparam TheVector   Vector type (e.g., gp_Vec, gp_Vec2d)
0046 //! @tparam TheSeqPOnC  Sequence of points on curve type
0047 template <typename TheCurve,
0048           typename TheCurveTool,
0049           typename ThePOnC,
0050           typename ThePoint,
0051           typename TheVector,
0052           typename TheSeqPOnC>
0053 class Extrema_GFuncExtPC : public math_FunctionWithDerivative
0054 {
0055 public:
0056   DEFINE_STANDARD_ALLOC
0057 
0058   //! Default constructor.
0059   Extrema_GFuncExtPC()
0060       : myC(nullptr),
0061         myU(0.0),
0062         myD1f(0.0),
0063         myPinit(false),
0064         myCinit(false),
0065         myD1Init(false),
0066         myTol(MinTol),
0067         myMaxDerivOrder(0),
0068         myUinfium(0.0),
0069         myUsupremum(0.0)
0070   {
0071   }
0072 
0073   //! Constructor with point and curve initialization.
0074   //! @param theP Point to compute distance from
0075   //! @param theC Curve to compute distance to
0076   Extrema_GFuncExtPC(const ThePoint& theP, const TheCurve& theC)
0077       : myP(theP),
0078         myC(const_cast<TheCurve*>(&theC)),
0079         myU(0.0),
0080         myD1f(0.0),
0081         myPinit(true),
0082         myCinit(true),
0083         myD1Init(false)
0084   {
0085     SubIntervalInitialize(TheCurveTool::FirstParameter(theC), TheCurveTool::LastParameter(theC));
0086 
0087     switch (TheCurveTool::GetType(theC))
0088     {
0089       case GeomAbs_BezierCurve:
0090       case GeomAbs_BSplineCurve:
0091       case GeomAbs_OffsetCurve:
0092       case GeomAbs_OtherCurve:
0093         myMaxDerivOrder = MaxOrder;
0094         myTol           = SearchOfTolerance();
0095         break;
0096       default:
0097         myMaxDerivOrder = 0;
0098         myTol           = MinTol;
0099         break;
0100     }
0101   }
0102 
0103   //! Sets the curve field.
0104   //! @param theC Curve to set
0105   void Initialize(const TheCurve& theC)
0106   {
0107     myC     = const_cast<TheCurve*>(&theC);
0108     myCinit = true;
0109     myPoint.Clear();
0110     mySqDist.Clear();
0111     myIsMin.Clear();
0112 
0113     SubIntervalInitialize(TheCurveTool::FirstParameter(theC), TheCurveTool::LastParameter(theC));
0114 
0115     switch (TheCurveTool::GetType(theC))
0116     {
0117       case GeomAbs_BezierCurve:
0118       case GeomAbs_BSplineCurve:
0119       case GeomAbs_OffsetCurve:
0120       case GeomAbs_OtherCurve:
0121         myMaxDerivOrder = MaxOrder;
0122         myTol           = SearchOfTolerance();
0123         break;
0124       default:
0125         myMaxDerivOrder = 0;
0126         myTol           = MinTol;
0127         break;
0128     }
0129   }
0130 
0131   //! Sets the point field.
0132   //! @param theP Point to set
0133   void SetPoint(const ThePoint& theP)
0134   {
0135     myP     = theP;
0136     myPinit = true;
0137     myPoint.Clear();
0138     mySqDist.Clear();
0139     myIsMin.Clear();
0140   }
0141 
0142   //! Calculation of F(u).
0143   //! @param theU Parameter value
0144   //! @param theF Output function value
0145   //! @return True if computation succeeded
0146   bool Value(const double theU, double& theF) override
0147   {
0148     if (!myPinit || !myCinit)
0149     {
0150       throw Standard_TypeMismatch("No init");
0151     }
0152 
0153     myU = theU;
0154     TheVector D1c;
0155     TheCurveTool::D1(*myC, myU, myPc, D1c);
0156 
0157     if (Precision::IsInfinite(D1c.X()) || Precision::IsInfinite(D1c.Y()))
0158     {
0159       theF = Precision::Infinite();
0160       return false;
0161     }
0162 
0163     double Ndu = D1c.Magnitude();
0164 
0165     if (myMaxDerivOrder != 0)
0166     {
0167       if (Ndu <= myTol) // Singular case
0168       {
0169         const double DivisionFactor = 1.e-3;
0170         double       du;
0171         if ((myUsupremum >= RealLast()) || (myUinfium <= RealFirst()))
0172           du = 0.0;
0173         else
0174           du = myUsupremum - myUinfium;
0175 
0176         const double aDelta = std::max(du * DivisionFactor, MinStep);
0177         // Derivative is approximated by Taylor-series
0178 
0179         int       n = 1; // Derivative order
0180         TheVector V;
0181         bool      IsDeriveFound;
0182 
0183         do
0184         {
0185           V             = TheCurveTool::DN(*myC, myU, ++n);
0186           Ndu           = V.Magnitude();
0187           IsDeriveFound = (Ndu > myTol);
0188         } while (!IsDeriveFound && n < myMaxDerivOrder);
0189 
0190         if (IsDeriveFound)
0191         {
0192           double u;
0193 
0194           if (myU - myUinfium < aDelta)
0195             u = myU + aDelta;
0196           else
0197             u = myU - aDelta;
0198 
0199           ThePoint P1, P2;
0200           TheCurveTool::D0(*myC, std::min(myU, u), P1);
0201           TheCurveTool::D0(*myC, std::max(myU, u), P2);
0202 
0203           TheVector V1(P1, P2);
0204           double    aDirFactor = V.Dot(V1);
0205 
0206           if (aDirFactor < 0.0)
0207             D1c = -V;
0208           else
0209             D1c = V;
0210         }
0211         else
0212         {
0213           // Derivative is approximated by three points
0214           ThePoint Ptemp;
0215           ThePoint P1, P2, P3;
0216           bool     IsParameterGrown;
0217 
0218           if (myU - myUinfium < 2 * aDelta)
0219           {
0220             TheCurveTool::D0(*myC, myU, P1);
0221             TheCurveTool::D0(*myC, myU + aDelta, P2);
0222             TheCurveTool::D0(*myC, myU + 2 * aDelta, P3);
0223             IsParameterGrown = true;
0224           }
0225           else
0226           {
0227             TheCurveTool::D0(*myC, myU - 2 * aDelta, P1);
0228             TheCurveTool::D0(*myC, myU - aDelta, P2);
0229             TheCurveTool::D0(*myC, myU, P3);
0230             IsParameterGrown = false;
0231           }
0232 
0233           TheVector V1(Ptemp, P1), V2(Ptemp, P2), V3(Ptemp, P3);
0234 
0235           if (IsParameterGrown)
0236             D1c = -3 * V1 + 4 * V2 - V3;
0237           else
0238             D1c = V1 - 4 * V2 + 3 * V3;
0239         }
0240         Ndu = D1c.Magnitude();
0241       }
0242     }
0243 
0244     if (Ndu <= MinTol)
0245     {
0246       // Warning: 1st derivative is equal to zero!
0247       return false;
0248     }
0249 
0250     TheVector PPc(myP, myPc);
0251     theF = PPc.Dot(D1c) / Ndu;
0252     return true;
0253   }
0254 
0255   //! Calculation of F'(u).
0256   //! @param theU Parameter value
0257   //! @param theDF Output derivative value
0258   //! @return True if computation succeeded
0259   bool Derivative(const double theU, double& theDF) override
0260   {
0261     if (!myPinit || !myCinit)
0262     {
0263       throw Standard_TypeMismatch();
0264     }
0265     double F;
0266     return Values(theU, F, theDF);
0267   }
0268 
0269   //! Calculation of F(u) and F'(u).
0270   //! @param theU Parameter value
0271   //! @param theF Output function value
0272   //! @param theDF Output derivative value
0273   //! @return True if computation succeeded
0274   bool Values(const double theU, double& theF, double& theDF) override
0275   {
0276     if (!myPinit || !myCinit)
0277     {
0278       throw Standard_TypeMismatch("No init");
0279     }
0280 
0281     ThePoint myPc_old = myPc, myP_old = myP;
0282 
0283     if (Value(theU, theF) == false)
0284     {
0285       myD1Init = false;
0286       return false;
0287     }
0288 
0289     myU  = theU;
0290     myPc = myPc_old;
0291     myP  = myP_old;
0292 
0293     TheVector D1c, D2c;
0294     TheCurveTool::D2(*myC, myU, myPc, D1c, D2c);
0295 
0296     double Ndu = D1c.Magnitude();
0297     if (Ndu <= myTol) // Singular case
0298     {
0299       // Derivative is approximated by three points
0300       const double DivisionFactor = 0.01;
0301       double       du;
0302       if ((myUsupremum >= RealLast()) || (myUinfium <= RealFirst()))
0303         du = 0.0;
0304       else
0305         du = myUsupremum - myUinfium;
0306 
0307       const double aDelta = std::max(du * DivisionFactor, MinStep);
0308 
0309       double F1, F2, F3;
0310 
0311       if (myU - myUinfium < 2 * aDelta)
0312       {
0313         F1              = theF;
0314         const double U2 = myU + aDelta;
0315         const double U3 = myU + aDelta * 2.0;
0316 
0317         if (!((Value(U2, F2)) && (Value(U3, F3))))
0318         {
0319           myD1Init = false;
0320           return false;
0321         }
0322 
0323         theDF = (-3 * F1 + 4 * F2 - F3) / (2.0 * aDelta);
0324       }
0325       else
0326       {
0327         F3              = theF;
0328         const double U1 = myU - aDelta * 2.0;
0329         const double U2 = myU - aDelta;
0330 
0331         if (!((Value(U2, F2)) && (Value(U1, F1))))
0332         {
0333           myD1Init = false;
0334           return false;
0335         }
0336 
0337         theDF = (F1 - 4 * F2 + 3 * F3) / (2.0 * aDelta);
0338       }
0339       myU  = theU;
0340       myPc = myPc_old;
0341       myP  = myP_old;
0342     }
0343     else
0344     {
0345       TheVector PPc(myP, myPc);
0346       theDF = Ndu + (PPc.Dot(D2c) / Ndu) - theF * (D1c.Dot(D2c)) / (Ndu * Ndu);
0347     }
0348 
0349     myD1f = theDF;
0350 
0351     myD1Init = true;
0352     return true;
0353   }
0354 
0355   //! Save the found extremum.
0356   //! @return State number
0357   int GetStateNumber() override
0358   {
0359     if (!myPinit || !myCinit)
0360     {
0361       throw Standard_TypeMismatch();
0362     }
0363     mySqDist.Append(myPc.SquareDistance(myP));
0364 
0365     // It is necessary to always compute myD1f.
0366     myD1Init = true;
0367     double FF, DD;
0368     Values(myU, FF, DD);
0369 
0370     int IntVal = 0;
0371     if (myD1f > 0.0)
0372     {
0373       IntVal = 1;
0374     }
0375 
0376     myIsMin.Append(IntVal);
0377     myPoint.Append(ThePOnC(myU, myPc));
0378     return 0;
0379   }
0380 
0381   //! Return the number of found extrema.
0382   int NbExt() const { return mySqDist.Length(); }
0383 
0384   //! Returns the Nth square distance.
0385   //! @param theN Index of extremum
0386   double SquareDistance(const int theN) const
0387   {
0388     if (!myPinit || !myCinit)
0389     {
0390       throw Standard_TypeMismatch();
0391     }
0392     return mySqDist.Value(theN);
0393   }
0394 
0395   //! Shows if the Nth distance is a minimum.
0396   //! @param theN Index of extremum
0397   bool IsMin(const int theN) const
0398   {
0399     if (!myPinit || !myCinit)
0400     {
0401       throw Standard_TypeMismatch();
0402     }
0403     return (myIsMin.Value(theN) == 1);
0404   }
0405 
0406   //! Returns the Nth extremum point.
0407   //! @param theN Index of extremum
0408   const ThePOnC& Point(const int theN) const
0409   {
0410     if (!myPinit || !myCinit)
0411     {
0412       throw Standard_TypeMismatch();
0413     }
0414     return myPoint.Value(theN);
0415   }
0416 
0417   //! Determines boundaries of subinterval for root finding.
0418   //! @param theUfirst First parameter bound
0419   //! @param theUlast Last parameter bound
0420   void SubIntervalInitialize(const double theUfirst, const double theUlast)
0421   {
0422     myUinfium   = theUfirst;
0423     myUsupremum = theUlast;
0424   }
0425 
0426   //! Computes a tolerance value. If 1st derivative of curve |D1| < Tol,
0427   //! it is considered D1=0.
0428   double SearchOfTolerance()
0429   {
0430     const int    NPoint = 10;
0431     const double aStep  = (myUsupremum - myUinfium) / static_cast<double>(NPoint);
0432 
0433     int    aNum = 0;
0434     double aMax = -Precision::Infinite();
0435 
0436     do
0437     {
0438       double u = myUinfium + aNum * aStep;
0439       if (u > myUsupremum)
0440         u = myUsupremum;
0441 
0442       ThePoint  Ptemp;
0443       TheVector VDer;
0444       TheCurveTool::D1(*myC, u, Ptemp, VDer);
0445 
0446       if (Precision::IsInfinite(VDer.X()) || Precision::IsInfinite(VDer.Y()))
0447       {
0448         continue;
0449       }
0450 
0451       double vm = VDer.Magnitude();
0452       if (vm > aMax)
0453         aMax = vm;
0454     } while (++aNum < NPoint + 1);
0455 
0456     return std::max(aMax * TolFactor, MinTol);
0457   }
0458 
0459 private:
0460   static constexpr double TolFactor = 1.e-12;
0461   static constexpr double MinTol    = 1.e-20;
0462   static constexpr double MinStep   = 1.e-7;
0463   static constexpr int    MaxOrder  = 3;
0464 
0465   ThePoint                     myP;
0466   TheCurve*                    myC;
0467   double                       myU;
0468   ThePoint                     myPc;
0469   double                       myD1f;
0470   NCollection_Sequence<double> mySqDist;
0471   NCollection_Sequence<int>    myIsMin;
0472   TheSeqPOnC                   myPoint;
0473   bool                         myPinit;
0474   bool                         myCinit;
0475   bool                         myD1Init;
0476   double                       myTol;
0477   int                          myMaxDerivOrder;
0478   double                       myUinfium;
0479   double                       myUsupremum;
0480 };
0481 
0482 #endif // _Extrema_GFuncExtPC_HeaderFile
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!