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

// Copyright (c) 1995-1999 Matra Datavision
// Copyright (c) 1999-2014 OPEN CASCADE SAS
//
// This file is part of Open CASCADE Technology software library.
//
// This library is free software; you can redistribute it and/or modify it under
// the terms of the GNU Lesser General Public License version 2.1 as published
// by the Free Software Foundation, with special exception defined in the file
// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
// distribution for complete text of the license and disclaimer of any warranty.
//
// Alternatively, this file may be used under the terms of Open CASCADE
// commercial license or contractual agreement.

//  Modified by skv - Thu Jul  7 14:37:05 2005 OCC9134

#include <ElCLib.hxx>
#include <ElSLib.hxx>
#include <Extrema_ExtElC.hxx>
#include <Extrema_ExtElCS.hxx>
#include <Extrema_ExtPElC.hxx>
#include <Extrema_ExtPElS.hxx>
#include <Extrema_POnCurv.hxx>
#include <Extrema_POnSurf.hxx>
#include <gp_Circ.hxx>
#include <gp_Cone.hxx>
#include <gp_Cylinder.hxx>
#include <gp_Hypr.hxx>
#include <gp_Lin.hxx>
#include <gp_Pln.hxx>
#include <gp_Sphere.hxx>
#include <gp_Torus.hxx>
#include <gp_Vec.hxx>
#include <IntAna_IntConicQuad.hxx>
#include <IntAna_Quadric.hxx>
#include <Precision.hxx>
#include <Standard_NotImplemented.hxx>
#include <Standard_OutOfRange.hxx>
#include <StdFail_InfiniteSolutions.hxx>
#include <StdFail_NotDone.hxx>

Extrema_ExtElCS::Extrema_ExtElCS() 
{
  myDone = Standard_False;
  myIsPar = Standard_False;
}


Extrema_ExtElCS::Extrema_ExtElCS(const gp_Lin& C,
				 const gp_Pln& S)
{
  Perform(C, S);
}


void Extrema_ExtElCS::Perform(const gp_Lin& C,
			      const gp_Pln& S)
{
  myDone = Standard_True;
  myIsPar = Standard_False;

  if (C.Direction().IsNormal(S.Axis().Direction(), 
			     Precision::Angular())) {
    mySqDist = new TColStd_HArray1OfReal(1, 1);
    mySqDist->SetValue(1, S.SquareDistance(C));
    myIsPar = Standard_True;
  }
  else {
    myNbExt = 0;
  }
  
}


Extrema_ExtElCS::Extrema_ExtElCS(const gp_Lin& C,
				 const gp_Cylinder& S)
{
  Perform(C, S);
}


void Extrema_ExtElCS::Perform(const gp_Lin& C,
                              const gp_Cylinder& S)
{
  myDone = Standard_False;
  myNbExt = 0;
  myIsPar = Standard_False;

  gp_Ax3 Pos = S.Position();
  gp_Pnt Origin = Pos.Location();
  gp_Pnt LineOrig = C.Location();

  Standard_Real radius = S.Radius();
  Extrema_ExtElC Extrem(gp_Lin(Pos.Axis()), C, Precision::Angular());
  if (Extrem.IsParallel()) {
    // Line direction is similar to cylinder axis of rotation.
    mySqDist = new TColStd_HArray1OfReal(1, 1);
    myPoint1 = new Extrema_HArray1OfPOnCurv(1, 1);
    myPoint2 = new Extrema_HArray1OfPOnSurf(1, 1);
    Standard_Real aDist = sqrt(Extrem.SquareDistance(1)) - radius;
    mySqDist->SetValue(1, aDist * aDist);
    Standard_Real u, v, w;
    gp_Vec aVec(LineOrig, Origin);
    gp_Vec aDirVec(C.Direction());
    w = aVec*aDirVec;
    gp_Pnt LinPoint = LineOrig.Translated(w * aDirVec);
    Extrema_POnCurv PonC(w, LinPoint);
    myPoint1->SetValue(1, PonC);
    gp_Pnt CylPoint;
    gp_Vec OrigToLine(Origin, LinPoint);
    if (OrigToLine.Magnitude() <= gp::Resolution())
    {
      u = 0.;
      v = 0.;
      CylPoint = ElSLib::Value(u, v, S);
    }
    else
    {
      OrigToLine.Normalize();
      CylPoint = Origin.Translated(radius * OrigToLine);
      ElSLib::CylinderParameters(Pos, radius, CylPoint, u, v);
    }
    Extrema_POnSurf PonS(u, v, CylPoint);
    myPoint2->SetValue(1, PonS);
    myDone = Standard_True;
    myIsPar = Standard_True;
  }
  else {
    Standard_Integer i, aStartIdx = 0;

    Extrema_POnCurv myPOnC1, myPOnC2;
    Extrem.Points(1, myPOnC1, myPOnC2);
    gp_Pnt PonAxis = myPOnC1.Value();
    gp_Pnt PC = myPOnC2.Value();

    // line intersects the cylinder
    if (radius - PonAxis.Distance(PC) > Precision::PConfusion())
    {
      IntAna_Quadric theQuadric(S);
      IntAna_IntConicQuad Inters(C, theQuadric);
      if (Inters.IsDone())
      {
        myNbExt = Inters.NbPoints();
        aStartIdx = myNbExt;
        if (myNbExt > 0)
        {
          // Not more than 2 additional points from perpendiculars.
          mySqDist = new TColStd_HArray1OfReal(1, myNbExt + 2);
          myPoint1 = new Extrema_HArray1OfPOnCurv(1, myNbExt + 2);
          myPoint2 = new Extrema_HArray1OfPOnSurf(1, myNbExt + 2);
          Standard_Real u, v, w;
          for (i = 1; i <= myNbExt; i++)
          {
            mySqDist->SetValue(i, 0.);
            gp_Pnt P_int = Inters.Point(i);
            w = Inters.ParamOnConic(i);
            Extrema_POnCurv PonC(w, P_int);
            myPoint1->SetValue(i, PonC);
            ElSLib::CylinderParameters(Pos, radius, P_int, u, v);
            Extrema_POnSurf PonS(u, v, P_int);
            myPoint2->SetValue(i, PonS);
          }
        }
      }
    }

    // line is tangent or outside of the cylinder
    Extrema_ExtPElS ExPS(PC, S, Precision::Confusion());
      if (ExPS.IsDone())
      {
        if (aStartIdx == 0)
        {
          myNbExt = ExPS.NbExt();
          mySqDist = new TColStd_HArray1OfReal(1, myNbExt);
          myPoint1 = new Extrema_HArray1OfPOnCurv(1, myNbExt);
          myPoint2 = new Extrema_HArray1OfPOnSurf(1, myNbExt);
        }
        else
          myNbExt += ExPS.NbExt();

        for (i = aStartIdx + 1; i <= myNbExt; i++) {
          myPoint1->SetValue(i, myPOnC2);
          myPoint2->SetValue(i, ExPS.Point(i - aStartIdx));
          mySqDist->SetValue(i,(myPOnC2.Value()).SquareDistance(ExPS.Point(i - aStartIdx).Value()));
        }
      }
    myDone = Standard_True;
  }

}


