16d8313e88687e06c6dd87c74bca7aaf0dc9844a
[occt.git] / 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 const 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 //=============================================================================