include/opencascade/ExtremaPC_Ellipse.hxx

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 _ExtremaPC_Ellipse_HeaderFile
0015 #define _ExtremaPC_Ellipse_HeaderFile
0016
0017 #include <ElCLib.hxx>
0018 #include <ExtremaPC.hxx>
0019 #include <gp_Elips.hxx>
0020 #include <gp_Pnt.hxx>
0021 #include <gp_Vec.hxx>
0022 #include <MathRoot_Trig.hxx>
0023 #include <Standard_DefineAlloc.hxx>
0024
0025 #include <cmath>
0026 #include <optional>
0027
0028 //! @brief Point-Ellipse extrema computation.
0029 //!
0030 //! Computes the extrema between a 3D point and an ellipse.
0031 //! Uses trigonometric equation solving via MathRoot::Trigonometric.
0032 //!
0033 //! The algorithm:
0034 //! 1. Projects point P onto the ellipse plane -> Pp
0035 //! 2. Solves: (B^2 - A^2)*cos(u)*sin(u) - B*Y*cos(u) + A*X*sin(u) = 0
0036 //!    where A = major radius, B = minor radius, and (X,Y) are coordinates
0037 //!    of Pp in the ellipse local coordinate system.
0038 //!
0039 //! @note Degenerate case: When P projects to the ellipse center and A = B
0040 //!       (i.e., the ellipse is a circle), returns Status::InfiniteSolutions.
0041 //!
0042 //! @note An ellipse can have up to 4 extrema.
0043 //!
0044 //! The domain is fixed at construction time for optimal performance.
0045 //! For full ellipse, construct without domain or with nullopt.
0046 class ExtremaPC_Ellipse
0047 {
0048 public:
0049   DEFINE_STANDARD_ALLOC
0050
0051   //! Constructor with ellipse geometry (full ellipse).
0052   //! @param[in] theEllipse the ellipse to compute extrema for
0053   explicit ExtremaPC_Ellipse(const gp_Elips& theEllipse)
0054       : myEllipse(theEllipse),
0055         myDomain(std::nullopt)
0056   {
0057   }
0058
0059   //! Constructor with ellipse geometry and parameter domain.
0060   //! @param[in] theEllipse the ellipse to compute extrema for
0061   //! @param[in] theDomain parameter domain in radians (fixed for all queries)
0062   ExtremaPC_Ellipse(const gp_Elips& theEllipse, const ExtremaPC::Domain1D& theDomain)
0063       : myEllipse(theEllipse),
0064         myDomain(theDomain.IsFullPeriod(2.0 * M_PI) ? std::nullopt
0065                                                     : std::optional<ExtremaPC::Domain1D>(theDomain))
0066   {
0067   }
0068
0069   //! Copy constructor is deleted.
0070   ExtremaPC_Ellipse(const ExtremaPC_Ellipse&) = delete;
0071
0072   //! Copy assignment operator is deleted.
0073   ExtremaPC_Ellipse& operator=(const ExtremaPC_Ellipse&) = delete;
0074
0075   //! Move constructor.
0076   ExtremaPC_Ellipse(ExtremaPC_Ellipse&&) = default;
0077
0078   //! Move assignment operator.
0079   ExtremaPC_Ellipse& operator=(ExtremaPC_Ellipse&&) = default;
0080
0081   //! Evaluates point on ellipse at parameter.
0082   //! @param theU parameter (radians)
0083   //! @return point on ellipse
0084   gp_Pnt Value(double theU) const { return ElCLib::Value(theU, myEllipse); }
0085
0086   //! Returns true if domain is bounded (partial arc).
0087   bool IsBounded() const { return myDomain.has_value(); }
0088
0089   //! Returns the domain (only valid if IsBounded() is true).
0090   const ExtremaPC::Domain1D& Domain() const { return *myDomain; }
0091
0092   //! Compute extrema between point P and the ellipse.
0093   //! Uses domain specified at construction time.
0094   //! @param theP query point
0095   //! @param theTol tolerance for degenerate case detection
0096   //! @param theMode search mode (MinMax, Min, or Max)
0097   //! @return const reference to result containing extrema or InfiniteSolutions status
0098   [[nodiscard]] const ExtremaPC::Result& Perform(
0099     const gp_Pnt&         theP,
0100     double                theTol,
0101     ExtremaPC::SearchMode theMode = ExtremaPC::SearchMode::MinMax) const
0102   {
0103     // Use stored domain or full parameter range [0, 2*PI]
0104     ExtremaPC::Domain1D aDomain = myDomain.value_or(ExtremaPC::Domain1D{0.0, 2.0 * M_PI});
0105     performCore(theP, aDomain, theTol, theMode);
0106     return myResult;
0107   }
0108
0109   //! Compute extrema between point P and the ellipse arc including endpoints.
0110   //! Uses domain specified at construction time.
0111   //! @param theP query point
0112   //! @param theTol tolerance for degenerate case detection
0113   //! @param theMode search mode (MinMax, Min, or Max)
0114   //! @return const reference to result containing interior + endpoint extrema or InfiniteSolutions
0115   //! status
0116   [[nodiscard]] const ExtremaPC::Result& PerformWithEndpoints(
0117     const gp_Pnt&         theP,
0118     double                theTol,
0119     ExtremaPC::SearchMode theMode = ExtremaPC::SearchMode::MinMax) const
0120   {
0121     (void)Perform(theP, theTol, theMode);
0122
0123     // Add endpoints if interior computation succeeded and domain is bounded
0124     if (myResult.Status == ExtremaPC::Status::OK && myDomain.has_value())
0125     {
0126       ExtremaPC::AddEndpointExtrema(myResult, theP, *myDomain, *this, theTol, theMode);
0127     }
0128
0129     return myResult;
0130   }
0131
0132   //! Returns the ellipse geometry.
0133   const gp_Elips& Ellipse() const { return myEllipse; }
0134
0135 private:
0136   //! Core algorithm - finds extrema with bounds checking.
0137   //! Stores results in myResult.
0138   void performCore(const gp_Pnt&              theP,
0139                    const ExtremaPC::Domain1D& theDomain,
0140                    double                     theTol,
0141                    ExtremaPC::SearchMode      theMode) const
0142   {
0143     myResult.Clear();
0144
0145     const double theUMin = theDomain.Min;
0146     const double theUMax = theDomain.Max;
0147
0148     // Step 1: Project point P onto the ellipse plane
0149     const gp_Pnt& aCenter = myEllipse.Location();
0150     const gp_Dir& aAxis   = myEllipse.Axis().Direction();
0151     gp_Vec        aToP(aCenter, theP);
0152     double        aHeight = aToP.Dot(gp_Vec(aAxis));
0153     gp_Vec        aTrsl   = gp_Vec(aAxis) * (-aHeight);
0154     gp_Pnt        aPp     = theP.Translated(aTrsl);
0155
0156     // Step 2: Get ellipse radii and compute local coordinates
0157     double aA = myEllipse.MajorRadius();
0158     double aB = myEllipse.MinorRadius();
0159
0160     gp_Vec aOPp(aCenter, aPp);
0161     double aOPpMag = aOPp.Magnitude();
0162
0163     // Check for degenerate case: point at center with circular ellipse
0164     if (aOPpMag < theTol)
0165     {
0166       if (std::abs(aA - aB) < theTol)
0167       {
0168         // Point at center of a circle - infinite solutions
0169         myResult.Status                 = ExtremaPC::Status::InfiniteSolutions;
0170         myResult.InfiniteSquareDistance = aA * aA + aHeight * aHeight;
0171         return;
0172       }
0173       // For non-circular ellipse at center, we still get valid extrema at semi-axes
0174     }
0175
0176     // Local coordinates
0177     double aX = aOPp.Dot(gp_Vec(myEllipse.XAxis().Direction()));
0178     double aY = aOPp.Dot(gp_Vec(myEllipse.YAxis().Direction()));
0179
0180     // Step 3: Solve trigonometric equation
0181     // (B^2 - A^2)*cos*sin - B*Y*cos + A*X*sin = 0
0182     // In MathRoot::Trigonometric form: a*cos^2 + 2*b*cos*sin + c*cos + d*sin + e = 0
0183     // a = 0, 2*b = (B^2 - A^2), c = -B*Y, d = A*X, e = 0
0184     double aKo2 = (aB * aB - aA * aA) / 2.0;
0185     double aKo3 = -aB * aY;
0186     double aKo4 = aA * aX;
0187
0188     // MathRoot::Trigonometric handles all special cases including Y ~= 0
0189     MathRoot::TrigResult aTrigRes =
0190       MathRoot::Trigonometric(0.0, aKo2, aKo3, aKo4, 0.0, theUMin, theUMax);
0191
0192     if (!aTrigRes.IsDone())
0193     {
0194       if (aTrigRes.InfiniteRoots)
0195       {
0196         myResult.Status = ExtremaPC::Status::InfiniteSolutions;
0197         // For infinite case, compute distance to any point (use U=0)
0198         gp_Pnt aPtOnCurve               = ElCLib::Value(0.0, myEllipse);
0199         myResult.InfiniteSquareDistance = theP.SquareDistance(aPtOnCurve);
0200       }
0201       else
0202       {
0203         myResult.Status = ExtremaPC::Status::NumericalError;
0204       }
0205       return;
0206     }
0207
0208     // Step 4: Collect extrema
0209     double aTol2 = theTol * theTol;
0210
0211     auto addExtremum = [&](double aU) {
0212       gp_Pnt aCurvePt = ElCLib::Value(aU, myEllipse);
0213
0214       // Check for duplicates using parameter proximity (more robust)
0215       for (int j = 0; j < myResult.Extrema.Length(); ++j)
0216       {
0217         if (std::abs(myResult.Extrema.Value(j).Parameter - aU) < theTol)
0218         {
0219           return;
0220         }
0221         if (aCurvePt.SquareDistance(myResult.Extrema.Value(j).Point) < aTol2)
0222         {
0223           return;
0224         }
0225       }
0226
0227       double aSqDist = theP.SquareDistance(aCurvePt);
0228
0229       // Determine if this is a minimum or maximum by checking neighboring points
0230       // Use step relative to parameter range
0231       double aStep      = std::max(ExtremaPC::THE_NEIGHBOR_STEP_RATIO,
0232                               (theUMax - theUMin) * ExtremaPC::THE_NEIGHBOR_STEP_RATIO);
0233       gp_Pnt aPtPlus    = ElCLib::Value(aU + aStep, myEllipse);
0234       gp_Pnt aPtMinus   = ElCLib::Value(aU - aStep, myEllipse);
0235       double aDistPlus  = theP.SquareDistance(aPtPlus);
0236       double aDistMinus = theP.SquareDistance(aPtMinus);
0237       bool   aIsMin     = (aSqDist <= aDistPlus) && (aSqDist <= aDistMinus);
0238
0239       // Filter by search mode
0240       if (theMode == ExtremaPC::SearchMode::Min && !aIsMin)
0241       {
0242         return;
0243       }
0244       if (theMode == ExtremaPC::SearchMode::Max && aIsMin)
0245       {
0246         return;
0247       }
0248
0249       ExtremaPC::ExtremumResult anExt;
0250       anExt.Parameter      = aU;
0251       anExt.Point          = aCurvePt;
0252       anExt.SquareDistance = aSqDist;
0253       anExt.IsMinimum      = aIsMin;
0254
0255       myResult.Extrema.Append(anExt);
0256     };
0257
0258     // Add solutions from solver
0259     for (int i = 0; i < aTrigRes.NbRoots; ++i)
0260     {
0261       addExtremum(aTrigRes.Roots[i]);
0262     }
0263
0264     myResult.Status = ExtremaPC::Status::OK;
0265   }
0266
0267   gp_Elips                           myEllipse; //!< Ellipse geometry
0268   std::optional<ExtremaPC::Domain1D> myDomain;  //!< Parameter domain (nullopt for full ellipse)
0269   mutable ExtremaPC::Result          myResult;  //!< Reusable result storage
0270 };
0271
0272 #endif // _ExtremaPC_Ellipse_HeaderFile