Extrema_ExtElCS::Extrema_ExtElCS(const gp_Lin& C,
				 const gp_Cone& S)
{  Perform(C, S);}


//void Extrema_ExtElCS::Perform(const gp_Lin& C,
//			      const gp_Cone& S)
void Extrema_ExtElCS::Perform(const gp_Lin& ,
			      const gp_Cone& )
{
  throw Standard_NotImplemented();

}


Extrema_ExtElCS::Extrema_ExtElCS(const gp_Lin& C,
				 const gp_Sphere& S)
{
  Perform(C, S);
}


void Extrema_ExtElCS::Perform(const gp_Lin& C,
                              const gp_Sphere& S)
{
  // In case of intersection - return four points:
  // 2 intersection points and 2 perpendicular.
  // No intersection - only min and max.

  myDone = Standard_False;
  myNbExt = 0;
  myIsPar = Standard_False;
  Standard_Integer aStartIdx = 0;

  gp_Pnt aCenter = S.Location();

  Extrema_ExtPElC Extrem(aCenter, C, Precision::Angular(), RealFirst(), RealLast());

  Standard_Integer i;
  if (Extrem.IsDone() &&
      Extrem.NbExt() > 0)
  {
    Extrema_POnCurv myPOnC1 =  Extrem.Point(1);
    if (myPOnC1.Value().Distance(aCenter) <= S.Radius())
    {
      IntAna_IntConicQuad aLinSphere(C, S);
      if (aLinSphere.IsDone())
      {
        myNbExt = aLinSphere.NbPoints();
        aStartIdx = myNbExt;
        // Not more than 2 additional points from perpendiculars.
        mySqDist = new TColStd_HArray1OfReal(1, myNbExt + 2);
        myPoint1 = new Extrema_HArray1OfPOnCurv(1, myNbExt + 2);
        myPoint2 = new Extrema_HArray1OfPOnSurf(1, myNbExt + 2);

        for (i = 1; i <= myNbExt; i++)
        {
          Extrema_POnCurv aCPnt(aLinSphere.ParamOnConic(i), aLinSphere.Point(i));

          Standard_Real u,v;
          ElSLib::Parameters(S, aLinSphere.Point(i), u, v);
          Extrema_POnSurf aSPnt(u, v, aLinSphere.Point(i));

          myPoint1->SetValue(i, aCPnt);
          myPoint2->SetValue(i, aSPnt);
          mySqDist->SetValue(i,(aCPnt.Value()).SquareDistance(aSPnt.Value()));
        }
      }
    }

    Extrema_ExtPElS ExPS(myPOnC1.Value(), S, Precision::Confusion());
    if (ExPS.IsDone())
    {
      if (aStartIdx == 0)
      {
        myNbExt = ExPS.NbExt();
        mySqDist = new TColStd_HArray1OfReal(1, myNbExt);
        myPoint1 = new Extrema_HArray1OfPOnCurv(1, myNbExt);
        myPoint2 = new Extrema_HArray1OfPOnSurf(1, myNbExt);
      }
      else
        myNbExt += ExPS.NbExt();

      for (i = aStartIdx + 1; i <= myNbExt; i++)
      {
        myPoint1->SetValue(i, myPOnC1);
        myPoint2->SetValue(i, ExPS.Point(i - aStartIdx));
        mySqDist->SetValue(i,(myPOnC1.Value()).SquareDistance(ExPS.Point(i - aStartIdx).Value()));
      }
    }
  }
  myDone = Standard_True;
}


Extrema_ExtElCS::Extrema_ExtElCS(const gp_Lin& C,
				 const gp_Torus& S)
{  Perform(C, S);}


//void Extrema_ExtElCS::Perform(const gp_Lin& C,
//			      const gp_Torus& S)
void Extrema_ExtElCS::Perform(const gp_Lin& ,
			      const gp_Torus& )
{
  throw Standard_NotImplemented();

}


//        Circle-?

Extrema_ExtElCS::Extrema_ExtElCS(const gp_Circ& C,
				 const gp_Pln& S)
{
  Perform(C, S);
}


