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>
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>
18 #include <gp_Trsf.hxx>
20 #include <Precision.hxx>
23 static gp_Ax2 GetPosition (const Adaptor3d_SurfaceOfRevolution& S)//const Handle(Adaptor_HCurve)& C)
25 Handle(Adaptor3d_HCurve) C = S.BasisCurve();
27 switch (C->GetType()) {
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
44 return gp_Ax2 (L.Location(), N, L.Direction());
47 return C->Circle().Position();
49 return C->Ellipse().Position();
50 case GeomAbs_Hyperbola:
51 return C->Hyperbola().Position();
52 case GeomAbs_Parabola:
53 return C->Parabola().Position();
59 //=======================================================================
60 //function : HasSingularity
62 //=======================================================================
64 static Standard_Boolean HasSingularity(const Adaptor3d_SurfaceOfRevolution& S)
67 const Handle(Adaptor3d_HCurve) C = S.BasisCurve();
68 gp_Dir N = S.AxeOfRevolution().Direction();
69 gp_Pnt P = S.AxeOfRevolution().Location();
73 P = C->Value(C->FirstParameter());
75 if(L.SquareDistance(P) < Precision::Confusion() * Precision::Confusion()) return Standard_True;
77 P = C->Value(C->LastParameter());
79 if(L.SquareDistance(P) < Precision::Confusion() * Precision::Confusion()) return Standard_True;
81 return Standard_False;
84 //=============================================================================
86 static void PerformExtPElC (Extrema_ExtPElC& E,
88 const Handle(Adaptor3d_HCurve)& C,
89 const Standard_Real Tol)
91 switch (C->GetType()) {
92 case GeomAbs_Hyperbola:
93 E.Perform(P, C->Hyperbola(), Tol, -Precision::Infinite(),Precision::Infinite());
96 E.Perform(P, C->Line(), Tol, -Precision::Infinite(),Precision::Infinite());
99 E.Perform(P, C->Circle(), Tol, 0.0, 2.0 * M_PI);
101 case GeomAbs_Ellipse:
102 E.Perform(P, C->Ellipse(), Tol, 0.0, 2.0 * M_PI);
104 case GeomAbs_Parabola:
105 E.Perform(P, C->Parabola(), Tol, -Precision::Infinite(),Precision::Infinite());
114 //=======================================================================
115 //function : IsCaseAnalyticallyComputable
117 //=======================================================================
119 static Standard_Boolean IsCaseAnalyticallyComputable
120 (const GeomAbs_CurveType& theType,
121 const gp_Ax2& theCurvePos,
122 const gp_Ax1& AxeOfRevolution)
128 case GeomAbs_Ellipse:
129 case GeomAbs_Hyperbola:
130 case GeomAbs_Parabola:
133 return Standard_False;
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();
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;
150 // return Standard_False;
153 //=======================================================================
154 //function : IsOriginalPnt
156 //=======================================================================
158 static Standard_Boolean IsOriginalPnt (const gp_Pnt& P,
159 const Extrema_POnSurf* Points,
160 const Standard_Integer NbPoints)
162 for (Standard_Integer i=1; i<=NbPoints; i++) {
163 if (Points[i].Value().IsEqual(P, Precision::Confusion())) {
164 return Standard_False;
167 return Standard_True;
169 //=======================================================================
170 //function : IsExtremum
172 //=======================================================================
174 static 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)
181 Dist2 = P.SquareDistance(E);
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)));
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)));
191 //=======================================================================
192 //function : Extrema_ExtPRevS
194 //=======================================================================
196 Extrema_ExtPRevS::Extrema_ExtPRevS()
199 myDone = Standard_False;
201 //=======================================================================
202 //function : Extrema_ExtPRevS
204 //=======================================================================
206 Extrema_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)
217 Umin, Usup, Vmin, Vsup,
221 //=======================================================================
222 //function : Extrema_ExtPRevS
224 //=======================================================================
226 Extrema_ExtPRevS::Extrema_ExtPRevS(const gp_Pnt& P,
227 const Adaptor3d_SurfaceOfRevolution& S,
228 const Standard_Real TolU,
229 const Standard_Real TolV)
233 S.FirstUParameter(), S.LastUParameter(),
234 S.FirstVParameter(), S.LastVParameter(),
238 //=======================================================================
239 //function : Initialize
241 //=======================================================================
243 void 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)
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());
262 if (!myIsAnalyticallyComputable) {
264 Standard_Integer nbu = 32, nbv = 32;
266 if(HasSingularity(S)) {nbv = 100;}
268 myExtPS.Initialize(S, nbu, nbv,
269 Umin, Usup, Vmin, Vsup,
273 //=======================================================================
276 //=======================================================================
278 void Extrema_ExtPRevS::Perform(const gp_Pnt& P)
280 myDone = Standard_False;
283 if (!myIsAnalyticallyComputable) {
286 myDone = myExtPS.IsDone();
287 myNbExt = myExtPS.NbExt();
291 Handle(Adaptor3d_HCurve) anACurve = myS->BasisCurve();
293 gp_Ax1 Ax = myS->AxeOfRevolution();
294 gp_Vec Dir = Ax.Direction(), Z = myPosition.Direction();
295 gp_Pnt O = Ax.Location();
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
304 Standard_Real OPpz = gp_Vec(O,Pp).Dot(Z);
305 if (Abs(OPpz) <= gp::Resolution()) {
310 Ppp = Pp.Translated(Z.Multiplied(-OPpz));
311 if (O.IsEqual(Ppp,Precision::Confusion()))
314 U = gp_Vec(O,Ppp).AngleWithRef(gp_Vec(O,Pp),Dir);
318 gp_Vec OPpp (O,Ppp), OPq (O, myS->Value(M_PI/2,0));
320 if (Abs(OPq.Magnitude()) <= gp::Resolution())
321 OPq = gp_Vec(O, myS->Value(M_PI/2,anACurve->LastParameter()/10));
322 if (OPpp.AngleWithRef(OPq,Dir) < 0)
327 T.SetRotation(Ax, -U);
328 P1 = P.Transformed(T);
334 Extrema_ExtPElC anExt;
335 PerformExtPElC(anExt, P1, anACurve, mytolv);
337 if (anExt.IsDone()) {
338 myDone = Standard_True;
339 for (i=1; i<=anExt.NbExt(); i++) {
340 Extrema_POnCurv POC=anExt.Point(i);
343 // if ( !IsExtremum (U, V = myvsup, P, myS, E, Dist2,
344 // Standard_True, anExt.IsMin(i))) continue;
345 Standard_Real newV = myvsup;
347 if((anACurve->GetType() == GeomAbs_Circle) ||
348 (anACurve->GetType() == GeomAbs_Ellipse)) {
349 newV = ElCLib::InPeriod(V, myvinf, myvinf + 2. * M_PI);
357 if ( !IsExtremum (U, V, P, myS, E, Dist2,
358 Standard_True, anExt.IsMin(i))) {
361 } else if (V < myvinf) {
362 // if ( !IsExtremum (U, V = myvinf, P, myS, E, Dist2,
363 // Standard_False, anExt.IsMin(i))) continue;
365 Standard_Real newV = myvinf;
367 if((anACurve->GetType() == GeomAbs_Circle) ||
368 (anACurve->GetType() == GeomAbs_Ellipse)) {
369 newV = ElCLib::InPeriod(V, myvsup - 2. * M_PI, myvsup);
376 if ( !IsExtremum (U, V, P, myS, E, Dist2,
377 Standard_False, anExt.IsMin(i))) continue;
380 Dist2 = P.SquareDistance(E);
382 if (IsOriginalPnt(E, myPoint, myNbExt)) {
383 myPoint[++myNbExt] = Extrema_POnSurf(U,V,E);
384 mySqDist[myNbExt] = Dist2;
388 T.SetRotation(Ax, M_PI);
391 PerformExtPElC(anExt, P1, anACurve, mytolv);
392 if (anExt.IsDone()) {
393 myDone = Standard_True;
397 for (i=1; i<=anExt.NbExt(); i++) {
398 Extrema_POnCurv POC=anExt.Point(i);
401 // if ( !IsExtremum (U, V = myvsup, P, myS, E, Dist2,
402 // Standard_True, anExt.IsMin(i))) continue;
404 Standard_Real newV = myvsup;
406 if((anACurve->GetType() == GeomAbs_Circle) ||
407 (anACurve->GetType() == GeomAbs_Ellipse)) {
408 newV = ElCLib::InPeriod(V, myvinf, myvinf + 2. * M_PI);
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;
423 if((anACurve->GetType() == GeomAbs_Circle) ||
424 (anACurve->GetType() == GeomAbs_Ellipse)) {
425 newV = ElCLib::InPeriod(V, myvsup - 2. * M_PI, myvsup);
432 if ( !IsExtremum (U, V, P, myS, E, Dist2,
433 Standard_False, anExt.IsMin(i))) continue;
436 Dist2 = P.SquareDistance(E);
438 if (IsOriginalPnt(E, myPoint, myNbExt)) {
439 myPoint[++myNbExt] = Extrema_POnSurf(U,V,E);
440 mySqDist[myNbExt] = Dist2;
448 //=======================================================================
451 //=======================================================================
453 Standard_Boolean Extrema_ExtPRevS::IsDone() const
459 //=======================================================================
462 //=======================================================================
464 Standard_Integer Extrema_ExtPRevS::NbExt() const
466 if (!IsDone()) { StdFail_NotDone::Raise(); }
471 //=======================================================================
474 //=======================================================================
476 Standard_Real Extrema_ExtPRevS::SquareDistance(const Standard_Integer N) const
478 if (!IsDone()) { StdFail_NotDone::Raise(); }
479 if ((N < 1) || (N > myNbExt)) { Standard_OutOfRange::Raise(); }
480 if (myIsAnalyticallyComputable)
483 return myExtPS.SquareDistance(N);
485 //=======================================================================
488 //=======================================================================
490 Extrema_POnSurf Extrema_ExtPRevS::Point(const Standard_Integer N) const
492 if (!IsDone()) { StdFail_NotDone::Raise(); }
493 if ((N < 1) || (N > myNbExt)) { Standard_OutOfRange::Raise(); }
494 if (myIsAnalyticallyComputable)
497 return myExtPS.Point(N);