[occt.git] / src / Extrema / Extrema_ExtPRevS.cxx

// Created on: 1999-09-21
// Created by: Edward AGAPOV
// Copyright (c) 1999-1999 Matra Datavision
// Copyright (c) 1999-2012 OPEN CASCADE SAS
//
// The content of this file is subject to the Open CASCADE Technology Public
// License Version 6.5 (the "License"). You may not use the content of this file
// except in compliance with the License. Please obtain a copy of the License
// at http://www.opencascade.org and read it completely before using this file.
//
// The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
// main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
//
// The Original Code and all software distributed under the License is
// distributed on an "AS IS" basis, without warranty of any kind, and the
// Initial Developer hereby disclaims all such warranties, including without
// limitation, any warranties of merchantability, fitness for a particular
// purpose or non-infringement. Please see the License for the specific terms
// and conditions governing the rights and limitations under the License.


#include <Extrema_ExtPRevS.ixx>
#include <Adaptor3d_HCurve.hxx>
#include <Extrema_ExtPElC.hxx>
#include <Extrema_GenExtPS.hxx>
#include <Extrema_POnCurv.hxx>
#include <Extrema_POnSurf.hxx>
#include <gp_Ax1.hxx>
#include <gp_Ax2.hxx>
#include <gp_Lin.hxx>
#include <gp_Pln.hxx>
#include <gp_Pnt.hxx>
#include <gp_Trsf.hxx>
#include <gp_Vec.hxx>
#include <Precision.hxx>
#include <ElCLib.hxx>

static gp_Ax2 GetPosition (const Adaptor3d_SurfaceOfRevolution& S)//const Handle(Adaptor_HCurve)& C)
{
  Handle(Adaptor3d_HCurve) C = S.BasisCurve();
  
  switch (C->GetType()) {
    
  case GeomAbs_Line: {
    gp_Lin L = C->Line();
    gp_Dir N = S.AxeOfRevolution().Direction();
    if (N.IsParallel(L.Direction(), Precision::Angular())) {
      gp_Vec OO (L.Location(), S.AxeOfRevolution().Location());
      if (OO.Magnitude() <= gp::Resolution()) {
	OO = gp_Vec(L.Location(), ElCLib::Value(100,L));
	if (N.IsParallel(OO, Precision::Angular()))
	  return gp_Ax2(); // Line and axe of revolution coinside
      }
      N ^= OO;
    }
    else {
      N ^= L.Direction();
    }
    return gp_Ax2 (L.Location(), N, L.Direction());
  }
  case GeomAbs_Circle:
    return C->Circle().Position();
  case GeomAbs_Ellipse:
    return C->Ellipse().Position();
  case GeomAbs_Hyperbola:
    return C->Hyperbola().Position();
  case GeomAbs_Parabola:
    return C->Parabola().Position();
  default:
    return gp_Ax2();
  }
}

//=======================================================================
//function : HasSingularity
//purpose  : 
//=======================================================================

static Standard_Boolean HasSingularity(const Adaptor3d_SurfaceOfRevolution& S) 
{

  const Handle(Adaptor3d_HCurve) C = S.BasisCurve();
  gp_Dir N = S.AxeOfRevolution().Direction();
  gp_Pnt P = S.AxeOfRevolution().Location();

  gp_Lin L(P, N);

  P = C->Value(C->FirstParameter());

  if(L.SquareDistance(P) < Precision::SquareConfusion()) return Standard_True;

  P = C->Value(C->LastParameter());

  if(L.SquareDistance(P) < Precision::SquareConfusion()) return Standard_True;
  
  return Standard_False;
}

//=============================================================================

static void PerformExtPElC (Extrema_ExtPElC& E,
			    const gp_Pnt& P,
			    const Handle(Adaptor3d_HCurve)& C,
			    const Standard_Real Tol)
{
  switch (C->GetType()) {
  case GeomAbs_Hyperbola:
    E.Perform(P, C->Hyperbola(), Tol, -Precision::Infinite(),Precision::Infinite());
    return;
  case GeomAbs_Line:
    E.Perform(P, C->Line(), Tol, -Precision::Infinite(),Precision::Infinite());
    return;
  case GeomAbs_Circle:
    E.Perform(P, C->Circle(), Tol, 0.0, 2.0 * M_PI);
    return;
  case GeomAbs_Ellipse:
    E.Perform(P, C->Ellipse(), Tol, 0.0, 2.0 * M_PI);
    return;
  case GeomAbs_Parabola:
    E.Perform(P, C->Parabola(), Tol, -Precision::Infinite(),Precision::Infinite());
    return;
  default:
    return ;
  }
}