void Extrema_ExtElCS::Perform(const gp_Circ& C,
                              const gp_Pln& S)
{
  myDone = Standard_True;
  myIsPar = Standard_False;

  gp_Ax2 Pos = C.Position();
  gp_Dir NCirc = Pos.Direction();
  gp_Dir NPln = S.Axis().Direction();

  if (NCirc.IsParallel(NPln, Precision::Angular())) {

    mySqDist = new TColStd_HArray1OfReal(1, 1);
    mySqDist->SetValue(1, S.SquareDistance(C.Location()));
    myIsPar = Standard_True;

  }
  else {

    gp_Dir ExtLine = NCirc ^ NPln;
    ExtLine = ExtLine ^ NCirc;
    //
    gp_Dir XDir = Pos.XDirection();
    Standard_Real T[2];
    T[0] = XDir.AngleWithRef(ExtLine, NCirc);
    if(T[0] < 0.)
    {
      //Put in period
      T[0] += M_PI;
    }
    T[1] = T[0] + M_PI;
    //
    myNbExt = 2;
    //Check intersection 
    IntAna_IntConicQuad anInter(C, S,
                                Precision::Angular(),
                                Precision::Confusion());
    if(anInter.IsDone())
    {
      if(anInter.NbPoints() > 1)
      {
        myNbExt += anInter.NbPoints();
      }
    }

    myPoint1 = new Extrema_HArray1OfPOnCurv(1, myNbExt);
    mySqDist = new TColStd_HArray1OfReal(1, myNbExt);
    myPoint2 = new Extrema_HArray1OfPOnSurf(1, myNbExt);

    Standard_Integer i;
    gp_Pnt PC, PP;
    Standard_Real U, V;
    Extrema_POnCurv POnC;
    Extrema_POnSurf POnS;
    for(i = 0; i < 2; ++i)
    {
      PC = ElCLib::CircleValue(T[i], C.Position(), C.Radius());
      POnC.SetValues(T[i], PC);
      myPoint1->SetValue(i+1, POnC);
      ElSLib::PlaneParameters(S.Position(), PC, U, V);
      PP = ElSLib::PlaneValue(U, V, S.Position());
      POnS.SetParameters(U, V, PP);
      myPoint2->SetValue(i+1, POnS);
      mySqDist->SetValue(i+1, PC.SquareDistance(PP));
    }
    //
    if(myNbExt > 2)
    {
      //Add intersection points
      for(i = 1; i <= anInter.NbPoints(); ++i)
      {
        Standard_Real t = anInter.ParamOnConic(i);
        PC = ElCLib::CircleValue(t, C.Position(), C.Radius());
        POnC.SetValues(t, PC);
        myPoint1->SetValue(i+2, POnC);
        ElSLib::PlaneParameters(S.Position(), PC, U, V);
        PP = ElSLib::PlaneValue(U, V, S.Position());
        POnS.SetParameters(U, V, PP);
        myPoint2->SetValue(i+2, POnS);
        mySqDist->SetValue(i+2, PC.SquareDistance(PP));
      }
    }
  }
  //  
}


Extrema_ExtElCS::Extrema_ExtElCS(const gp_Circ& C,
				 const gp_Cylinder& S)
{
  Perform(C, S);
}


//  Modified by skv - Thu Jul  7 14:37:05 2005 OCC9134 Begin
// Implementation of the method.
void Extrema_ExtElCS::Perform(const gp_Circ& C,
			      const gp_Cylinder& S)
{
  myDone  = Standard_False;
  myIsPar = Standard_False;
  myNbExt = 0;

  // Get an axis line of the cylinder.
  gp_Lin anAxis(S.Axis());

  // Compute extrema between the circle and the line.
  Extrema_ExtElC anExtC(anAxis, C, 0.);

  if (anExtC.IsDone()) {
    if (anExtC.IsParallel()) {
      myIsPar =     Standard_True;
      mySqDist = new TColStd_HArray1OfReal(1, 1);
      Standard_Real aDist = sqrt (anExtC.SquareDistance(1)) - S.Radius();
      mySqDist->SetValue(1, aDist * aDist);
    } else {
      Standard_Integer aNbExt   = anExtC.NbExt();
      Standard_Integer i;
      Standard_Integer aCurI    = 1;
      Standard_Real    aTolConf = Precision::Confusion();
      Standard_Real    aCylRad  = S.Radius();

      // Check whether two objects have intersection points
      IntAna_Quadric aCylQuad(S);
      IntAna_IntConicQuad aCircCylInter(C, aCylQuad);
      Standard_Integer aNbInter = 0;
      if (aCircCylInter.IsDone())
        aNbInter = aCircCylInter.NbPoints();

      // Compute the extremas.
      myNbExt  =     2*aNbExt + aNbInter;
      mySqDist  = new TColStd_HArray1OfReal(1, myNbExt);
      myPoint1 = new Extrema_HArray1OfPOnCurv(1, myNbExt);
      myPoint2 = new Extrema_HArray1OfPOnSurf(1, myNbExt);

      for (i = 1; i <= aNbExt; i++) {
        Extrema_POnCurv aPOnAxis;
        Extrema_POnCurv aPOnCirc;
        Standard_Real   aSqDist = anExtC.SquareDistance(i);
        Standard_Real   aDist = sqrt (aSqDist);

        anExtC.Points(i, aPOnAxis, aPOnCirc);

        if (aSqDist <= (aTolConf * aTolConf)) {
          myNbExt -= 2;
          continue;
        }

        gp_Dir aDir(aPOnAxis.Value().XYZ().Subtracted(aPOnCirc.Value().XYZ()));
        Standard_Real aShift[2] = { aDist + aCylRad, aDist - aCylRad };
        Standard_Integer j;

        for (j = 0; j < 2; j++) {
          gp_Vec aVec(aDir);
          gp_Pnt aPntOnCyl;

          aVec.Multiply(aShift[j]);
          aPntOnCyl = aPOnCirc.Value().Translated(aVec);

          Standard_Real aU;
          Standard_Real aV;

          ElSLib::Parameters(S, aPntOnCyl, aU, aV);

          Extrema_POnSurf aPOnSurf(aU, aV, aPntOnCyl);

          myPoint1->SetValue(aCurI, aPOnCirc);
          myPoint2->SetValue(aCurI, aPOnSurf);
          mySqDist->SetValue(aCurI++, aShift[j] * aShift[j]);
        }
      }

      // Adding intersection points to the list of extremas
      for (i=1; i<=aNbInter; i++)
      {
        Standard_Real aU;
        Standard_Real aV;

        gp_Pnt aInterPnt = aCircCylInter.Point(i);

        aU = ElCLib::Parameter(C, aInterPnt);
        Extrema_POnCurv aPOnCirc(aU, aInterPnt);

        ElSLib::Parameters(S, aInterPnt, aU, aV);
        Extrema_POnSurf aPOnCyl(aU, aV, aInterPnt);
        myPoint1->SetValue(aCurI, aPOnCirc);
        myPoint2->SetValue(aCurI, aPOnCyl);
        mySqDist->SetValue(aCurI++, 0.0);
      }
    }

    myDone = Standard_True;
  }
}
//  Modified by skv - Thu Jul  7 14:37:05 2005 OCC9134 End


Extrema_ExtElCS::Extrema_ExtElCS(const gp_Circ& C,
				 const gp_Cone& S)
{  Perform(C, S);}


//void Extrema_ExtElCS::Perform(const gp_Circ& C,
//			 const gp_Cone& S)
void Extrema_ExtElCS::Perform(const gp_Circ& ,
			 const gp_Cone& )
{
  throw Standard_NotImplemented();

}


Extrema_ExtElCS::Extrema_ExtElCS(const gp_Circ& C,
				 const gp_Sphere& S)
{  Perform(C, S);}


