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