//=======================================================================
//function : IsCaseAnalyticallyComputable
//purpose  : 
//=======================================================================

static Standard_Boolean IsCaseAnalyticallyComputable
  (const GeomAbs_CurveType& theType,
   const gp_Ax2& theCurvePos,
   const gp_Ax1& AxeOfRevolution) 
{
  // check type
  switch (theType) {
  case GeomAbs_Line:
  case GeomAbs_Circle:
  case GeomAbs_Ellipse:
  case GeomAbs_Hyperbola:
  case GeomAbs_Parabola:
    break;
  default:
    return  Standard_False;
  }
//  the axe of revolution must be in the plane of the curve.
  gp_Pln pl(theCurvePos.Location(), theCurvePos.Direction());
  gp_Pnt p1 = AxeOfRevolution.Location();
  Standard_Real dist = 100., dist2 = dist * dist;
  Standard_Real aThreshold = Precision::Angular() * Precision::Angular() * dist2;
  gp_Pnt p2 = AxeOfRevolution.Location().XYZ() + dist * AxeOfRevolution.Direction().XYZ();

  if((pl.SquareDistance(p1) < aThreshold) && 
     (pl.SquareDistance(p2) < aThreshold))
    return Standard_True;
  return Standard_False;
  //   gp_Vec V (AxeOfRevolution.Location(),theCurvePos.Location());
  //   if (Abs( V * theCurvePos.Direction()) <= gp::Resolution())
  //     return Standard_True;
  //   else
  //     return Standard_False;
}

//=======================================================================
//function : IsOriginalPnt
//purpose  : 
//=======================================================================

static Standard_Boolean IsOriginalPnt (const gp_Pnt& P,
				       const Extrema_POnSurf* Points,
				       const Standard_Integer NbPoints)
{
  for (Standard_Integer i=1; i<=NbPoints; i++) {
    if (Points[i].Value().IsEqual(P, Precision::Confusion())) {
      return Standard_False;
    }
  }
  return Standard_True;
}
//=======================================================================
//function : IsExtremum
//purpose  : 
//=======================================================================

static Standard_Boolean IsExtremum (const Standard_Real U, const Standard_Real V,
				    const gp_Pnt& P, const Adaptor3d_SurfacePtr& S,
				    gp_Pnt& E,        Standard_Real& Dist2,
				    const Standard_Boolean IsVSup,
				    const Standard_Boolean IsMin)
{
  E = S->Value(U,V);
  Dist2 = P.SquareDistance(E);
  if (IsMin) 
    return (Dist2 < P.SquareDistance(S->Value(U+1,V)) &&
	    Dist2 < P.SquareDistance(S->Value(U-1,V)) &&
	    Dist2 < P.SquareDistance(S->Value(U, IsVSup ? V-1 : V+1)));
  else
    return (Dist2 > P.SquareDistance(S->Value(U+1,V)) &&
	    Dist2 > P.SquareDistance(S->Value(U-1,V)) &&
	    Dist2 > P.SquareDistance(S->Value(U, IsVSup ? V-1 : V+1)));
}  
//=======================================================================
//function : Extrema_ExtPRevS
//purpose  : 
//=======================================================================

Extrema_ExtPRevS::Extrema_ExtPRevS() 
{
  myS=NULL;
  myDone = Standard_False;
}
//=======================================================================
//function : Extrema_ExtPRevS
//purpose  : 
//=======================================================================

Extrema_ExtPRevS::Extrema_ExtPRevS(const gp_Pnt& P,
				   const Adaptor3d_SurfaceOfRevolution& S,
				   const Standard_Real Umin,
				   const Standard_Real Usup,
				   const Standard_Real Vmin,
				   const Standard_Real Vsup,
				   const Standard_Real TolU,
				   const Standard_Real TolV) 
{
  myS=NULL;
  Initialize (S,
	      Umin, Usup, Vmin, Vsup,
	      TolU, TolV);
  Perform(P);
}
//=======================================================================
//function : Extrema_ExtPRevS
//purpose  : 
//=======================================================================