//void Extrema_ExtElCS::Perform(const gp_Circ& C,
//			      const gp_Sphere& S)
void Extrema_ExtElCS::Perform(const gp_Circ& ,
			      const gp_Sphere& )
{
  throw Standard_NotImplemented();

}

Extrema_ExtElCS::Extrema_ExtElCS(const gp_Circ& C,
				 const gp_Torus& S)
{  Perform(C, S);}


//void Extrema_ExtElCS::Perform(const gp_Circ& C,
//			      const gp_Torus& S)
void Extrema_ExtElCS::Perform(const gp_Circ& ,
			      const gp_Torus& )
{
  throw Standard_NotImplemented();

}

Extrema_ExtElCS::Extrema_ExtElCS(const gp_Hypr& C,
				 const gp_Pln& S)
{
  Perform(C, S);
}


void Extrema_ExtElCS::Perform(const gp_Hypr& C,
			      const gp_Pln& S)
{
  myDone = Standard_True;
  myIsPar = Standard_False;

  gp_Ax2 Pos = C.Position();
  gp_Dir NHypr = Pos.Direction();
  gp_Dir NPln = S.Axis().Direction();

  if (NHypr.IsParallel(NPln, Precision::Angular())) {

    mySqDist = new TColStd_HArray1OfReal(1, 1);
    mySqDist->SetValue(1, S.SquareDistance(C.Location()));
    myIsPar = Standard_True;

  }
  else {

    gp_Dir XDir = Pos.XDirection();
    gp_Dir YDir = Pos.YDirection();

    Standard_Real A = C.MinorRadius()*(NPln.Dot(YDir)); 
    Standard_Real B = C.MajorRadius()*(NPln.Dot(XDir)); 

    if(Abs(B) > Abs(A)) {
      Standard_Real T = -0.5 * Log((A+B)/(B-A));
      gp_Pnt Ph = ElCLib::HyperbolaValue(T, Pos, C.MajorRadius(), C.MinorRadius());
      Extrema_POnCurv PC(T, Ph);
      myPoint1 = new Extrema_HArray1OfPOnCurv(1,1);
      myPoint1->SetValue(1, PC);

      mySqDist = new TColStd_HArray1OfReal(1, 1);
      mySqDist->SetValue(1, S.SquareDistance(Ph));

      Standard_Real U, V;
      ElSLib::PlaneParameters(S.Position(), Ph, U, V);
      gp_Pnt Pp = ElSLib::PlaneValue(U, V, S.Position());
      Extrema_POnSurf PS(U, V, Pp);
      myPoint2 = new Extrema_HArray1OfPOnSurf(1,1);
      myPoint2->SetValue(1, PS);

      myNbExt = 1;
    }
    else {
      myNbExt = 0;
    }

  }
  
}


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


Standard_Integer Extrema_ExtElCS::NbExt() const
{
  if (myIsPar) throw StdFail_InfiniteSolutions();
  return myNbExt;
}

Standard_Real Extrema_ExtElCS::SquareDistance(const Standard_Integer N) const
{
  if (myIsPar && N != 1) throw StdFail_InfiniteSolutions();
  return mySqDist->Value(N);
}


void Extrema_ExtElCS::Points(const Standard_Integer N,
			     Extrema_POnCurv& P1,
			     Extrema_POnSurf& P2) const
{
  if (myIsPar) throw StdFail_InfiniteSolutions();
  P1 = myPoint1->Value(N);
  P2 = myPoint2->Value(N);
}


Standard_Boolean Extrema_ExtElCS::IsParallel() const
{
  return myIsPar;
}
Commit	Line	Data
b311480e	1	// Copyright (c) 1995-1999 Matra Datavision
973c2be1	2	// Copyright (c) 1999-2014 OPEN CASCADE SAS
b311480e	3	//
973c2be1	4	// This file is part of Open CASCADE Technology software library.
b311480e	5	//
d5f74e42	6	// This library is free software; you can redistribute it and/or modify it under
d5f74e42	7	// the terms of the GNU Lesser General Public License version 2.1 as published
973c2be1	8	// by the Free Software Foundation, with special exception defined in the file
	9	// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
	10	// distribution for complete text of the license and disclaimer of any warranty.
b311480e	11	//
973c2be1	12	// Alternatively, this file may be used under the terms of Open CASCADE
973c2be1	13	// commercial license or contractual agreement.
7fd59977	14
	15	// Modified by skv - Thu Jul 7 14:37:05 2005 OCC9134
	16
42cf5bc1	17	#include <ElCLib.hxx>
42cf5bc1	18	#include <ElSLib.hxx>
7fd59977	19	#include <Extrema_ExtElC.hxx>
42cf5bc1	20	#include <Extrema_ExtElCS.hxx>
	21	#include <Extrema_ExtPElC.hxx>
	22	#include <Extrema_ExtPElS.hxx>
7fd59977	23	#include <Extrema_POnCurv.hxx>
42cf5bc1	24	#include <Extrema_POnSurf.hxx>
	25	#include <gp_Circ.hxx>
	26	#include <gp_Cone.hxx>
	27	#include <gp_Cylinder.hxx>
	28	#include <gp_Hypr.hxx>
	29	#include <gp_Lin.hxx>
	30	#include <gp_Pln.hxx>
	31	#include <gp_Sphere.hxx>
	32	#include <gp_Torus.hxx>
7fd59977	33	#include <gp_Vec.hxx>
f34cd0d1	34	#include <IntAna_IntConicQuad.hxx>
42cf5bc1	35	#include <IntAna_Quadric.hxx>
	36	#include <Precision.hxx>
	37	#include <Standard_NotImplemented.hxx>
	38	#include <Standard_OutOfRange.hxx>
	39	#include <StdFail_InfiniteSolutions.hxx>
	40	#include <StdFail_NotDone.hxx>
