projects / occt.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
0022241: The bug is appendix to the Salome Bug 0021148
[occt.git] / src / Extrema / Extrema_ExtPRevS.cxx
CommitLineData
7fd59977 1// File: Extrema_ExtPRevS.cxx
2// Created: Tue Sep 21 15:50:01 1999
3// Author: Edward AGAPOV
4// <eap@strelox.nnov.matra-dtv.fr>
5
6
7#include <Extrema_ExtPRevS.ixx>
8#include <Adaptor3d_HCurve.hxx>
9#include <Extrema_ExtPElC.hxx>
10#include <Extrema_GenExtPS.hxx>
11#include <Extrema_POnCurv.hxx>
12#include <Extrema_POnSurf.hxx>
13#include <gp_Ax1.hxx>
14#include <gp_Ax2.hxx>
15#include <gp_Lin.hxx>
16#include <gp_Pln.hxx>
17#include <gp_Pnt.hxx>
18#include <gp_Trsf.hxx>
19#include <gp_Vec.hxx>
20#include <Precision.hxx>
21#include <ElCLib.hxx>
22
23static gp_Ax2 GetPosition (const Adaptor3d_SurfaceOfRevolution& S)//const Handle(Adaptor_HCurve)& C)
24{
25 Handle(Adaptor3d_HCurve) C = S.BasisCurve();
26
27 switch (C->GetType()) {
28
29 case GeomAbs_Line: {
30 gp_Lin L = C->Line();
31 gp_Dir N = S.AxeOfRevolution().Direction();
32 if (N.IsParallel(L.Direction(), Precision::Angular())) {
33 gp_Vec OO (L.Location(), S.AxeOfRevolution().Location());
34 if (OO.Magnitude() <= gp::Resolution()) {
35 OO = gp_Vec(L.Location(), ElCLib::Value(100,L));
36 if (N.IsParallel(OO, Precision::Angular()))
37 return gp_Ax2(); // Line and axe of revolution coinside
38 }
39 N ^= OO;
40 }
41 else {
42 N ^= L.Direction();
43 }
44 return gp_Ax2 (L.Location(), N, L.Direction());
45 }
46 case GeomAbs_Circle:
47 return C->Circle().Position();
48 case GeomAbs_Ellipse:
49 return C->Ellipse().Position();
50 case GeomAbs_Hyperbola:
51 return C->Hyperbola().Position();
52 case GeomAbs_Parabola:
53 return C->Parabola().Position();
54 default:
55 return gp_Ax2();
56 }
57}
58
59//=======================================================================
60//function : HasSingularity
61//purpose :
62//=======================================================================
63
64static Standard_Boolean HasSingularity(const Adaptor3d_SurfaceOfRevolution& S)
65{
66
67 const Handle(Adaptor3d_HCurve) C = S.BasisCurve();
68 gp_Dir N = S.AxeOfRevolution().Direction();
69 gp_Pnt P = S.AxeOfRevolution().Location();
70
71 gp_Lin L(P, N);
72
73 P = C->Value(C->FirstParameter());
74
75 if(L.SquareDistance(P) < Precision::Confusion() * Precision::Confusion()) return Standard_True;
76
77 P = C->Value(C->LastParameter());
78
79 if(L.SquareDistance(P) < Precision::Confusion() * Precision::Confusion()) return Standard_True;
80
81 return Standard_False;
82}
83
84//=============================================================================
85
86static void PerformExtPElC (Extrema_ExtPElC& E,
87 const gp_Pnt& P,
88 const Handle(Adaptor3d_HCurve)& C,
89 const Standard_Real Tol)
90{
91 switch (C->GetType()) {
92 case GeomAbs_Hyperbola:
93 E.Perform(P, C->Hyperbola(), Tol, -Precision::Infinite(),Precision::Infinite());
94 return;
95 case GeomAbs_Line:
96 E.Perform(P, C->Line(), Tol, -Precision::Infinite(),Precision::Infinite());
97 return;
98 case GeomAbs_Circle:
99 E.Perform(P, C->Circle(), Tol, 0.0, 2.0 * PI);
100 return;
101 case GeomAbs_Ellipse:
102 E.Perform(P, C->Ellipse(), Tol, 0.0, 2.0 * PI);
103 return;
104 case GeomAbs_Parabola:
105 E.Perform(P, C->Parabola(), Tol, -Precision::Infinite(),Precision::Infinite());
106 return;
107#ifndef DEB
108 default:
109 return ;
110#endif
111 }
112}
113
114//=======================================================================
115//function : IsCaseAnalyticallyComputable
116//purpose :
117//=======================================================================
118
119static Standard_Boolean IsCaseAnalyticallyComputable
120 (const GeomAbs_CurveType& theType,
121 const gp_Ax2& theCurvePos,
122 const gp_Ax1& AxeOfRevolution)
123{
124 // check type
125 switch (theType) {
126 case GeomAbs_Line:
127 case GeomAbs_Circle:
128 case GeomAbs_Ellipse:
129 case GeomAbs_Hyperbola:
130 case GeomAbs_Parabola:
131 break;
132 default:
133 return Standard_False;
134 }
135// the axe of revolution must be in the plane of the curve.
136 gp_Pln pl(theCurvePos.Location(), theCurvePos.Direction());
137 gp_Pnt p1 = AxeOfRevolution.Location();
138 Standard_Real dist = 100., dist2 = dist * dist;
139 Standard_Real aThreshold = Precision::Angular() * Precision::Angular() * dist2;
140 gp_Pnt p2 = AxeOfRevolution.Location().XYZ() + dist * AxeOfRevolution.Direction().XYZ();
141
142 if((pl.SquareDistance(p1) < aThreshold) &&
143 (pl.SquareDistance(p2) < aThreshold))
144 return Standard_True;
145 return Standard_False;
146 // gp_Vec V (AxeOfRevolution.Location(),theCurvePos.Location());
147 // if (Abs( V * theCurvePos.Direction()) <= gp::Resolution())
148 // return Standard_True;
149 // else
150 // return Standard_False;
151}
152
153//=======================================================================
154//function : IsOriginalPnt
155//purpose :
156//=======================================================================
157
158static Standard_Boolean IsOriginalPnt (const gp_Pnt& P,
159 const Extrema_POnSurf* Points,
160 const Standard_Integer NbPoints)
161{
162 for (Standard_Integer i=1; i<=NbPoints; i++) {
163 if (Points[i].Value().IsEqual(P, Precision::Confusion())) {
164 return Standard_False;
165 }
166 }
167 return Standard_True;
168}
169//=======================================================================
170//function : IsExtremum
171//purpose :
172//=======================================================================
173
174static Standard_Boolean IsExtremum (const Standard_Real U, const Standard_Real V,
175 const gp_Pnt& P, const Adaptor3d_SurfacePtr& S,
176 gp_Pnt& E, Standard_Real& Dist2,
177 const Standard_Boolean IsVSup,
178 const Standard_Boolean IsMin)
179{
180 E = S->Value(U,V);
181 Dist2 = P.SquareDistance(E);
182 if (IsMin)
183 return (Dist2 < P.SquareDistance(S->Value(U+1,V)) &&
184 Dist2 < P.SquareDistance(S->Value(U-1,V)) &&
185 Dist2 < P.SquareDistance(S->Value(U, IsVSup ? V-1 : V+1)));
186 else
187 return (Dist2 > P.SquareDistance(S->Value(U+1,V)) &&
188 Dist2 > P.SquareDistance(S->Value(U-1,V)) &&
189 Dist2 > P.SquareDistance(S->Value(U, IsVSup ? V-1 : V+1)));
190}
191//=======================================================================
192//function : Extrema_ExtPRevS
193//purpose :
194//=======================================================================
195
196Extrema_ExtPRevS::Extrema_ExtPRevS()
197{
198 myS=NULL;
199 myDone = Standard_False;
200}
201//=======================================================================
202//function : Extrema_ExtPRevS
203//purpose :
204//=======================================================================
205
206Extrema_ExtPRevS::Extrema_ExtPRevS(const gp_Pnt& P,
207 const Adaptor3d_SurfaceOfRevolution& S,
208 const Standard_Real Umin,
209 const Standard_Real Usup,
210 const Standard_Real Vmin,
211 const Standard_Real Vsup,
212 const Standard_Real TolU,
213 const Standard_Real TolV)
214{
215 myS=NULL;
216 Initialize (S,
217 Umin, Usup, Vmin, Vsup,
218 TolU, TolV);
219 Perform(P);
220}
221//=======================================================================
222//function : Extrema_ExtPRevS
223//purpose :
224//=======================================================================
225
226Extrema_ExtPRevS::Extrema_ExtPRevS(const gp_Pnt& P,
227 const Adaptor3d_SurfaceOfRevolution& S,
228 const Standard_Real TolU,
229 const Standard_Real TolV)
230{
231 myS=NULL;
232 Initialize (S,
233 S.FirstUParameter(), S.LastUParameter(),
234 S.FirstVParameter(), S.LastVParameter(),
235 TolU, TolV);
236 Perform(P);
237}
238//=======================================================================
239//function : Initialize
240//purpose :
241//=======================================================================
242
243void Extrema_ExtPRevS::Initialize(const Adaptor3d_SurfaceOfRevolution& S,
244 const Standard_Real Umin,
245 const Standard_Real Usup,
246 const Standard_Real Vmin,
247 const Standard_Real Vsup,
248 const Standard_Real TolU,
249 const Standard_Real TolV)
250{
251 myvinf=Vmin;
252 myvsup=Vsup;
253 mytolv=TolV;
254
255 Handle(Adaptor3d_HCurve) anACurve = S.BasisCurve();
256 if (!myS || myS != (Adaptor3d_SurfacePtr)&S) {
257 myS = Adaptor3d_SurfacePtr(&S);
258 myPosition = GetPosition(S);
259 myIsAnalyticallyComputable =
260 IsCaseAnalyticallyComputable (anACurve->GetType(),myPosition,S.AxeOfRevolution());
261 }
262 if (!myIsAnalyticallyComputable) {
263
264 Standard_Integer nbu = 32, nbv = 32;
265
266 if(HasSingularity(S)) {nbv = 100;}
267
268 myExtPS.Initialize(S, nbu, nbv,
269 Umin, Usup, Vmin, Vsup,
270 TolU, TolV);
271 }
272}
273//=======================================================================
274//function : Perform
275//purpose :
276//=======================================================================
277
278void Extrema_ExtPRevS::Perform(const gp_Pnt& P)
279{
280 myDone = Standard_False;
281 myNbExt = 0;
282
283 if (!myIsAnalyticallyComputable) {
284
285 myExtPS.Perform(P);
286 myDone = myExtPS.IsDone();
287 myNbExt = myExtPS.NbExt();
288 return;
289 }
290
291 Handle(Adaptor3d_HCurve) anACurve = myS->BasisCurve();
292
293 gp_Ax1 Ax = myS->AxeOfRevolution();
294 gp_Vec Dir = Ax.Direction(), Z = myPosition.Direction();
295 gp_Pnt O = Ax.Location();
296
297 Standard_Real OPdir = gp_Vec(O,P).Dot(Dir);
298 gp_Pnt Pp = P.Translated(Dir.Multiplied(-OPdir));
299 if (O.IsEqual(Pp,Precision::Confusion())) // P is on the AxeOfRevolution
300 return;
301
302 Standard_Real U,V;
303 gp_Pnt P1, Ppp;
304 Standard_Real OPpz = gp_Vec(O,Pp).Dot(Z);
305 if (Abs(OPpz) <= gp::Resolution()) {
306 Ppp = Pp;
307 U = 0;
308 }
309 else {
310 Ppp = Pp.Translated(Z.Multiplied(-OPpz));
311 if (O.IsEqual(Ppp,Precision::Confusion()))
312 U = PI/2;
313 else {
314 U = gp_Vec(O,Ppp).AngleWithRef(gp_Vec(O,Pp),Dir);
315 }
316 }
317
318 gp_Vec OPpp (O,Ppp), OPq (O, myS->Value(PI/2,0));
319 if (U != PI/2) {
320 if (Abs(OPq.Magnitude()) <= gp::Resolution())
321 OPq = gp_Vec(O, myS->Value(PI/2,anACurve->LastParameter()/10));
322 if (OPpp.AngleWithRef(OPq,Dir) < 0)
323 U += PI;
324 }
325
326 gp_Trsf T;
327 T.SetRotation(Ax, -U);
328 P1 = P.Transformed(T);
329
330 gp_Pnt E;
331 Standard_Real Dist2;
332 Standard_Integer i;
333
334 Extrema_ExtPElC anExt;
335 PerformExtPElC(anExt, P1, anACurve, mytolv);
336
337 if (anExt.IsDone()) {
338 myDone = Standard_True;
339 for (i=1; i<=anExt.NbExt(); i++) {
340 Extrema_POnCurv POC=anExt.Point(i);
341 V = POC.Parameter();
342 if (V > myvsup) {
343 // if ( !IsExtremum (U, V = myvsup, P, myS, E, Dist2,
344 // Standard_True, anExt.IsMin(i))) continue;
345 Standard_Real newV = myvsup;
346
347 if((anACurve->GetType() == GeomAbs_Circle) ||
348 (anACurve->GetType() == GeomAbs_Ellipse)) {
349 newV = ElCLib::InPeriod(V, myvinf, myvinf + 2. * PI);
350
351 if (newV > myvsup) {
352 newV = myvsup;
353 }
354 }
355 V = newV;
356
357 if ( !IsExtremum (U, V, P, myS, E, Dist2,
358 Standard_True, anExt.IsMin(i))) {
359 continue;
360 }
361 } else if (V < myvinf) {
362 // if ( !IsExtremum (U, V = myvinf, P, myS, E, Dist2,
363 // Standard_False, anExt.IsMin(i))) continue;
364
365 Standard_Real newV = myvinf;
366
367 if((anACurve->GetType() == GeomAbs_Circle) ||
368 (anACurve->GetType() == GeomAbs_Ellipse)) {
369 newV = ElCLib::InPeriod(V, myvsup - 2. * PI, myvsup);
370
371 if(newV < myvinf)
372 newV = myvinf;
373 }
374 V = newV;
375
376 if ( !IsExtremum (U, V, P, myS, E, Dist2,
377 Standard_False, anExt.IsMin(i))) continue;
378 } else {
379 E = myS->Value(U,V);
380 Dist2 = P.SquareDistance(E);
381 }
382 if (IsOriginalPnt(E, myPoint, myNbExt)) {
383 myPoint[++myNbExt] = Extrema_POnSurf(U,V,E);
384 mySqDist[myNbExt] = Dist2;
385 }
386 }
387 }
388 T.SetRotation(Ax, PI);
389 P1.Transform(T);
390
391 PerformExtPElC(anExt, P1, anACurve, mytolv);
392 if (anExt.IsDone()) {
393 myDone = Standard_True;
394
395 U += PI;
396
397 for (i=1; i<=anExt.NbExt(); i++) {
398 Extrema_POnCurv POC=anExt.Point(i);
399 V = POC.Parameter();
400 if (V > myvsup) {
401 // if ( !IsExtremum (U, V = myvsup, P, myS, E, Dist2,
402 // Standard_True, anExt.IsMin(i))) continue;
403
404 Standard_Real newV = myvsup;
405
406 if((anACurve->GetType() == GeomAbs_Circle) ||
407 (anACurve->GetType() == GeomAbs_Ellipse)) {
408 newV = ElCLib::InPeriod(V, myvinf, myvinf + 2. * PI);
409
410 if (newV > myvsup) {
411 newV = myvsup;
412 }
413 }
414 V = newV;
415
416 if ( !IsExtremum (U, V, P, myS, E, Dist2,
417 Standard_True, anExt.IsMin(i))) continue;
418 } else if (V < myvinf) {
419 // if ( !IsExtremum (U, V = myvinf, P, myS, E, Dist2,
420 // Standard_False, anExt.IsMin(i))) continue;
421 Standard_Real newV = myvinf;
422
423 if((anACurve->GetType() == GeomAbs_Circle) ||
424 (anACurve->GetType() == GeomAbs_Ellipse)) {
425 newV = ElCLib::InPeriod(V, myvsup - 2. * PI, myvsup);
426
427 if(newV < myvinf)
428 newV = myvinf;
429 }
430 V = newV;
431
432 if ( !IsExtremum (U, V, P, myS, E, Dist2,
433 Standard_False, anExt.IsMin(i))) continue;
434 } else {
435 E = myS->Value(U,V);
436 Dist2 = P.SquareDistance(E);
437 }
438 if (IsOriginalPnt(E, myPoint, myNbExt)) {
439 myPoint[++myNbExt] = Extrema_POnSurf(U,V,E);
440 mySqDist[myNbExt] = Dist2;
441
442 }
443 }
444 }
445}
446
447
448//=======================================================================
449//function : IsDone
450//purpose :
451//=======================================================================
452
453Standard_Boolean Extrema_ExtPRevS::IsDone() const
454{
455 return myDone;
456}
457
458
459//=======================================================================
460//function : NbExt
461//purpose :
462//=======================================================================
463
464Standard_Integer Extrema_ExtPRevS::NbExt() const
465{
466 if (!IsDone()) { StdFail_NotDone::Raise(); }
467 return myNbExt;
468}
469
470
471//=======================================================================
472//function : Value
473//purpose :
474//=======================================================================
475
476Standard_Real Extrema_ExtPRevS::SquareDistance(const Standard_Integer N) const
477{
478 if (!IsDone()) { StdFail_NotDone::Raise(); }
479 if ((N < 1) || (N > myNbExt)) { Standard_OutOfRange::Raise(); }
480 if (myIsAnalyticallyComputable)
481 return mySqDist[N];
482 else
483 return myExtPS.SquareDistance(N);
484}
485//=======================================================================
486//function : Point
487//purpose :
488//=======================================================================
489
490Extrema_POnSurf Extrema_ExtPRevS::Point(const Standard_Integer N) const
491{
492 if (!IsDone()) { StdFail_NotDone::Raise(); }
493 if ((N < 1) || (N > myNbExt)) { Standard_OutOfRange::Raise(); }
494 if (myIsAnalyticallyComputable)
495 return myPoint[N];
496 else
497 return myExtPS.Point(N);
498}
499
500
Open CASCADE Technology repository