Extrema_ExtPRevS::Extrema_ExtPRevS(const gp_Pnt& P,
				   const Adaptor3d_SurfaceOfRevolution& S,
				   const Standard_Real TolU,
				   const Standard_Real TolV) 
{
  myS=NULL;
  Initialize (S,
	      S.FirstUParameter(), S.LastUParameter(),
	      S.FirstVParameter(), S.LastVParameter(),
	      TolU, TolV);
  Perform(P);
}
//=======================================================================
//function : Initialize
//purpose  : 
//=======================================================================

void Extrema_ExtPRevS::Initialize(const Adaptor3d_SurfaceOfRevolution& S,
				  const Standard_Real Umin,
				  const Standard_Real Usup,
				  const Standard_Real Vmin,
				  const Standard_Real Vsup,
				  const Standard_Real TolU,
				  const Standard_Real TolV)
{
  myvinf=Vmin;
  myvsup=Vsup;
  mytolv=TolV;

  Handle(Adaptor3d_HCurve) anACurve = S.BasisCurve();
  if (!myS || myS != (Adaptor3d_SurfacePtr)&S) {
    myS = Adaptor3d_SurfacePtr(&S);
    myPosition = GetPosition(S);
    myIsAnalyticallyComputable =
      IsCaseAnalyticallyComputable (anACurve->GetType(),myPosition,S.AxeOfRevolution());
  }
  if (!myIsAnalyticallyComputable) {
    
    Standard_Integer nbu = 32, nbv = 32;

    if(HasSingularity(S)) {nbv = 100;}

    myExtPS.Initialize(S, nbu, nbv,
		       Umin, Usup, Vmin, Vsup,
		       TolU, TolV);
  }
}
//=======================================================================
//function : Perform
//purpose  : 
//=======================================================================