7fd59977	41
	42	Extrema_ExtElCS::Extrema_ExtElCS()
	43	{
	44	myDone = Standard_False;
	45	myIsPar = Standard_False;
	46	}
	47
	48
	49	Extrema_ExtElCS::Extrema_ExtElCS(const gp_Lin& C,
	50	const gp_Pln& S)
	51	{
	52	Perform(C, S);
	53	}
	54
	55
	56
	57	void Extrema_ExtElCS::Perform(const gp_Lin& C,
	58	const gp_Pln& S)
	59	{
	60	myDone = Standard_True;
	61	myIsPar = Standard_False;
	62
	63	if (C.Direction().IsNormal(S.Axis().Direction(),
	64	Precision::Angular())) {
	65	mySqDist = new TColStd_HArray1OfReal(1, 1);
	66	mySqDist->SetValue(1, S.SquareDistance(C));
	67	myIsPar = Standard_True;
	68	}
	69	else {
	70	myNbExt = 0;
	71	}
	72
	73	}
	74
	75
	76	Extrema_ExtElCS::Extrema_ExtElCS(const gp_Lin& C,
	77	const gp_Cylinder& S)
	78	{
	79	Perform(C, S);
	80	}
	81
	82
	83
	84	void Extrema_ExtElCS::Perform(const gp_Lin& C,
db914841	85	const gp_Cylinder& S)
7fd59977	86	{
	87	myDone = Standard_False;
	88	myNbExt = 0;
	89	myIsPar = Standard_False;
	90
	91	gp_Ax3 Pos = S.Position();
f34cd0d1	92	gp_Pnt Origin = Pos.Location();
	93	gp_Pnt LineOrig = C.Location();
	94
7fd59977	95	Standard_Real radius = S.Radius();
	96	Extrema_ExtElC Extrem(gp_Lin(Pos.Axis()), C, Precision::Angular());
	97	if (Extrem.IsParallel()) {
db914841	98	// Line direction is similar to cylinder axis of rotation.
7fd59977	99	mySqDist = new TColStd_HArray1OfReal(1, 1);
f34cd0d1	100	myPoint1 = new Extrema_HArray1OfPOnCurv(1, 1);
	101	myPoint2 = new Extrema_HArray1OfPOnSurf(1, 1);
	102	Standard_Real aDist = sqrt(Extrem.SquareDistance(1)) - radius;
7fd59977	103	mySqDist->SetValue(1, aDist * aDist);
f34cd0d1	104	Standard_Real u, v, w;
	105	gp_Vec aVec(LineOrig, Origin);
	106	gp_Vec aDirVec(C.Direction());
	107	w = aVec*aDirVec;
	108	gp_Pnt LinPoint = LineOrig.Translated(w * aDirVec);
	109	Extrema_POnCurv PonC(w, LinPoint);
	110	myPoint1->SetValue(1, PonC);
	111	gp_Pnt CylPoint;
	112	gp_Vec OrigToLine(Origin, LinPoint);
	113	if (OrigToLine.Magnitude() <= gp::Resolution())
	114	{
	115	u = 0.;
	116	v = 0.;
	117	CylPoint = ElSLib::Value(u, v, S);
	118	}
	119	else
	120	{
	121	OrigToLine.Normalize();
	122	CylPoint = Origin.Translated(radius * OrigToLine);
	123	ElSLib::CylinderParameters(Pos, radius, CylPoint, u, v);
	124	}
	125	Extrema_POnSurf PonS(u, v, CylPoint);
	126	myPoint2->SetValue(1, PonS);
7fd59977	127	myDone = Standard_True;
	128	myIsPar = Standard_True;
	129	}
	130	else {
db914841	131	Standard_Integer i, aStartIdx = 0;
db914841	132
f34cd0d1	133	Extrema_POnCurv myPOnC1, myPOnC2;
	134	Extrem.Points(1, myPOnC1, myPOnC2);
	135	gp_Pnt PonAxis = myPOnC1.Value();
	136	gp_Pnt PC = myPOnC2.Value();
	137
f34cd0d1	138	// line intersects the cylinder
db914841	139	if (radius - PonAxis.Distance(PC) > Precision::PConfusion())
f34cd0d1	140	{
	141	IntAna_Quadric theQuadric(S);
	142	IntAna_IntConicQuad Inters(C, theQuadric);
	143	if (Inters.IsDone())
	144	{
	145	myNbExt = Inters.NbPoints();
db914841	146	aStartIdx = myNbExt;
f34cd0d1	147	if (myNbExt > 0)
f34cd0d1	148	{
db914841	149	// Not more than 2 additional points from perpendiculars.
	150	mySqDist = new TColStd_HArray1OfReal(1, myNbExt + 2);
	151	myPoint1 = new Extrema_HArray1OfPOnCurv(1, myNbExt + 2);
	152	myPoint2 = new Extrema_HArray1OfPOnSurf(1, myNbExt + 2);
f34cd0d1	153	Standard_Real u, v, w;
	154	for (i = 1; i <= myNbExt; i++)
	155	{
	156	mySqDist->SetValue(i, 0.);
	157	gp_Pnt P_int = Inters.Point(i);
	158	w = Inters.ParamOnConic(i);
	159	Extrema_POnCurv PonC(w, P_int);
	160	myPoint1->SetValue(i, PonC);
	161	ElSLib::CylinderParameters(Pos, radius, P_int, u, v);
	162	Extrema_POnSurf PonS(u, v, P_int);
	163	myPoint2->SetValue(i, PonS);
	164	}
	165	}
7fd59977	166	}
7fd59977	167	}
db914841	168
	169	// line is tangent or outside of the cylinder
	170	Extrema_ExtPElS ExPS(PC, S, Precision::Confusion());
	171	if (ExPS.IsDone())
	172	{
	173	if (aStartIdx == 0)
	174	{
	175	myNbExt = ExPS.NbExt();
	176	mySqDist = new TColStd_HArray1OfReal(1, myNbExt);
	177	myPoint1 = new Extrema_HArray1OfPOnCurv(1, myNbExt);
	178	myPoint2 = new Extrema_HArray1OfPOnSurf(1, myNbExt);
	179	}
	180	else
	181	myNbExt += ExPS.NbExt();
	182
	183	for (i = aStartIdx + 1; i <= myNbExt; i++) {
	184	myPoint1->SetValue(i, myPOnC2);
	185	myPoint2->SetValue(i, ExPS.Point(i - aStartIdx));
	186	mySqDist->SetValue(i,(myPOnC2.Value()).SquareDistance(ExPS.Point(i - aStartIdx).Value()));
	187	}
	188	}
f34cd0d1	189	myDone = Standard_True;
7fd59977	190	}
	191
	192	}
	193
	194
	195
	196	Extrema_ExtElCS::Extrema_ExtElCS(const gp_Lin& C,
	197	const gp_Cone& S)
