// 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.

#include <algorithm>
#include <memory>
#include <TopoDS_Edge.hxx>
#include <Geom_Curve.hxx>
#include <BRepAdaptor_Curve.hxx>
#include <Adaptor3d_Surface.hxx>
#include <Adaptor3d_CurveOnSurface.hxx>
#include <Adaptor3d_CurveOnSurface.hxx>
#include <GeomAbs_SurfaceType.hxx>
#include <BRep_Tool.hxx>
#include <Geom_Line.hxx>
#include <Geom_Plane.hxx>
#include <Geom_CylindricalSurface.hxx>
#include <Geom_ConicalSurface.hxx>
#include <Geom_SphericalSurface.hxx>
#include <Geom_ToroidalSurface.hxx>
#include <gp_Lin.hxx>
#include <gp_Vec.hxx>
#include <gp_Dir.hxx>
#include <gp_Cylinder.hxx>
#include <gp_Ax1.hxx>
#include <gp_Lin.hxx>

#include <GeomAdaptor_Curve.hxx>
#include <GeomAdaptor_Surface.hxx>
#include <Precision.hxx>
#include <Extrema_ExtCC.hxx>
//#include <Extrema_ExtCS.hxx>
#include <Extrema_POnCurv.hxx>
#include <IntCurveSurface_HInter.hxx>

#include <math_FunctionSample.hxx>
#include <math_FunctionAllRoots.hxx>
#include <TColgp_SequenceOfPnt.hxx>

//  Modified by skv - Tue Aug 31 12:13:51 2004 OCC569

#include <Precision.hxx>
#include <IntSurf_Quadric.hxx>
#include <math_Function.hxx>
#include <math_BrentMinimum.hxx>
#include <math_Matrix.hxx>
#include <math_Vector.hxx>
#include <NCollection_Array1.hxx>

#ifdef OCCT_DEBUG
#include <Geom_Circle.hxx>
#include <Geom_Ellipse.hxx>
#include <Geom_Hyperbola.hxx>
#include <Geom_Parabola.hxx>
#include <Geom_BezierCurve.hxx>
#include <Geom_BSplineCurve.hxx>
#include <GeomLib.hxx>
#endif


static Standard_Boolean IsDegenerated(const Handle(Adaptor3d_CurveOnSurface)& theCurve);
static Standard_Boolean IsDegenerated(const IntSurf_Quadric& theQuadric);

static void FindVertex (const TheArc&,
                        const Handle(TheTopolTool)&,
                        TheFunction&,
                        IntStart_SequenceOfPathPoint&,
                        const Standard_Real);

                        
static void BoundedArc (const TheArc& A,
                        const Handle(TheTopolTool)& Domain,
                        const Standard_Real Pdeb,
                        const Standard_Real Pfin,
                        TheFunction& Func,
                        IntStart_SequenceOfPathPoint& pnt,
                        IntStart_SequenceOfSegment& seg,
                        const Standard_Real TolBoundary,
                        const Standard_Real TolTangency,
                        Standard_Boolean& Arcsol,
                        const Standard_Boolean RecheckOnRegularity);
                 
static void PointProcess (const gp_Pnt&,
                          const Standard_Real,
                          const TheArc&,
                          const Handle(TheTopolTool)&,
                          IntStart_SequenceOfPathPoint&,
                          const Standard_Real,
                          Standard_Integer&);

static Standard_Integer TreatLC (const TheArc& A,
                                 const Handle(TheTopolTool)& aDomain,
                                 const IntSurf_Quadric& aQuadric,
                                 const Standard_Real TolBoundary,
                                 IntStart_SequenceOfPathPoint& pnt);

static Standard_Boolean IsRegularity(const TheArc& A,
                                     const Handle(TheTopolTool)& aDomain);

class MinFunction : public math_Function
{
public:
  MinFunction(TheFunction &theFunc) : myFunc(&theFunc) {};

  //returns value of the one-dimension-function when parameter
  //is equal to theX
  virtual Standard_Boolean Value(const Standard_Real theX,
                                 Standard_Real& theFVal)
  {
    if(!myFunc->Value(theX, theFVal))
      return Standard_False;

    theFVal *= theFVal;
    return Standard_True;
  }

  //see analogical method for abstract owner class math_Function
  virtual Standard_Integer GetStateNumber()
  {
    return 0;
  }

private:
  TheFunction *myFunc;
};


//=======================================================================
//function : FindVertex
//purpose  : 
//=======================================================================
void FindVertex (const TheArc& A,
                 const Handle(TheTopolTool)& Domain,
                 TheFunction& Func,
                 IntStart_SequenceOfPathPoint& pnt,
                 const Standard_Real Toler) 
{

// Find the vertex of the arc A restriction solutions. It stores
// Vertex in the list solutions pnt.


  TheVertex vtx;
  Standard_Real param,valf;
  Standard_Integer itemp;

  Domain->Initialize(A);
  Domain->InitVertexIterator();
  while (Domain->MoreVertex()) {
    vtx = Domain->Vertex();
    param = TheSOBTool::Parameter(vtx,A);

    // Evaluate the function and look compared to tolerance of the
    // Vertex. If distance <= tolerance then add a vertex to the list of solutions.
    // The arc is already assumed in the load function.

    Func.Value(param,valf);
    if (Abs(valf) <= Toler) {
      itemp = Func.GetStateNumber();
      pnt.Append(IntStart_ThePathPoint(Func.Valpoint(itemp),Toler, vtx,A,param));
      // Solution is added
    }
    Domain->NextVertex();
  }
}

