src/Extrema/Extrema_ExtPElC.cxx

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