void Extrema_ExtPRevS::Perform(const gp_Pnt& P)
{
  myDone = Standard_False;
  myNbExt = 0;
  
  if (!myIsAnalyticallyComputable) {
    
    myExtPS.Perform(P);
    myDone = myExtPS.IsDone();
    myNbExt = myExtPS.NbExt();
    return;
  }
  
  Handle(Adaptor3d_HCurve) anACurve = myS->BasisCurve();

  gp_Ax1 Ax = myS->AxeOfRevolution();
  gp_Vec Dir = Ax.Direction(), Z = myPosition.Direction();
  gp_Pnt O = Ax.Location();

  Standard_Real OPdir = gp_Vec(O,P).Dot(Dir);
  gp_Pnt Pp = P.Translated(Dir.Multiplied(-OPdir));
  if (O.IsEqual(Pp,Precision::Confusion())) // P is on the AxeOfRevolution
    return;
  
  Standard_Real U,V;
  gp_Pnt P1, Ppp;
  Standard_Real OPpz = gp_Vec(O,Pp).Dot(Z);
  if (Abs(OPpz) <= gp::Resolution()) {
    Ppp = Pp;
    U = 0;
  }
  else {
    Ppp = Pp.Translated(Z.Multiplied(-OPpz));
    if (O.IsEqual(Ppp,Precision::Confusion())) 
      U = M_PI/2;
    else {
      U = gp_Vec(O,Ppp).AngleWithRef(gp_Vec(O,Pp),Dir);
    }
  }

  gp_Vec OPpp (O,Ppp), OPq (O, myS->Value(M_PI/2,0));
  if (U != M_PI/2) {
    if (Abs(OPq.Magnitude()) <= gp::Resolution()) 
      OPq = gp_Vec(O, myS->Value(M_PI/2,anACurve->LastParameter()/10));
    if (OPpp.AngleWithRef(OPq,Dir) < 0)
      U += M_PI;
  }
  
  gp_Trsf T;
  T.SetRotation(Ax, -U);
  P1 = P.Transformed(T);
  
  gp_Pnt E;
  Standard_Real Dist2;
  Standard_Integer i;
  
  Extrema_ExtPElC anExt;
  PerformExtPElC(anExt, P1, anACurve, mytolv);
  
  if (anExt.IsDone()) {
    myDone = Standard_True;
    for (i=1; i<=anExt.NbExt(); i++) {
      Extrema_POnCurv POC=anExt.Point(i);
      V = POC.Parameter();
      if (V > myvsup) {
	// 	 if ( !IsExtremum (U, V = myvsup, P, myS, E, Dist2,
	// 			   Standard_True, anExt.IsMin(i))) continue;
	Standard_Real newV = myvsup;

	if((anACurve->GetType() == GeomAbs_Circle) || 
	   (anACurve->GetType() == GeomAbs_Ellipse)) {
	  newV = ElCLib::InPeriod(V, myvinf, myvinf + 2. * M_PI);

	  if (newV > myvsup) {
            newV -= 2. * M_PI;

            if (newV + mytolv < myvinf) {
              newV = myvsup;
            } else if (newV < myvinf) {
              newV = myvinf;
            }
	  }
	}
	V = newV;

	if ( !IsExtremum (U, V, P, myS, E, Dist2,
			   Standard_True, anExt.IsMin(i))) {
	  continue;
	}
      } else if (V < myvinf) {
	// 	if ( !IsExtremum (U, V = myvinf, P, myS, E, Dist2,
	// 			  Standard_False, anExt.IsMin(i))) continue;

	Standard_Real newV = myvinf;

	if((anACurve->GetType() == GeomAbs_Circle) || 
	   (anACurve->GetType() == GeomAbs_Ellipse)) {
	  newV = ElCLib::InPeriod(V, myvsup - 2. * M_PI, myvsup);
	  
          if(newV < myvinf) {
            newV += 2. * M_PI;
 
            if (newV - mytolv > myvsup) {
              newV = myvinf;
            } else if (newV > myvsup) {
              newV = myvsup;
            }
          }
	}
	V = newV;

	if ( !IsExtremum (U, V, P, myS, E, Dist2,
			  Standard_False, anExt.IsMin(i))) continue;
      } else {
	E = myS->Value(U,V);
	Dist2 = P.SquareDistance(E);
      }
      if (IsOriginalPnt(E, myPoint, myNbExt)) {
	myPoint[++myNbExt] = Extrema_POnSurf(U,V,E);
	mySqDist[myNbExt] = Dist2;
      }
    }
  }
  T.SetRotation(Ax, M_PI);
  P1.Transform(T);
  
  PerformExtPElC(anExt, P1, anACurve, mytolv);
  if (anExt.IsDone()) {
    myDone = Standard_True;

    U += M_PI;
    
    for (i=1; i<=anExt.NbExt(); i++) {
      Extrema_POnCurv POC=anExt.Point(i);
      V = POC.Parameter();
      if (V > myvsup) {
	// 	if ( !IsExtremum (U, V = myvsup, P, myS, E, Dist2,
	// 			   Standard_True, anExt.IsMin(i))) continue;

	Standard_Real newV = myvsup;

	if((anACurve->GetType() == GeomAbs_Circle) || 
	   (anACurve->GetType() == GeomAbs_Ellipse)) {
	  newV = ElCLib::InPeriod(V, myvinf, myvinf + 2. * M_PI);

	  if (newV > myvsup) {
            newV -= 2. * M_PI;

            if (newV + mytolv < myvinf) {
              newV = myvsup;
            } else if (newV < myvinf) {
              newV = myvinf;
            }
	  }
	}
	V = newV;
	
	if ( !IsExtremum (U, V, P, myS, E, Dist2,
			   Standard_True, anExt.IsMin(i))) continue;
      } else if (V < myvinf) {
	// 	if ( !IsExtremum (U, V = myvinf, P, myS, E, Dist2,
	// 			  Standard_False, anExt.IsMin(i))) continue;
	Standard_Real newV = myvinf;

	if((anACurve->GetType() == GeomAbs_Circle) || 
	   (anACurve->GetType() == GeomAbs_Ellipse)) {
	  newV = ElCLib::InPeriod(V, myvsup - 2. * M_PI, myvsup);
	  
          if(newV < myvinf) {
            newV += 2. * M_PI;
 
            if (newV - mytolv > myvsup) {
              newV = myvinf;
            } else if (newV > myvsup) {
              newV = myvsup;
            }
          }
	}
	V = newV;

	if ( !IsExtremum (U, V, P, myS, E, Dist2,
			  Standard_False, anExt.IsMin(i))) continue;
      } else {
	E = myS->Value(U,V);
	Dist2 = P.SquareDistance(E);
      }
      if (IsOriginalPnt(E, myPoint, myNbExt)) {
	myPoint[++myNbExt] = Extrema_POnSurf(U,V,E);
	mySqDist[myNbExt] = Dist2;
      
      }
    }
  }
}


//=======================================================================
//function : IsDone
//purpose  : 
//=======================================================================

Standard_Boolean Extrema_ExtPRevS::IsDone() const
{
  return myDone; 
}


//=======================================================================
//function : NbExt
//purpose  : 
//=======================================================================