Standard_Boolean IsDegenerated(const Handle(Adaptor3d_CurveOnSurface)& theCurve)
{
  if (theCurve->GetType() == GeomAbs_Circle)
  {
    gp_Circ aCirc = theCurve->Circle();
    if (aCirc.Radius() <= Precision::Confusion())
      return Standard_True;
  }
  return Standard_False;
}

Standard_Boolean IsDegenerated(const IntSurf_Quadric& theQuadric)
{
  GeomAbs_SurfaceType TypeQuad = theQuadric.TypeQuadric();
  if (TypeQuad == GeomAbs_Cone)
  {
    gp_Cone aCone = theQuadric.Cone();
    Standard_Real aSemiAngle = Abs(aCone.SemiAngle());
    if (aSemiAngle < 0.02 || aSemiAngle > 1.55)
      return Standard_True;
  }
  return Standard_False;
}

class SolInfo
{
public:
  SolInfo() : myMathIndex(-1), myValue(RealLast())
  {
  }

  void Init(const math_FunctionAllRoots& theSolution, const Standard_Integer theIndex)
  {
    myMathIndex = theIndex;
    myValue = theSolution.GetPoint(theIndex);
  }

  void Init(const IntCurveSurface_HInter& theSolution, const Standard_Integer theIndex)
  {
    myMathIndex = theIndex;
    myValue = theSolution.Point(theIndex).W();
  }

  Standard_Real Value() const
  {
    return myValue;
  }

  Standard_Integer Index() const
  {
    return myMathIndex;
  }

  bool operator>(const SolInfo& theOther) const
  {
    return myValue > theOther.myValue;
  }

  bool operator<(const SolInfo& theOther) const
  {
    return myValue < theOther.myValue;
  }

  bool operator==(const SolInfo& theOther) const
  {
    return myValue == theOther.myValue;
  }

  Standard_Real& ChangeValue()
  {
    return myValue;
  }

private:
  Standard_Integer myMathIndex;
  Standard_Real myValue;
};

