src/Extrema/Extrema_ExtPElS.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
  16 #include <ElSLib.hxx>
  17 #include <Extrema_ExtPElS.hxx>
  18 #include <Extrema_POnSurf.hxx>
  19 #include <gp_Cone.hxx>
  20 #include <gp_Cylinder.hxx>
  21 #include <gp_Pln.hxx>
  22 #include <gp_Pnt.hxx>
  23 #include <gp_Sphere.hxx>
  24 #include <gp_Torus.hxx>
  25 #include <Standard_NotImplemented.hxx>
  26 #include <Standard_OutOfRange.hxx>
  27 #include <StdFail_NotDone.hxx>
  28
  29 static const Standard_Real ExtPElS_MyEps = Epsilon(2. * M_PI);
  30 //=============================================================================
  31
  32 Extrema_ExtPElS::Extrema_ExtPElS () { myDone = Standard_False; }
  33 //=============================================================================
  34
  35 Extrema_ExtPElS::Extrema_ExtPElS (const gp_Pnt& P,
  36                                   const gp_Cylinder& S,
  37                                   const Standard_Real Tol)
  38 {
  39
  40  Perform(P, S, Tol);
  41 }
  42 /*-----------------------------------------------------------------------------
  43 Function:
  44 Find 2 extreme distances between point P and cylinder S.
  45
  46 Method:
  47   Let Pp be the projection of P in plane XOY of the cylinder;
  48   2 cases are considered:
  49   1- distance(Pp,O) < Tol:
  50      There are infinite solutions; IsDone() = Standard_False.
  51   2- distance(Pp,O) > Tol:
  52      let V = OP.OZ,
  53           U1 = angle(OX,OPp) with 0 < U1 < 2.*M_PI
  54           U2 = U1 + M_PI with 0 < U2 < 2.*M_PI;
  55      then (U1,V) corresponds to the min distance.
  56      and  (U2,V) corresponds to the max distance.
  57 -----------------------------------------------------------------------------*/
  58
  59 void Extrema_ExtPElS::Perform(const gp_Pnt&       P,
  60                               const gp_Cylinder&  S,
  61                               const Standard_Real Tol)
  62 {
  63   myDone = Standard_False;
  64   myNbExt = 0;
  65
  66 // Projection of point P in plane XOY of the cylinder ...
  67   gp_Ax3 Pos = S.Position();
  68   gp_Pnt O = Pos.Location();
  69   gp_Vec OZ (Pos.Direction());
  70   Standard_Real V = gp_Vec(O,P).Dot(OZ);
  71   gp_Pnt Pp = P.Translated(OZ.Multiplied(-V));
  72
  73 // Calculation of extrema
  74   gp_Vec OPp (O,Pp);
  75   if (OPp.Magnitude() < Tol) { return; }
  76   gp_Vec myZ = Pos.XDirection()^Pos.YDirection();
  77   Standard_Real U1 = gp_Vec(Pos.XDirection()).AngleWithRef(OPp,myZ); //-M_PI<U1<M_PI
  78   if (U1 > -ExtPElS_MyEps && U1 < ExtPElS_MyEps) { U1 = 0.; }
  79   Standard_Real U2 = U1 + M_PI;
  80   if (U1 < 0.) { U1 += 2. * M_PI; }
  81
  82   gp_Pnt Ps;
  83   Ps = ElSLib::Value(U1,V,S);
  84   mySqDist[0] = Ps.SquareDistance(P);
  85   myPoint[0] = Extrema_POnSurf(U1,V,Ps);
  86   Ps = ElSLib::Value(U2,V,S);
  87   mySqDist[1] = Ps.SquareDistance(P);
  88   myPoint[1] = Extrema_POnSurf(U2,V,Ps);
  89
  90   myNbExt = 2;
  91   myDone = Standard_True;
  92 }
  93 //=============================================================================
  94
  95 Extrema_ExtPElS::Extrema_ExtPElS (const gp_Pnt&       P,
  96                                   const gp_Cone&      S,
  97                                   const Standard_Real Tol)
  98 {
  99   Perform(P, S, Tol);
 100 }
 101 /*-----------------------------------------------------------------------------
 102 Function:
 103    Find 2 extreme distances between point P and cone S.
 104
 105 Method:
 106   Let M the top of the cone.
 107   2 cases are considered:
 108   1- distance(P,M) < Tol:
 109      there is a minimum in M.
 110   2- distance(P,M) > Tol:
 111      Let Pp the projection of P in the plane XOY of the cone;
 112      2 cases are considered:
 113      1- distance(Pp,O) < Tol:
 114      There is an infinite number of solutions; IsDone() = Standard_False.
 115      2- distance(Pp,O) > Tol:
 116         There exist 2 extrema:
 117         Let Vm = value of v for point M,
 118              Vp = value of v for point P,
 119              U1 = angle(OX,OPp) if Vp > Vm )
 120                  -angle(OX,OPp) otherwise      ) with 0. < U1 < 2*M_PI,
 121              U2 = U1 + M_PI with 0. < U2 < 2*M_PI;
 122         We are in plane PpOZ.
 123        Let A the angle of the cone,
 124              B = angle(MP,MO) with 0. < B < M_PI,
 125              L = longueur(MP),
 126              V1 = (L * cos(B-A)) + Vm,
 127              V2 = (L * cos(B+A)) + Vm;
 128        then (U1,V1) and (U2,V2) correspond to min distances.
 129 -----------------------------------------------------------------------------*/
 130
 131 void Extrema_ExtPElS::Perform(const gp_Pnt&       P,
 132                               const gp_Cone&      S,
 133                               const Standard_Real Tol)
 134 {
 135   myDone = Standard_False;
 136   myNbExt = 0;
 137
 138   gp_Pnt M = S.Apex();
 139   gp_Ax3 Pos = S.Position();
 140   gp_Pnt O = Pos.Location();
 141   Standard_Real A = S.SemiAngle();
 142   gp_Vec OZ (Pos.Direction());
 143   gp_Vec myZ = Pos.XDirection()^Pos.YDirection();
 144   gp_Vec MP (M,P);
 145
 146   Standard_Real L2 = MP.SquareMagnitude();
 147   Standard_Real Vm = -(S.RefRadius() / Sin(A));
 148
 149 // Case when P is mixed with S ...
 150   if (L2 < Tol * Tol) {
 151     mySqDist[0] = L2;
 152     myPoint[0] = Extrema_POnSurf(0.,Vm,M);
 153     myNbExt = 1;
 154     myDone = Standard_True;
 155     return;
 156   }
 157     gp_Vec DirZ;
 158     if (M.SquareDistance(O)<Tol * Tol)
 159     { DirZ=OZ;
 160      if( A<0) DirZ.Multiplied(-1.);
 161     }
 162     else
 163      DirZ=gp_Vec(M,O);
 164 // Projection of P in the reference plane of the cone ...
 165   Standard_Real Zp = gp_Vec(O, P).Dot(OZ);
 166
 167   gp_Pnt Pp = P.Translated(OZ.Multiplied(-Zp));
 168   gp_Vec OPp(O, Pp);
 169   if (OPp.SquareMagnitude() < Tol * Tol) return;
 170   Standard_Real B, U1, V1, U2, V2;
 171   Standard_Boolean Same = DirZ.Dot(MP) >= 0.0;
 172   U1 = gp_Vec(Pos.XDirection()).AngleWithRef(OPp,myZ); //-M_PI<U1<M_PI
 173   if (U1 > -ExtPElS_MyEps && U1 < ExtPElS_MyEps) { U1 = 0.; }
 174   B = MP.Angle(DirZ);
 175   if (!Same) { U1 += M_PI; }
 176   U2 = U1 + M_PI;
 177   if (U1 < 0.) { U1 += 2. * M_PI; }
 178   if (U2 > 2.*M_PI) { U2 -= 2. * M_PI; }
 179   B = MP.Angle(DirZ);
 180   A = Abs(A);
 181   Standard_Real L = sqrt(L2);
 182   if (!Same) {
 183     B = M_PI-B;
 184     V1 = -L*cos(B-A);
 185     V2 = -L*cos(B+A);
 186   }
 187   else {
 188     V1 = L * cos(B-A);
 189     V2 = L * cos(B+A);
 190   }
 191   Standard_Real Sense = OZ.Dot(gp_Dir(DirZ));
 192   V1 *= Sense;   V2 *= Sense;
 193   V1 += Vm;   V2 += Vm;
 194
 195   gp_Pnt Ps;
 196   Ps = ElSLib::Value(U1,V1,S);
 197   mySqDist[0] = Ps.SquareDistance(P);
 198   myPoint[0] = Extrema_POnSurf(U1,V1,Ps);
 199   Ps = ElSLib::Value(U2,V2,S);
 200   mySqDist[1] = Ps.SquareDistance(P);
 201   myPoint[1] = Extrema_POnSurf(U2,V2,Ps);
 202
 203   myNbExt = 2;
 204   myDone = Standard_True;
 205 }
 206 //=============================================================================
 207
 208 Extrema_ExtPElS::Extrema_ExtPElS (const gp_Pnt&       P,
 209                                   const gp_Sphere&    S,
 210                                   const Standard_Real Tol)
 211 {
 212   Perform(P, S, Tol);
 213 }
 214 /*-----------------------------------------------------------------------------
 215 Function:
 216   Find 2 extreme distances between point P and sphere S.
 217
 218 Method:
 219    Let O be the origin of the sphere.
 220   2 cases are considered:
 221   1- distance(P,O) < Tol:
 222      There is an infinite number of solutions; IsDone() = Standard_False
 223   2- distance(P,O) > Tol:
 224      Let Pp be the projection of point P in the plane XOY of the sphere;
 225      2 cases are considered:
 226      1- distance(Pp,O) < Tol:
 227         2 solutions are: (0,-M_PI/2.) and (0.,M_PI/2.)
 228      2- distance(Pp,O) > Tol:
 229         Let U1 = angle(OX,OPp) with 0. < U1 < 2.*M_PI,
 230              U2 = U1 + M_PI avec 0. < U2 < 2*M_PI,
 231              V1 = angle(OPp,OP) with -M_PI/2. < V1 < M_PI/2. ,
 232         then (U1, V1) corresponds to the min distance
 233         and  (U2,-V1) corresponds to the max distance.
 234 -----------------------------------------------------------------------------*/
 235
 236 void Extrema_ExtPElS::Perform(const gp_Pnt&       P,
 237                               const gp_Sphere&    S,
 238                               const Standard_Real Tol)
 239 {
 240   myDone = Standard_False;
 241   myNbExt = 0;
 242
 243   gp_Ax3 Pos = S.Position();
 244   gp_Vec OP (Pos.Location(),P);
 245
 246 // Case when P is mixed with O ...
 247   if (OP.SquareMagnitude() < Tol * Tol) { return; }
 248
 249 // Projection if P in plane XOY of the sphere ...
 250   gp_Pnt O = Pos.Location();
 251   gp_Vec OZ (Pos.Direction());
 252   Standard_Real Zp = OP.Dot(OZ);
 253   gp_Pnt Pp = P.Translated(OZ.Multiplied(-Zp));
 254
 255 // Calculation of extrema ...
 256   gp_Vec OPp (O,Pp);
 257   Standard_Real U1, U2, V;
 258   if (OPp.SquareMagnitude() < Tol * Tol) {
 259     U1 = 0.;
 260     U2 = 0.;
 261     if (Zp < 0.) { V = -M_PI / 2.; }
 262     else { V = M_PI / 2.; }
 263   }
 264   else {
 265     gp_Vec myZ = Pos.XDirection()^Pos.YDirection();
 266     U1 = gp_Vec(Pos.XDirection()).AngleWithRef(OPp,myZ);
 267     if (U1 > -ExtPElS_MyEps && U1 < ExtPElS_MyEps) { U1 = 0.; }
 268     U2 = U1 + M_PI;
 269     if (U1 < 0.) { U1 += 2. * M_PI; }
 270     V = OP.Angle(OPp);
 271     if (Zp < 0.) { V = -V; }
 272   }
 273
 274   gp_Pnt Ps;
 275   Ps = ElSLib::Value(U1,V,S);
 276   mySqDist[0] = Ps.SquareDistance(P);
 277   myPoint[0] = Extrema_POnSurf(U1,V,Ps);
 278   Ps = ElSLib::Value(U2,-V,S);
 279   mySqDist[1] = Ps.SquareDistance(P);
 280   myPoint[1] = Extrema_POnSurf(U2,-V,Ps);
 281
 282   myNbExt = 2;
 283   myDone = Standard_True;
 284 }
 285 //=============================================================================
 286
 287 Extrema_ExtPElS::Extrema_ExtPElS (const gp_Pnt&       P,
 288                                   const gp_Torus&     S,
 289                                   const Standard_Real Tol)
 290 {
 291   Perform(P, S, Tol);
 292 }
 293 /*-----------------------------------------------------------------------------
 294 Function:
 295   Find 2 extreme distances between point P and torus S.
 296
 297   Method:
 298   Let Pp be the projection of point P in plane XOY of the torus;
 299   2 cases are consideres:
 300   1- distance(Pp,O) < Tol:
 301      There is an infinite number of solutions; IsDone() = Standard_False.
 302   2- distance(Pp,O) > Tol:
 303      One is located in plane PpOZ;
 304      Let V1 = angle(OX,OPp) with 0. < V1 < 2.*M_PI,
 305          V2 = V1 + M_PI with 0. < V2 < 2.*M_PI,
 306          O1 and O2 centers of circles (O1 on coord. posit.)
 307          U1 = angle(OPp,O1P),
 308          U2 = angle(OPp,PO2);
 309      then (U1,V1) corresponds to the min distance
 310      and  (U2,V2) corresponds to the max distance.
 311 -----------------------------------------------------------------------------*/
 312 void Extrema_ExtPElS::Perform(const gp_Pnt&       P,
 313                               const gp_Torus&     S,
 314                               const Standard_Real Tol)
 315 {
 316   const Standard_Real tol2 = Tol*Tol;
 317   myDone = Standard_False;
 318   myNbExt = 0;
 319
 320 // Projection of P in plane XOY ...
 321   gp_Ax3 Pos = S.Position();
 322   gp_Pnt O = Pos.Location();
 323   gp_Vec OZ (Pos.Direction());
 324   gp_Pnt Pp = P.Translated(OZ.Multiplied(-(gp_Vec(O,P).Dot(Pos.Direction()))));
 325
 326 // Calculation of extrema ...
 327   gp_Vec OPp (O,Pp);
 328   Standard_Real R2 = OPp.SquareMagnitude();
 329   if (R2 < tol2) { return; }
 330
 331   gp_Vec myZ = Pos.XDirection()^Pos.YDirection();
 332   Standard_Real U1 = gp_Vec(Pos.XDirection()).AngleWithRef(OPp,myZ);
 333   if (U1 > -ExtPElS_MyEps && U1 < ExtPElS_MyEps) { U1 = 0.; }
 334   Standard_Real U2 = U1 + M_PI;
 335   if (U1 < 0.) { U1 += 2. * M_PI; }
 336   Standard_Real R = sqrt(R2);
 337   gp_Vec OO1 = OPp.Divided(R).Multiplied(S.MajorRadius());
 338   gp_Vec OO2 = OO1.Multiplied(-1.);
 339   gp_Pnt O1 = O.Translated(OO1);
 340   gp_Pnt O2 = O.Translated(OO2);
 341
 342   if(O1.SquareDistance(P) < tol2) { return; }
 343   if(O2.SquareDistance(P) < tol2) { return; }
 344
 345   Standard_Real V1 = OPp.AngleWithRef(gp_Vec(O1,P),OPp.Crossed(OZ));
 346   if (V1 > -ExtPElS_MyEps && V1 < ExtPElS_MyEps) { V1 = 0.; }
 347   OPp.Reverse();
 348   Standard_Real V2 = OPp.AngleWithRef(gp_Vec(P,O2),OPp.Crossed(OZ));
 349   if (V2 > -ExtPElS_MyEps && V2 < ExtPElS_MyEps) { V2 = 0.; }
 350
 351   if (V1 < 0.) { V1 += 2. * M_PI; }
 352   if (V2 < 0.) { V2 += 2. * M_PI; }
 353
 354   gp_Pnt Ps;
 355   Ps = ElSLib::Value(U1,V1,S);
 356   mySqDist[0] = Ps.SquareDistance(P);
 357   myPoint[0] = Extrema_POnSurf(U1,V1,Ps);
 358
 359   Ps = ElSLib::Value(U1,V1+M_PI,S);
 360   mySqDist[1] = Ps.SquareDistance(P);
 361   myPoint[1] = Extrema_POnSurf(U1,V1+M_PI,Ps);
 362
 363   Ps = ElSLib::Value(U2,V2,S);
 364   mySqDist[2] = Ps.SquareDistance(P);
 365   myPoint[2] = Extrema_POnSurf(U2,V2,Ps);
 366
 367   Ps = ElSLib::Value(U2,V2+M_PI,S);
 368   mySqDist[3] = Ps.SquareDistance(P);
 369   myPoint[3] = Extrema_POnSurf(U2,V2+M_PI,Ps);
 370
 371   myNbExt = 4;
 372   myDone = Standard_True;
 373 }
 374
 375
 376 Extrema_ExtPElS::Extrema_ExtPElS (const gp_Pnt&       P,
 377                                   const gp_Pln&       S,
 378                                   const Standard_Real Tol)
 379 {
 380   Perform(P, S, Tol);
 381 }
 382
 383 void Extrema_ExtPElS::Perform (const gp_Pnt&       P,
 384                                const gp_Pln&       S,
 385 //                             const Standard_Real Tol)
 386                                const Standard_Real )
 387 {
 388   myDone = Standard_False;
 389   myNbExt = 0;
 390
 391 // Projection of point P in plane XOY of the cylinder ...
 392   gp_Pnt O = S.Location();
 393   gp_Vec OZ (S.Axis().Direction());
 394   Standard_Real U, V = gp_Vec(O,P).Dot(OZ);
 395   gp_Pnt Pp = P.Translated(OZ.Multiplied(-V));
 396
 397   ElSLib::Parameters(S, P, U, V);
 398   mySqDist[0] = Pp.SquareDistance(P);
 399   myPoint[0] = Extrema_POnSurf(U,V,Pp);
 400   myNbExt = 1;
 401   myDone = Standard_True;
 402 }
 403
 404
 405 //=============================================================================
 406
 407 Standard_Boolean Extrema_ExtPElS::IsDone () const { return myDone; }
 408 //=============================================================================
 409
 410 Standard_Integer Extrema_ExtPElS::NbExt () const
 411 {
 412   if (!IsDone()) { throw StdFail_NotDone(); }
 413   return myNbExt;
 414 }
 415 //=============================================================================
 416
 417 Standard_Real Extrema_ExtPElS::SquareDistance (const Standard_Integer N) const
 418 {
 419   if (!IsDone()) { throw StdFail_NotDone(); }
 420   if ((N < 1) || (N > myNbExt)) { throw Standard_OutOfRange(); }
 421   return mySqDist[N-1];
 422 }
 423 //=============================================================================
 424
 425 const Extrema_POnSurf& Extrema_ExtPElS::Point (const Standard_Integer N) const
 426 {
 427   if (!IsDone()) { throw StdFail_NotDone(); }
 428   if ((N < 1) || (N > myNbExt)) { throw Standard_OutOfRange(); }
 429   return myPoint[N-1];
 430 }
 431 //=============================================================================