Commit | Line | Data |
---|---|---|
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 | ||
13 | Extrema_ExtPElC::Extrema_ExtPElC () { myDone = Standard_False; } | |
14 | //============================================================================= | |
15 | ||
16 | Extrema_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 | ||
25 | void 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 | ||
53 | Extrema_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 | ||
62 | void 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 | /*----------------------------------------------------------------------------- | |
68 | Fonction: | |
69 | Recherche des valeurs de parametre u telle que: | |
70 | - dist(P,C(u)) passe par un extremum, | |
71 | - Uinf <= u <= Usup. | |
72 | ||
73 | Methode: | |
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 | ||
142 | Extrema_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 | ||
153 | void 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 | /*----------------------------------------------------------------------------- | |
159 | Fonction: | |
160 | Recherche des valeurs de parametre u telle que: | |
161 | - dist(P,C(u)) passe par un extremum, | |
162 | - Uinf <= u <= Usup. | |
163 | ||
164 | Methode: | |
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 | ||
222 | Extrema_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 | ||
232 | void 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 | /*----------------------------------------------------------------------------- | |
238 | Fonction: | |
239 | Recherche des valeurs de parametre u telle que: | |
240 | - dist(P,C(u)) passe par un extremum, | |
241 | - Uinf <= u <= Usup. | |
242 | ||
243 | Methode: | |
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 | ||
327 | Extrema_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 | ||
337 | void 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 | /*----------------------------------------------------------------------------- | |
344 | Fonction: | |
345 | Recherche des valeurs de parametre u telle que: | |
346 | - dist(P,C(u)) passe par un extremum, | |
347 | - Uinf <= u <= Usup. | |
348 | ||
349 | Methode: | |
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 | ||
417 | Standard_Boolean Extrema_ExtPElC::IsDone () const { return myDone; } | |
418 | //============================================================================= | |
419 | ||
420 | Standard_Integer Extrema_ExtPElC::NbExt () const | |
421 | { | |
422 | if (!IsDone()) { StdFail_NotDone::Raise(); } | |
423 | return myNbExt; | |
424 | } | |
425 | //============================================================================= | |
426 | ||
427 | Standard_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 | ||
434 | Standard_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 | ||
441 | Extrema_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 | //============================================================================= |