static
void BoundedArc (const TheArc& A,
                 const Handle(TheTopolTool)& Domain,
                 const Standard_Real Pdeb,
                 const Standard_Real Pfin,
                 TheFunction& Func,
                 IntStart_SequenceOfPathPoint& pnt,
                 IntStart_SequenceOfSegment& seg,
                 const Standard_Real TolBoundary,
                 const Standard_Real TolTangency,
                 Standard_Boolean& Arcsol,
                 const Standard_Boolean RecheckOnRegularity)
{
  // Recherche des points solutions et des bouts d arc solution sur un arc donne.
  // On utilise la fonction math_FunctionAllRoots. Ne convient donc que pour
  // des arcs ayant un point debut et un point de fin (intervalle ferme de
  // parametrage).

  Standard_Integer i, Nbi = 0, Nbp = 0;

  gp_Pnt ptdeb,ptfin;
  Standard_Real pardeb = 0., parfin = 0.;
  Standard_Integer ideb,ifin,range,ranged,rangef;

  // Creer l echantillonage (math_FunctionSample ou classe heritant)
  // Appel a math_FunctionAllRoots

  //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
  //@@@ La Tolerance est asociee a l arc  ( Incoherence avec le cheminement )
  //@@@   ( EpsX ~ 1e-5   et ResolutionU et V ~ 1e-9 )
  //@@@   le vertex trouve ici n'est pas retrouve comme point d arret d une 
  //@@@   ligne de cheminement
  //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
  Standard_Real EpsX = 1.e-10;
  //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
  //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
  //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

  //  Standard_Integer NbEchant = TheSOBTool::NbSamplesOnArc(A); 
  Standard_Integer NbEchant = Func.NbSamples(); 
  if(NbEchant<100) NbEchant = 100; //-- lbr le 22 Avril 96 
  //-- Toujours des pbs 
  
  //-- Modif 24  Aout 93 -----------------------------
  Standard_Real nTolTangency = TolTangency;
  if((Pfin - Pdeb) < (TolTangency*10.0)) { 
    nTolTangency=(Pfin-Pdeb)*0.1;
  }   
  if(EpsX>(nTolTangency+nTolTangency)) { 
    EpsX = nTolTangency * 0.1; 
  }

  //--------------------------------------------------
  //-- Plante avec un edge avec 2 Samples  
  //-- dont les extremites son solutions (f=0) 
  //-- et ou la derivee est nulle 
  //-- Exemple : un segment diametre d une sphere
  //-- if(NbEchant<3) NbEchant = 3; //-- lbr le 19 Avril 95
  //--------------------------------------------------
  Standard_Real para=0,dist,maxdist;
  
  //-------------------------------------------------------------- REJECTIONS le 15 oct 98 
  Standard_Boolean Rejection=Standard_True;  
  Standard_Real maxdr,maxr,minr,ur,dur;
  minr=RealLast();
  maxr=-minr;
  maxdr=-minr;
  dur=(Pfin-Pdeb)*0.2;
  for(i=1,ur=Pdeb;i<=6;i++) { 
    Standard_Real F,D;
    if(Func.Values(ur,F,D)) { 
      Standard_Real lminr,lmaxr;
      if(D<0.0) D=-D;
      D*=dur+dur;
      if(D>maxdr) maxdr=D;
      lminr=F-D;
      lmaxr=F+D;
      if(lminr<minr) minr=lminr;
      if(lmaxr>maxr) maxr=lmaxr;
      if(minr<0.0 && maxr>0.0)  {
        Rejection=Standard_False;
        break;
      }
    }
    ur+=dur;
  }
  if(Rejection)
  {
    dur=0.001+maxdr+(maxr-minr)*0.1;
    minr-=dur;
    maxr+=dur;
    if(minr<0.0 && maxr>0.0)  { 	
      Rejection=Standard_False;
    }
  }

  Arcsol=Standard_False;

  if(Rejection==Standard_False)
  {
    const IntSurf_Quadric& aQuadric = Func.Quadric();
    GeomAbs_SurfaceType TypeQuad = aQuadric.TypeQuadric();
    GeomAbs_CurveType TypeConS = GeomAbs_OtherCurve;
    
    IntCurveSurface_HInter IntCS;
    Standard_Boolean IsIntCSdone = Standard_False;
    TColStd_SequenceOfReal Params;
    
    std::unique_ptr<math_FunctionAllRoots> pSol;
    
    math_FunctionSample Echant(Pdeb,Pfin,NbEchant);

    Standard_Boolean aelargir=Standard_True;
    //modified by NIZNHY-PKV Thu Apr 12 09:25:19 2001 f
    //
    //maxdist = 100.0*TolBoundary;
    maxdist = TolBoundary+TolTangency;
    //
    //modified by NIZNHY-PKV Thu Apr 12 09:25:23 2001 t
    for(i=1; i<=NbEchant && aelargir;i++) { 
      Standard_Real u = Echant.GetParameter(i);
      if(Func.Value(u,dist)) { 
        if(dist>maxdist || -dist>maxdist) {
          aelargir=Standard_False;
        }
      }
    }
    if(!(aelargir && maxdist<0.01)) { 
      maxdist = TolBoundary;
    }

    if (TypeQuad != GeomAbs_OtherSurface) //intersection of boundary curve and quadric surface
    {
      //Exact solution
      Handle(Adaptor3d_Surface) aSurf = Func.Surface();
      Adaptor3d_CurveOnSurface ConS(A, aSurf);
      TypeConS = ConS.GetType();
#ifdef OCCT_DEBUG
      Handle(Geom_Curve) CurveConS;
      switch(TypeConS)
      {
      case GeomAbs_Line:
        {
          CurveConS = new Geom_Line(ConS.Line());
          break;
        }
      case GeomAbs_Circle:
        {
          CurveConS = new Geom_Circle(ConS.Circle());
          break;
        }
      case GeomAbs_Ellipse:
        {
          CurveConS = new Geom_Ellipse(ConS.Ellipse());
          break;
        }
      case GeomAbs_Hyperbola:
        {
          CurveConS = new Geom_Hyperbola(ConS.Hyperbola());
          break;
        }
      case GeomAbs_Parabola:
        {
          CurveConS = new Geom_Parabola(ConS.Parabola());
          break;
        }
      case GeomAbs_BezierCurve:
        {
          CurveConS = ConS.Bezier();
          break;
        }
      case GeomAbs_BSplineCurve:
        {
          CurveConS = ConS.BSpline();
          break;
        }
      default:
        {
          Standard_Real MaxDeviation, AverageDeviation;
          GeomLib::BuildCurve3d(1.e-5, ConS, ConS.FirstParameter(), ConS.LastParameter(),
                                CurveConS, MaxDeviation, AverageDeviation);
          break;
        }
      }
#endif
      Handle(Adaptor3d_CurveOnSurface) HConS = new Adaptor3d_CurveOnSurface(ConS);
      Handle(Geom_Surface) QuadSurf;
      switch (TypeQuad)
      {
      case GeomAbs_Plane:
        {
          QuadSurf = new Geom_Plane(aQuadric.Plane());
          break;
        }
      case GeomAbs_Cylinder:
        {
          QuadSurf = new Geom_CylindricalSurface(aQuadric.Cylinder());
          break;
        }
      case GeomAbs_Cone:
        {
          QuadSurf = new Geom_ConicalSurface(aQuadric.Cone());
          break;
        }
      case GeomAbs_Sphere:
        {
          QuadSurf = new Geom_SphericalSurface(aQuadric.Sphere());
          break;
        }
      case GeomAbs_Torus:
        {
          QuadSurf = new Geom_ToroidalSurface(aQuadric.Torus());
          break;
        }
      default:
        break;
      }
      Handle(GeomAdaptor_Surface) GAHsurf = new GeomAdaptor_Surface(QuadSurf);
      
      if ((TypeConS == GeomAbs_Line ||
           TypeConS == GeomAbs_Circle ||
           TypeConS == GeomAbs_Ellipse ||
           TypeConS == GeomAbs_Parabola ||
           TypeConS == GeomAbs_Hyperbola) &&
          TypeQuad != GeomAbs_Torus &&
          !IsDegenerated(HConS) &&
          !IsDegenerated(aQuadric))
      {
        //exact intersection for only canonic curves and real quadric surfaces
        IntCS.Perform(HConS, GAHsurf);
      }
      
      IsIntCSdone = IntCS.IsDone();
      if (IsIntCSdone)
      {
        Nbp = IntCS.NbPoints();
        Nbi = IntCS.NbSegments();
      }
      //If we have not got intersection, it may be touch with some tolerance,
      //need to be checked
      if (Nbp == 0 && Nbi == 0)
        IsIntCSdone = Standard_False;

    } //if (TypeQuad != GeomAbs_OtherSurface) - intersection of boundary curve and quadric surface
    
    if (!IsIntCSdone)
    {
      pSol.reset(new math_FunctionAllRoots(Func,Echant,EpsX,maxdist,maxdist)); //-- TolBoundary,nTolTangency);
      
      if (!pSol->IsDone()) {throw Standard_Failure();}
      
      Nbp=pSol->NbPoints();
    }
    //
    //jgv: build solution on the whole boundary
    if (RecheckOnRegularity && Nbp > 0 && IsRegularity(A, Domain))
    {
      //Standard_Real theTol = Domain->MaxTolerance(A);
      //theTol += theTol;
      Standard_Real theTol = 5.e-4;
      math_FunctionAllRoots SolAgain(Func,Echant,EpsX,theTol,theTol); //-- TolBoundary,nTolTangency);

      if (!SolAgain.IsDone()) {throw Standard_Failure();}

      Standard_Integer Nbi_again = SolAgain.NbIntervals();

      if (Nbi_again > 0)
      {
        Standard_Integer NbSamples = 10;
        Standard_Real delta = (Pfin - Pdeb)/NbSamples;
        Standard_Real GlobalTol = theTol*10;
        Standard_Boolean SolOnBoundary = Standard_True;
        for (i = 0; i <= NbSamples; i++)
        {
          Standard_Real aParam = Pdeb + i*delta;
          Standard_Real aValue;
          Func.Value(aParam, aValue);
          if (Abs(aValue) > GlobalTol)
          {
            SolOnBoundary = Standard_False;
            break;
          }
        }

        if (SolOnBoundary)
        {
          for (i = 1; i <= Nbi_again; i++)
          {
            IntStart_TheSegment newseg;
            newseg.SetValue(A);
            // Recuperer point debut et fin, et leur parametre.
            SolAgain.GetInterval(i,pardeb,parfin);

            if (Abs(pardeb - Pdeb) <= Precision::PConfusion())
              pardeb = Pdeb;
            if (Abs(parfin - Pfin) <= Precision::PConfusion())
              parfin = Pfin;

            SolAgain.GetIntervalState(i,ideb,ifin);

            //-- cout<<" Debug : IntStart_SearchOnBoundaries_1.gxx : i= "<<i<<" ParDeb:"<<pardeb<<"  ParFin:"<<parfin<<endl;

            ptdeb=Func.Valpoint(ideb);
            ptfin=Func.Valpoint(ifin);

            PointProcess(ptdeb,pardeb,A,Domain,pnt,theTol,ranged);
            newseg.SetLimitPoint(pnt.Value(ranged),Standard_True);
            PointProcess(ptfin,parfin,A,Domain,pnt,theTol,rangef);
            newseg.SetLimitPoint(pnt.Value(rangef),Standard_False);
            seg.Append(newseg);
          }
          Arcsol=Standard_True;
          return;
        }
      }
    } //if (RecheckOnRegularity && Nbp > 0 && IsRegularity(A, Domain))
    ////////////////////////////////////////////

    //-- detection du cas ou la fonction est quasi tangente et que les 
    //-- zeros sont quasi confondus. 
    //-- Dans ce cas on prend le point "milieu"
    //-- On suppose que les solutions sont triees. 

    if(Nbp) { 
      NCollection_Array1<SolInfo> aSI(1, Nbp);

      for(i=1;i<=Nbp;i++)
      {
        if (IsIntCSdone)
          aSI(i).Init(IntCS, i);
        else
          aSI(i).Init(*pSol, i);
      }

      std::sort(aSI.begin(), aSI.end());

      //modified by NIZNHY-PKV Wed Mar 21 18:34:18 2001 f
      //////////////////////////////////////////////////////////
      // The treatment of the situation when line(arc) that is 
      // tangent to cylinder(domain). 
      // We should have only one solution i.e Nbp=1. Ok?
      // But we have 2,3,.. solutions.     That is wrong ersult.
      // The TreatLC(...) function is dedicated to solve the pb.
      //                           PKV Fri Mar 23 12:17:29 2001

      Standard_Integer ip = TreatLC (A, Domain, aQuadric, TolBoundary, pnt);
      if (ip) {
        //////////////////////////////////////////////////////////
        //modified by NIZNHY-PKV Wed Mar 21 18:34:23 2001 t
        // 
        // Using of old usual way proposed by Laurent 
        //
        for(i=1;i<Nbp;i++) { 
          Standard_Real parap1 = aSI(i + 1).Value();
          para = aSI(i).Value();

          Standard_Real param=(para+parap1)*0.5;
          Standard_Real yf = 0.0;
          Standard_Real ym = 0.0;
          Standard_Real yl = 0.0;
          if(Func.Value(param,ym) && Abs(ym) < maxdist) {
            Standard_Real sm = Sign(1., ym);
            Standard_Boolean aTang = Func.Value(para,yf) && Func.Value(parap1,yl);
            if (aTang) {
              //Line can be tangent surface if all distances less then maxdist
              aTang = aTang && Abs(yf) < maxdist && Abs(yl) < maxdist;
            }
            if (aTang && IsIntCSdone && TypeConS == GeomAbs_Line) {
              //Interval is got by exact intersection
              //Line can be tangent if all points are on the same side of surface
              //it means that signs of all distances are the same
              Standard_Real sf = Sign(1., yf), sl = Sign(1., yl);
              aTang = aTang && (sm == sf) && (sm == sl);
            }
            if(aTang) { 
              //  Modified by skv - Tue Aug 31 12:13:51 2004 OCC569 Begin
              // Consider this interval as tangent one. Treat it to find
              // parameter with the lowest function value.
              // Compute the number of nodes.
              Standard_Real    aTol = TolBoundary*1000.0;
              if(aTol > 0.001)
                aTol = 0.001;

              // fix floating point exception 569, chl-922-e9
              parap1 = (Abs(parap1) < 1.e9) ? parap1 : ((parap1 >= 0.) ? 1.e9 : -1.e9);
              para   = (Abs(para) < 1.e9) ? para : ((para >= 0.) ? 1.e9 : -1.e9);

              Standard_Integer aNbNodes = RealToInt(Ceiling((parap1 - para)/aTol));

              Standard_Real    aVal     = RealLast();
              Standard_Real aValMax = 0.;
              //Standard_Integer aNbNodes = 23;
              Standard_Real    aDelta   = (parap1 - para)/(aNbNodes + 1.);
              Standard_Integer ii;
              Standard_Real    aCurPar;
              Standard_Real    aCurVal;

              for (ii = 0; ii <= aNbNodes + 1; ii++) {
                aCurPar = (ii < aNbNodes + 1) ? para + ii*aDelta : parap1;

                if (Func.Value(aCurPar, aCurVal)) {
                  Standard_Real anAbsVal = Abs(aCurVal);
                  if (anAbsVal < aVal) {
                    aVal  = anAbsVal;
                    param = aCurPar;
                  }
                  if (anAbsVal > aValMax)
                  {
                    aValMax = anAbsVal;
                  }
                }
              }
              // At last, interval got by exact intersection can be considered as tangent if
              // minimal distance is inside interval and
              // minimal and maximal values are almost the same
              if (IsIntCSdone && aNbNodes > 1) {
                aTang = Abs(param - para) > EpsX && Abs(parap1 - param) > EpsX &&
                  0.01*aValMax <= aVal;
              }
              if (aTang)
              {
                aSI(i).ChangeValue() = Pdeb - 1;
                aSI(i + 1).ChangeValue() = param;
              }
            }
          }
        }

        for (i=1; i<=Nbp; i++) {
          para = aSI(i).Value();
          if((para-Pdeb)<EpsX || (Pfin-para)<EpsX)
            continue;

          if(!Func.Value(para,dist))
            continue;

          dist = Abs(dist);

          Standard_Integer anIndx = -1;
          //const Standard_Real aParam = Sol->GetPoint(aSI(i).Index());
          const Standard_Real aParam = aSI(i).Value();
          if (dist < maxdist)
          {
            if (!IsIntCSdone &&
                (Abs(aParam - Pdeb) <= Precision::PConfusion() || Abs(aParam - Pfin) <= Precision::PConfusion()))
            {
              anIndx = pSol->GetPointState(aSI(i).Index());
            }
          }

          gp_Pnt aPnt(anIndx < 0 ? Func.LastComputedPoint() : Func.Valpoint(anIndx));

          if (dist > 0.1*Precision::Confusion())
          {
            //Precise found points. It results in following:
            //  1. Make the vertex nearer to the intersection line
            //    (see description to issue #27252 in order to 
            //    understand necessity).
            //  2. Merge two near vertices to single point.

            //All members in TabSol array has already been sorted in increase order.
            //Now, we limit precise boundaries in order to avoid changing this order.
            const Standard_Real aFPar = (i == 1) ? Pdeb : (para + aSI(i - 1).Value()) / 2.0;
            const Standard_Real aLPar = (i == Nbp) ? Pfin : (para + aSI(i + 1).Value()) / 2.0;

            MinFunction aNewFunc(Func);
            math_BrentMinimum aMin(Precision::Confusion());

            aMin.Perform(aNewFunc, aFPar, para, aLPar);
            if(aMin.IsDone())
            {
              para = aMin.Location();
              const gp_Pnt2d aP2d(A->Value(para));
              aPnt = Func.Surface()->Value(aP2d.X(), aP2d.Y());
            }
          }

          PointProcess(aPnt, para, A, Domain, pnt, TolBoundary, range);
        }
      }// end of if(ip)
    } // end of if(Nbp)  

    // Pour chaque intervalle trouve faire
    //   Traiter les extremites comme des points
    //   Ajouter intervalle dans la liste des segments

    if (!IsIntCSdone)
      Nbi = pSol->NbIntervals();

    if (!RecheckOnRegularity && Nbp) { 
      //--cout<<" Debug : IntStart_SearchOnBoundaries_1.gxx :Nbp>0  0 <- Nbi "<<Nbi<<endl;
      Nbi=0; 
    }

    //-- cout<<" Debug : IntStart_SearchOnBoundaries_1.gxx : Nbi : "<<Nbi<<endl;

    for (i=1; i<=Nbi; i++) {
      IntStart_TheSegment newseg;
      newseg.SetValue(A);
      // Recuperer point debut et fin, et leur parametre.
      if (IsIntCSdone)
      {
        IntCurveSurface_IntersectionSegment IntSeg = IntCS.Segment(i);
        IntCurveSurface_IntersectionPoint End1 = IntSeg.FirstPoint();
        IntCurveSurface_IntersectionPoint End2 = IntSeg.SecondPoint();
        pardeb = End1.W();
        parfin = End2.W();
        ptdeb  = End1.Pnt();
        ptfin  = End2.Pnt();
      }
      else
      {
        pSol->GetInterval(i,pardeb,parfin);
        pSol->GetIntervalState(i,ideb,ifin);

        //-- cout<<" Debug : IntStart_SearchOnBoundaries_1.gxx : i= "<<i<<" ParDeb:"<<pardeb<<"  ParFin:"<<parfin<<endl;
        
        ptdeb=Func.Valpoint(ideb);
        ptfin=Func.Valpoint(ifin);
      }

      PointProcess(ptdeb,pardeb,A,Domain,pnt,TolBoundary,ranged);
      newseg.SetLimitPoint(pnt.Value(ranged),Standard_True);
      PointProcess(ptfin,parfin,A,Domain,pnt,TolBoundary,rangef);
      newseg.SetLimitPoint(pnt.Value(rangef),Standard_False);
      seg.Append(newseg);
    }

    if (Nbi==1) {
      if((Abs(pardeb - Pdeb) < Precision::PConfusion()) &&
         (Abs(parfin - Pfin) < Precision::PConfusion()))
      {
        Arcsol=Standard_True;
      }
    }
  }
}

