projects / occt.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
0022692: A wrong function is called inside the Extrema_CurveTool::IsRational
[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/*-----------------------------------------------------------------------------
68Fonction:
69 Recherche des valeurs de parametre u telle que:
70 - dist(P,C(u)) passe par un extremum,
71 - Uinf <= u <= Usup.
72
73Methode:
74 On procede en 3 temps:
75 1- Projection du point P dans le plan du cercle,
76 2- Calculs des u solutions dans [0.,2.*PI]:
77 Soit Pp, le point projete et
78 O, le centre du cercle;
79 2 cas:
80 - si Pp est confondu avec O, il y a une infinite de solutions;
81 IsDone() renvoie Standard_False.
82 - sinon, 2 points sont solutions pour le cercle complet:
83 . Us1 = angle(OPp,OX) correspond au minimum,
84 . soit Us2 = ( Us1 + PI si Us1 < PI,
85 ( Us1 - PI sinon;
86 Us2 correspond au maximum.
87 3- Calcul des extrema dans [Uinf,Usup].
88-----------------------------------------------------------------------------*/
89{
90 myDone = Standard_False;
91 myNbExt = 0;
92
93// 1- Projection du point P dans le plan du cercle -> Pp ...
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
100// 2- Calcul des u solutions dans [0.,2.*PI] ...
101
102 gp_Vec OPp (O,Pp);
103 if (OPp.Magnitude() < Tol) { return; }
104 Standard_Real Usol[2];
105 Usol[0] = C.XAxis().Direction().AngleWithRef(OPp,Axe); // -PI<U1<PI
106 Usol[1] = Usol[0] + PI;
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
416cd0b7
S		 118 ElCLib::AdjustPeriodic(Uinf, Uinf+2*PI, TolU, myuinf, Usol[0]);
119 ElCLib::AdjustPeriodic(Uinf, Uinf+2*PI, TolU, myuinf, Usol[1]);
7fd59977 120 if (((Usol[0]-2*PI-Uinf) < TolU) && ((Usol[0]-2*PI-Uinf) > -TolU)) Usol[0] = Uinf;
121 if (((Usol[1]-2*PI-Uinf) < TolU) && ((Usol[1]-2*PI-Uinf) > -TolU)) Usol[1] = Uinf;
122
123
124// 3- Calcul des extrema dans [Umin,Umax] ...
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/*-----------------------------------------------------------------------------
159Fonction:
160 Recherche des valeurs de parametre u telle que:
161 - dist(P,C(u)) passe par un extremum,
162 - Uinf <= u <= Usup.
163
164Methode:
165 On procede en 2 temps:
166 1- Projection du point P dans le plan de l'ellipse,
167 2- Calculs des solutions:
168 Soit Pp, le point projete; on recherche les valeurs u telles que:
169 (C(u)-Pp).C'(u) = 0. (1)
170 Soit Cos = cos(u) et Sin = sin(u),
171 C(u) = (A*Cos,B*Sin) et Pp = (X,Y);
172 Alors, (1) <=> (A*Cos-X,B*Sin-Y).(-A*Sin,B*Cos) = 0.
173 (B**2-A**2)*Cos*Sin - B*Y*Cos + A*X*Sin = 0.
174 On utilise l'algorithme math_TrigonometricFunctionRoots pour resoudre
175 cette equation.
176-----------------------------------------------------------------------------*/
177{
178 myDone = Standard_False;
179 myNbExt = 0;
180
181// 1- Projection du point P dans le plan de l ellipse -> Pp ...
182
183 gp_Pnt O = C.Location();
184 gp_Vec Axe (C.Axis().Direction());
185 gp_Vec Trsl = Axe.Multiplied(-(gp_Vec(O,P).Dot(Axe)));
186 gp_Pnt Pp = P.Translated(Trsl);
187
188// 2- Calculs des solutions ...
189
190 Standard_Integer NoSol, NbSol;
191 Standard_Real A = C.MajorRadius();
192 Standard_Real B = C.MinorRadius();
193 gp_Vec OPp (O,Pp);
194 Standard_Real OPpMagn = OPp.Magnitude();
195 if (OPpMagn < Tol) { if (Abs(A-B) < Tol) { return; } }
196 Standard_Real X = OPp.Dot(gp_Vec(C.XAxis().Direction()));
197 Standard_Real Y = OPp.Dot(gp_Vec(C.YAxis().Direction()));
198 // Standard_Real Y = Sqrt(OPpMagn*OPpMagn-X*X);
199
200 Standard_Real ko2 = (B*B-A*A)/2., ko3 = -B*Y, ko4 = A*X;
201 if(Abs(ko3) < 1.e-16*Max(Abs(ko2), Abs(ko3))) ko3 = 0.0;
202
203// math_TrigonometricFunctionRoots Sol(0.,(B*B-A*A)/2.,-B*Y,A*X,0.,Uinf,Usup);
204 math_TrigonometricFunctionRoots Sol(0.,ko2, ko3, ko4, 0.,Uinf,Usup);
205
206 if (!Sol.IsDone()) { return; }
207 gp_Pnt Cu;
208 Standard_Real Us;
209 NbSol = Sol.NbSolutions();
210 for (NoSol = 1; NoSol <= NbSol; NoSol++) {
211 Us = Sol.Value(NoSol);
212 Cu = ElCLib::Value(Us,C);
213 mySqDist[myNbExt] = Cu.SquareDistance(P);
214 myIsMin[myNbExt] = (NoSol == 1);
215 myPoint[myNbExt] = Extrema_POnCurv(Us,Cu);
216 myNbExt++;
217 }
218 myDone = Standard_True;
219}
220//=============================================================================
221
222Extrema_ExtPElC::Extrema_ExtPElC (const gp_Pnt& P,
223 const gp_Hypr& C,
224 const Standard_Real Tol,
225 const Standard_Real Uinf,
226 const Standard_Real Usup)
227{
228 Perform(P, C, Tol, Uinf, Usup);
229}
230
231
232void Extrema_ExtPElC::Perform(const gp_Pnt& P,
233 const gp_Hypr& C,
234 const Standard_Real Tol,
235 const Standard_Real Uinf,
236 const Standard_Real Usup)
237/*-----------------------------------------------------------------------------
238Fonction:
239 Recherche des valeurs de parametre u telle que:
240 - dist(P,C(u)) passe par un extremum,
241 - Uinf <= u <= Usup.
242
243Methode:
244 On procede en 2 temps:
245 1- Projection du point P dans le plan de l'hyperbole,
246 2- Calculs des solutions:
247 Soit Pp, le point projete; on recherche les valeurs u telles que:
248 (C(u)-Pp).C'(u) = 0. (1)
249 Soit R et r les rayons de l'hyperbole,
250 Chu = Cosh(u) et Shu = Sinh(u),
251 C(u) = (R*Chu,r*Shu) et Pp = (X,Y);
252 Alors, (1) <=> (R*Chu-X,r*Shu-Y).(R*Shu,r*Chu) = 0.
253 (R**2+r**2)*Chu*Shu - X*R*Shu - Y*r*Chu = 0. (2)
254 Soit v = e**u;
255 Alors, en utilisant Chu = (e**u+e**(-u))/2. et Sh = (e**u-e**(-u)))/2.
256 (2) <=> ((R**2+r**2)/4.) * (v**2-v**(-2)) -
257 ((X*R+Y*r)/2.) * v +
258 ((X*R-Y*r)/2.) * v**(-1) = 0.
259 (2)* v**2 <=> ((R**2+r**2)/4.) * v**4 -
260 ((X*R+Y*r)/2.) * v**3 +
261 ((X*R-Y*r)/2.) * v -
262 ((R**2+r**2)/4.) = 0.
263 On utilise l'algorithme math_DirectPolynomialRoots pour resoudre
264 cette equation en v.
265-----------------------------------------------------------------------------*/
266{
267 myDone = Standard_False;
268 myNbExt = 0;
269
270// 1- Projection du point P dans le plan de l hyperbole -> Pp ...
271
272 gp_Pnt O = C.Location();
273 gp_Vec Axe (C.Axis().Direction());
274 gp_Vec Trsl = Axe.Multiplied(-(gp_Vec(O,P).Dot(Axe)));
275 gp_Pnt Pp = P.Translated(Trsl);
276
277// 2- Calculs des solutions ...
278
279 Standard_Real Tol2 = Tol * Tol;
280 Standard_Real R = C.MajorRadius();
281 Standard_Real r = C.MinorRadius();
282 gp_Vec OPp (O,Pp);
283#ifdef DEB
284 Standard_Real OPpMagn = OPp.Magnitude();
285#else
286 OPp.Magnitude();
287#endif
288 Standard_Real X = OPp.Dot(gp_Vec(C.XAxis().Direction()));
289 Standard_Real Y = OPp.Dot(gp_Vec(C.YAxis().Direction()));
290
291 Standard_Real C1 = (R*R+r*r)/4.;
292 math_DirectPolynomialRoots Sol(C1,-(X*R+Y*r)/2.,0.,(X*R-Y*r)/2.,-C1);
293 if (!Sol.IsDone()) { return; }
294 gp_Pnt Cu;
295 Standard_Real Us, Vs;
296 Standard_Integer NbSol = Sol.NbSolutions();
297 Standard_Boolean DejaEnr;
298 Standard_Integer NoExt;
299 gp_Pnt TbExt[4];
300 for (Standard_Integer NoSol = 1; NoSol <= NbSol; NoSol++) {
301 Vs = Sol.Value(NoSol);
302 if (Vs > 0.) {
303 Us = Log(Vs);
304 if ((Us >= Uinf) && (Us <= Usup)) {
305 Cu = ElCLib::Value(Us,C);
306 DejaEnr = Standard_False;
307 for (NoExt = 0; NoExt < myNbExt; NoExt++) {
308 if (TbExt[NoExt].SquareDistance(Cu) < Tol2) {
309 DejaEnr = Standard_True;
310 break;
311 }
312 }
313 if (!DejaEnr) {
314 TbExt[myNbExt] = Cu;
315 mySqDist[myNbExt] = Cu.SquareDistance(P);
316 myIsMin[myNbExt] = mySqDist[myNbExt] < P.SquareDistance(ElCLib::Value(Us+1,C));
317 myPoint[myNbExt] = Extrema_POnCurv(Us,Cu);
318 myNbExt++;
319 }
320 } // if ((Us >= Uinf) && (Us <= Usup))
321 } // if (Vs > 0.)
322 } // for (Standard_Integer NoSol = 1; ...
323 myDone = Standard_True;
324}
325//=============================================================================
326
327Extrema_ExtPElC::Extrema_ExtPElC (const gp_Pnt& P,
328 const gp_Parab& C,
329 const Standard_Real Tol,
330 const Standard_Real Uinf,
331 const Standard_Real Usup)
332{
333 Perform(P, C, Tol, Uinf, Usup);
334}
335
336
337void Extrema_ExtPElC::Perform(const gp_Pnt& P,
338 const gp_Parab& C,
339// const Standard_Real Tol,
340 const Standard_Real ,
341 const Standard_Real Uinf,
342 const Standard_Real Usup)
343/*-----------------------------------------------------------------------------
344Fonction:
345 Recherche des valeurs de parametre u telle que:
346 - dist(P,C(u)) passe par un extremum,
347 - Uinf <= u <= Usup.
348
349Methode:
350 On procede en 2 temps:
351 1- Projection du point P dans le plan de la parabole,
352 2- Calculs des solutions:
353 Soit Pp, le point projete; on recherche les valeurs u telles que:
354 (C(u)-Pp).C'(u) = 0. (1)
355 Soit F la focale de la parabole,
356 C(u) = ((u*u)/(4.*F),u) et Pp = (X,Y);
357 Alors, (1) <=> ((u*u)/(4.*F)-X,u-Y).(u/(2.*F),1) = 0.
358 (1./(4.*F)) * U**3 + (2.*F-X) * U - 2*F*Y = 0.
359 On utilise l'algorithme math_DirectPolynomialRoots pour resoudre
360 cette equation en U.
361-----------------------------------------------------------------------------*/
362{
363 myDone = Standard_False;
364 myNbExt = 0;
365
366// 1- Projection du point P dans le plan de la parabole -> Pp ...
367
368 gp_Pnt O = C.Location();
369 gp_Vec Axe (C.Axis().Direction());
370 gp_Vec Trsl = Axe.Multiplied(-(gp_Vec(O,P).Dot(Axe)));
371 gp_Pnt Pp = P.Translated(Trsl);
372
373// 2- Calculs des solutions ...
374
375 Standard_Real F = C.Focal();
376 gp_Vec OPp (O,Pp);
377#ifdef DEB
378 Standard_Real OPpMagn = OPp.Magnitude();
379#else
380 OPp.Magnitude();
381#endif
382 Standard_Real X = OPp.Dot(gp_Vec(C.XAxis().Direction()));
383// Standard_Real Y = Sqrt(OPpMagn*OPpMagn-X*X);
384 Standard_Real Y = OPp.Dot(gp_Vec(C.YAxis().Direction()));
385 math_DirectPolynomialRoots Sol(1./(4.*F),0.,2.*F-X,-2.*F*Y);
386 if (!Sol.IsDone()) { return; }
387 gp_Pnt Cu;
388 Standard_Real Us;
389 Standard_Integer NbSol = Sol.NbSolutions();
390 Standard_Boolean DejaEnr;
391 Standard_Integer NoExt;
392 gp_Pnt TbExt[3];
393 for (Standard_Integer NoSol = 1; NoSol <= NbSol; NoSol++) {
394 Us = Sol.Value(NoSol);
395 if ((Us >= Uinf) && (Us <= Usup)) {
396 Cu = ElCLib::Value(Us,C);
397 DejaEnr = Standard_False;
398 for (NoExt = 0; NoExt < myNbExt; NoExt++) {
399 if (TbExt[NoExt].SquareDistance(Cu) < Precision::Confusion() * Precision::Confusion()) {
400 DejaEnr = Standard_True;
401 break;
402 }
403 }
404 if (!DejaEnr) {
405 TbExt[myNbExt] = Cu;
406 mySqDist[myNbExt] = Cu.SquareDistance(P);
407 myIsMin[myNbExt] = mySqDist[myNbExt] < P.SquareDistance(ElCLib::Value(Us+1,C));
408 myPoint[myNbExt] = Extrema_POnCurv(Us,Cu);
409 myNbExt++;
410 }
411 } // if ((Us >= Uinf) && (Us <= Usup))
412 } // for (Standard_Integer NoSol = 1; ...
413 myDone = Standard_True;
414}
415//=============================================================================
416
417Standard_Boolean Extrema_ExtPElC::IsDone () const { return myDone; }
418//=============================================================================
419
420Standard_Integer Extrema_ExtPElC::NbExt () const
421{
422 if (!IsDone()) { StdFail_NotDone::Raise(); }
423 return myNbExt;
424}
425//=============================================================================
426
427Standard_Real Extrema_ExtPElC::SquareDistance (const Standard_Integer N) const
428{
429 if ((N < 1) || (N > NbExt())) { Standard_OutOfRange::Raise(); }
430 return mySqDist[N-1];
431}
432//=============================================================================
433
434Standard_Boolean Extrema_ExtPElC::IsMin (const Standard_Integer N) const
435{
436 if ((N < 1) || (N > NbExt())) { Standard_OutOfRange::Raise(); }
437 return myIsMin[N-1];
438}
439//=============================================================================
440
441Extrema_POnCurv Extrema_ExtPElC::Point (const Standard_Integer N) const
442{
443 if ((N < 1) || (N > NbExt())) { Standard_OutOfRange::Raise(); }
444 return myPoint[N-1];
445}
446//=============================================================================
Open CASCADE Technology repository