EIC code displayed by LXR master/​include/​opencascade/​ExtremaPC_GridEvaluator.hxx	Source navigation Diff markup Identifier search General search
	Version: master test Revert   Change

File indexing completed on 2026-09-28 09:19:57

 // Copyright (c) 2025 OPEN CASCADE SAS
 //
 // This file is part of Open CASCADE Technology software library.
 //
 // This library is free software; you can redistribute it and/or modify it under
 // the terms of the GNU Lesser General Public License version 2.1 as published
 // by the Free Software Foundation, with special exception defined in the file
 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
 // distribution for complete text of the license and disclaimer of any warranty.
 //
 // Alternatively, this file may be used under the terms of Open CASCADE
 // commercial license or contractual agreement.
 
 #ifndef _ExtremaPC_GridEvaluator_HeaderFile
 #define _ExtremaPC_GridEvaluator_HeaderFile
 
 #include <Adaptor3d_Curve.hxx>
 #include <ExtremaPC.hxx>
 #include <ExtremaPC_DistanceFunction.hxx>
 #include <GeomGridEval.hxx>
 #include <math_Vector.hxx>
 #include <MathRoot_Newton.hxx>
 #include <MathUtils_Config.hxx>
 #include <NCollection_Array1.hxx>
 #include <NCollection_DynamicArray.hxx>
 
 #include <algorithm>
 #include <cmath>
 #include <limits>
 #include <utility>
 
 //! @brief Grid-based point-curve extrema computation class.
 //!
 //! Provides grid-based extrema finding algorithm with cached state for
 //! optimal performance on repeated queries. Used by BSpline, Bezier,
 //! Offset, and other curve evaluators.
 //!
 //! Algorithm:
 //! 1. Build grid of (parameter, point, D1) from GeomGridEval
 //! 2. Linear scan of grid to find candidate intervals (sign changes in F(u))
 //! 3. Newton refinement on each candidate
 //! 4. Optional endpoint handling
 //!
 //! All temporary vectors are stored as mutable fields and reused via Clear()
 //! to avoid repeated heap allocations.
 class ExtremaPC_GridEvaluator
 {
 public:
   //! Cached grid point with pre-computed data.
   struct GridPoint
   {
     double Param; //!< Parameter value
     gp_Pnt Point; //!< Curve point C(u)
     gp_Vec D1;    //!< First derivative C'(u)
   };
 
   //! Type of candidate extremum detected during grid scan.
   enum class CandidateType
   {
     SignChange, //!< F(u) changes sign between grid points
     NearZero    //!< F(u) is very small at grid point
   };
 
   //! Candidate interval for Newton refinement.
   struct Candidate
   {
     CandidateType Type;   //!< Type of candidate
     int           IdxLo;  //!< Lower grid index
     int           IdxHi;  //!< Upper grid index (same as IdxLo for NearZero)
     double        StartU; //!< Starting point for Newton
   };
 
   //! Default constructor.
   ExtremaPC_GridEvaluator() = default;
 
   //! @brief Build grid from GeomGridEval D1 results.
   //!
   //! @tparam GridEval type with EvaluateGridD1(params) method
   //! @param theEval grid evaluator
   //! @param theParams parameter values (math_Vector with Array1() accessor)
   template <typename GridEval>
   void BuildGrid(GridEval& theEval, const math_Vector& theParams)
   {
     // Use Array1() accessor to pass to GeomGridEval which expects NCollection_Array1
     NCollection_Array1<GeomGridEval::CurveD1> aD1Grid = theEval.EvaluateGridD1(theParams.Array1());
 
     const int aNbParams = theParams.Length();
 
     // Resize grid if needed
     if (myGrid.Length() != aNbParams)
     {
       myGrid = NCollection_Array1<GridPoint>(0, aNbParams - 1);
     }
 
     for (int i = 0; i < aNbParams; ++i)
     {
       const int aD1Idx = aD1Grid.Lower() + i;
 
       myGrid[i].Param = theParams(theParams.Lower() + i);
       myGrid[i].Point = aD1Grid.Value(aD1Idx).Point;
       myGrid[i].D1    = aD1Grid.Value(aD1Idx).D1;
     }
   }
 
   //! Returns the cached grid.
   const NCollection_Array1<GridPoint>& Grid() const { return myGrid; }
 
   //! Returns mutable reference to the result for post-processing.
   ExtremaPC::Result& Result() const { return myResult; }
 
   //! @brief Perform extrema computation using cached grid (interior only).
   //!
   //! @param theCurve curve adaptor
   //! @param theP query point
   //! @param theDomain parameter domain
   //! @param theTol tolerance
   //! @param theMode search mode
   //! @return const reference to result with interior extrema only
   [[nodiscard]] const ExtremaPC::Result& Perform(const Adaptor3d_Curve&     theCurve,
                                                  const gp_Pnt&              theP,
                                                  const ExtremaPC::Domain1D& theDomain,
                                                  double                     theTol,
                                                  ExtremaPC::SearchMode      theMode) const
   {
     myResult.Clear();
     scanGrid(theP, theTol, theMode);
     refineCandidates(theCurve, theP, theDomain, theTol, theMode);
 
     if (!myResult.Extrema.IsEmpty())
     {
       myResult.Status = ExtremaPC::Status::OK;
     }
     return myResult;
   }
 
   //! @brief Build uniform parameter grid.
   //! @return math_Vector with 1-based indexing
   static math_Vector BuildUniformParams(double theUMin, double theUMax, int theNbSamples)
   {
     math_Vector  aParams(1, theNbSamples);
     const double aStep = (theUMax - theUMin) / (theNbSamples - 1);
 
     for (int i = 1; i <= theNbSamples; ++i)
     {
       aParams(i) = theUMin + (i - 1) * aStep;
     }
     aParams(theNbSamples) = theUMax; // Ensure exact endpoint
 
     return aParams;
   }
 
 private:
   //! @brief Scan grid to find candidate intervals for extrema.
   void scanGrid(const gp_Pnt& theP, double theTol, ExtremaPC::SearchMode theMode) const
   {
     myCandidates.Clear();
     const int aNbGrid = myGrid.Length();
 
     if (aNbGrid < 2)
     {
       return;
     }
 
     // Resize processed array if needed
     if (myProcessed.Length() != aNbGrid)
     {
       myProcessed = NCollection_Array1<bool>(0, aNbGrid - 1);
     }
     myProcessed.Init(false);
 
     double aPrevF     = 0.0;
     double aPrevDist  = 0.0;
     bool   aPrevValid = false;
 
     for (int i = 0; i < aNbGrid; ++i)
     {
       const GridPoint& aGP = myGrid[i];
 
       // Compute distance function value: F(u) = (C(u) - P) . C'(u)
       gp_Vec aVec(theP, aGP.Point);
       double aF    = aVec.Dot(aGP.D1);
       double aDist = aVec.SquareMagnitude();
 
       // Check for sign change with previous point
       if (aPrevValid && aPrevF * aF < 0.0 && !myProcessed[i - 1])
       {
         Candidate aCand;
         aCand.Type  = CandidateType::SignChange;
         aCand.IdxLo = i - 1;
         aCand.IdxHi = i;
         // Use linear interpolation for better starting point (secant method)
         double aFLo  = aPrevF;
         double aFHi  = aF;
         double aULo  = myGrid[i - 1].Param;
         double aUHi  = aGP.Param;
         aCand.StartU = aULo - aFLo * (aUHi - aULo) / (aFHi - aFLo);
         myCandidates.Append(aCand);
         myProcessed[i - 1] = true;
         myProcessed[i]     = true;
       }
 
       // Check for near-zero F (direct hit on extremum)
       if (std::abs(aF) < theTol * ExtremaPC::THE_NEAR_ZERO_F_FACTOR && !myProcessed[i])
       {
         Candidate aCand;
         aCand.Type   = CandidateType::NearZero;
         aCand.IdxLo  = i;
         aCand.IdxHi  = i;
         aCand.StartU = aGP.Param;
         myCandidates.Append(aCand);
         myProcessed[i] = true;
       }
 
       // Check for local extremum by distance comparison (3-point test)
       if (i > 0 && i < aNbGrid - 1 && !myProcessed[i])
       {
         double aNextDist = theP.SquareDistance(myGrid[i + 1].Point);
 
         // Local minimum: distance decreases then increases
         bool aIsLocalMin = (aDist <= aPrevDist && aDist <= aNextDist);
         // Local maximum: distance increases then decreases
         bool aIsLocalMax = (aDist >= aPrevDist && aDist >= aNextDist);
 
         if ((theMode == ExtremaPC::SearchMode::Min || theMode == ExtremaPC::SearchMode::MinMax)
             && aIsLocalMin)
         {
           Candidate aCand;
           aCand.Type   = CandidateType::NearZero;
           aCand.IdxLo  = i;
           aCand.IdxHi  = i;
           aCand.StartU = aGP.Param;
           myCandidates.Append(aCand);
           myProcessed[i] = true;
         }
         else if ((theMode == ExtremaPC::SearchMode::Max || theMode == ExtremaPC::SearchMode::MinMax)
                  && aIsLocalMax && !aIsLocalMin)
         {
           Candidate aCand;
           aCand.Type   = CandidateType::NearZero;
           aCand.IdxLo  = i;
           aCand.IdxHi  = i;
           aCand.StartU = aGP.Param;
           myCandidates.Append(aCand);
           myProcessed[i] = true;
         }
       }
 
       aPrevF     = aF;
       aPrevDist  = aDist;
       aPrevValid = true;
     }
   }
 
   //! @brief Refine candidates using Newton's method.
   void refineCandidates(const Adaptor3d_Curve&     theCurve,
                         const gp_Pnt&              theP,
                         const ExtremaPC::Domain1D& theDomain,
                         double                     theTol,
                         ExtremaPC::SearchMode      theMode) const
   {
     myResult.Status = ExtremaPC::Status::OK;
     myFoundRoots.Clear();
     mySortedIndices.Clear();
 
     ExtremaPC_DistanceFunction aFunc(theCurve, theP);
 
     // Newton configuration
     MathUtils::Config aConfig;
     aConfig.XTolerance    = theTol * ExtremaPC::THE_NEWTON_XTOL_FACTOR;
     aConfig.FTolerance    = theTol * ExtremaPC::THE_NEWTON_FTOL_FACTOR;
     aConfig.MaxIterations = ExtremaPC::THE_MAX_NEWTON_ITERATIONS;
 
     // Build sorted indices by estimated distance
     for (int c = 0; c < myCandidates.Length(); ++c)
     {
       const Candidate& aCand     = myCandidates.Value(c);
       double           anEstDist = theP.SquareDistance(myGrid[aCand.IdxLo].Point);
       mySortedIndices.Append(std::make_pair(c, anEstDist));
     }
 
     // Sort by estimated distance for Min mode (ascending), Max mode (descending)
     if (theMode == ExtremaPC::SearchMode::Min)
     {
       std::sort(mySortedIndices.begin(),
                 mySortedIndices.end(),
                 [](const std::pair<int, double>& a, const std::pair<int, double>& b) {
                   return a.second < b.second;
                 });
     }
     else if (theMode == ExtremaPC::SearchMode::Max)
     {
       std::sort(mySortedIndices.begin(),
                 mySortedIndices.end(),
                 [](const std::pair<int, double>& a, const std::pair<int, double>& b) {
                   return a.second > b.second;
                 });
     }
 
     // Best distance found so far (for early termination)
     double aBestSqDist = (theMode == ExtremaPC::SearchMode::Min)
                            ? std::numeric_limits<double>::max()
                            : -std::numeric_limits<double>::max();
 
     for (int s = 0; s < mySortedIndices.Length(); ++s)
     {
       int              c         = mySortedIndices.Value(s).first;
       double           anEstDist = mySortedIndices.Value(s).second;
       const Candidate& aCand     = myCandidates.Value(c);
 
       // Early termination: skip candidates that are clearly worse than the best found.
       // For Min mode: skip if estimated distance > best * (2.0 - threshold), i.e., ~1.1x best.
       // For Max mode: skip if estimated distance < best * threshold, i.e., ~0.9x best.
       constexpr double aMinSkipThreshold = 2.0 - ExtremaPC::THE_MAX_SKIP_THRESHOLD;
       if (theMode == ExtremaPC::SearchMode::Min && anEstDist > aBestSqDist * aMinSkipThreshold)
       {
         break;
       }
       if (theMode == ExtremaPC::SearchMode::Max
           && anEstDist < aBestSqDist * ExtremaPC::THE_MAX_SKIP_THRESHOLD)
       {
         break;
       }
 
       // Skip if too close to already found root
       bool aSkip = false;
       for (int r = 0; r < myFoundRoots.Length(); ++r)
       {
         if (std::abs(aCand.StartU - myFoundRoots.Value(r)) < theTol)
         {
           aSkip = true;
           break;
         }
       }
       if (aSkip)
       {
         continue;
       }
 
       // Determine Newton bounds
       double aULo, aUHi;
       if (aCand.Type == CandidateType::SignChange)
       {
         aULo = myGrid[aCand.IdxLo].Param;
         aUHi = myGrid[aCand.IdxHi].Param;
       }
       else
       {
         double aExpand = (theDomain.Max - theDomain.Min) * ExtremaPC::THE_INTERVAL_EXPAND_RATIO;
         aULo           = std::max(theDomain.Min, aCand.StartU - aExpand);
         aUHi           = std::min(theDomain.Max, aCand.StartU + aExpand);
       }
 
       // Try Newton refinement
       MathUtils::ScalarResult aNewtonRes =
         MathRoot::NewtonBounded(aFunc, aCand.StartU, aULo, aUHi, aConfig);
 
       double aRootU     = 0.0;
       bool   aConverged = false;
 
       if (aNewtonRes.IsDone())
       {
         aRootU     = std::max(theDomain.Min, std::min(theDomain.Max, *aNewtonRes.Root));
         aConverged = true;
       }
       else
       {
         // Try iterative grid refinement as fallback
         double aBestU    = aCand.StartU;
         double aBestDist = std::numeric_limits<double>::max();
         double aRefUMin  = aULo;
         double aRefUMax  = aUHi;
 
         for (int aPass = 0; aPass < ExtremaPC::THE_REFINEMENT_NB_PASSES; ++aPass)
         {
           const int    aNbSamples = ExtremaPC::THE_REFINEMENT_NB_SAMPLES;
           const double aStep      = (aRefUMax - aRefUMin) / (aNbSamples - 1);
 
           for (int i = 0; i < aNbSamples; ++i)
           {
             double aU    = aRefUMin + i * aStep;
             gp_Pnt aPt   = theCurve.Value(aU);
             double aDist = theP.SquareDistance(aPt);
 
             if (aDist < aBestDist)
             {
               aBestDist = aDist;
               aBestU    = aU;
             }
           }
 
           // Narrow range
           double aRangeHalf = (aRefUMax - aRefUMin) * ExtremaPC::THE_RANGE_NARROWING_FACTOR * 0.5;
           aRefUMin          = std::max(theDomain.Min, aBestU - aRangeHalf);
           aRefUMax          = std::min(theDomain.Max, aBestU + aRangeHalf);
 
           // Try Newton with refined point
           MathUtils::ScalarResult aRetryRes =
             MathRoot::NewtonBounded(aFunc, aBestU, aRefUMin, aRefUMax, aConfig);
           if (aRetryRes.IsDone())
           {
             aRootU     = std::max(theDomain.Min, std::min(theDomain.Max, *aRetryRes.Root));
             aConverged = true;
             break;
           }
         }
 
         // Use best grid point as fallback
         if (!aConverged)
         {
           gp_Pnt aPt;
           gp_Vec aD1;
           theCurve.D1(aBestU, aPt, aD1);
           gp_Vec aVec(theP, aPt);
           double aF = aVec.Dot(aD1);
 
           if (std::abs(aF) < theTol * ExtremaPC::THE_FALLBACK_F_FACTOR)
           {
             aRootU     = aBestU;
             aConverged = true;
           }
         }
       }
 
       if (!aConverged)
         continue;
 
       // Check for duplicate
       bool aDuplicate = false;
       for (int r = 0; r < myFoundRoots.Length(); ++r)
       {
         if (std::abs(aRootU - myFoundRoots.Value(r)) < theTol)
         {
           aDuplicate = true;
           break;
         }
       }
       if (aDuplicate)
         continue;
 
       gp_Pnt aPt     = theCurve.Value(aRootU);
       double aSqDist = theP.SquareDistance(aPt);
 
       // Classify as min/max using neighbor sampling
       double aStep = (theDomain.Max - theDomain.Min) * ExtremaPC::THE_REFINEMENT_STEP_RATIO;
       double aDistPlus =
         theP.SquareDistance(theCurve.Value(std::min(theDomain.Max, aRootU + aStep)));
       double aDistMinus =
         theP.SquareDistance(theCurve.Value(std::max(theDomain.Min, aRootU - aStep)));
       bool aIsMin = (aSqDist <= aDistPlus) && (aSqDist <= aDistMinus);
 
       // Filter by mode
       bool aKeep = false;
       if (theMode == ExtremaPC::SearchMode::MinMax)
       {
         aKeep = true;
       }
       else if (theMode == ExtremaPC::SearchMode::Min && aIsMin)
       {
         aKeep = true;
       }
       else if (theMode == ExtremaPC::SearchMode::Max && !aIsMin)
       {
         aKeep = true;
       }
 
       if (aKeep)
       {
         ExtremaPC::ExtremumResult anExt;
         anExt.Parameter      = aRootU;
         anExt.Point          = aPt;
         anExt.SquareDistance = aSqDist;
         anExt.IsMinimum      = aIsMin;
         myResult.Extrema.Append(anExt);
 
         myFoundRoots.Append(aRootU);
 
         // Update best distance for early termination
         if (theMode == ExtremaPC::SearchMode::Min && aSqDist < aBestSqDist)
         {
           aBestSqDist = aSqDist;
         }
         else if (theMode == ExtremaPC::SearchMode::Max && aSqDist > aBestSqDist)
         {
           aBestSqDist = aSqDist;
         }
       }
     }
 
     if (myResult.Extrema.IsEmpty() && myCandidates.IsEmpty())
     {
       myResult.Status = ExtremaPC::Status::NoSolution;
     }
   }
 
 private:
   NCollection_Array1<GridPoint> myGrid; //!< Cached grid
 
   // Mutable cached temporaries (reused via Clear())
   mutable ExtremaPC::Result                   myResult;     //!< Reusable result
   mutable NCollection_DynamicArray<Candidate> myCandidates; //!< Candidates from grid scan
   mutable NCollection_DynamicArray<double>    myFoundRoots; //!< Found roots for dedup
   mutable NCollection_DynamicArray<std::pair<int, double>>
                                    mySortedIndices; //!< Sorted candidate indices
   mutable NCollection_Array1<bool> myProcessed;     //!< Processed flags for grid scan
 };
 
 #endif // _ExtremaPC_GridEvaluator_HeaderFile

This page was automatically generated by the 2.3.7 LXR engine. The LXR team