//=======================================================================
//function : ComputeBoundsfromInfinite
//purpose  : 
//=======================================================================
// - PROVISIONAL - TEMPORARY - NOT GOOD - NYI - TO DO
// - Temporary - temporary - not good - nyi - to do
void ComputeBoundsfromInfinite(TheFunction& Func,
                               Standard_Real& PDeb,
                               Standard_Real& PFin,
                               Standard_Integer& NbEchant) 
{ 
  
  // - We are looking for parameters for start and end of the arc (2d curve)
  // - Infinity, a way to intersect the quadric with a portion of arc
  // - Finished.
  //
  // - The quadric is a plane, a cylinder, a cone and a sphere.
  // - Idea: We take any point on the arc and the fact grow
  // - Terminals to the signed distance function values or is likely
  // - S cancel.
  //
  // - WARNING: The following calculations provide a very estimated coarse parameters.
  // - This avoids the raises and allows a case of Boxes
  // - Inifinies walk. It will take this code
  // - With curve surface intersections.

  NbEchant = 100;

  Standard_Real U0    = 0.0;
  Standard_Real dU    = 0.001;
  Standard_Real Dist0,Dist1;

  Func.Value(U0   , Dist0);
  Func.Value(U0+dU, Dist1);
  Standard_Real dDist = Dist1 - Dist0;
  if(dDist) { 
    U0  -=  dU*Dist0 / dDist;
    PDeb = PFin = U0;
    Standard_Real Umin = U0 - 1e5;
    Func.Value(Umin   , Dist0);
    Func.Value(Umin+dU, Dist1);
    dDist = Dist1-Dist0;
    if(dDist) { 
      Umin  -=  dU*Dist0 / dDist;
    }
    else { 
      Umin-=10.0; 
    }
    Standard_Real Umax = U0 + 1e8;
    Func.Value(Umax   , Dist0);
    Func.Value(Umax+dU, Dist1);
    dDist = Dist1-Dist0;
    if(dDist) { 
      Umax  -=  dU*Dist0 / dDist;
    }
    else { 
      Umax+=10.0; 
    }
    if(Umin>U0) { Umin=U0-10.0; } 
    if(Umax<U0) { Umax=U0+10.0; } 
    
    PFin = Umax + 10. * (Umax - Umin);
    PDeb = Umin - 10. * (Umax - Umin);
  }
  else { 
    //-- Possibilite de Arc totalement inclu ds Quad
    PDeb = 1e10;
    PFin = -1e10;
  }
} 

