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