Standard_Integer Extrema_ExtPRevS::NbExt() const
{
  if (!IsDone()) { StdFail_NotDone::Raise(); }
  return myNbExt;
}


//=======================================================================
//function : Value
//purpose  : 
//=======================================================================

Standard_Real Extrema_ExtPRevS::SquareDistance(const Standard_Integer N) const
{
  if (!IsDone()) { StdFail_NotDone::Raise(); }
  if ((N < 1) || (N > myNbExt)) { Standard_OutOfRange::Raise(); }
  if (myIsAnalyticallyComputable)
    return mySqDist[N];
  else
    return myExtPS.SquareDistance(N);
}
//=======================================================================
//function : Point
//purpose  : 
//=======================================================================

Extrema_POnSurf Extrema_ExtPRevS::Point(const Standard_Integer N) const
{
  if (!IsDone()) { StdFail_NotDone::Raise(); }
  if ((N < 1) || (N > myNbExt)) { Standard_OutOfRange::Raise(); }
  if (myIsAnalyticallyComputable)
    return myPoint[N];
  else
    return myExtPS.Point(N);
}
Commit	Line	Data
b311480e	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
	5	//
	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.
	10	//
	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.
	13	//
	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.
	20
7fd59977	21
	22
	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>
	29	#include <gp_Ax1.hxx>
	30	#include <gp_Ax2.hxx>
	31	#include <gp_Lin.hxx>
	32	#include <gp_Pln.hxx>
	33	#include <gp_Pnt.hxx>
	34	#include <gp_Trsf.hxx>
	35	#include <gp_Vec.hxx>
	36	#include <Precision.hxx>
	37	#include <ElCLib.hxx>
	38
	39	static gp_Ax2 GetPosition (const Adaptor3d_SurfaceOfRevolution& S)//const Handle(Adaptor_HCurve)& C)
	40	{
	41	Handle(Adaptor3d_HCurve) C = S.BasisCurve();
	42
	43	switch (C->GetType()) {
	44
	45	case GeomAbs_Line: {
	46	gp_Lin L = C->Line();
	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
	54	}
	55	N ^= OO;
	56	}
	57	else {
	58	N ^= L.Direction();
	59	}
	60	return gp_Ax2 (L.Location(), N, L.Direction());
	61	}
	62	case GeomAbs_Circle:
	63	return C->Circle().Position();
	64	case GeomAbs_Ellipse:
	65	return C->Ellipse().Position();
	66	case GeomAbs_Hyperbola:
	67	return C->Hyperbola().Position();
	68	case GeomAbs_Parabola:
	69	return C->Parabola().Position();
	70	default:
	71	return gp_Ax2();
	72	}
	73	}
	74
	75	//=======================================================================
	76	//function : HasSingularity
	77	//purpose :
	78	//=======================================================================
	79
	80	static Standard_Boolean HasSingularity(const Adaptor3d_SurfaceOfRevolution& S)
	81	{
	82
	83	const Handle(Adaptor3d_HCurve) C = S.BasisCurve();
	84	gp_Dir N = S.AxeOfRevolution().Direction();
85	gp_Pnt P = S.AxeOfRevolution().Location();
86
87	gp_Lin L(P, N);
88
89	P = C->Value(C->FirstParameter());
90
08cd2f6b	91	if(L.SquareDistance(P) < Precision::SquareConfusion()) return Standard_True;
7fd59977	92
	93	P = C->Value(C->LastParameter());
	94
08cd2f6b	95	if(L.SquareDistance(P) < Precision::SquareConfusion()) return Standard_True;
7fd59977	96
	97	return Standard_False;
	98	}
	99
	100	//=============================================================================
	101
	102	static void PerformExtPElC (Extrema_ExtPElC& E,
	103	const gp_Pnt& P,
	104	const Handle(Adaptor3d_HCurve)& C,
	105	const Standard_Real Tol)
	106	{
	107	switch (C->GetType()) {
	108	case GeomAbs_Hyperbola:
	109	E.Perform(P, C->Hyperbola(), Tol, -Precision::Infinite(),Precision::Infinite());
	110	return;
	111	case GeomAbs_Line:
	112	E.Perform(P, C->Line(), Tol, -Precision::Infinite(),Precision::Infinite());
	113	return;
	114	case GeomAbs_Circle:
c6541a0c	115	E.Perform(P, C->Circle(), Tol, 0.0, 2.0 * M_PI);
7fd59977	116	return;
7fd59977	117	case GeomAbs_Ellipse:
c6541a0c	118	E.Perform(P, C->Ellipse(), Tol, 0.0, 2.0 * M_PI);
7fd59977	119	return;
	120	case GeomAbs_Parabola:
	121	E.Perform(P, C->Parabola(), Tol, -Precision::Infinite(),Precision::Infinite());
	122	return;
7fd59977	123	default:
7fd59977	124	return ;
7fd59977	125	}
	126	}
	127
	128	//=======================================================================
	129	//function : IsCaseAnalyticallyComputable
	130	//purpose :
	131	//=======================================================================
	132
	133	static Standard_Boolean IsCaseAnalyticallyComputable
	134	(const GeomAbs_CurveType& theType,
	135	const gp_Ax2& theCurvePos,
	136	const gp_Ax1& AxeOfRevolution)
	137	{
	138	// check type
	139	switch (theType) {
	140	case GeomAbs_Line:
	141	case GeomAbs_Circle:
	142	case GeomAbs_Ellipse:
	143	case GeomAbs_Hyperbola:
	144	case GeomAbs_Parabola:
	145	break;
	146	default:
	147	return Standard_False;
	148	}
	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();
	155
	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;
	163	// else
	164	// return Standard_False;
	165	}
	166
	167	//=======================================================================
	168	//function : IsOriginalPnt
	169	//purpose :
	170	//=======================================================================
	171
	172	static Standard_Boolean IsOriginalPnt (const gp_Pnt& P,
	173	const Extrema_POnSurf* Points,
	174	const Standard_Integer NbPoints)
	175	{
	176	for (Standard_Integer i=1; i<=NbPoints; i++) {
	177	if (Points[i].Value().IsEqual(P, Precision::Confusion())) {
	178	return Standard_False;
	179	}
	180	}
	181	return Standard_True;
	182	}
	183	//=======================================================================
	184	//function : IsExtremum
	185	//purpose :
	186	//=======================================================================
	187
	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)