//=======================================================================
//function : PointProcess
//purpose  : 
//=======================================================================
void PointProcess (const gp_Pnt& Pt,
                   const Standard_Real Para,
                   const TheArc& A,
                   const Handle(TheTopolTool)& Domain,
                   IntStart_SequenceOfPathPoint& pnt,
                   const Standard_Real Tol,
                   Standard_Integer& Range) 
{

// Check to see if a solution point is coincident with a vertex.
// If confused, you should find this vertex in the list of
// Start. It then returns the position of this point in the list pnt.
// Otherwise, add the point in the list.
  
  Standard_Integer k;
  Standard_Boolean found,goon;
  Standard_Real dist,toler;

  Standard_Integer Nbsol = pnt.Length();
  TheVertex vtx;
  IntStart_ThePathPoint ptsol;

  Domain->Initialize(A);
  Domain->InitVertexIterator();
  found = Standard_False;
  goon = Domain->MoreVertex();
  while (goon) {
    vtx = Domain->Vertex();
    dist= Abs(Para-TheSOBTool::Parameter(vtx,A));
    toler = TheSOBTool::Tolerance(vtx,A);
#ifdef OCCT_DEBUG
    if(toler>0.1) { 
      std::cout<<"IntStart_SearchOnBoundaries_1.gxx  : ** WARNING ** Tol Vertex="<<toler<<std::endl;
      std::cout<<"                                     Ou Edge degenere Ou Kro pointu"<<std::endl;
      if(toler>10000) toler=1e-7;
    }
#endif

    if (dist <= toler) {
      // Locate the vertex in the list of solutions
      k=1;
      found = (k>Nbsol);
      while (!found) {
        ptsol = pnt.Value(k);
        if (!ptsol.IsNew()) {
        //jag 940608  if (ptsol.Vertex() == vtx && ptsol.Arc()    == A) {
          if (Domain->Identical(ptsol.Vertex(),vtx) &&
                    ptsol.Arc()    == A &&
                    Abs(ptsol.Parameter()-Para) <= toler) {
            found=Standard_True;
          }
          else {
            k=k+1;
            found=(k>Nbsol);
          }
        }
        else {
          k=k+1;
          found=(k>Nbsol);
        }
      }
      if (k<=Nbsol) {     // We find the vertex
        Range = k;
      }
      else {              // Otherwise
        ptsol.SetValue(Pt,Tol,vtx,A,Para);
        pnt.Append(ptsol);
        Range = pnt.Length();
      }
      found = Standard_True;
      goon = Standard_False;
    }
    else {
      Domain->NextVertex();
      goon = Domain->MoreVertex();
    }
  }

  if (!found) {   // No one is falling on a vertex
    //jgv: do not add segment's extremities if they already exist
    Standard_Boolean found_internal = Standard_False;
    for (k = 1; k <= pnt.Length(); k++)
    {
      ptsol = pnt.Value(k);
      if (ptsol.Arc() != A ||
          !ptsol.IsNew()) //vertex
        continue;
      if (Abs(ptsol.Parameter()-Para) <= Precision::PConfusion())
      {
        found_internal = Standard_True;
        Range = k;
      }
    }
    /////////////////////////////////////////////////////////////

    if (!found_internal)
    {
      Standard_Real TOL=Tol;
      TOL*=1000.0; 
      //if(TOL>0.001) TOL=0.001;
      if(TOL>0.005) TOL=0.005; //#24643
      
      ptsol.SetValue(Pt,TOL,A,Para);
      pnt.Append(ptsol);
      Range = pnt.Length();
    }
  }
}