9775fa61	198	{ Perform(C, S);}
7fd59977	199
	200
	201
	202	//void Extrema_ExtElCS::Perform(const gp_Lin& C,
	203	// const gp_Cone& S)
	204	void Extrema_ExtElCS::Perform(const gp_Lin& ,
	205	const gp_Cone& )
	206	{
9775fa61	207	throw Standard_NotImplemented();
7fd59977	208
	209	}
	210
	211
	212
	213	Extrema_ExtElCS::Extrema_ExtElCS(const gp_Lin& C,
	214	const gp_Sphere& S)
	215	{
	216	Perform(C, S);
	217	}
	218
	219
	220
	221	void Extrema_ExtElCS::Perform(const gp_Lin& C,
db914841	222	const gp_Sphere& S)
7fd59977	223	{
db914841	224	// In case of intersection - return four points:
	225	// 2 intersection points and 2 perpendicular.
	226	// No intersection - only min and max.
	227
7fd59977	228	myDone = Standard_False;
	229	myNbExt = 0;
	230	myIsPar = Standard_False;
db914841	231	Standard_Integer aStartIdx = 0;
7fd59977	232
db914841	233	gp_Pnt aCenter = S.Location();
7fd59977	234
db914841	235	Extrema_ExtPElC Extrem(aCenter, C, Precision::Angular(), RealFirst(), RealLast());
7fd59977	236
7fd59977	237	Standard_Integer i;
db914841	238	if (Extrem.IsDone() &&
	239	Extrem.NbExt() > 0)
	240	{
7fd59977	241	Extrema_POnCurv myPOnC1 = Extrem.Point(1);
db914841	242	if (myPOnC1.Value().Distance(aCenter) <= S.Radius())
	243	{
	244	IntAna_IntConicQuad aLinSphere(C, S);
	245	if (aLinSphere.IsDone())
	246	{
	247	myNbExt = aLinSphere.NbPoints();
	248	aStartIdx = myNbExt;
	249	// Not more than 2 additional points from perpendiculars.
	250	mySqDist = new TColStd_HArray1OfReal(1, myNbExt + 2);
	251	myPoint1 = new Extrema_HArray1OfPOnCurv(1, myNbExt + 2);
	252	myPoint2 = new Extrema_HArray1OfPOnSurf(1, myNbExt + 2);
	253
	254	for (i = 1; i <= myNbExt; i++)
	255	{
	256	Extrema_POnCurv aCPnt(aLinSphere.ParamOnConic(i), aLinSphere.Point(i));
	257
	258	Standard_Real u,v;
	259	ElSLib::Parameters(S, aLinSphere.Point(i), u, v);
	260	Extrema_POnSurf aSPnt(u, v, aLinSphere.Point(i));
	261
	262	myPoint1->SetValue(i, aCPnt);
	263	myPoint2->SetValue(i, aSPnt);
	264	mySqDist->SetValue(i,(aCPnt.Value()).SquareDistance(aSPnt.Value()));
	265	}
	266	}
	267	}
	268
7fd59977	269	Extrema_ExtPElS ExPS(myPOnC1.Value(), S, Precision::Confusion());
db914841	270	if (ExPS.IsDone())
	271	{
	272	if (aStartIdx == 0)
	273	{
	274	myNbExt = ExPS.NbExt();
	275	mySqDist = new TColStd_HArray1OfReal(1, myNbExt);
	276	myPoint1 = new Extrema_HArray1OfPOnCurv(1, myNbExt);
	277	myPoint2 = new Extrema_HArray1OfPOnSurf(1, myNbExt);
	278	}
	279	else
	280	myNbExt += ExPS.NbExt();
	281
	282	for (i = aStartIdx + 1; i <= myNbExt; i++)
	283	{
	284	myPoint1->SetValue(i, myPOnC1);
	285	myPoint2->SetValue(i, ExPS.Point(i - aStartIdx));
	286	mySqDist->SetValue(i,(myPOnC1.Value()).SquareDistance(ExPS.Point(i - aStartIdx).Value()));
7fd59977	287	}
	288	}
	289	}
db914841	290	myDone = Standard_True;
7fd59977	291	}
	292
	293
	294	Extrema_ExtElCS::Extrema_ExtElCS(const gp_Lin& C,
	295	const gp_Torus& S)
9775fa61	296	{ Perform(C, S);}
7fd59977	297
	298
	299
	300	//void Extrema_ExtElCS::Perform(const gp_Lin& C,
	301	// const gp_Torus& S)
	302	void Extrema_ExtElCS::Perform(const gp_Lin& ,
	303	const gp_Torus& )
	304	{
9775fa61	305	throw Standard_NotImplemented();
7fd59977	306
	307	}
	308
	309
	310	// Circle-?
	311
	312	Extrema_ExtElCS::Extrema_ExtElCS(const gp_Circ& C,
	313	const gp_Pln& S)
	314	{
	315	Perform(C, S);
	316	}
	317
	318
	319
