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