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