9dfbbfe6	320	void Extrema_ExtElCS::Perform(const gp_Circ& C,
9dfbbfe6	321	const gp_Pln& S)
7fd59977	322	{
9dfbbfe6	323	myDone = Standard_True;
	324	myIsPar = Standard_False;
	325
	326	gp_Ax2 Pos = C.Position();
	327	gp_Dir NCirc = Pos.Direction();
	328	gp_Dir NPln = S.Axis().Direction();
	329
	330	if (NCirc.IsParallel(NPln, Precision::Angular())) {
	331
	332	mySqDist = new TColStd_HArray1OfReal(1, 1);
	333	mySqDist->SetValue(1, S.SquareDistance(C.Location()));
	334	myIsPar = Standard_True;
	335
	336	}
	337	else {
	338
	339	gp_Dir ExtLine = NCirc ^ NPln;
	340	ExtLine = ExtLine ^ NCirc;
	341	//
	342	gp_Dir XDir = Pos.XDirection();
	343	Standard_Real T[2];
	344	T[0] = XDir.AngleWithRef(ExtLine, NCirc);
	345	if(T[0] < 0.)
	346	{
	347	//Put in period
	348	T[0] += M_PI;
	349	}
	350	T[1] = T[0] + M_PI;
	351	//
	352	myNbExt = 2;
	353	//Check intersection
	354	IntAna_IntConicQuad anInter(C, S,
	355	Precision::Angular(),
	356	Precision::Confusion());
	357	if(anInter.IsDone())
	358	{
	359	if(anInter.NbPoints() > 1)
	360	{
	361	myNbExt += anInter.NbPoints();
	362	}
	363	}
	364
	365	myPoint1 = new Extrema_HArray1OfPOnCurv(1, myNbExt);
	366	mySqDist = new TColStd_HArray1OfReal(1, myNbExt);
	367	myPoint2 = new Extrema_HArray1OfPOnSurf(1, myNbExt);
7fd59977	368
9dfbbfe6	369	Standard_Integer i;
	370	gp_Pnt PC, PP;
	371	Standard_Real U, V;
	372	Extrema_POnCurv POnC;
	373	Extrema_POnSurf POnS;
	374	for(i = 0; i < 2; ++i)
	375	{
	376	PC = ElCLib::CircleValue(T[i], C.Position(), C.Radius());
	377	POnC.SetValues(T[i], PC);
	378	myPoint1->SetValue(i+1, POnC);
	379	ElSLib::PlaneParameters(S.Position(), PC, U, V);
	380	PP = ElSLib::PlaneValue(U, V, S.Position());
	381	POnS.SetParameters(U, V, PP);
	382	myPoint2->SetValue(i+1, POnS);
	383	mySqDist->SetValue(i+1, PC.SquareDistance(PP));
	384	}
	385	//
	386	if(myNbExt > 2)
	387	{
	388	//Add intersection points
	389	for(i = 1; i <= anInter.NbPoints(); ++i)
	390	{
	391	Standard_Real t = anInter.ParamOnConic(i);
	392	PC = ElCLib::CircleValue(t, C.Position(), C.Radius());
	393	POnC.SetValues(t, PC);
	394	myPoint1->SetValue(i+2, POnC);
	395	ElSLib::PlaneParameters(S.Position(), PC, U, V);
	396	PP = ElSLib::PlaneValue(U, V, S.Position());
	397	POnS.SetParameters(U, V, PP);
	398	myPoint2->SetValue(i+2, POnS);
	399	mySqDist->SetValue(i+2, PC.SquareDistance(PP));
	400	}
	401	}
	402	}
	403	//
7fd59977	404	}
	405
	406
	407
	408	Extrema_ExtElCS::Extrema_ExtElCS(const gp_Circ& C,
	409	const gp_Cylinder& S)
	410	{
	411	Perform(C, S);
	412	}
	413
	414
	415
	416	// Modified by skv - Thu Jul 7 14:37:05 2005 OCC9134 Begin
	417	// Implementation of the method.
	418	void Extrema_ExtElCS::Perform(const gp_Circ& C,
	419	const gp_Cylinder& S)
	420	{
	421	myDone = Standard_False;
	422	myIsPar = Standard_False;
	423	myNbExt = 0;
	424
	425	// Get an axis line of the cylinder.
	426	gp_Lin anAxis(S.Axis());
	427
	428	// Compute extrema between the circle and the line.
	429	Extrema_ExtElC anExtC(anAxis, C, 0.);
	430
	431	if (anExtC.IsDone()) {
	432	if (anExtC.IsParallel()) {
	433	myIsPar = Standard_True;
	434	mySqDist = new TColStd_HArray1OfReal(1, 1);
	435	Standard_Real aDist = sqrt (anExtC.SquareDistance(1)) - S.Radius();
	436	mySqDist->SetValue(1, aDist * aDist);
	437	} else {
	438	Standard_Integer aNbExt = anExtC.NbExt();
7fd59977	439	Standard_Integer i;
	440	Standard_Integer aCurI = 1;
	441	Standard_Real aTolConf = Precision::Confusion();
	442	Standard_Real aCylRad = S.Radius();
	443
19875353	444	// Check whether two objects have intersection points
	445	IntAna_Quadric aCylQuad(S);
	446	IntAna_IntConicQuad aCircCylInter(C, aCylQuad);
fc0180e5	447	Standard_Integer aNbInter = 0;
	448	if (aCircCylInter.IsDone())
	449	aNbInter = aCircCylInter.NbPoints();
19875353	450
7fd59977	451	// Compute the extremas.
19875353	452	myNbExt = 2*aNbExt + aNbInter;
7fd59977	453	mySqDist = new TColStd_HArray1OfReal(1, myNbExt);
	454	myPoint1 = new Extrema_HArray1OfPOnCurv(1, myNbExt);
	455	myPoint2 = new Extrema_HArray1OfPOnSurf(1, myNbExt);
	456
	457	for (i = 1; i <= aNbExt; i++) {
19875353	458	Extrema_POnCurv aPOnAxis;
	459	Extrema_POnCurv aPOnCirc;
	460	Standard_Real aSqDist = anExtC.SquareDistance(i);
	461	Standard_Real aDist = sqrt (aSqDist);
	462
	463	anExtC.Points(i, aPOnAxis, aPOnCirc);
	464
	465	if (aSqDist <= (aTolConf * aTolConf)) {
	466	myNbExt -= 2;
	467	continue;
	468	}
	469
	470	gp_Dir aDir(aPOnAxis.Value().XYZ().Subtracted(aPOnCirc.Value().XYZ()));
	471	Standard_Real aShift[2] = { aDist + aCylRad, aDist - aCylRad };
	472	Standard_Integer j;
7fd59977	473
19875353	474	for (j = 0; j < 2; j++) {
	475	gp_Vec aVec(aDir);
	476	gp_Pnt aPntOnCyl;
7fd59977	477
19875353	478	aVec.Multiply(aShift[j]);
19875353	479	aPntOnCyl = aPOnCirc.Value().Translated(aVec);
7fd59977	480
19875353	481	Standard_Real aU;
19875353	482	Standard_Real aV;
7fd59977	483
19875353	484	ElSLib::Parameters(S, aPntOnCyl, aU, aV);
7fd59977	485
19875353	486	Extrema_POnSurf aPOnSurf(aU, aV, aPntOnCyl);
7fd59977	487
19875353	488	myPoint1->SetValue(aCurI, aPOnCirc);
	489	myPoint2->SetValue(aCurI, aPOnSurf);
	490	mySqDist->SetValue(aCurI++, aShift[j] * aShift[j]);
	491	}
	492	}
	493
	494	// Adding intersection points to the list of extremas
	495	for (i=1; i<=aNbInter; i++)
	496	{
	497	Standard_Real aU;
	498	Standard_Real aV;
7fd59977	499
19875353	500	gp_Pnt aInterPnt = aCircCylInter.Point(i);
7fd59977	501
19875353	502	aU = ElCLib::Parameter(C, aInterPnt);
19875353	503	Extrema_POnCurv aPOnCirc(aU, aInterPnt);
7fd59977	504
19875353	505	ElSLib::Parameters(S, aInterPnt, aU, aV);
	506	Extrema_POnSurf aPOnCyl(aU, aV, aInterPnt);
	507	myPoint1->SetValue(aCurI, aPOnCirc);
	508	myPoint2->SetValue(aCurI, aPOnCyl);
	509	mySqDist->SetValue(aCurI++, 0.0);
7fd59977	510	}
	511	}
	512
	513	myDone = Standard_True;
	514	}
	515	}
	516	// Modified by skv - Thu Jul 7 14:37:05 2005 OCC9134 End
	517
	518
	519
	520	Extrema_ExtElCS::Extrema_ExtElCS(const gp_Circ& C,
	521	const gp_Cone& S)