//=======================================================================
//function : IsRegularity
//purpose  : 
//=======================================================================
Standard_Boolean IsRegularity(const TheArc& /*A*/,
                              const Handle(TheTopolTool)& aDomain)
{
  Standard_Address anEAddress=aDomain->Edge();
  if (anEAddress==NULL) {
    return Standard_False;
  }
  
  TopoDS_Edge* anE=(TopoDS_Edge*)anEAddress;
  
  return (BRep_Tool::HasContinuity(*anE));
}

//=======================================================================
//function : TreatLC
//purpose  : 
//=======================================================================
Standard_Integer TreatLC (const TheArc& A,
                          const Handle(TheTopolTool)& aDomain,
                          const IntSurf_Quadric& aQuadric,
                          const Standard_Real TolBoundary,
                          IntStart_SequenceOfPathPoint& pnt)
{
  Standard_Integer anExitCode=1, aNbExt;
  
  Standard_Address anEAddress=aDomain->Edge();
  if (anEAddress==NULL) {
    return anExitCode;
  }
  
  TopoDS_Edge* anE=(TopoDS_Edge*)anEAddress;

  if (BRep_Tool::Degenerated(*anE)) {
    return anExitCode;
  }
  
  GeomAbs_CurveType   aTypeE;
  BRepAdaptor_Curve aBAC(*anE);
  aTypeE=aBAC.GetType();
  
  if (aTypeE!=GeomAbs_Line) {
    return anExitCode;
  }
  
  GeomAbs_SurfaceType aTypeS;
  aTypeS=aQuadric.TypeQuadric();
  
  if (aTypeS!=GeomAbs_Cylinder) {
    return anExitCode;
  }
  
  Standard_Real f, l, U1f, U1l, U2f, U2l, UEgde, TOL, aDist, aR, aRRel, Tol;
  Handle(Geom_Curve) aCEdge=BRep_Tool::Curve(*anE, f, l);
  
  gp_Cylinder aCyl=aQuadric.Cylinder();
  const gp_Ax1& anAx1=aCyl.Axis();
  gp_Lin aLin(anAx1);
  Handle(Geom_Line) aCAxis=new Geom_Line (aLin);
  aR=aCyl.Radius();
  
  U1f = aCAxis->FirstParameter();
  U1l = aCAxis->LastParameter();
  
  U2f = aCEdge->FirstParameter();
  U2l = aCEdge->LastParameter();
  

  GeomAdaptor_Curve C1, C2;
  
  C1.Load(aCAxis);
  C2.Load(aCEdge);
  
  Tol = Precision::PConfusion();

  Extrema_ExtCC anExtCC(C1, C2, U1f, U1l, U2f, U2l, Tol, Tol); 

  aNbExt=anExtCC.NbExt();
  if (aNbExt!=1) {
    return anExitCode;
  }

  gp_Pnt P1,PEdge;
  Extrema_POnCurv PC1, PC2;
  
  anExtCC.Points(1, PC1, PC2);
  
  P1   =PC1.Value();
  PEdge=PC2.Value();
  
  UEgde=PC2.Parameter();
  
  aDist=PEdge.Distance(P1);
  aRRel=fabs(aDist-aR)/aR;
  if (aRRel > TolBoundary) {
    return anExitCode;
  }

  if (UEgde < (f+TolBoundary) || UEgde > (l-TolBoundary)) {
    return anExitCode;
  }
  //
  // Do not wonder !
  // It was done as into PointProcess(...) function 
  //printf("TreatLC()=> tangent line is found\n");
  TOL=1000.*TolBoundary;
  if(TOL>0.001) TOL=0.001;
  
  IntStart_ThePathPoint ptsol;
  ptsol.SetValue(PEdge, TOL, A, UEgde);
  pnt.Append(ptsol);

  anExitCode=0;
  return anExitCode;

}


