src/Extrema/Extrema_ExtPElC.cxx

   1 #include <Extrema_ExtPElC.ixx>
   2 #include <StdFail_NotDone.hxx>
   3 #include <math_DirectPolynomialRoots.hxx>
   4 #include <math_TrigonometricFunctionRoots.hxx>
   5 #include <ElCLib.hxx>
   6 #include <Standard_OutOfRange.hxx>
   7 #include <Standard_NotImplemented.hxx>
   8 #include <Precision.hxx>
   9
  10
  11 //=============================================================================
  12
  13 Extrema_ExtPElC::Extrema_ExtPElC () { myDone = Standard_False; }
  14 //=============================================================================
  15
  16 Extrema_ExtPElC::Extrema_ExtPElC (const gp_Pnt&       P,
  17                                   const gp_Lin&       L,
  18                                   const Standard_Real Tol,
  19                                   const Standard_Real Uinf,
  20                                   const Standard_Real Usup)
  21 {
  22   Perform(P, L, Tol, Uinf, Usup);
  23 }
  24
  25 void Extrema_ExtPElC::Perform(const gp_Pnt&       P,
  26                               const gp_Lin&       L,
  27                               const Standard_Real Tol,
  28                               const Standard_Real Uinf,
  29                               const Standard_Real Usup)
  30 {
  31   myDone = Standard_False;
  32   myNbExt = 0;
  33   gp_Vec V1 = gp_Vec(L.Direction());
  34   gp_Pnt OR = L.Location();
  35   gp_Vec V(OR, P);
  36   Standard_Real Mydist = V1.Dot(V);
  37   if ((Mydist >= Uinf-Tol) &&
  38       (Mydist <= Usup+Tol)){
  39
  40     gp_Pnt MyP = OR.Translated(Mydist*V1);
  41     Extrema_POnCurv MyPOnCurve(Mydist, MyP);
  42     mySqDist[0] = P.SquareDistance(MyP);
  43     myPoint[0] = MyPOnCurve;
  44     myIsMin[0] = Standard_True;
  45     myNbExt = 1;
  46     myDone = Standard_True;
  47   }
  48
  49 }
  50
  51
  52
  53 Extrema_ExtPElC::Extrema_ExtPElC (const gp_Pnt&       P,
  54                                   const gp_Circ&      C,
  55                                   const Standard_Real Tol,
  56                                   const Standard_Real Uinf,
  57                                   const Standard_Real Usup)
  58 {
  59   Perform(P, C, Tol, Uinf, Usup);
  60 }
  61
  62 void Extrema_ExtPElC::Perform(const gp_Pnt&       P,
  63                               const gp_Circ&      C,
  64                               const Standard_Real Tol,
  65                               const Standard_Real Uinf,
  66                               const Standard_Real Usup)
  67 /*-----------------------------------------------------------------------------
  68 Function:
  69   Find values of parameter u such as:
  70    - dist(P,C(u)) pass by an extrema,
  71    - Uinf <= u <= Usup.
  72
  73 Method:
  74   Pass 3 stages:
  75   1- Projection of point P in the plane of the circle,
  76   2- Calculation of u solutions in [0.,2.*M_PI]:
  77       Let Pp, the projected point and
  78            O, the center of the circle;
  79      2 cases:
  80      - if Pp is mixed with 0, there is an infinite number of solutions;
  81        IsDone() renvoie Standard_False.
  82      - otherwise, 2 points are solutions for the complete circle:
  83        . Us1 = angle(OPp,OX) corresponds to the minimum,
  84        . let Us2 = ( Us1 + M_PI if Us1 < M_PI,
  85                     ( Us1 - M_PI otherwise;
  86          Us2 corresponds to the maximum.
  87   3- Calculate the extrema in [Uinf,Usup].
  88 -----------------------------------------------------------------------------*/
  89 {
  90   myDone = Standard_False;
  91   myNbExt = 0;
  92
  93 // 1- Projection of the point P in the plane of circle -> Pp ...
  94
  95   gp_Pnt O = C.Location();
  96   gp_Vec Axe (C.Axis().Direction());
  97   gp_Vec Trsl = Axe.Multiplied(-(gp_Vec(O,P).Dot(Axe)));
  98   gp_Pnt Pp = P.Translated(Trsl);
  99
 100 // 2- Calculate u solutions in [0.,2.*PI] ...
 101
 102   gp_Vec OPp (O,Pp);
 103   if (OPp.Magnitude() < Tol) { return; }
 104   Standard_Real Usol[2];
 105   Usol[0] = C.XAxis().Direction().AngleWithRef(OPp,Axe); // -M_PI<U1<M_PI
 106   Usol[1] = Usol[0] + M_PI;
 107
 108   Standard_Real myuinf = Uinf;
 109   //modified by NIZNHY-PKV Fri Apr 20 15:03:28 2001 f
 110   //Standard_Real TolU = Tol*C.Radius();
 111   Standard_Real TolU, aR;
 112   aR=C.Radius();
 113   TolU=Precision::Infinite();
 114   if (aR > gp::Resolution()) {
 115     TolU= Tol/aR;
 116   }
 117   //modified by NIZNHY-PKV Fri Apr 20 15:03:32 2001 t
 118   ElCLib::AdjustPeriodic(Uinf, Uinf+2*M_PI, TolU, myuinf, Usol[0]);
 119   ElCLib::AdjustPeriodic(Uinf, Uinf+2*M_PI, TolU, myuinf, Usol[1]);
 120   if (((Usol[0]-2*M_PI-Uinf) < TolU) && ((Usol[0]-2*M_PI-Uinf) > -TolU)) Usol[0] = Uinf;
 121   if (((Usol[1]-2*M_PI-Uinf) < TolU) && ((Usol[1]-2*M_PI-Uinf) > -TolU)) Usol[1] = Uinf;
 122
 123
 124 // 3- Calculate extrema in [Umin,Umax] ...
 125
 126   gp_Pnt Cu;
 127   Standard_Real Us;
 128   for (Standard_Integer NoSol = 0; NoSol <= 1; NoSol++) {
 129     Us = Usol[NoSol];
 130     if (((Uinf-Us) < TolU) && ((Us-Usup) < TolU)) {
 131       Cu = ElCLib::Value(Us,C);
 132       mySqDist[myNbExt] = Cu.SquareDistance(P);
 133       myIsMin[myNbExt] = (NoSol == 0);
 134       myPoint[myNbExt] = Extrema_POnCurv(Us,Cu);
 135       myNbExt++;
 136     }
 137   }
 138   myDone = Standard_True;
 139 }
 140 //=============================================================================
 141
 142 Extrema_ExtPElC::Extrema_ExtPElC (const gp_Pnt&       P,
 143                                   const gp_Elips&     C,
 144                                   const Standard_Real Tol,
 145                                   const Standard_Real Uinf,
 146                                   const Standard_Real Usup)
 147 {
 148   Perform(P, C, Tol, Uinf, Usup);
 149 }
 150
 151
 152
 153 void Extrema_ExtPElC::Perform (const gp_Pnt&       P,
 154                                const gp_Elips&     C,
 155                                const Standard_Real Tol,
 156                                const Standard_Real Uinf,
 157                                const Standard_Real Usup)
 158 /*-----------------------------------------------------------------------------
 159 Function:
 160   Find values of parameter u so that:
 161    - dist(P,C(u)) passes by an extremum,
 162    - Uinf <= u <= Usup.
 163
 164 Method:
 165   Takes 2 steps:
 166   1- Projection of point P in the plane of the ellipse,
 167   2- Calculation of the solutions:
 168      Let Pp, the projected point; find values u so that:
 169        (C(u)-Pp).C'(u) = 0. (1)
 170      Let Cos = cos(u) and Sin = sin(u),
 171           C(u) = (A*Cos,B*Sin) and Pp = (X,Y);
 172      Then, (1) <=> (A*Cos-X,B*Sin-Y).(-A*Sin,B*Cos) = 0.
 173                     (B**2-A**2)*Cos*Sin - B*Y*Cos + A*X*Sin = 0.
 174      Use algorithm math_TrigonometricFunctionRoots to solve this equation.
 175 -----------------------------------------------------------------------------*/
 176 {
 177   myDone = Standard_False;
 178   myNbExt = 0;
 179
 180 // 1- Projection of point P in the plane of the ellipse -> Pp ...
 181
 182   gp_Pnt O = C.Location();
 183   gp_Vec Axe (C.Axis().Direction());
 184   gp_Vec Trsl = Axe.Multiplied(-(gp_Vec(O,P).Dot(Axe)));
 185   gp_Pnt Pp = P.Translated(Trsl);
 186
 187 // 2- Calculation of solutions ...
 188
 189   Standard_Integer NoSol, NbSol;
 190   Standard_Real A = C.MajorRadius();
 191   Standard_Real B = C.MinorRadius();
 192   gp_Vec OPp (O,Pp);
 193   Standard_Real OPpMagn = OPp.Magnitude();
 194   if (OPpMagn < Tol) { if (Abs(A-B) < Tol) { return; } }
 195   Standard_Real X = OPp.Dot(gp_Vec(C.XAxis().Direction()));
 196   Standard_Real Y = OPp.Dot(gp_Vec(C.YAxis().Direction()));
 197   //  Standard_Real Y = Sqrt(OPpMagn*OPpMagn-X*X);
 198
 199   Standard_Real ko2 = (B*B-A*A)/2., ko3 = -B*Y, ko4 = A*X;
 200   if(Abs(ko3) < 1.e-16*Max(Abs(ko2), Abs(ko3))) ko3 = 0.0;
 201
 202 //  math_TrigonometricFunctionRoots Sol(0.,(B*B-A*A)/2.,-B*Y,A*X,0.,Uinf,Usup);
 203   math_TrigonometricFunctionRoots Sol(0.,ko2, ko3, ko4, 0.,Uinf,Usup);
 204
 205   if (!Sol.IsDone()) { return; }
 206   gp_Pnt Cu;
 207   Standard_Real Us;
 208   NbSol = Sol.NbSolutions();
 209   for (NoSol = 1; NoSol <= NbSol; NoSol++) {
 210     Us = Sol.Value(NoSol);
 211     Cu = ElCLib::Value(Us,C);
 212     mySqDist[myNbExt] = Cu.SquareDistance(P);
 213     myIsMin[myNbExt] = (NoSol == 1);
 214     myPoint[myNbExt] = Extrema_POnCurv(Us,Cu);
 215     myNbExt++;
 216   }
 217   myDone = Standard_True;
 218 }
 219 //=============================================================================
 220
 221 Extrema_ExtPElC::Extrema_ExtPElC (const gp_Pnt&       P,
 222                                   const gp_Hypr&      C,
 223                                   const Standard_Real Tol,
 224                                   const Standard_Real Uinf,
 225                                   const Standard_Real Usup)
 226 {
 227   Perform(P, C, Tol, Uinf, Usup);
 228 }
 229
 230
 231 void Extrema_ExtPElC::Perform(const gp_Pnt&       P,
 232                               const gp_Hypr&      C,
 233                               const Standard_Real Tol,
 234                               const Standard_Real Uinf,
 235                               const Standard_Real Usup)
 236 /*-----------------------------------------------------------------------------
 237 Function:
 238   Find values of parameter u so that:
 239    - dist(P,C(u)) passes by an extremum,
 240    - Uinf <= u <= Usup.
 241
 242 Method:
 243   Takes 2 steps:
 244   1- Projection of point P in the plane of the hyperbola,
 245   2- Calculation of solutions:
 246       Let Pp, le point projete; on recherche les valeurs u telles que:
 247        (C(u)-Pp).C'(u) = 0. (1)
 248       Let R and r be the radiuses of the hyperbola,
 249           Chu = Cosh(u) and Shu = Sinh(u),
 250           C(u) = (R*Chu,r*Shu) and Pp = (X,Y);
 251      Then, (1) <=> (R*Chu-X,r*Shu-Y).(R*Shu,r*Chu) = 0.
 252                     (R**2+r**2)*Chu*Shu - X*R*Shu - Y*r*Chu = 0. (2)
 253      Let v = e**u;
 254      Then, by using Chu = (e**u+e**(-u))/2. and Sh = (e**u-e**(-u)))/2.
 255             (2) <=> ((R**2+r**2)/4.) * (v**2-v**(-2)) -
 256                     ((X*R+Y*r)/2.) * v +
 257                     ((X*R-Y*r)/2.) * v**(-1) = 0.
 258       (2)* v**2 <=> ((R**2+r**2)/4.) * v**4 -
 259                      ((X*R+Y*r)/2.)  * v**3 +
 260                      ((X*R-Y*r)/2.)  * v    -
 261                     ((R**2+r**2)/4.)        = 0.
 262   Use algorithm math_DirectPolynomialRoots to solve this equation by v.
 263 -----------------------------------------------------------------------------*/
 264 {
 265   myDone = Standard_False;
 266   myNbExt = 0;
 267
 268 // 1- Projection of point P in the plane of hyperbola -> Pp ...
 269
 270   gp_Pnt O = C.Location();
 271   gp_Vec Axe (C.Axis().Direction());
 272   gp_Vec Trsl = Axe.Multiplied(-(gp_Vec(O,P).Dot(Axe)));
 273   gp_Pnt Pp = P.Translated(Trsl);
 274
 275 // 2- Calculation of solutions ...
 276
 277   Standard_Real Tol2 = Tol * Tol;
 278   Standard_Real R = C.MajorRadius();
 279   Standard_Real r = C.MinorRadius();
 280   gp_Vec OPp (O,Pp);
 281 #ifdef DEB
 282   Standard_Real OPpMagn = OPp.Magnitude();
 283 #else
 284   OPp.Magnitude();
 285 #endif
 286   Standard_Real X = OPp.Dot(gp_Vec(C.XAxis().Direction()));
 287   Standard_Real Y = OPp.Dot(gp_Vec(C.YAxis().Direction()));
 288
 289   Standard_Real C1 = (R*R+r*r)/4.;
 290   math_DirectPolynomialRoots Sol(C1,-(X*R+Y*r)/2.,0.,(X*R-Y*r)/2.,-C1);
 291   if (!Sol.IsDone()) { return; }
 292   gp_Pnt Cu;
 293   Standard_Real Us, Vs;
 294   Standard_Integer NbSol = Sol.NbSolutions();
 295   Standard_Boolean DejaEnr;
 296   Standard_Integer NoExt;
 297   gp_Pnt TbExt[4];
 298   for (Standard_Integer NoSol = 1; NoSol <= NbSol; NoSol++) {
 299     Vs = Sol.Value(NoSol);
 300     if (Vs > 0.) {
 301       Us = Log(Vs);
 302       if ((Us >= Uinf) && (Us <= Usup)) {
 303         Cu = ElCLib::Value(Us,C);
 304         DejaEnr = Standard_False;
 305         for (NoExt = 0; NoExt < myNbExt; NoExt++) {
 306           if (TbExt[NoExt].SquareDistance(Cu) < Tol2) {
 307             DejaEnr = Standard_True;
 308             break;
 309           }
 310         }
 311         if (!DejaEnr) {
 312           TbExt[myNbExt] = Cu;
 313           mySqDist[myNbExt] = Cu.SquareDistance(P);
 314           myIsMin[myNbExt] = mySqDist[myNbExt] < P.SquareDistance(ElCLib::Value(Us+1,C));
 315           myPoint[myNbExt] = Extrema_POnCurv(Us,Cu);
 316           myNbExt++;
 317         }
 318       } // if ((Us >= Uinf) && (Us <= Usup))
 319     } // if (Vs > 0.)
 320   } // for (Standard_Integer NoSol = 1; ...
 321   myDone = Standard_True;
 322 }
 323 //=============================================================================
 324
 325 Extrema_ExtPElC::Extrema_ExtPElC (const gp_Pnt&       P,
 326                                   const gp_Parab&     C,
 327                                   const Standard_Real Tol,
 328                                   const Standard_Real Uinf,
 329                                   const Standard_Real Usup)
 330 {
 331   Perform(P, C, Tol, Uinf, Usup);
 332 }
 333
 334
 335 void Extrema_ExtPElC::Perform(const gp_Pnt&       P,
 336                               const gp_Parab&     C,
 337 //                            const Standard_Real Tol,
 338                               const Standard_Real ,
 339                               const Standard_Real Uinf,
 340                               const Standard_Real Usup)
 341 /*-----------------------------------------------------------------------------
 342 Function:
 343   Find values of parameter u so that:
 344    - dist(P,C(u)) pass by an extremum,
 345    - Uinf <= u <= Usup.
 346
 347 Method:
 348   Takes 2 steps:
 349   1- Projection of point P in the plane of the parabola,
 350   2- Calculation of solutions:
 351      Let Pp, the projected point; find values u so that:
 352        (C(u)-Pp).C'(u) = 0. (1)
 353      Let F the focus of the parabola,
 354           C(u) = ((u*u)/(4.*F),u) and Pp = (X,Y);
 355      Alors, (1) <=> ((u*u)/(4.*F)-X,u-Y).(u/(2.*F),1) = 0.
 356                     (1./(4.*F)) * U**3 + (2.*F-X) * U - 2*F*Y = 0.
 357     Use algorithm math_DirectPolynomialRoots to solve this equation by U.
 358 -----------------------------------------------------------------------------*/
 359 {
 360   myDone = Standard_False;
 361   myNbExt = 0;
 362
 363 // 1- Projection of point P in the plane of the parabola -> Pp ...
 364
 365   gp_Pnt O = C.Location();
 366   gp_Vec Axe (C.Axis().Direction());
 367   gp_Vec Trsl = Axe.Multiplied(-(gp_Vec(O,P).Dot(Axe)));
 368   gp_Pnt Pp = P.Translated(Trsl);
 369
 370 // 2- Calculation of solutions ...
 371
 372   Standard_Real F = C.Focal();
 373   gp_Vec OPp (O,Pp);
 374 #ifdef DEB
 375   Standard_Real OPpMagn = OPp.Magnitude();
 376 #else
 377   OPp.Magnitude();
 378 #endif
 379   Standard_Real X = OPp.Dot(gp_Vec(C.XAxis().Direction()));
 380 //  Standard_Real Y = Sqrt(OPpMagn*OPpMagn-X*X);
 381   Standard_Real Y = OPp.Dot(gp_Vec(C.YAxis().Direction()));
 382   math_DirectPolynomialRoots Sol(1./(4.*F),0.,2.*F-X,-2.*F*Y);
 383   if (!Sol.IsDone()) { return; }
 384   gp_Pnt Cu;
 385   Standard_Real Us;
 386   Standard_Integer NbSol = Sol.NbSolutions();
 387   Standard_Boolean DejaEnr;
 388   Standard_Integer NoExt;
 389   gp_Pnt TbExt[3];
 390   for (Standard_Integer NoSol = 1; NoSol <= NbSol; NoSol++) {
 391     Us = Sol.Value(NoSol);
 392     if ((Us >= Uinf) && (Us <= Usup)) {
 393       Cu = ElCLib::Value(Us,C);
 394       DejaEnr = Standard_False;
 395       for (NoExt = 0; NoExt < myNbExt; NoExt++) {
 396         if (TbExt[NoExt].SquareDistance(Cu) < Precision::Confusion() * Precision::Confusion()) {
 397           DejaEnr = Standard_True;
 398           break;
 399         }
 400       }
 401       if (!DejaEnr) {
 402         TbExt[myNbExt] = Cu;
 403         mySqDist[myNbExt] = Cu.SquareDistance(P);
 404         myIsMin[myNbExt] = mySqDist[myNbExt] < P.SquareDistance(ElCLib::Value(Us+1,C));
 405         myPoint[myNbExt] = Extrema_POnCurv(Us,Cu);
 406         myNbExt++;
 407       }
 408     } // if ((Us >= Uinf) && (Us <= Usup))
 409   } // for (Standard_Integer NoSol = 1; ...
 410   myDone = Standard_True;
 411 }
 412 //=============================================================================
 413
 414 Standard_Boolean Extrema_ExtPElC::IsDone () const { return myDone; }
 415 //=============================================================================
 416
 417 Standard_Integer Extrema_ExtPElC::NbExt () const
 418 {
 419   if (!IsDone()) { StdFail_NotDone::Raise(); }
 420   return myNbExt;
 421 }
 422 //=============================================================================
 423
 424 Standard_Real Extrema_ExtPElC::SquareDistance (const Standard_Integer N) const
 425 {
 426   if ((N < 1) || (N > NbExt())) { Standard_OutOfRange::Raise(); }
 427   return mySqDist[N-1];
 428 }
 429 //=============================================================================
 430
 431 Standard_Boolean Extrema_ExtPElC::IsMin (const Standard_Integer N) const
 432 {
 433   if ((N < 1) || (N > NbExt())) { Standard_OutOfRange::Raise(); }
 434   return myIsMin[N-1];
 435 }
 436 //=============================================================================
 437
 438 Extrema_POnCurv Extrema_ExtPElC::Point (const Standard_Integer N) const
 439 {
 440   if ((N < 1) || (N > NbExt())) { Standard_OutOfRange::Raise(); }
 441   return myPoint[N-1];
 442 }
 443 //=============================================================================