9775fa61	522	{ Perform(C, S);}
7fd59977	523
	524
	525
	526	//void Extrema_ExtElCS::Perform(const gp_Circ& C,
	527	// const gp_Cone& S)
	528	void Extrema_ExtElCS::Perform(const gp_Circ& ,
	529	const gp_Cone& )
	530	{
9775fa61	531	throw Standard_NotImplemented();
7fd59977	532
	533	}
	534
	535
	536
	537	Extrema_ExtElCS::Extrema_ExtElCS(const gp_Circ& C,
	538	const gp_Sphere& S)
9775fa61	539	{ Perform(C, S);}
7fd59977	540
	541
	542
	543	//void Extrema_ExtElCS::Perform(const gp_Circ& C,
	544	// const gp_Sphere& S)
	545	void Extrema_ExtElCS::Perform(const gp_Circ& ,
	546	const gp_Sphere& )
	547	{
9775fa61	548	throw Standard_NotImplemented();
7fd59977	549
	550	}
	551
	552	Extrema_ExtElCS::Extrema_ExtElCS(const gp_Circ& C,
	553	const gp_Torus& S)
9775fa61	554	{ Perform(C, S);}
7fd59977	555
	556
	557
	558	//void Extrema_ExtElCS::Perform(const gp_Circ& C,
	559	// const gp_Torus& S)
	560	void Extrema_ExtElCS::Perform(const gp_Circ& ,
	561	const gp_Torus& )
	562	{
9775fa61	563	throw Standard_NotImplemented();
7fd59977	564
	565	}
	566
	567	Extrema_ExtElCS::Extrema_ExtElCS(const gp_Hypr& C,
	568	const gp_Pln& S)
	569	{
	570	Perform(C, S);
	571	}
	572
	573
	574
	575	void Extrema_ExtElCS::Perform(const gp_Hypr& C,
	576	const gp_Pln& S)
	577	{
	578	myDone = Standard_True;
	579	myIsPar = Standard_False;
	580
	581	gp_Ax2 Pos = C.Position();
	582	gp_Dir NHypr = Pos.Direction();
	583	gp_Dir NPln = S.Axis().Direction();
	584
	585	if (NHypr.IsParallel(NPln, Precision::Angular())) {
	586
	587	mySqDist = new TColStd_HArray1OfReal(1, 1);
	588	mySqDist->SetValue(1, S.SquareDistance(C.Location()));
	589	myIsPar = Standard_True;
	590
	591	}
	592	else {
	593
	594	gp_Dir XDir = Pos.XDirection();
	595	gp_Dir YDir = Pos.YDirection();
	596
	597	Standard_Real A = C.MinorRadius()*(NPln.Dot(YDir));
	598	Standard_Real B = C.MajorRadius()*(NPln.Dot(XDir));
	599
	600	if(Abs(B) > Abs(A)) {
	601	Standard_Real T = -0.5 * Log((A+B)/(B-A));
	602	gp_Pnt Ph = ElCLib::HyperbolaValue(T, Pos, C.MajorRadius(), C.MinorRadius());
	603	Extrema_POnCurv PC(T, Ph);
	604	myPoint1 = new Extrema_HArray1OfPOnCurv(1,1);
	605	myPoint1->SetValue(1, PC);
	606
	607	mySqDist = new TColStd_HArray1OfReal(1, 1);
	608	mySqDist->SetValue(1, S.SquareDistance(Ph));
	609
	610	Standard_Real U, V;
	611	ElSLib::PlaneParameters(S.Position(), Ph, U, V);
	612	gp_Pnt Pp = ElSLib::PlaneValue(U, V, S.Position());
	613	Extrema_POnSurf PS(U, V, Pp);
	614	myPoint2 = new Extrema_HArray1OfPOnSurf(1,1);
	615	myPoint2->SetValue(1, PS);
	616
	617	myNbExt = 1;
	618	}
	619	else {
	620	myNbExt = 0;
	621	}
	622
	623	}
	624
	625	}
	626
	627
628	Standard_Boolean Extrema_ExtElCS::IsDone() const
629	{
630	return myDone;
631	}
632
633
634	Standard_Integer Extrema_ExtElCS::NbExt() const
635	{
9775fa61	636	if (myIsPar) throw StdFail_InfiniteSolutions();
7fd59977	637	return myNbExt;
	638	}
	639
	640	Standard_Real Extrema_ExtElCS::SquareDistance(const Standard_Integer N) const
	641	{
9775fa61	642	if (myIsPar && N != 1) throw StdFail_InfiniteSolutions();
7fd59977	643	return mySqDist->Value(N);
	644	}
	645
	646
	647	void Extrema_ExtElCS::Points(const Standard_Integer N,
	648	Extrema_POnCurv& P1,
	649	Extrema_POnSurf& P2) const
	650	{
9775fa61	651	if (myIsPar) throw StdFail_InfiniteSolutions();
7fd59977	652	P1 = myPoint1->Value(N);
	653	P2 = myPoint2->Value(N);
	654	}
	655
	656
	657	Standard_Boolean Extrema_ExtElCS::IsParallel() const
	658	{
	659	return myIsPar;
	660	}