src/Extrema/Extrema_ExtPElC.cxx

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