193	{
194	E = S->Value(U,V);
195	Dist2 = P.SquareDistance(E);
196	if (IsMin)
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)));
200	else
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)));
204	}
205	//=======================================================================
206	//function : Extrema_ExtPRevS
207	//purpose :
208	//=======================================================================
209
210	Extrema_ExtPRevS::Extrema_ExtPRevS()
211	{
212	myS=NULL;
213	myDone = Standard_False;
214	}
215	//=======================================================================
216	//function : Extrema_ExtPRevS
217	//purpose :
218	//=======================================================================
219
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)
228	{
229	myS=NULL;
230	Initialize (S,
231	Umin, Usup, Vmin, Vsup,
232	TolU, TolV);
233	Perform(P);
234	}
235	//=======================================================================
236	//function : Extrema_ExtPRevS
237	//purpose :
238	//=======================================================================
239
240	Extrema_ExtPRevS::Extrema_ExtPRevS(const gp_Pnt& P,
241	const Adaptor3d_SurfaceOfRevolution& S,
242	const Standard_Real TolU,
243	const Standard_Real TolV)
244	{
245	myS=NULL;
246	Initialize (S,
247	S.FirstUParameter(), S.LastUParameter(),
248	S.FirstVParameter(), S.LastVParameter(),
249	TolU, TolV);
250	Perform(P);
251	}
252	//=======================================================================
253	//function : Initialize
254	//purpose :
255	//=======================================================================
256
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)
264	{
265	myvinf=Vmin;
266	myvsup=Vsup;
267	mytolv=TolV;
268
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());
275	}
276	if (!myIsAnalyticallyComputable) {
277
278	Standard_Integer nbu = 32, nbv = 32;
279
280	if(HasSingularity(S)) {nbv = 100;}
281
282	myExtPS.Initialize(S, nbu, nbv,
283	Umin, Usup, Vmin, Vsup,
284	TolU, TolV);
285	}
286	}
287	//=======================================================================
288	//function : Perform
289	//purpose :
290	//=======================================================================
291
292	void Extrema_ExtPRevS::Perform(const gp_Pnt& P)
293	{
294	myDone = Standard_False;
295	myNbExt = 0;
296
297	if (!myIsAnalyticallyComputable) {
298
299	myExtPS.Perform(P);
300	myDone = myExtPS.IsDone();
301	myNbExt = myExtPS.NbExt();
302	return;
303	}
304
305	Handle(Adaptor3d_HCurve) anACurve = myS->BasisCurve();
306
307	gp_Ax1 Ax = myS->AxeOfRevolution();
308	gp_Vec Dir = Ax.Direction(), Z = myPosition.Direction();
309	gp_Pnt O = Ax.Location();
310
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
314	return;
315
316	Standard_Real U,V;
317	gp_Pnt P1, Ppp;
318	Standard_Real OPpz = gp_Vec(O,Pp).Dot(Z);
319	if (Abs(OPpz) <= gp::Resolution()) {
320	Ppp = Pp;
321	U = 0;
322	}
323	else {
324	Ppp = Pp.Translated(Z.Multiplied(-OPpz));
325	if (O.IsEqual(Ppp,Precision::Confusion()))
c6541a0c	326	U = M_PI/2;
7fd59977	327	else {
	328	U = gp_Vec(O,Ppp).AngleWithRef(gp_Vec(O,Pp),Dir);
	329	}
	330	}
	331
c6541a0c D	332	gp_Vec OPpp (O,Ppp), OPq (O, myS->Value(M_PI/2,0));
c6541a0c D	333	if (U != M_PI/2) {
7fd59977	334	if (Abs(OPq.Magnitude()) <= gp::Resolution())
c6541a0c	335	OPq = gp_Vec(O, myS->Value(M_PI/2,anACurve->LastParameter()/10));
7fd59977	336	if (OPpp.AngleWithRef(OPq,Dir) < 0)
c6541a0c	337	U += M_PI;
7fd59977	338	}
	339
	340	gp_Trsf T;
	341	T.SetRotation(Ax, -U);
	342	P1 = P.Transformed(T);
	343
	344	gp_Pnt E;
	345	Standard_Real Dist2;
	346	Standard_Integer i;
	347
	348	Extrema_ExtPElC anExt;
	349	PerformExtPElC(anExt, P1, anACurve, mytolv);
	350
	351	if (anExt.IsDone()) {
	352	myDone = Standard_True;
	353	for (i=1; i<=anExt.NbExt(); i++) {
	354	Extrema_POnCurv POC=anExt.Point(i);
	355	V = POC.Parameter();
	356	if (V > myvsup) {
	357	// if ( !IsExtremum (U, V = myvsup, P, myS, E, Dist2,
	358	// Standard_True, anExt.IsMin(i))) continue;
	359	Standard_Real newV = myvsup;
	360
	361	if((anACurve->GetType() == GeomAbs_Circle) \|\|
	362	(anACurve->GetType() == GeomAbs_Ellipse)) {
c6541a0c	363	newV = ElCLib::InPeriod(V, myvinf, myvinf + 2. * M_PI);
7fd59977	364
7fd59977	365	if (newV > myvsup) {
a20400b3	366	newV -= 2. * M_PI;
	367
	368	if (newV + mytolv < myvinf) {
	369	newV = myvsup;
	370	} else if (newV < myvinf) {
	371	newV = myvinf;
	372	}
7fd59977	373	}
	374	}
	375	V = newV;
	376
	377	if ( !IsExtremum (U, V, P, myS, E, Dist2,
	378	Standard_True, anExt.IsMin(i))) {
	379	continue;
	380	}
	381	} else if (V < myvinf) {
	382	// if ( !IsExtremum (U, V = myvinf, P, myS, E, Dist2,
	383	// Standard_False, anExt.IsMin(i))) continue;
	384
	385	Standard_Real newV = myvinf;
	386
	387	if((anACurve->GetType() == GeomAbs_Circle) \|\|
	388	(anACurve->GetType() == GeomAbs_Ellipse)) {
c6541a0c	389	newV = ElCLib::InPeriod(V, myvsup - 2. * M_PI, myvsup);
7fd59977	390
a20400b3	391	if(newV < myvinf) {
	392	newV += 2. * M_PI;
	393
	394	if (newV - mytolv > myvsup) {
	395	newV = myvinf;
	396	} else if (newV > myvsup) {
	397	newV = myvsup;
	398	}
	399	}
7fd59977	400	}
	401	V = newV;
	402
	403	if ( !IsExtremum (U, V, P, myS, E, Dist2,
	404	Standard_False, anExt.IsMin(i))) continue;
	405	} else {
	406	E = myS->Value(U,V);
	407	Dist2 = P.SquareDistance(E);
	408	}
	409	if (IsOriginalPnt(E, myPoint, myNbExt)) {
	410	myPoint[++myNbExt] = Extrema_POnSurf(U,V,E);
	411	mySqDist[myNbExt] = Dist2;
	412	}
	413	}
	414	}
c6541a0c	415	T.SetRotation(Ax, M_PI);
7fd59977	416	P1.Transform(T);
	417
	418	PerformExtPElC(anExt, P1, anACurve, mytolv);
	419	if (anExt.IsDone()) {
	420	myDone = Standard_True;
	421
c6541a0c	422	U += M_PI;
7fd59977	423
	424	for (i=1; i<=anExt.NbExt(); i++) {
	425	Extrema_POnCurv POC=anExt.Point(i);
	426	V = POC.Parameter();
	427	if (V > myvsup) {
	428	// if ( !IsExtremum (U, V = myvsup, P, myS, E, Dist2,
	429	// Standard_True, anExt.IsMin(i))) continue;
	430
	431	Standard_Real newV = myvsup;
	432
	433	if((anACurve->GetType() == GeomAbs_Circle) \|\|
	434	(anACurve->GetType() == GeomAbs_Ellipse)) {
c6541a0c	435	newV = ElCLib::InPeriod(V, myvinf, myvinf + 2. * M_PI);
7fd59977	436
7fd59977	437	if (newV > myvsup) {
a20400b3	438	newV -= 2. * M_PI;
	439
	440	if (newV + mytolv < myvinf) {
	441	newV = myvsup;
	442	} else if (newV < myvinf) {
	443	newV = myvinf;
	444	}
7fd59977	445	}
	446	}
	447	V = newV;
	448
	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;
	455
	456	if((anACurve->GetType() == GeomAbs_Circle) \|\|
	457	(anACurve->GetType() == GeomAbs_Ellipse)) {
c6541a0c	458	newV = ElCLib::InPeriod(V, myvsup - 2. * M_PI, myvsup);
7fd59977	459
a20400b3	460	if(newV < myvinf) {
	461	newV += 2. * M_PI;
	462
	463	if (newV - mytolv > myvsup) {
	464	newV = myvinf;
	465	} else if (newV > myvsup) {
	466	newV = myvsup;
	467	}
	468	}
7fd59977	469	}
	470	V = newV;
	471
	472	if ( !IsExtremum (U, V, P, myS, E, Dist2,
	473	Standard_False, anExt.IsMin(i))) continue;
	474	} else {
	475	E = myS->Value(U,V);
	476	Dist2 = P.SquareDistance(E);
	477	}
	478	if (IsOriginalPnt(E, myPoint, myNbExt)) {
	479	myPoint[++myNbExt] = Extrema_POnSurf(U,V,E);
	480	mySqDist[myNbExt] = Dist2;
	481
	482	}
	483	}
	484	}
	485	}
	486
	487
	488	//=======================================================================
	489	//function : IsDone
	490	//purpose :
	491	//=======================================================================
	492
	493	Standard_Boolean Extrema_ExtPRevS::IsDone() const
	494	{
	495	return myDone;
	496	}
	497
	498
	499	//=======================================================================
	500	//function : NbExt
	501	//purpose :
	502	//=======================================================================
	503
	504	Standard_Integer Extrema_ExtPRevS::NbExt() const
	505	{
	506	if (!IsDone()) { StdFail_NotDone::Raise(); }
	507	return myNbExt;
	508	}
	509
	510
	511	//=======================================================================
	512	//function : Value
	513	//purpose :
	514	//=======================================================================
	515
	516	Standard_Real Extrema_ExtPRevS::SquareDistance(const Standard_Integer N) const
	517	{
	518	if (!IsDone()) { StdFail_NotDone::Raise(); }
	519	if ((N < 1) \|\| (N > myNbExt)) { Standard_OutOfRange::Raise(); }
	520	if (myIsAnalyticallyComputable)
	521	return mySqDist[N];
	522	else
	523	return myExtPS.SquareDistance(N);
	524	}
	525	//=======================================================================
	526	//function : Point
	527	//purpose :
	528	//=======================================================================
	529
	530	Extrema_POnSurf Extrema_ExtPRevS::Point(const Standard_Integer N) const
	531	{
	532	if (!IsDone()) { StdFail_NotDone::Raise(); }
533	if ((N < 1) \|\| (N > myNbExt)) { Standard_OutOfRange::Raise(); }
534	if (myIsAnalyticallyComputable)
535	return myPoint[N];
536	else
537	return myExtPS.Point(N);
538	}
539
540