include/opencascade/Extrema_GGExtPC.hxx

0001 // Created on: 1992-10-19
0002 // Created by: Laurent PAINNOT
0003 // Copyright (c) 1992-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_GGExtPC_HeaderFile
0018 #define _Extrema_GGExtPC_HeaderFile
0019
0020 #include <Standard.hxx>
0021 #include <Standard_DefineAlloc.hxx>
0022
0023 #include <ElCLib.hxx>
0024 #include <GeomAbs_CurveType.hxx>
0025 #include <NCollection_Array1.hxx>
0026 #include <Precision.hxx>
0027 #include <Standard_Failure.hxx>
0028 #include <Standard_OutOfRange.hxx>
0029 #include <StdFail_NotDone.hxx>
0030 #include <NCollection_HArray1.hxx>
0031 #include <NCollection_Sequence.hxx>
0032
0033 //! Generic class for computing extremal distances between a point and a curve.
0034 //!
0035 //! This template class provides comprehensive extremum search functionality,
0036 //! handling different curve types (lines, circles, ellipses, parabolas,
0037 //! hyperbolas, Bezier, BSpline, and general curves) with optimized algorithms.
0038 //!
0039 //! @tparam TheCurve       The curve type (e.g., Adaptor3d_Curve, Adaptor2d_Curve2d)
0040 //! @tparam TheCurveTool   The curve tool providing static methods
0041 //! @tparam TheExtPElC     The elementary curve extremum class
0042 //! @tparam ThePoint       The point type (e.g., gp_Pnt, gp_Pnt2d)
0043 //! @tparam TheVector      The vector type (e.g., gp_Vec, gp_Vec2d)
0044 //! @tparam ThePOnC        The point-on-curve type
0045 //! @tparam TheSequenceOfPOnC The sequence of points on curve
0046 //! @tparam TheEPC         The general extremum point-curve class
0047 template <typename TheCurve,
0048           typename TheCurveTool,
0049           typename TheExtPElC,
0050           typename ThePoint,
0051           typename TheVector,
0052           typename ThePOnC,
0053           typename TheSequenceOfPOnC,
0054           typename TheEPC>
0055 class Extrema_GGExtPC
0056 {
0057 public:
0058   DEFINE_STANDARD_ALLOC
0059
0060   //! Default constructor.
0061   Extrema_GGExtPC()
0062       : myC(nullptr),
0063         mydone(false),
0064         mydist1(RealLast()),
0065         mydist2(RealLast()),
0066         mytolu(0.0),
0067         mytolf(0.0),
0068         mysample(17),
0069         myintuinf(Precision::Infinite()),
0070         myintusup(Precision::Infinite()),
0071         myuinf(Precision::Infinite()),
0072         myusup(Precision::Infinite()),
0073         type(GeomAbs_OtherCurve)
0074   {
0075   }
0076
0077   //! Calculates all extremum distances between point P and curve C
0078   //! within the specified parameter range.
0079   //! @param theP    The point
0080   //! @param theC    The curve
0081   //! @param theUinf Lower bound of parameter range
0082   //! @param theUsup Upper bound of parameter range
0083   //! @param theTolF Tolerance on function value (default 1.0e-10)
0084   Extrema_GGExtPC(const ThePoint& theP,
0085                   const TheCurve& theC,
0086                   const double    theUinf,
0087                   const double    theUsup,
0088                   const double    theTolF = 1.0e-10)
0089   {
0090     Initialize(theC, theUinf, theUsup, theTolF);
0091     Perform(theP);
0092   }
0093
0094   //! Calculates all extremum distances between point P and curve C
0095   //! using the full curve parameter range.
0096   //! @param theP    The point
0097   //! @param theC    The curve
0098   //! @param theTolF Tolerance on function value (default 1.0e-10)
0099   Extrema_GGExtPC(const ThePoint& theP, const TheCurve& theC, const double theTolF = 1.0e-10)
0100   {
0101     Initialize(theC,
0102                TheCurveTool::FirstParameter(theC),
0103                TheCurveTool::LastParameter(theC),
0104                theTolF);
0105     Perform(theP);
0106   }
0107
0108   //! Initializes the algorithm with curve and parameter range.
0109   //! @param theC    The curve
0110   //! @param theUinf Lower bound of parameter range
0111   //! @param theUsup Upper bound of parameter range
0112   //! @param theTolF Tolerance on function value (default 1.0e-10)
0113   void Initialize(const TheCurve& theC,
0114                   const double    theUinf,
0115                   const double    theUsup,
0116                   const double    theTolF = 1.0e-10)
0117   {
0118     myC       = const_cast<TheCurve*>(&theC);
0119     myintuinf = myuinf = theUinf;
0120     myintusup = myusup = theUsup;
0121     mytolf             = theTolF;
0122     mytolu             = TheCurveTool::Resolution(*myC, Precision::Confusion());
0123     type               = TheCurveTool::GetType(theC);
0124     mydone             = false;
0125     mydist1            = RealLast();
0126     mydist2            = RealLast();
0127     mysample           = 17;
0128   }
0129
0130   //! Performs the extremum computation for the given point.
0131   //! @param theP The point to find extrema from
0132   void Perform(const ThePoint& theP)
0133   {
0134     mySqDist.Clear();
0135     mypoint.Clear();
0136     myismin.Clear();
0137     int    i, NbExt, n;
0138     double U;
0139     mysample             = 17;
0140     constexpr double t3d = Precision::Confusion();
0141
0142     if (Precision::IsInfinite(myuinf))
0143       mydist1 = RealLast();
0144     else
0145     {
0146       Pf      = TheCurveTool::Value(*myC, myuinf);
0147       mydist1 = theP.SquareDistance(Pf);
0148     }
0149
0150     if (Precision::IsInfinite(myusup))
0151       mydist2 = RealLast();
0152     else
0153     {
0154       Pl      = TheCurveTool::Value(*myC, myusup);
0155       mydist2 = theP.SquareDistance(Pl);
0156     }
0157
0158     TheCurve& aCurve = *myC;
0159
0160     switch (type)
0161     {
0162       case GeomAbs_Circle: {
0163         myExtPElC.Perform(theP, TheCurveTool::Circle(aCurve), t3d, myuinf, myusup);
0164         break;
0165       }
0166       case GeomAbs_Ellipse: {
0167         myExtPElC.Perform(theP, TheCurveTool::Ellipse(aCurve), t3d, myuinf, myusup);
0168         break;
0169       }
0170       case GeomAbs_Parabola: {
0171         myExtPElC.Perform(theP, TheCurveTool::Parabola(aCurve), t3d, myuinf, myusup);
0172         break;
0173       }
0174       case GeomAbs_Hyperbola: {
0175         myExtPElC.Perform(theP, TheCurveTool::Hyperbola(aCurve), t3d, myuinf, myusup);
0176         break;
0177       }
0178       case GeomAbs_Line: {
0179         myExtPElC.Perform(theP, TheCurveTool::Line(aCurve), t3d, myuinf, myusup);
0180         break;
0181       }
0182       case GeomAbs_BezierCurve: {
0183         myintuinf = myuinf;
0184         myintusup = myusup;
0185         mysample  = (TheCurveTool::Bezier(aCurve))->NbPoles() * 2;
0186         myExtPC.Initialize(aCurve);
0187         IntervalPerform(theP);
0188         return;
0189       }
0190       case GeomAbs_BSplineCurve: {
0191         auto                              aBSpline  = TheCurveTool::BSpline(aCurve);
0192         const int                         aFirstIdx = aBSpline->FirstUKnotIndex();
0193         const int                         aLastIdx  = aBSpline->LastUKnotIndex();
0194         const NCollection_Array1<double>& aKnots    = aBSpline->Knots();
0195
0196         double       aPeriodJump = 0.0;
0197         const double aTolCoeff   = (myusup - myuinf) * Precision::PConfusion();
0198         if (TheCurveTool::IsPeriodic(aCurve))
0199         {
0200           int aPeriodShift = int((myuinf - aKnots(aFirstIdx)) / TheCurveTool::Period(aCurve));
0201           if (myuinf < aKnots(aFirstIdx) - aTolCoeff)
0202             aPeriodShift--;
0203           aPeriodJump = TheCurveTool::Period(aCurve) * aPeriodShift;
0204         }
0205
0206         int anIdx;
0207         int aFirstUsedKnot = aFirstIdx, aLastUsedKnot = aLastIdx;
0208         for (anIdx = aFirstIdx; anIdx <= aLastIdx; anIdx++)
0209         {
0210           double aKnot = aKnots(anIdx) + aPeriodJump;
0211           if (myuinf >= aKnot - aTolCoeff)
0212             aFirstUsedKnot = anIdx;
0213           else
0214             break;
0215         }
0216         for (anIdx = aLastIdx; anIdx >= aFirstIdx; anIdx--)
0217         {
0218           double aKnot = aKnots(anIdx) + aPeriodJump;
0219           if (myusup <= aKnot + aTolCoeff)
0220             aLastUsedKnot = anIdx;
0221           else
0222             break;
0223         }
0224
0225         if (aFirstUsedKnot == aLastUsedKnot)
0226         {
0227           aFirstUsedKnot = aFirstIdx;
0228           aLastUsedKnot  = aFirstIdx + 1;
0229         }
0230
0231         mysample = aBSpline->Degree() + 1;
0232
0233         if (mysample == 2)
0234         {
0235           ThePoint aPmin;
0236           double   tmin = 0., distmin = RealLast();
0237           double   aMin1 = 0., aMin2 = 0.;
0238           myExtPC.Initialize(aCurve);
0239           for (anIdx = aFirstUsedKnot; anIdx < aLastUsedKnot; anIdx++)
0240           {
0241             double aF = aKnots(anIdx) + aPeriodJump, aL = aKnots(anIdx + 1) + aPeriodJump;
0242
0243             if (anIdx == aFirstUsedKnot)
0244               aF = myuinf;
0245             else if (anIdx == aLastUsedKnot - 1)
0246               aL = myusup;
0247
0248             ThePoint aP1, aP2;
0249             TheCurveTool::D0(aCurve, aF, aP1);
0250             TheCurveTool::D0(aCurve, aL, aP2);
0251             TheVector aBase1(theP, aP1), aBase2(theP, aP2);
0252             TheVector aV(aP2, aP1);
0253             double    aVal1 = aV.Dot(aBase1);
0254             double    aVal2 = aV.Dot(aBase2);
0255             if (anIdx == aFirstUsedKnot)
0256               aMin1 = theP.SquareDistance(aP1);
0257             else
0258             {
0259               aMin1 = aMin2;
0260               if (distmin > aMin1)
0261               {
0262                 distmin = aMin1;
0263                 tmin    = aF;
0264                 aPmin   = aP1;
0265               }
0266             }
0267             aMin2             = theP.SquareDistance(aP2);
0268             double aMinSqDist = std::min(aMin1, aMin2);
0269             double aMinDer    = std::min(std::abs(aVal1), std::abs(aVal2));
0270             if (!(Precision::IsInfinite(aVal1) || Precision::IsInfinite(aVal2)))
0271             {
0272               if (aVal1 * aVal2 <= 0.0 || aMinSqDist < 100. * Precision::SquareConfusion()
0273                   || 2. * aMinDer < Precision::Confusion())
0274               {
0275                 myintuinf = aF;
0276                 myintusup = aL;
0277                 IntervalPerform(theP);
0278               }
0279             }
0280           }
0281           if (!Precision::IsInfinite(distmin))
0282           {
0283             bool isToAdd = true;
0284             NbExt        = mypoint.Length();
0285             for (i = 1; i <= NbExt && isToAdd; i++)
0286             {
0287               double t = mypoint.Value(i).Parameter();
0288               isToAdd  = (distmin < mySqDist(i)) && (std::abs(t - tmin) > mytolu);
0289             }
0290             if (isToAdd)
0291             {
0292               ThePOnC PC(tmin, aPmin);
0293               mySqDist.Append(distmin);
0294               myismin.Append(true);
0295               mypoint.Append(PC);
0296             }
0297           }
0298         }
0299         else
0300         {
0301           int                        aValIdx = 1;
0302           NCollection_Array1<double> aVal(1, (mysample) * (aLastUsedKnot - aFirstUsedKnot) + 1);
0303           NCollection_Array1<double> aParam(1, (mysample) * (aLastUsedKnot - aFirstUsedKnot) + 1);
0304           for (anIdx = aFirstUsedKnot; anIdx < aLastUsedKnot; anIdx++)
0305           {
0306             double aF = aKnots(anIdx) + aPeriodJump, aL = aKnots(anIdx + 1) + aPeriodJump;
0307
0308             if (anIdx == aFirstUsedKnot)
0309               aF = myuinf;
0310             if (anIdx == aLastUsedKnot - 1)
0311               aL = myusup;
0312
0313             double aStep = (aL - aF) / mysample;
0314             for (int aPntIdx = 0; aPntIdx < mysample; aPntIdx++)
0315             {
0316               double aCurrentParam = aF + aStep * aPntIdx;
0317               aVal(aValIdx)   = TheCurveTool::Value(aCurve, aCurrentParam).SquareDistance(theP);
0318               aParam(aValIdx) = aCurrentParam;
0319               aValIdx++;
0320             }
0321           }
0322           aVal(aValIdx)   = TheCurveTool::Value(aCurve, myusup).SquareDistance(theP);
0323           aParam(aValIdx) = myusup;
0324
0325           myExtPC.Initialize(aCurve);
0326
0327           for (anIdx = aVal.Lower() + 1; anIdx < aVal.Upper(); anIdx++)
0328           {
0329             if (aVal(anIdx) <= Precision::SquareConfusion())
0330             {
0331               mySqDist.Append(aVal(anIdx));
0332               myismin.Append(true);
0333               mypoint.Append(ThePOnC(aParam(anIdx), TheCurveTool::Value(aCurve, aParam(anIdx))));
0334             }
0335             if ((aVal(anIdx) >= aVal(anIdx + 1) && aVal(anIdx) >= aVal(anIdx - 1))
0336                 || (aVal(anIdx) <= aVal(anIdx + 1) && aVal(anIdx) <= aVal(anIdx - 1)))
0337             {
0338               myintuinf = aParam(anIdx - 1);
0339               myintusup = aParam(anIdx + 1);
0340               IntervalPerform(theP);
0341             }
0342           }
0343
0344           if (mydist1 > Precision::SquareConfusion() && !Precision::IsPositiveInfinite(mydist1))
0345           {
0346             ThePoint  aP1, aP2;
0347             TheVector aV1, aV2;
0348             TheCurveTool::D1(aCurve, aParam.Value(aParam.Lower()), aP1, aV1);
0349             TheCurveTool::D1(aCurve, aParam.Value(aParam.Lower() + 1), aP2, aV2);
0350             TheVector aBase1(theP, aP1), aBase2(theP, aP2);
0351             double    aVal1 = aV1.Dot(aBase1);
0352             double    aVal2 = aV2.Dot(aBase2);
0353             if (!(Precision::IsInfinite(aVal1) || Precision::IsInfinite(aVal2)))
0354             {
0355               if (aVal1 * aVal2 <= 0.0 || aBase1.Dot(aBase2) <= 0.0
0356                   || 2.0 * std::abs(aVal1) < Precision::Confusion())
0357               {
0358                 myintuinf = aParam(aVal.Lower());
0359                 myintusup = aParam(aVal.Lower() + 1);
0360                 IntervalPerform(theP);
0361               }
0362             }
0363           }
0364
0365           if (mydist2 > Precision::SquareConfusion() && !Precision::IsPositiveInfinite(mydist2))
0366           {
0367             ThePoint  aP1, aP2;
0368             TheVector aV1, aV2;
0369             TheCurveTool::D1(aCurve, aParam.Value(aParam.Upper() - 1), aP1, aV1);
0370             TheCurveTool::D1(aCurve, aParam.Value(aParam.Upper()), aP2, aV2);
0371             TheVector aBase1(theP, aP1), aBase2(theP, aP2);
0372             double    aVal1 = aV1.Dot(aBase1);
0373             double    aVal2 = aV2.Dot(aBase2);
0374
0375             if (!(Precision::IsInfinite(aVal1) || Precision::IsInfinite(aVal2)))
0376             {
0377               if (aVal1 * aVal2 <= 0.0 || aBase1.Dot(aBase2) <= 0.0
0378                   || 2.0 * std::abs(aVal2) < Precision::Confusion())
0379               {
0380                 myintuinf = aParam(aVal.Upper() - 1);
0381                 myintusup = aParam(aVal.Upper());
0382                 IntervalPerform(theP);
0383               }
0384             }
0385           }
0386         }
0387         mydone = true;
0388         break;
0389       }
0390       default: {
0391         constexpr int                            aMaxSample   = 17;
0392         bool                                     IntExtIsDone = false;
0393         bool                                     IntIsNotValid;
0394         occ::handle<NCollection_HArray1<double>> theHInter;
0395         n = TheCurveTool::NbIntervals(aCurve, GeomAbs_C2);
0396         if (n > 1)
0397         {
0398           theHInter = new NCollection_HArray1<double>(1, n + 1);
0399           TheCurveTool::Intervals(aCurve, theHInter->ChangeArray1(), GeomAbs_C2);
0400         }
0401         else
0402         {
0403           theHInter = TheCurveTool::DeflCurvIntervals(aCurve);
0404           n         = theHInter->Length() - 1;
0405         }
0406         mysample      = std::max(mysample / n, aMaxSample);
0407         double maxint = 0.;
0408         for (i = 1; i <= n; ++i)
0409         {
0410           double dt = theHInter->Value(i + 1) - theHInter->Value(i);
0411           if (maxint < dt)
0412             maxint = dt;
0413         }
0414         bool      isPeriodic = TheCurveTool::IsPeriodic(aCurve);
0415         TheVector V1;
0416         ThePoint  PP;
0417         double    s1 = 0.0;
0418         double    s2 = 0.0;
0419         myExtPC.Initialize(aCurve);
0420         for (i = 1; i <= n; i++)
0421         {
0422           myintuinf = theHInter->Value(i);
0423           myintusup = theHInter->Value(i + 1);
0424           mysample  = std::max(RealToInt(aMaxSample * (myintusup - myintuinf) / maxint), 3);
0425
0426           double anInfToCheck = myintuinf;
0427           double aSupToCheck  = myintusup;
0428
0429           if (isPeriodic)
0430           {
0431             double aPeriod = TheCurveTool::Period(aCurve);
0432             anInfToCheck   = ElCLib::InPeriod(myintuinf, myuinf, myuinf + aPeriod);
0433             aSupToCheck    = myintusup + (anInfToCheck - myintuinf);
0434           }
0435           IntIsNotValid = (myuinf > aSupToCheck) || (myusup < anInfToCheck);
0436
0437           if (IntIsNotValid)
0438             continue;
0439
0440           if (myuinf >= anInfToCheck)
0441             anInfToCheck = myuinf;
0442           if (myusup <= aSupToCheck)
0443             aSupToCheck = myusup;
0444           if ((aSupToCheck - anInfToCheck) <= mytolu)
0445             continue;
0446
0447           if (i != 1)
0448           {
0449             TheCurveTool::D1(aCurve, myintuinf, PP, V1);
0450             s1 = (TheVector(theP, PP)) * V1;
0451             if (s1 * s2 < 0.0)
0452             {
0453               mySqDist.Append(PP.SquareDistance(theP));
0454               myismin.Append((s1 < 0.0));
0455               mypoint.Append(ThePOnC(myintuinf, PP));
0456             }
0457           }
0458           if (i != n)
0459           {
0460             TheCurveTool::D1(aCurve, myintusup, PP, V1);
0461             s2 = (TheVector(theP, PP)) * V1;
0462           }
0463
0464           IntervalPerform(theP);
0465           IntExtIsDone = IntExtIsDone || mydone;
0466         }
0467
0468         mydone = IntExtIsDone;
0469         break;
0470       }
0471     }
0472
0473     // Postprocessing
0474     if (type == GeomAbs_BSplineCurve || type == GeomAbs_OffsetCurve || type == GeomAbs_OtherCurve)
0475     {
0476       if (mydist1 < Precision::SquareConfusion() || mydist2 < Precision::SquareConfusion())
0477       {
0478         bool isFirstAdded = false;
0479         bool isLastAdded  = false;
0480         int  aNbPoints    = mypoint.Length();
0481         for (i = 1; i <= aNbPoints; i++)
0482         {
0483           U = mypoint.Value(i).Parameter();
0484           if (std::abs(U - myuinf) < mytolu)
0485             isFirstAdded = true;
0486           else if (std::abs(myusup - U) < mytolu)
0487             isLastAdded = true;
0488         }
0489         if (!isFirstAdded && mydist1 < Precision::SquareConfusion())
0490         {
0491           mySqDist.Prepend(mydist1);
0492           myismin.Prepend(true);
0493           mypoint.Prepend(ThePOnC(myuinf, Pf));
0494         }
0495         if (!isLastAdded && mydist2 < Precision::SquareConfusion())
0496         {
0497           mySqDist.Append(mydist2);
0498           myismin.Append(true);
0499           mypoint.Append(ThePOnC(myusup, Pl));
0500         }
0501         mydone = true;
0502       }
0503     }
0504     else
0505     {
0506       mydone = myExtPElC.IsDone();
0507       if (mydone)
0508       {
0509         NbExt = myExtPElC.NbExt();
0510         for (i = 1; i <= NbExt; i++)
0511         {
0512           ThePOnC PC = myExtPElC.Point(i);
0513           U          = PC.Parameter();
0514           if (TheCurveTool::IsPeriodic(aCurve))
0515           {
0516             U = ElCLib::InPeriod(U, myuinf, myuinf + TheCurveTool::Period(aCurve));
0517           }
0518           if ((U >= myuinf - mytolu) && (U <= myusup + mytolu))
0519           {
0520             PC.SetValues(U, myExtPElC.Point(i).Value());
0521             mySqDist.Append(myExtPElC.SquareDistance(i));
0522             myismin.Append(myExtPElC.IsMin(i));
0523             mypoint.Append(PC);
0524           }
0525         }
0526       }
0527     }
0528   }
0529
0530   //! Returns true if the distances are found.
0531   bool IsDone() const { return mydone; }
0532
0533   //! Returns the Nth extremum square distance.
0534   //! @param theN Index of the extremum (1-based)
0535   //! @return Square distance value
0536   double SquareDistance(const int theN) const
0537   {
0538     if ((theN < 1) || (theN > NbExt()))
0539       throw Standard_OutOfRange();
0540     return mySqDist.Value(theN);
0541   }
0542
0543   //! Returns the number of extremum distances.
0544   //! @return Number of extrema found
0545   int NbExt() const
0546   {
0547     if (!IsDone())
0548       throw StdFail_NotDone();
0549     return mySqDist.Length();
0550   }
0551
0552   //! Returns true if the Nth extremum distance is a minimum.
0553   //! @param theN Index of the extremum (1-based)
0554   //! @return true if minimum, false if maximum
0555   bool IsMin(const int theN) const
0556   {
0557     if ((theN < 1) || (theN > NbExt()))
0558       throw Standard_OutOfRange();
0559     return myismin.Value(theN);
0560   }
0561
0562   //! Returns the point of the Nth extremum distance.
0563   //! @param theN Index of the extremum (1-based)
0564   //! @return The point on curve
0565   const ThePOnC& Point(const int theN) const
0566   {
0567     if ((theN < 1) || (theN > NbExt()))
0568       throw Standard_OutOfRange();
0569     return mypoint.Value(theN);
0570   }
0571
0572   //! Returns the distances at curve endpoints.
0573   //! @param[out] theDist1 Square distance to first point
0574   //! @param[out] theDist2 Square distance to last point
0575   //! @param[out] theP1 First point on curve
0576   //! @param[out] theP2 Last point on curve
0577   void TrimmedSquareDistances(double&   theDist1,
0578                               double&   theDist2,
0579                               ThePoint& theP1,
0580                               ThePoint& theP2) const
0581   {
0582     theDist1 = mydist1;
0583     theDist2 = mydist2;
0584     theP1    = Pf;
0585     theP2    = Pl;
0586   }
0587
0588 protected:
0589   //! Performs extremum search on the current interval.
0590   void IntervalPerform(const ThePoint& theP)
0591   {
0592     int    i;
0593     double U;
0594     myExtPC.Initialize(mysample, myintuinf, myintusup, mytolu, mytolf);
0595     myExtPC.Perform(theP);
0596     mydone = myExtPC.IsDone();
0597     if (mydone)
0598     {
0599       int NbExt = myExtPC.NbExt();
0600       for (i = 1; i <= NbExt; i++)
0601       {
0602         ThePOnC PC = myExtPC.Point(i);
0603         U          = PC.Parameter();
0604         if (TheCurveTool::IsPeriodic(*myC))
0605         {
0606           U = ElCLib::InPeriod(U, myuinf, myuinf + TheCurveTool::Period(*myC));
0607         }
0608         if ((U >= myuinf - mytolu) && (U <= myusup + mytolu))
0609         {
0610           AddSol(U, PC.Value(), myExtPC.SquareDistance(i), myExtPC.IsMin(i));
0611         }
0612       }
0613     }
0614   }
0615
0616   //! Adds a solution if not already present.
0617   void AddSol(const double theU, const ThePoint& theP, const double theSqDist, const bool isMin)
0618   {
0619     int i, NbExt = mypoint.Length();
0620     for (i = 1; i <= NbExt; i++)
0621     {
0622       double t = mypoint.Value(i).Parameter();
0623       if (std::abs(t - theU) <= mytolu)
0624       {
0625         return;
0626       }
0627     }
0628     ThePOnC PC(theU, theP);
0629     mySqDist.Append(theSqDist);
0630     myismin.Append(isMin);
0631     mypoint.Append(PC);
0632   }
0633
0634 private:
0635   TheCurve*                    myC;
0636   ThePoint                     Pf;
0637   ThePoint                     Pl;
0638   TheExtPElC                   myExtPElC;
0639   TheSequenceOfPOnC            mypoint;
0640   bool                         mydone;
0641   double                       mydist1;
0642   double                       mydist2;
0643   TheEPC                       myExtPC;
0644   double                       mytolu;
0645   double                       mytolf;
0646   int                          mysample;
0647   double                       myintuinf;
0648   double                       myintusup;
0649   double                       myuinf;
0650   double                       myusup;
0651   GeomAbs_CurveType            type;
0652   NCollection_Sequence<bool>   myismin;
0653   NCollection_Sequence<double> mySqDist;
0654 };
0655
0656 #endif // _Extrema_GGExtPC_HeaderFile