1 // Created on: 1999-09-21
2 // Created by: Edward AGAPOV
3 // Copyright (c) 1999-1999 Matra Datavision
4 // Copyright (c) 1999-2012 OPEN CASCADE SAS
6 // The content of this file is subject to the Open CASCADE Technology Public
7 // License Version 6.5 (the "License"). You may not use the content of this file
8 // except in compliance with the License. Please obtain a copy of the License
9 // at http://www.opencascade.org and read it completely before using this file.
11 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
12 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
14 // The Original Code and all software distributed under the License is
15 // distributed on an "AS IS" basis, without warranty of any kind, and the
16 // Initial Developer hereby disclaims all such warranties, including without
17 // limitation, any warranties of merchantability, fitness for a particular
18 // purpose or non-infringement. Please see the License for the specific terms
19 // and conditions governing the rights and limitations under the License.
23 #include <Extrema_ExtPRevS.ixx>
24 #include <Adaptor3d_HCurve.hxx>
25 #include <Extrema_ExtPElC.hxx>
26 #include <Extrema_GenExtPS.hxx>
27 #include <Extrema_POnCurv.hxx>
28 #include <Extrema_POnSurf.hxx>
34 #include <gp_Trsf.hxx>
36 #include <Precision.hxx>
39 static gp_Ax2 GetPosition (const Adaptor3d_SurfaceOfRevolution& S)//const Handle(Adaptor_HCurve)& C)
41 Handle(Adaptor3d_HCurve) C = S.BasisCurve();
43 switch (C->GetType()) {
47 gp_Dir N = S.AxeOfRevolution().Direction();
48 if (N.IsParallel(L.Direction(), Precision::Angular())) {
49 gp_Vec OO (L.Location(), S.AxeOfRevolution().Location());
50 if (OO.Magnitude() <= gp::Resolution()) {
51 OO = gp_Vec(L.Location(), ElCLib::Value(100,L));
52 if (N.IsParallel(OO, Precision::Angular()))
53 return gp_Ax2(); // Line and axe of revolution coinside
60 return gp_Ax2 (L.Location(), N, L.Direction());
63 return C->Circle().Position();
65 return C->Ellipse().Position();
66 case GeomAbs_Hyperbola:
67 return C->Hyperbola().Position();
68 case GeomAbs_Parabola:
69 return C->Parabola().Position();
75 //=======================================================================
76 //function : HasSingularity
78 //=======================================================================
80 static Standard_Boolean HasSingularity(const Adaptor3d_SurfaceOfRevolution& S)
83 const Handle(Adaptor3d_HCurve) C = S.BasisCurve();
84 gp_Dir N = S.AxeOfRevolution().Direction();
85 gp_Pnt P = S.AxeOfRevolution().Location();
89 P = C->Value(C->FirstParameter());
91 if(L.SquareDistance(P) < Precision::SquareConfusion()) return Standard_True;
93 P = C->Value(C->LastParameter());
95 if(L.SquareDistance(P) < Precision::SquareConfusion()) return Standard_True;
97 return Standard_False;
100 //=============================================================================
102 static void PerformExtPElC (Extrema_ExtPElC& E,
104 const Handle(Adaptor3d_HCurve)& C,
105 const Standard_Real Tol)
107 switch (C->GetType()) {
108 case GeomAbs_Hyperbola:
109 E.Perform(P, C->Hyperbola(), Tol, -Precision::Infinite(),Precision::Infinite());
112 E.Perform(P, C->Line(), Tol, -Precision::Infinite(),Precision::Infinite());
115 E.Perform(P, C->Circle(), Tol, 0.0, 2.0 * M_PI);
117 case GeomAbs_Ellipse:
118 E.Perform(P, C->Ellipse(), Tol, 0.0, 2.0 * M_PI);
120 case GeomAbs_Parabola:
121 E.Perform(P, C->Parabola(), Tol, -Precision::Infinite(),Precision::Infinite());
128 //=======================================================================
129 //function : IsCaseAnalyticallyComputable
131 //=======================================================================
133 static Standard_Boolean IsCaseAnalyticallyComputable
134 (const GeomAbs_CurveType& theType,
135 const gp_Ax2& theCurvePos,
136 const gp_Ax1& AxeOfRevolution)
142 case GeomAbs_Ellipse:
143 case GeomAbs_Hyperbola:
144 case GeomAbs_Parabola:
147 return Standard_False;
149 // the axe of revolution must be in the plane of the curve.
150 gp_Pln pl(theCurvePos.Location(), theCurvePos.Direction());
151 gp_Pnt p1 = AxeOfRevolution.Location();
152 Standard_Real dist = 100., dist2 = dist * dist;
153 Standard_Real aThreshold = Precision::Angular() * Precision::Angular() * dist2;
154 gp_Pnt p2 = AxeOfRevolution.Location().XYZ() + dist * AxeOfRevolution.Direction().XYZ();
156 if((pl.SquareDistance(p1) < aThreshold) &&
157 (pl.SquareDistance(p2) < aThreshold))
158 return Standard_True;
159 return Standard_False;
160 // gp_Vec V (AxeOfRevolution.Location(),theCurvePos.Location());
161 // if (Abs( V * theCurvePos.Direction()) <= gp::Resolution())
162 // return Standard_True;
164 // return Standard_False;
167 //=======================================================================
168 //function : IsOriginalPnt
170 //=======================================================================
172 static Standard_Boolean IsOriginalPnt (const gp_Pnt& P,
173 const Extrema_POnSurf* Points,
174 const Standard_Integer NbPoints)
176 for (Standard_Integer i=1; i<=NbPoints; i++) {
177 if (Points[i].Value().IsEqual(P, Precision::Confusion())) {
178 return Standard_False;
181 return Standard_True;
183 //=======================================================================
184 //function : IsExtremum
186 //=======================================================================
188 static Standard_Boolean IsExtremum (const Standard_Real U, const Standard_Real V,
189 const gp_Pnt& P, const Adaptor3d_SurfacePtr& S,
190 gp_Pnt& E, Standard_Real& Dist2,
191 const Standard_Boolean IsVSup,
192 const Standard_Boolean IsMin)
195 Dist2 = P.SquareDistance(E);
197 return (Dist2 < P.SquareDistance(S->Value(U+1,V)) &&
198 Dist2 < P.SquareDistance(S->Value(U-1,V)) &&
199 Dist2 < P.SquareDistance(S->Value(U, IsVSup ? V-1 : V+1)));
201 return (Dist2 > P.SquareDistance(S->Value(U+1,V)) &&
202 Dist2 > P.SquareDistance(S->Value(U-1,V)) &&
203 Dist2 > P.SquareDistance(S->Value(U, IsVSup ? V-1 : V+1)));
205 //=======================================================================
206 //function : Extrema_ExtPRevS
208 //=======================================================================
210 Extrema_ExtPRevS::Extrema_ExtPRevS()
213 myDone = Standard_False;
215 //=======================================================================
216 //function : Extrema_ExtPRevS
218 //=======================================================================
220 Extrema_ExtPRevS::Extrema_ExtPRevS(const gp_Pnt& P,
221 const Adaptor3d_SurfaceOfRevolution& S,
222 const Standard_Real Umin,
223 const Standard_Real Usup,
224 const Standard_Real Vmin,
225 const Standard_Real Vsup,
226 const Standard_Real TolU,
227 const Standard_Real TolV)
231 Umin, Usup, Vmin, Vsup,
235 //=======================================================================
236 //function : Extrema_ExtPRevS
238 //=======================================================================
240 Extrema_ExtPRevS::Extrema_ExtPRevS(const gp_Pnt& P,
241 const Adaptor3d_SurfaceOfRevolution& S,
242 const Standard_Real TolU,
243 const Standard_Real TolV)
247 S.FirstUParameter(), S.LastUParameter(),
248 S.FirstVParameter(), S.LastVParameter(),
252 //=======================================================================
253 //function : Initialize
255 //=======================================================================
257 void Extrema_ExtPRevS::Initialize(const Adaptor3d_SurfaceOfRevolution& S,
258 const Standard_Real Umin,
259 const Standard_Real Usup,
260 const Standard_Real Vmin,
261 const Standard_Real Vsup,
262 const Standard_Real TolU,
263 const Standard_Real TolV)
269 Handle(Adaptor3d_HCurve) anACurve = S.BasisCurve();
270 if (!myS || myS != (Adaptor3d_SurfacePtr)&S) {
271 myS = Adaptor3d_SurfacePtr(&S);
272 myPosition = GetPosition(S);
273 myIsAnalyticallyComputable =
274 IsCaseAnalyticallyComputable (anACurve->GetType(),myPosition,S.AxeOfRevolution());
276 if (!myIsAnalyticallyComputable) {
278 Standard_Integer nbu = 32, nbv = 32;
280 if(HasSingularity(S)) {nbv = 100;}
282 myExtPS.Initialize(S, nbu, nbv,
283 Umin, Usup, Vmin, Vsup,
287 //=======================================================================
290 //=======================================================================
292 void Extrema_ExtPRevS::Perform(const gp_Pnt& P)
294 myDone = Standard_False;
297 if (!myIsAnalyticallyComputable) {
300 myDone = myExtPS.IsDone();
301 myNbExt = myExtPS.NbExt();
305 Handle(Adaptor3d_HCurve) anACurve = myS->BasisCurve();
307 gp_Ax1 Ax = myS->AxeOfRevolution();
308 gp_Vec Dir = Ax.Direction(), Z = myPosition.Direction();
309 gp_Pnt O = Ax.Location();
311 Standard_Real OPdir = gp_Vec(O,P).Dot(Dir);
312 gp_Pnt Pp = P.Translated(Dir.Multiplied(-OPdir));
313 if (O.IsEqual(Pp,Precision::Confusion())) // P is on the AxeOfRevolution
318 Standard_Real OPpz = gp_Vec(O,Pp).Dot(Z);
319 if (Abs(OPpz) <= gp::Resolution()) {
324 Ppp = Pp.Translated(Z.Multiplied(-OPpz));
325 if (O.IsEqual(Ppp,Precision::Confusion()))
328 U = gp_Vec(O,Ppp).AngleWithRef(gp_Vec(O,Pp),Dir);
332 gp_Vec OPpp (O,Ppp), OPq (O, myS->Value(M_PI/2,0));
334 if (Abs(OPq.Magnitude()) <= gp::Resolution())
335 OPq = gp_Vec(O, myS->Value(M_PI/2,anACurve->LastParameter()/10));
336 if (OPpp.AngleWithRef(OPq,Dir) < 0)
341 T.SetRotation(Ax, -U);
342 P1 = P.Transformed(T);
348 Extrema_ExtPElC anExt;
349 PerformExtPElC(anExt, P1, anACurve, mytolv);
351 if (anExt.IsDone()) {
352 myDone = Standard_True;
353 for (i=1; i<=anExt.NbExt(); i++) {
354 Extrema_POnCurv POC=anExt.Point(i);
357 // if ( !IsExtremum (U, V = myvsup, P, myS, E, Dist2,
358 // Standard_True, anExt.IsMin(i))) continue;
359 Standard_Real newV = myvsup;
361 if((anACurve->GetType() == GeomAbs_Circle) ||
362 (anACurve->GetType() == GeomAbs_Ellipse)) {
363 newV = ElCLib::InPeriod(V, myvinf, myvinf + 2. * M_PI);
368 if (newV + mytolv < myvinf) {
370 } else if (newV < myvinf) {
377 if ( !IsExtremum (U, V, P, myS, E, Dist2,
378 Standard_True, anExt.IsMin(i))) {
381 } else if (V < myvinf) {
382 // if ( !IsExtremum (U, V = myvinf, P, myS, E, Dist2,
383 // Standard_False, anExt.IsMin(i))) continue;
385 Standard_Real newV = myvinf;
387 if((anACurve->GetType() == GeomAbs_Circle) ||
388 (anACurve->GetType() == GeomAbs_Ellipse)) {
389 newV = ElCLib::InPeriod(V, myvsup - 2. * M_PI, myvsup);
394 if (newV - mytolv > myvsup) {
396 } else if (newV > myvsup) {
403 if ( !IsExtremum (U, V, P, myS, E, Dist2,
404 Standard_False, anExt.IsMin(i))) continue;
407 Dist2 = P.SquareDistance(E);
409 if (IsOriginalPnt(E, myPoint, myNbExt)) {
410 myPoint[++myNbExt] = Extrema_POnSurf(U,V,E);
411 mySqDist[myNbExt] = Dist2;
415 T.SetRotation(Ax, M_PI);
418 PerformExtPElC(anExt, P1, anACurve, mytolv);
419 if (anExt.IsDone()) {
420 myDone = Standard_True;
424 for (i=1; i<=anExt.NbExt(); i++) {
425 Extrema_POnCurv POC=anExt.Point(i);
428 // if ( !IsExtremum (U, V = myvsup, P, myS, E, Dist2,
429 // Standard_True, anExt.IsMin(i))) continue;
431 Standard_Real newV = myvsup;
433 if((anACurve->GetType() == GeomAbs_Circle) ||
434 (anACurve->GetType() == GeomAbs_Ellipse)) {
435 newV = ElCLib::InPeriod(V, myvinf, myvinf + 2. * M_PI);
440 if (newV + mytolv < myvinf) {
442 } else if (newV < myvinf) {
449 if ( !IsExtremum (U, V, P, myS, E, Dist2,
450 Standard_True, anExt.IsMin(i))) continue;
451 } else if (V < myvinf) {
452 // if ( !IsExtremum (U, V = myvinf, P, myS, E, Dist2,
453 // Standard_False, anExt.IsMin(i))) continue;
454 Standard_Real newV = myvinf;
456 if((anACurve->GetType() == GeomAbs_Circle) ||
457 (anACurve->GetType() == GeomAbs_Ellipse)) {
458 newV = ElCLib::InPeriod(V, myvsup - 2. * M_PI, myvsup);
463 if (newV - mytolv > myvsup) {
465 } else if (newV > myvsup) {
472 if ( !IsExtremum (U, V, P, myS, E, Dist2,
473 Standard_False, anExt.IsMin(i))) continue;
476 Dist2 = P.SquareDistance(E);
478 if (IsOriginalPnt(E, myPoint, myNbExt)) {
479 myPoint[++myNbExt] = Extrema_POnSurf(U,V,E);
480 mySqDist[myNbExt] = Dist2;
488 //=======================================================================
491 //=======================================================================
493 Standard_Boolean Extrema_ExtPRevS::IsDone() const
499 //=======================================================================
502 //=======================================================================
504 Standard_Integer Extrema_ExtPRevS::NbExt() const
506 if (!IsDone()) { StdFail_NotDone::Raise(); }
511 //=======================================================================
514 //=======================================================================
516 Standard_Real Extrema_ExtPRevS::SquareDistance(const Standard_Integer N) const
518 if (!IsDone()) { StdFail_NotDone::Raise(); }
519 if ((N < 1) || (N > myNbExt)) { Standard_OutOfRange::Raise(); }
520 if (myIsAnalyticallyComputable)
523 return myExtPS.SquareDistance(N);
525 //=======================================================================
528 //=======================================================================
530 Extrema_POnSurf Extrema_ExtPRevS::Point(const Standard_Integer N) const
532 if (!IsDone()) { StdFail_NotDone::Raise(); }
533 if ((N < 1) || (N > myNbExt)) { Standard_OutOfRange::Raise(); }
534 if (myIsAnalyticallyComputable)
537 return myExtPS.Point(N);