//=======================================================================
//function : IntStart_SearchOnBoundaries::IntStart_SearchOnBoundaries
//purpose  : 
//=======================================================================
IntStart_SearchOnBoundaries::IntStart_SearchOnBoundaries ()
:  done(Standard_False),
   all(Standard_False)  
{
}  

//=======================================================================
//function : Perform
//purpose  : 
//=======================================================================
  void IntStart_SearchOnBoundaries::Perform (TheFunction& Func,
                                             const Handle(TheTopolTool)& Domain,
                                             const Standard_Real TolBoundary,
                                             const Standard_Real TolTangency,
                                             const Standard_Boolean RecheckOnRegularity)
{
  
  done = Standard_False;
  spnt.Clear();
  sseg.Clear();

  Standard_Boolean Arcsol;
  Standard_Real PDeb,PFin, prm, tol;
  Standard_Integer i, nbknown, nbfound,index;
  gp_Pnt pt;
  
  Domain->Init();

  if (Domain->More()) {
    all  = Standard_True;
  }
  else {
    all = Standard_False;
  }

  while (Domain->More()) {
    TheArc A = Domain->Value();
    if (!TheSOBTool::HasBeenSeen(A)) {
      Func.Set(A);
      FindVertex(A,Domain,Func,spnt,TolBoundary);
      TheSOBTool::Bounds(A,PDeb,PFin);
      if(Precision::IsNegativeInfinite(PDeb) || 
         Precision::IsPositiveInfinite(PFin)) { 
        Standard_Integer NbEchant;
        ComputeBoundsfromInfinite(Func,PDeb,PFin,NbEchant);
      }
      BoundedArc(A,Domain,PDeb,PFin,Func,spnt,sseg,TolBoundary,TolTangency,Arcsol,RecheckOnRegularity);
      all = (all && Arcsol);
    }
    
    else {
      // as it seems we'll never be here, because 
      // TheSOBTool::HasBeenSeen(A) always returns FALSE
      nbfound = spnt.Length();

      // On recupere les points connus
      nbknown = TheSOBTool::NbPoints(A);
      for (i=1; i<=nbknown; i++) {
        TheSOBTool::Value(A,i,pt,tol,prm);
        if (TheSOBTool::IsVertex(A,i)) {
          TheVertex vtx;
          TheSOBTool::Vertex(A,i,vtx);
          spnt.Append(IntStart_ThePathPoint(pt,tol,vtx,A,prm));
        }
        else {
          spnt.Append(IntStart_ThePathPoint(pt,tol,A,prm));
        }
      }
      // On recupere les arcs solutions
      nbknown = TheSOBTool::NbSegments(A);
      for (i=1; i<=nbknown; i++) {
        IntStart_TheSegment newseg;
        newseg.SetValue(A);
        if (TheSOBTool::HasFirstPoint(A,i,index)) {
          newseg.SetLimitPoint(spnt.Value(nbfound+index),Standard_True);
        }
        if (TheSOBTool::HasLastPoint(A,i,index)) {
          newseg.SetLimitPoint(spnt.Value(nbfound+index),Standard_False);
        }
        sseg.Append(newseg);
      }
      all = (all& TheSOBTool::IsAllSolution(A));
    }
    Domain->Next();
  }
  done = Standard_True;
}