Logo Search packages:      
Sourcecode: opencascade version File versions  Download package

ChFi3d_Builder_0.cxx

// File:    ChFi3d_Builder_0.cxx
// Created: Thu Dec 16 15:54:39 1993
// Author:  Isabelle GRIGNON
//          <isg@zerox>

//  modified by ofv - Thu Feb 26 11:18:16 2004 OCC5246
//  Modified by skv - Fri Oct 24 14:24:47 2003 OCC4077
//  Modified by skv - Mon Jun 16 15:50:44 2003 OCC615

#include <ChFi3d.hxx>
#include <Precision.hxx>

#include <Standard_NotImplemented.hxx>
#include <Standard_ConstructionError.hxx>

#include <gp.hxx>
#include <gp_Circ.hxx>
#include <gp_Elips.hxx>
#include <gp_Lin.hxx>
#include <gp_Pnt.hxx>
#include <gp_Pnt2d.hxx>
#include <gp_Lin2d.hxx>
#include <ElCLib.hxx>
#include <ElSLib.hxx>
#include <BSplCLib.hxx>
#include <GeomLib.hxx>

#include <TColgp_Array1OfPnt2d.hxx>
#include <TColgp_Array1OfPnt.hxx>
#include <TColgp_Array1OfXYZ.hxx>
#include <TColStd_Array1OfInteger.hxx>
#include <TColStd_Array1OfReal.hxx>

#include <Geom_TrimmedCurve.hxx>
#include <Geom_BSplineCurve.hxx>
#include <Geom_Surface.hxx>
#include <Geom_CylindricalSurface.hxx>
#include <Geom_RectangularTrimmedSurface.hxx>
#include <Geom_Plane.hxx>
#include <Geom_Line.hxx>
#include <Geom_Circle.hxx>
#include <Geom_Ellipse.hxx>
#include <Geom2d_BezierCurve.hxx>
#include <Geom2d_BSplineCurve.hxx>
#include <Geom2d_Line.hxx>
#include <Geom2d_Circle.hxx>
#include <Geom2d_Ellipse.hxx>
#include <Geom2d_Hyperbola.hxx>
#include <Geom2d_Parabola.hxx>
#include <Geom2d_TrimmedCurve.hxx>
#include <Geom2d_Line.hxx>
#include <Geom2d_OffsetCurve.hxx>
#include <Geom2dAdaptor_Curve.hxx>
#include <Geom2dAdaptor_HCurve.hxx>
#include <Adaptor3d_TopolTool.hxx>
#include <Adaptor3d_CurveOnSurface.hxx>
#include <Adaptor3d_HCurveOnSurface.hxx>
#include <GeomAdaptor_HSurface.hxx>

#include <FairCurve_Batten.hxx>
#include <FairCurve_AnalysisCode.hxx>
#include <Convert_ParameterisationType.hxx>
#include <GeomConvert_CompCurveToBSplineCurve.hxx>
#include <GeomConvert.hxx>
#include <GeomLib_Interpolate.hxx>
#include <GeomAPI_ProjectPointOnSurf.hxx>
#include <GeomAPI_ProjectPointOnCurve.hxx>
#include <GC_MakeCircle.hxx>
#include <BRepAdaptor_Curve.hxx>
#include <BRepAdaptor_HCurve.hxx>
#include <BRepAdaptor_HCurve2d.hxx>
#include <BRepAdaptor_Surface.hxx>
#include <BRepTopAdaptor_HVertex.hxx>

#include <BRep_Tool.hxx>
#include <BRep_Builder.hxx>
#include <BRepTools.hxx>
#include <BRepTools_WireExplorer.hxx>
#include <BRepLib.hxx>
#include <BRepLib_MakeEdge.hxx>
#include <BRepLib_MakeWire.hxx>
#include <BRepLib_MakeFace.hxx>

#include <TopAbs.hxx>
#include <TopoDS_Shape.hxx>
#include <TopoDS_Edge.hxx>
#include <TopoDS_Vertex.hxx>
#include <TopoDS_Wire.hxx>
#include <TopoDS_Face.hxx>
#include <TopExp.hxx>
#include <TopExp_Explorer.hxx>
#include <TopTools_Array1OfShape.hxx>


#include <GeomAbs_Shape.hxx>
#include <Bnd_Box2d.hxx>

//#include <math_FunctionSample.hxx>
//#include <math_FunctionAllRoots.hxx>
#include <GCPnts_AbscissaPoint.hxx>

#include <IntCurveSurface_TheQuadCurvFuncOfTheQuadCurvExactHInter.hxx>
#include <IntCurveSurface_HInter.hxx>
#include <IntCurveSurface_IntersectionPoint.hxx>
#include <IntSurf_Quadric.hxx>
#include <IntSurf_PntOn2S.hxx>
#include <IntSurf_LineOn2S.hxx>
#include <IntAna_QuadQuadGeo.hxx>
#include <IntAna2d_AnaIntersection.hxx>
#include <IntRes2d_IntersectionPoint.hxx>
#include <IntPatch_ThePWalkingInterOfThePPIntOfIntersection.hxx>
#include <IntPatch_TheWLineOfIntersection.hxx>
#include <Geom2dInt_GInter.hxx>
#include <GeomInt_WLApprox.hxx>
#include <GeomInt_IntSS.hxx>
#include <AppParCurves_MultiBSpCurve.hxx>
#include <Approx_SameParameter.hxx>

#include <TopAbs.hxx>
#include <TopoDS_Shape.hxx>
#include <TopoDS_Edge.hxx>
#include <TopExp.hxx>

#include <TopOpeBRepDS.hxx>
#include <TopOpeBRepDS_Surface.hxx>
#include <TopOpeBRepDS_Point.hxx>
#include <TopOpeBRepDS_SolidSurfaceInterference.hxx>
#include <TopOpeBRepDS_CurvePointInterference.hxx>
#include <TopOpeBRepDS_ListOfInterference.hxx>
#include <TopOpeBRepDS_InterferenceIterator.hxx>
#include <ProjLib_ProjectedCurve.hxx>

#include <BRepBlend_PointOnRst.hxx>

#include <ChFiDS_HData.hxx>
#include <ChFiDS_SurfData.hxx>
#include <ChFiDS_FaceInterference.hxx>
#include <ChFiDS_Spine.hxx>
#include <ChFiDS_FilSpine.hxx>
#include <ChFiDS_SequenceOfSurfData.hxx>
#include <ChFiDS_Regul.hxx>
#include <Law_Function.hxx>
#include <Law_Composite.hxx>
#include <GeomAPI_PointsToBSpline.hxx>
#include <GeomLProp_CLProps.hxx>

#include <ChFi3d_Builder_0.hxx>

#ifdef DRAW
#include <Draw_Appli.hxx>
#include <Draw.hxx>
#include <Draw_Marker3D.hxx>
#include <Draw_Marker2D.hxx>
#include <DrawTrSurf.hxx>
#include <DBRep.hxx>
static Standard_Integer NbSp = 1;
static Standard_Integer NbelFa = 1;
#endif

#ifdef DEB
#include <OSD_Chronometer.hxx>
extern Standard_Boolean ChFi3d_GetcontextFORCEBLEND(); 
extern Standard_Boolean ChFi3d_GettraceDRAWINT();
extern Standard_Boolean ChFi3d_GettraceDRAWENLARGE();
extern Standard_Boolean ChFi3d_GettraceDRAWSPINE();
extern Standard_Real  t_sameparam, t_batten;
extern void ChFi3d_InitChron(OSD_Chronometer& ch);
extern void ChFi3d_ResultChron(OSD_Chronometer & ch,Standard_Real& time);
#endif

#include <stdio.h>

#include <GeomAdaptor_HCurve.hxx>
#include <BRepAdaptor_HSurface.hxx>

Standard_Real ChFi3d_InPeriod(const Standard_Real U, 
                        const Standard_Real UFirst, 
                        const Standard_Real ULast,
                        const Standard_Real Eps)
{
  Standard_Real u = U, period = ULast - UFirst;
  while (Eps < (UFirst-u)) u += period;
  while (Eps > (ULast -u)) u -= period;
  if ( u < UFirst) u = UFirst;
  return u;
}

//=======================================================================
//function : Boite 
//purpose  : Calcul des min/max uv du conge a intersecter.
//=======================================================================

void ChFi3d_Boite(const gp_Pnt2d& p1,const gp_Pnt2d& p2,
              Standard_Real& mu,Standard_Real& Mu,
              Standard_Real& mv,Standard_Real& Mv)
{
  mu = Min(p1.X(),p2.X()); Mu = Max(p1.X(),p2.X());
  mv = Min(p1.Y(),p2.Y()); Mv = Max(p1.Y(),p2.Y());
}

//=======================================================================
//function : Boite 
//purpose  : Calcul des min/max uv du conge a intersecter.
//=======================================================================

void ChFi3d_Boite(const gp_Pnt2d& p1,const gp_Pnt2d& p2,
              const gp_Pnt2d& p3,const gp_Pnt2d& p4,
              Standard_Real& Du,Standard_Real& Dv,
              Standard_Real& mu,Standard_Real& Mu,
              Standard_Real& mv,Standard_Real& Mv)
{
  Standard_Real a,b;
  a = Min(p1.X(),p2.X());  b = Min(p3.X(),p4.X()); mu = Min(a,b);
  a = Max(p1.X(),p2.X());  b = Max(p3.X(),p4.X()); Mu = Max(a,b);
  a = Min(p1.Y(),p2.Y());  b = Min(p3.Y(),p4.Y()); mv = Min(a,b);
  a = Max(p1.Y(),p2.Y());  b = Max(p3.Y(),p4.Y()); Mv = Max(a,b);
  Du = Mu - mu;
  Dv = Mv - mv;
}

//=======================================================================
//function : EnlargeBox et ses amis.
//purpose  : 
//=======================================================================

static Handle(Adaptor3d_HSurface) Geometry(TopOpeBRepDS_DataStructure& DStr,
                               const Standard_Integer      ind)
{
  if(ind == 0) return Handle(Adaptor3d_HSurface)();
  if(ind > 0) {
    TopoDS_Face F = TopoDS::Face(DStr.Shape(ind));
    if(F.IsNull()) return Handle(Adaptor3d_HSurface)();
    Handle(BRepAdaptor_HSurface) HS = new BRepAdaptor_HSurface();
    HS->ChangeSurface().Initialize(F,0);
    return HS;
  }
  else{
    Handle(Geom_Surface) S  = DStr.Surface(-ind).Surface();
    if(S.IsNull()) return Handle(Adaptor3d_HSurface)();
    return new GeomAdaptor_HSurface(S);
  }
}

void ChFi3d_SetPointTolerance(TopOpeBRepDS_DataStructure& DStr,
                        const Bnd_Box&              box,
                        const Standard_Integer      IP)
{
  Standard_Real a,b,c,d,e,f,vtol;
  box.Get(a,b,c,d,e,f); 
  d-=a; e-=b; f-=c; 
  d*=d; e*=e; f*=f;
  vtol = sqrt(d + e + f) * 1.5;// on prend un petit rab.
  DStr.ChangePoint(IP).Tolerance(vtol);
}

void ChFi3d_EnlargeBox(const Handle(Geom_Curve)& C,
                   const Standard_Real       wd,
                   const Standard_Real       wf,
                   Bnd_Box&                  box1,
                   Bnd_Box&                  box2)
{
  box1.Add(C->Value(wd));
  box2.Add(C->Value(wf));
}

void ChFi3d_EnlargeBox(const Handle(Adaptor3d_HSurface)& S,
                   const Handle(Geom2d_Curve)&     PC,
                   const Standard_Real             wd,
                   const Standard_Real             wf,
                   Bnd_Box&                        box1,
                   Bnd_Box&                        box2)
{
  Standard_Real u,v;
  PC->Value(wd).Coord(u,v);
  box1.Add(S->Value(u,v));
  PC->Value(wf).Coord(u,v);
  box2.Add(S->Value(u,v));
}

void ChFi3d_EnlargeBox(const TopoDS_Edge&           E,
                   const TopTools_ListOfShape&  LF,
                   const Standard_Real          w,
                   Bnd_Box&                     box)

{
  BRepAdaptor_Curve BC(E);
  box.Add(BC.Value(w));
  TopTools_ListIteratorOfListOfShape It;
  for(It.Initialize(LF); It.More(); It.Next()) {
    TopoDS_Face F = TopoDS::Face(It.Value());
    if(!F.IsNull()) {
      BC.Initialize(E,F);
      box.Add(BC.Value(w));
    }
  }
}

void ChFi3d_EnlargeBox(TopOpeBRepDS_DataStructure&    DStr,
                   const Handle(ChFiDS_Stripe)&   st, 
                   const Handle(ChFiDS_SurfData)& sd,
                   Bnd_Box&                       b1,
                   Bnd_Box&                       b2,
                   const Standard_Boolean         isfirst)
{
  Standard_Real u,v;
  const ChFiDS_CommonPoint& cp1 = sd->Vertex(isfirst,1);
  const ChFiDS_CommonPoint& cp2 = sd->Vertex(isfirst,2);
  b1.Add(cp1.Point());
  b2.Add(cp2.Point());
  const ChFiDS_FaceInterference& fi1 = sd->InterferenceOnS1();
  const ChFiDS_FaceInterference& fi2 = sd->InterferenceOnS2();
  const Handle(Geom_Surface)& S = DStr.Surface(sd->Surf()).Surface();
  const Handle(Geom2d_Curve)& pcs1 = fi1.PCurveOnSurf();
  const Handle(Geom2d_Curve)& pcs2 = fi2.PCurveOnSurf();
  const Handle(Geom_Curve)& c3d1 = DStr.Curve(fi1.LineIndex()).Curve();
  const Handle(Geom_Curve)& c3d2 = DStr.Curve(fi2.LineIndex()).Curve();
  Handle(Adaptor3d_HSurface) F1 = Geometry(DStr,sd->IndexOfS1());
  Handle(Adaptor3d_HSurface) F2 = Geometry(DStr,sd->IndexOfS2());
  Standard_Real p1 = fi1.Parameter(isfirst);
  if(!c3d1.IsNull()) b1.Add(c3d1->Value(p1));
  if(!pcs1.IsNull()) {
    pcs1->Value(p1).Coord(u,v);
    b1.Add(S->Value(u,v));
  }
  if(!F1.IsNull()) {
    const Handle(Geom2d_Curve)& pcf1 = fi1.PCurveOnFace();
    if(!pcf1.IsNull()) {
      pcf1->Value(p1).Coord(u,v);
      b1.Add(F1->Value(u,v));
    }
  }
  Standard_Real p2 = fi2.Parameter(isfirst);
  if(!c3d2.IsNull()) b2.Add(c3d2->Value(p2));
  if(!pcs2.IsNull()) {
    pcs2->Value(p2).Coord(u,v);
    b2.Add(S->Value(u,v));
  }
  if(!F2.IsNull()) {
    const Handle(Geom2d_Curve)& pcf2 = fi2.PCurveOnFace();
    if(!pcf2.IsNull()) {
      pcf2->Value(p2).Coord(u,v);
      b2.Add(F2->Value(u,v));
    }
  }
  if(!st.IsNull()) {
    const Handle(Geom_Curve)& c3d = DStr.Curve(st->Curve(isfirst)).Curve();
    const Handle(Geom2d_Curve)& c2d = st->PCurve(isfirst);
    if(st->Orientation(isfirst) == TopAbs_FORWARD) st->Parameters(isfirst,p1,p2);
    else st->Parameters(isfirst,p2,p1);
    if(!c3d.IsNull()) {
      b1.Add(c3d->Value(p1));
      b2.Add(c3d->Value(p2));
    }
    if(!c2d.IsNull()) {
      c2d->Value(p1).Coord(u,v);
      b1.Add(S->Value(u,v)); 
      c2d->Value(p2).Coord(u,v);
      b2.Add(S->Value(u,v));
    }
  }
}

//=======================================================================
//function : conexfaces
//purpose  : 
//=======================================================================

void ChFi3d_conexfaces(const TopoDS_Edge& E,
                   TopoDS_Face&       F1,
                   TopoDS_Face&       F2,
                   const ChFiDS_Map&  EFMap)
{
  TopTools_ListIteratorOfListOfShape It;
  F1.Nullify();
  F2.Nullify();
  for(It.Initialize(EFMap(E));It.More();It.Next()) {  
    if (F1.IsNull()) {
      F1 = TopoDS::Face(It.Value());
    }
    else {
      F2 = TopoDS::Face(It.Value());
      if(!F2.IsSame(F1) || BRep_Tool::IsClosed(E,F1)) {
      break;
      }
      else F2.Nullify();
    }
  }  
}

//=======================================================================
//function : EdgeState
//purpose  : examun des concavites pour les sommets a 3 aretes.
//=======================================================================

ChFiDS_State ChFi3d_EdgeState(TopoDS_Edge* E,
                        const ChFiDS_Map&  EFMap)
{
  ChFiDS_State sst;
  Standard_Integer i,j;
  TopoDS_Face F[3];
  TopoDS_Face F1,F2,F3,F4,F5,F6;
  ChFi3d_conexfaces(E[0],F1,F2,EFMap);
  ChFi3d_conexfaces(E[1],F3,F4,EFMap);
  ChFi3d_conexfaces(E[2],F5,F6,EFMap);
  if(F1.IsSame(F2)) {
    F[0] = F[1] = F1;
    if(F1.IsSame(F3)) F[2] = F4;
    else F[2] = F3;
  }
  else if(F3.IsSame(F4)) {
    F[0] = F[2] = F3;
    if(F3.IsSame(F1)) F[1] = F2;
    else F[1] = F1;
  }
  else if(F5.IsSame(F6)) {
    F[1] = F[2] = F5;
    if(F5.IsSame(F1)) F[0] = F2;
    else F[0] = F1;
  }
  else{
    if(F1.IsSame(F3) || F1.IsSame(F4)) F[0] = F1;
    else F[0] = F2;
    if(F3.IsSame(F[0])) F[2] = F4;
    else F[2] = F3;
    if(F5.IsSame(F[2])) F[1] = F6;
    else F[1] = F5;
    
  }
  
  if(F[0].IsNull() || F[1].IsNull() || F[2].IsNull()) sst = ChFiDS_FreeBoundary;
  else{
    TopAbs_Orientation o01,o02,o11,o12,o21,o22;
    i=ChFi3d::ConcaveSide(F[0],F[1],E[0],o01,o02);
    i=ChFi3d::ConcaveSide(F[0],F[2],E[1],o11,o12);
    j=ChFi3d::ConcaveSide(F[1],F[2],E[2],o21,o22);
    if(o01==o11 && o02==o21 && o12==o22) sst = ChFiDS_AllSame;
    else if(o12==o22 || i ==10 || j ==10) sst = ChFiDS_OnDiff;
    else sst = ChFiDS_OnSame;
  }
  return sst;
}

//=======================================================================
//function : evalconti
//purpose  : Methode tres rapide, a la limite de l imposture pour
//           coder les regularites CN. Il faudra affiner le traitement
//           quand le bug reviendra.
//=======================================================================

GeomAbs_Shape ChFi3d_evalconti(const TopoDS_Edge& /*E*/,
                         const TopoDS_Face& F1,
                         const TopoDS_Face& F2)
{
  GeomAbs_Shape cont = GeomAbs_G1;
  if(!F1.IsSame(F2)) return cont;
  TopoDS_Face F = F1;
  F.Orientation(TopAbs_FORWARD);
  BRepAdaptor_Surface S(F,Standard_False);
  GeomAbs_SurfaceType typ = S.GetType();
  if(typ != GeomAbs_Cone && 
     typ != GeomAbs_Sphere && 
     typ != GeomAbs_Torus) return cont;
  return GeomAbs_CN;
}

//=======================================================================
//function : KParticular
//purpose  : 
//=======================================================================

Standard_Boolean ChFi3d_KParticular
(const Handle(ChFiDS_Spine)& Spine,
 const Standard_Integer      IE,
 const BRepAdaptor_Surface&  S1,
 const BRepAdaptor_Surface&  S2)
{
#ifdef DEB
  if(ChFi3d_GetcontextFORCEBLEND()) return Standard_False; 
#endif
  Handle(ChFiDS_FilSpine) fs = Handle(ChFiDS_FilSpine)::DownCast(Spine);
  if(!fs.IsNull() && !fs->IsConstant(IE)) return Standard_False;
  const BRepAdaptor_Curve& bc = Spine->CurrentElementarySpine(IE);
  GeomAbs_CurveType ctyp = bc.GetType();
  
  switch (S1.GetType()) {
    
  case GeomAbs_Plane :
    
    switch(S2.GetType()) {
      
    case GeomAbs_Plane :
      
      if (ctyp == GeomAbs_Line) return Standard_True;
      else return Standard_False;
      
    case GeomAbs_Cylinder :
      
      if (ctyp == GeomAbs_Line && 
        (S1.Plane().Axis().Direction().
         IsNormal(S2.Cylinder().Axis().Direction(),
                Precision::Angular())))
      return Standard_True;
      else if (ctyp == GeomAbs_Circle && 
             (S1.Plane().Axis().Direction().
            IsParallel(S2.Cylinder().Axis().Direction(),
                     Precision::Angular())))
      return Standard_True;
      else return Standard_False;
      
    case GeomAbs_Cone :
      
      if (ctyp != GeomAbs_Circle) return Standard_False;
      else if (S1.Plane().Axis().Direction().
             IsParallel(S2.Cone().Axis().Direction(),
                    Precision::Angular()))
      return Standard_True;
      else return Standard_False;
      
    case GeomAbs_Sphere :
      /*
       if (!Constant) return Standard_False;
       else if (S1.Plane().Axis().Direction().
       IsParallel(S2.Sphere().Position().Axis().Direction(),
       Precision::Angular()))
       return Standard_True;
       else return Standard_False;
       */
      return Standard_False;
    case GeomAbs_Torus :
      /*
       if (!Constant) return Standard_False;
       else if (S1.Plane().Axis().Direction().
       IsParallel(S2.Torus().Axis().Direction(),
       Precision::Angular()))
       return Standard_True;
       else return Standard_False;;
       */
      return Standard_False;;
    default:
      
      return Standard_False;
    }
  case GeomAbs_Cylinder :
    
    if (S2.GetType() == GeomAbs_Plane) {
      if (ctyp == GeomAbs_Line && 
        (S2.Plane().Axis().Direction().
         IsNormal(S1.Cylinder().Axis().Direction(),
                Precision::Angular())))
      return Standard_True;
      else if (ctyp == GeomAbs_Circle && 
             (S2.Plane().Axis().Direction().
            IsParallel(S1.Cylinder().Axis().Direction(),
                     Precision::Angular())))
      return Standard_True;
      else return Standard_False;
    }
    else return Standard_False;
    
  case GeomAbs_Cone :
    
    if (S2.GetType() == GeomAbs_Plane) {
      if (ctyp != GeomAbs_Circle) return Standard_False;
      else if (S2.Plane().Axis().Direction().
             IsParallel(S1.Cone().Axis().Direction(),
                    Precision::Angular()))
      return Standard_True;
      else return Standard_False;
      
    }
    else return Standard_False;
    
  case GeomAbs_Sphere :
    /*
       if (S2.GetType() == GeomAbs_Plane) {
       if (!Constant) return Standard_False;
       else if (S2.Plane().Axis().Direction().
       IsParallel(S1.Sphere().Position().Axis().Direction(),
       Precision::Angular()))
       return Standard_True;
       else return Standard_False;
       }
       else return Standard_False;
       */
    return Standard_False;
  case GeomAbs_Torus :
    /*
       if (S2.GetType() == GeomAbs_Plane) {
       if (!Constant) return Standard_False;
       else if (S2.Plane().Axis().Direction().
       IsParallel(S1.Torus().Axis().Direction(),
       Precision::Angular()))
       return Standard_True;
       else return Standard_False;
       }
       else return Standard_False;
       */
    return Standard_False;
  default:
    return Standard_False;
  }
  return Standard_False;
}


//=======================================================================
//function : BoundFac
//purpose  : Resize les bornes d une surface au voisinage de la boite 
//           donnee. Utile pour les intersections dont on connait les
//           extremites. 
//=======================================================================

void ChFi3d_BoundFac(BRepAdaptor_Surface& S,
                 const Standard_Real uumin,
                 const Standard_Real uumax,
                 const Standard_Real vvmin,
                 const Standard_Real vvmax,
                 const Standard_Boolean checknaturalbounds)
{
  ChFi3d_BoundSrf(S.ChangeSurface(), uumin,uumax,vvmin,vvmax,checknaturalbounds);
}
//=======================================================================
//function : ChFi3d_BoundSrf
//purpose  : Resize les bornes d une surface au voisinage de la boite 
//           donnee. Utile pour les intersections dont on connait les
//           extremites. 
//=======================================================================

void ChFi3d_BoundSrf(GeomAdaptor_Surface& S,
                 const Standard_Real uumin,
                 const Standard_Real uumax,
                 const Standard_Real vvmin,
                 const Standard_Real vvmax,
                 const Standard_Boolean checknaturalbounds)
{
  Standard_Real umin = uumin, umax = uumax, vmin = vvmin, vmax = vvmax; 
  Handle(Geom_Surface) surface = S.Surface();
  Handle(Geom_RectangularTrimmedSurface) 
    trs = Handle(Geom_RectangularTrimmedSurface)::DownCast(surface);
  if(!trs.IsNull()) surface = trs->BasisSurface();
  Standard_Real u1,u2,v1,v2;
  surface->Bounds(u1,u2,v1,v2);
  Standard_Real peru=0, perv=0;
  if(surface->IsUPeriodic()) {
    peru = surface->UPeriod();
//     if((umax - umin) > (0.5 * peru)) {
//       umin = uumax - peru;
//       umax = peru - (uumax - uumin);
//     }
  }
  if(surface->IsVPeriodic()) {
    perv = surface->VPeriod();
//     if((vmax - vmin) > (0.5 * perv)) {
//       vmin = vvmax - perv;
//       vmax = perv - (vvmax - vvmin);
//     }
  }
  Standard_Real Stepu = umax - umin;
  Standard_Real Stepv = vmax - vmin;

  //On table sur le fait que la boite uv est non nulle dans
  //une des directions au moins.
  Standard_Real scalu = S.UResolution(1.);
  Standard_Real scalv = S.VResolution(1.);

  Standard_Real step3du = Stepu/scalu; 
  Standard_Real step3dv = Stepv/scalv;
  
  if(step3du > step3dv) Stepv = step3du*scalv;
  if(step3dv > step3du) Stepu = step3dv*scalu;

  if (peru > 0) Stepu = 0.1 * (peru - (umax - umin));
  if (perv > 0) Stepv = 0.1 * (perv - (vmax - vmin));

  Standard_Real uu1 = umin - Stepu;
  Standard_Real uu2 = umax + Stepu;
  Standard_Real vv1 = vmin - Stepv;
  Standard_Real vv2 = vmax + Stepv;
  if(checknaturalbounds) {
    if(!S.IsUPeriodic()) {uu1 = Max(uu1,u1);  uu2 = Min(uu2,u2);}
    if(!S.IsVPeriodic()) {vv1 = Max(vv1,v1);  vv2 = Min(vv2,v2);}
  }
  S.Load(surface,uu1,uu2,vv1,vv2);
}

//=======================================================================
//function : ChFi3d_InterPlaneEdge
//purpose  : 
//=======================================================================

Standard_Boolean  ChFi3d_InterPlaneEdge (Handle(Adaptor3d_HSurface)& Plan,
                               Handle(Adaptor3d_HCurve)&  C,
                               Standard_Real& W,
                               const Standard_Boolean Sens,
                               const Standard_Real tolc)
{ 
  IntCurveSurface_HInter Intersection;
  Standard_Integer isol = 0, nbp ,iip;
  Standard_Real uf = C->FirstParameter(),ul = C->LastParameter();
  Standard_Real CW;
  
  Intersection.Perform(C,Plan);
  
  if(Intersection.IsDone()) {
    nbp = Intersection.NbPoints();
    for (iip = 1; iip <= nbp; iip++) {
      CW = Intersection.Point(iip).W();
      if(C->IsPeriodic())
      CW = ElCLib::InPeriod(CW,uf-tolc,uf-tolc+C->Period());
      if(uf - tolc <= CW && ul + tolc >= CW) {
      if (isol == 0) {
        isol = iip; W = CW;
      }
      else {
        if      ( Sens && CW < W) {W = CW; isol = iip;}
        else if (!Sens && CW > W) {W = CW; isol = iip;}
      }
      }
    }
  }
  if(isol == 0) return Standard_False;
  return Standard_True;
}

//=======================================================================
//function : ExtrSpineCarac
//purpose  : 
//=======================================================================

void ChFi3d_ExtrSpineCarac(const TopOpeBRepDS_DataStructure& DStr,
                     const Handle(ChFiDS_Stripe)& cd,
                     const Standard_Integer i,
                     const Standard_Real p,
                     const Standard_Integer jf,
                     const Standard_Integer sens,
                     gp_Pnt& P,
                     gp_Vec& V,
                     Standard_Real& R) //voir s il ne faudrait pas rajouter D1,D2 et DR
{
  // Attention pour les surfaces approximees on assume que
  // le parametrage de la pcurve est le meme que celui de 
  // l elspine qui a servi a la construction.
  const Handle(Geom_Surface)& fffil = 
    DStr.Surface(cd->SetOfSurfData()->Value(i)->Surf()).Surface();
  gp_Pnt2d pp = cd->SetOfSurfData()->Value(i)->Interference(jf).
    PCurveOnSurf()->Value(p);
  GeomAdaptor_Surface gs(fffil);
  P = fffil->Value(pp.X(),pp.Y());
  gp_Pnt Pbid; gp_Vec Vbid;
  switch (gs.GetType()) {
  case GeomAbs_Cylinder :
    {
      gp_Cylinder cyl = gs.Cylinder();
      R = cyl.Radius();
      ElSLib::D1(pp.X(),pp.Y(),cyl,Pbid,Vbid,V);
    }
    break;
  case GeomAbs_Torus :
    {
      gp_Torus tor = gs.Torus();
      R = tor.MinorRadius();
      ElSLib::D1(pp.X(),pp.Y(),tor,Pbid,V,Vbid);
    }
    break;
  default:
    { Standard_Integer nbelspine;
      const Handle(ChFiDS_Spine)& sp = cd->Spine();
      Handle(ChFiDS_FilSpine) fsp = Handle(ChFiDS_FilSpine)::DownCast(sp);
      nbelspine=sp->NbEdges();
      Handle(ChFiDS_HElSpine) hels;
      if   (nbelspine==1) hels = sp->ElSpine(1);
      else hels = sp->ElSpine(p);
      if(fsp->IsConstant()) { R = fsp->Radius(); }
      else { R = fsp->Law(hels)->Value(p); }
      hels->D1(p,Pbid,V);
    }
    break;
  }
  V.Normalize();
  if(sens == 1) V.Reverse();
}

//=======================================================================
//function : ChFi3d_CircularSpine
//purpose  : Calcule une ligne guide ciculaire pour le du coin a partir 
//           des points et vecteurs tangents calcules aux extremites des
//           lignes guides des conges deb et fin.
//=======================================================================
Handle(Geom_Circle) ChFi3d_CircularSpine(Standard_Real&      WFirst,
                               Standard_Real&      WLast,
                               const gp_Pnt&       Pdeb,
                               const gp_Vec&       Vdeb,
                               const gp_Pnt&       Pfin,
                               const gp_Vec&       Vfin,
                               const Standard_Real rad)
{
  gp_Circ ccc;
  gp_Pln Pl1(Pdeb,gp_Dir(Vdeb)),Pl2(Pfin,gp_Dir(Vfin));
  IntAna_QuadQuadGeo LInt (Pl1,Pl2,Precision::Angular(),
                     Precision::Confusion());
  gp_Lin li;
  if (LInt.IsDone()) {
    li = LInt.Line(1);
    gp_Pnt cendeb = ElCLib::Value(ElCLib::Parameter(li,Pdeb),li);
    gp_Pnt cenfin = ElCLib::Value(ElCLib::Parameter(li,Pfin),li);
    gp_Vec vvdeb(cendeb,Pdeb);
    gp_Vec vvfin(cenfin,Pfin);
    gp_Dir dddeb(vvdeb);
    gp_Dir ddfin(vvfin);
    if(Vdeb.Crossed(vvdeb).Dot(Vfin.Crossed(vvfin)) > 0.) {
      return Handle(Geom_Circle)();
    }
    gp_Ax2 circax2(cendeb,dddeb^ddfin,dddeb);
    ccc.SetPosition(circax2);
    ccc.SetRadius(rad);
    WFirst = 0.;
    WLast = dddeb.Angle(ddfin);
    return new Geom_Circle(ccc);
  }
#ifdef DEB
  else {
    cout <<"Construction de la ligne guide du coin impossible"<<endl;
  }
#endif
  return Handle(Geom_Circle)();
}

//=======================================================================
//function : ChFi3d_Spine
//purpose  : Calcule les poles de la ligne guide du coin a partir des
//           points et vecteurs tangents calcules aux extremites des
//           lignes guides des conges deb et fin.
//=======================================================================

Handle(Geom_BezierCurve) ChFi3d_Spine(const gp_Pnt&       pd,
                              gp_Vec&             vd,
                              const gp_Pnt&       pf,
                              gp_Vec&             vf,
                              const Standard_Real R)
{     
  TColgp_Array1OfPnt pol(1,4);
  Standard_Real fac = 0.5 * tan((PI-vd.Angle(vf))/2.);
  pol(1) = pd;
  vd.Multiply(fac*R);
  pol(2).SetCoord(pd.X()+vd.X(),pd.Y()+vd.Y(),pd.Z()+vd.Z());
  pol(4) = pf;
  vf.Multiply(fac*R);
  pol(3).SetCoord(pf.X()+vf.X(),pf.Y()+vf.Y(),pf.Z()+vf.Z());
  return new Geom_BezierCurve(pol);
}

//=======================================================================
//function : IsInFront
//purpose  : regarde si les surfdata i1 et i2 sont en vis a vis
//=======================================================================

Standard_Boolean ChFi3d_IsInFront(TopOpeBRepDS_DataStructure& DStr,
                          const Handle(ChFiDS_Stripe)& cd1, 
                          const Handle(ChFiDS_Stripe)& cd2,
                          const Standard_Integer i1,
                          const Standard_Integer i2,
                          const Standard_Integer sens1,
                          const Standard_Integer sens2,
                          Standard_Real& p1,
                          Standard_Real& p2,
                          TopoDS_Face& face,
                          Standard_Boolean& sameside,
                          Standard_Integer& jf1,
                          Standard_Integer& jf2,
                          Standard_Boolean& visavis,
                          const TopoDS_Vertex& Vtx,
                          const Standard_Boolean Check2dDistance,
                          const Standard_Boolean enlarge)
{
  Standard_Boolean isf1 = (sens1 == 1), isf2 = (sens2 == 1);
  const Handle(ChFiDS_SurfData)& fd1 = cd1->SetOfSurfData()->Value(i1);
  const Handle(ChFiDS_SurfData)& fd2 = cd2->SetOfSurfData()->Value(i2);

  TopAbs_Orientation Or,OrSave1,OrSave2,OrFace1,OrFace2;
  visavis = Standard_False;
#ifndef DEB
  Standard_Real u1 = 0.,u2 = 0.; 
  Standard_Boolean ss = 0,ok = 0;
  Standard_Integer j1 = 0,j2 = 0;
#else
  Standard_Real u1,u2; 
  Standard_Boolean ss,ok = 0;
  Standard_Integer j1,j2;
#endif
  TopoDS_Face ff;
  if(fd1->IndexOfS1() == fd2->IndexOfS1()) { 
    jf1 = 1; jf2 = 1; 
    face = TopoDS::Face(DStr.Shape(fd1->Index(jf1)));
    OrSave1 = cd1->Orientation(jf1);
    Or = OrFace1 = face.Orientation();
    OrSave2 = cd2->Orientation(jf2);
    OrFace2 = DStr.Shape(fd2->Index(jf2)).Orientation();
    visavis = Standard_True;
    sameside = ChFi3d::SameSide(Or,OrSave1,OrSave2,OrFace1,OrFace2);
    // On ne detrompe plus a l'aide des parametres de l'autre cote ca posait des problemes.
    Standard_Integer kf1 = jf1, kf2 = jf2;
    Standard_Real pref1 = fd1->Interference(kf1).Parameter(isf1);
    Standard_Real pref2 = fd2->Interference(kf2).Parameter(isf2);
    gp_Pnt2d P2d;
    if (Check2dDistance)
      P2d = BRep_Tool::Parameters( Vtx, face );
    if(ChFi3d_IntTraces(fd1,pref1,p1,jf1,sens1,fd2,pref2,p2,jf2,sens2,P2d,Check2dDistance,enlarge)) {
      u1 = p1; u2 = p2; ss = sameside; j1 = jf1; j2 = jf2; ff = face;
      ok = 1;
    }
  }
  if(fd1->IndexOfS2() == fd2->IndexOfS1()) { 
    jf1 = 2; jf2 = 1; 
    face = TopoDS::Face(DStr.Shape(fd1->Index(jf1)));
    OrSave1 = cd1->Orientation(jf1);
    Or = OrFace1 = face.Orientation();
    OrSave2 = cd2->Orientation(jf2);
    OrFace2 = DStr.Shape(fd2->Index(jf2)).Orientation();
    visavis = Standard_True;
    sameside = ChFi3d::SameSide(Or,OrSave1,OrSave2,OrFace1,OrFace2);
    // On ne detrompe plus a l'aide des parametres de l'autre cote ca posait des problemes.
    Standard_Integer kf1 = jf1, kf2 = jf2;
    Standard_Real pref1 = fd1->Interference(kf1).Parameter(isf1);
    Standard_Real pref2 = fd2->Interference(kf2).Parameter(isf2);
    gp_Pnt2d P2d;
    if (Check2dDistance)
      P2d = BRep_Tool::Parameters( Vtx, face );
    if(ChFi3d_IntTraces(fd1,pref1,p1,jf1,sens1,fd2,pref2,p2,jf2,sens2,P2d,Check2dDistance,enlarge)) {
      Standard_Boolean restore = 
      ok && ((j1 == jf1 && sens1*(p1 - u1) > 0.) || 
             (j2 == jf2 && sens2*(p2 - u2) > 0.));
      ok = 1;
      if(restore) {
      p1 = u1; p2 = u2; sameside = ss; jf1 = j1; jf2 = j2; face = ff;
      }
      else {
      u1 = p1; u2 = p2; ss = sameside; j1 = jf1; j2 = jf2; ff = face;
      }
    }
    //on rajoute ces re-initialisations au cas ou p1,... auraient pris des valeurs fausses
    else if (ok) {
      p1 = u1; p2 = u2; sameside = ss; jf1 = j1; jf2 = j2; face = ff;
    }
  }
  if(fd1->IndexOfS1() == fd2->IndexOfS2()) { 
    jf1 = 1; jf2 = 2; 
    face = TopoDS::Face(DStr.Shape(fd1->Index(jf1)));
    OrSave1 = cd1->Orientation(jf1);
    Or = OrFace1 = face.Orientation();
    OrSave2 = cd2->Orientation(jf2);
    OrFace2 = DStr.Shape(fd2->Index(jf2)).Orientation();
    visavis = Standard_True;
    sameside = ChFi3d::SameSide(Or,OrSave1,OrSave2,OrFace1,OrFace2);
    // On ne detrompe plus a l'aide des parametres de l'autre cote.
    Standard_Integer kf1 = jf1, kf2 = jf2;
    Standard_Real pref1 = fd1->Interference(kf1).Parameter(isf1);
    Standard_Real pref2 = fd2->Interference(kf2).Parameter(isf2);
    gp_Pnt2d P2d;
    if (Check2dDistance)
      P2d = BRep_Tool::Parameters( Vtx, face );
    if(ChFi3d_IntTraces(fd1,pref1,p1,jf1,sens1,fd2,pref2,p2,jf2,sens2,P2d,Check2dDistance,enlarge)) {
      Standard_Boolean restore = 
      ok && ((j1 == jf1 && sens1*(p1 - u1) > 0.) || 
             (j2 == jf2 && sens2*(p2 - u2) > 0.));
      ok = 1;
      if(restore) {
      p1 = u1; p2 = u2; sameside = ss; jf1 = j1; jf2 = j2; face = ff;
      }
      else {
      u1 = p1; u2 = p2; ss = sameside; j1 = jf1; j2 = jf2; ff = face;
      }
    }
    //on rajoute ces re-initialisations au cas ou p1,... auraient pris des valeurs fausses
    else if (ok) {
      p1 = u1; p2 = u2; sameside = ss; jf1 = j1; jf2 = j2; face = ff;
    }
  }
  if(fd1->IndexOfS2() == fd2->IndexOfS2()) { 
    jf1 = 2; jf2 = 2; 
    face = TopoDS::Face(DStr.Shape(fd1->Index(jf1)));
    OrSave1 = cd1->Orientation(jf1);
    Or = OrFace1 = face.Orientation();
    OrSave2 = cd2->Orientation(jf2);
    OrFace2 = DStr.Shape(fd2->Index(jf2)).Orientation();
    visavis = Standard_True;
    sameside = ChFi3d::SameSide(Or,OrSave1,OrSave2,OrFace1,OrFace2);
    // On ne detrompe plus a l'aide des parametres de l'autre cote.
    Standard_Integer kf1 = jf1, kf2 = jf2;
    Standard_Real pref1 = fd1->Interference(kf1).Parameter(isf1);
    Standard_Real pref2 = fd2->Interference(kf2).Parameter(isf2);
    gp_Pnt2d P2d;
    if (Check2dDistance)
      P2d = BRep_Tool::Parameters( Vtx, face );
    if(ChFi3d_IntTraces(fd1,pref1,p1,jf1,sens1,fd2,pref2,p2,jf2,sens2,P2d,Check2dDistance,enlarge)) {
      Standard_Boolean restore = 
      ok && ((j1 == jf1 && sens1*(p1 - u1) > 0.) || 
             (j2 == jf2 && sens2*(p2 - u2) > 0.));
      ok = 1;
      if(restore) {
      p1 = u1; p2 = u2; sameside = ss; jf1 = j1; jf2 = j2; face = ff;
      }
      else {
      u1 = p1; u2 = p2; ss = sameside; j1 = jf1; j2 = jf2; ff = face;
      }
    }
    //on rajoute ces re-initialisations au cas ou p1,... auraient pris des valeurs fausses
    else if (ok) {
      p1 = u1; p2 = u2; sameside = ss; jf1 = j1; jf2 = j2; face = ff;
    }
  }
  return ok;

}

static Standard_Real recadre(const Standard_Real p,
                       const Standard_Real ref,
                       const Standard_Integer sens,
                       const Standard_Real first,
                       const Standard_Real last)
{
  Standard_Real pp = p;
  if(sens > 0) pp -= (last - first);
  else pp += (last - first);
  if(Abs(pp - ref) < Abs(p - ref)) return pp;
  return p;
}

Standard_Boolean ChFi3d_IntTraces(const Handle(ChFiDS_SurfData)& fd1,
                          const Standard_Real            pref1,
                          Standard_Real&                 p1,
                          const Standard_Integer         jf1,
                          const Standard_Integer         sens1,
                          const Handle(ChFiDS_SurfData)& fd2,
                          const Standard_Real            pref2,
                          Standard_Real&                 p2,
                          const Standard_Integer         jf2,
                          const Standard_Integer         sens2,
                          const gp_Pnt2d&                RefP2d,
                          const Standard_Boolean         Check2dDistance,
                          const Standard_Boolean         enlarge)
{
  Geom2dAdaptor_Curve C1;
  Geom2dAdaptor_Curve C2;
  // on aggrandit les pcurves pour etre sur qu'il y a intersection
  // par ailleurs on prend la totalite des courbes periodiques, on
  // triera les points sur celles-ci avec un critere specifique.

  Standard_Real first,last,delta = 0.;
  first = fd1->Interference(jf1).FirstParameter();
  last = fd1->Interference(jf1).LastParameter();
  if ((last-first) < Precision::PConfusion())
    return Standard_False;
  if(enlarge) delta = Min(0.1,0.05*(last-first));
  Handle(Geom2d_Curve) pcf1 = fd1->Interference(jf1).PCurveOnFace();
  if(pcf1.IsNull()) return Standard_False;
  Standard_Boolean isper1 = pcf1->IsPeriodic();
  if(isper1) {
    Handle(Geom2d_TrimmedCurve) tr1 = Handle(Geom2d_TrimmedCurve)::DownCast(pcf1);
    if(!tr1.IsNull()) pcf1 = tr1->BasisCurve();
    C1.Load(pcf1);
  }
  else C1.Load(pcf1,first-delta,last+delta);
  Standard_Real first1 = pcf1->FirstParameter(), last1 = pcf1->LastParameter();
  
  first = fd2->Interference(jf2).FirstParameter();
  last = fd2->Interference(jf2).LastParameter();
  if ((last-first) < Precision::PConfusion())
    return Standard_False;
  if(enlarge) delta = Min(0.1,0.05*(last-first));
  Handle(Geom2d_Curve) pcf2 = fd2->Interference(jf2).PCurveOnFace();
  if(pcf2.IsNull()) return Standard_False;
  Standard_Boolean isper2 = pcf2->IsPeriodic();
  if(isper2) {
    Handle(Geom2d_TrimmedCurve) tr2 = Handle(Geom2d_TrimmedCurve)::DownCast(pcf2);
    if(!tr2.IsNull()) pcf2 = tr2->BasisCurve();
    C2.Load(pcf2);
  }
  else C2.Load(fd2->Interference(jf2).PCurveOnFace(),first-delta,last+delta);
  Standard_Real first2 = pcf2->FirstParameter(), last2 = pcf2->LastParameter();
    
  IntRes2d_IntersectionPoint int2d;
  Geom2dInt_GInter Intersection;
  Standard_Integer nbpt,nbseg;
  gp_Pnt2d p2d;
  if(fd1->Interference(jf1).PCurveOnFace() == fd2->Interference(jf2).PCurveOnFace()) {
    Intersection.Perform(C1,
                   Precision::PIntersection(),
                   Precision::PIntersection());
  }
  else{
    Intersection.Perform(C1,C2,
                   Precision::PIntersection(),
                   Precision::PIntersection());
  }
  if (Intersection.IsDone()) {
    if (!Intersection.IsEmpty()) {
      nbseg = Intersection.NbSegments();
      if ( nbseg > 0 ) { 
#ifdef DEB
      cout<< "segments d intersection sur les restrictions"<<endl; 
#endif
      }
      nbpt = Intersection.NbPoints();
      if ( nbpt >= 1 ) {
      // On manque de billes pour trier les points trouves de facon
        // rigoureuse. On adopte donc deux criteres differents un peu
        // arbitraires :
      // - courbes periodiques : le plus proche de la borne.
      // - courbes non periodiques : le plus a gauche sur les 2 courbes
      //                             modulo sens1 et sens2
      int2d = Intersection.Point(1);
      p2d = int2d.Value();
      p1 = int2d.ParamOnFirst();
      p2 = int2d.ParamOnSecond();
      if(isper1) p1 = recadre(p1,pref1,sens1,first1,last1);
      if(isper2) p2 = recadre(p2,pref2,sens2,first2,last2);
      for(Standard_Integer i = 2; i<=nbpt; i++) {
        int2d = Intersection.Point(i);
        if(isper1) {
          Standard_Real pp1 = int2d.ParamOnFirst();
          pp1 = recadre(pp1,pref1,sens1,first1,last1);
          if((Abs(pp1 - pref1) < Abs(p1 - pref1))) {
            p1 = pp1;
            p2 = int2d.ParamOnSecond();
            p2d = int2d.Value();
          }
//  Modified by skv - Mon Jun 16 15:51:21 2003 OCC615 Begin
          else if (Check2dDistance &&
                 RefP2d.Distance(int2d.Value()) < RefP2d.Distance(p2d)) {
            Standard_Real pp2 = int2d.ParamOnSecond();

            if(isper2)
            pp2 = recadre(pp2,pref2,sens2,first2,last2);

            p1  = pp1;
            p2  = pp2;
            p2d = int2d.Value();
          }
//  Modified by skv - Mon Jun 16 15:51:22 2003 OCC615 End
        }
        else if(isper2) {
          Standard_Real pp2 = int2d.ParamOnSecond();
          pp2 = recadre(pp2,pref2,sens2,first2,last2);
          if((Abs(pp2 - pref2) < Abs(p2 - pref2))) {
            p2 = pp2;
            p1 = int2d.ParamOnFirst();
            p2d = int2d.Value();
          }
//  Modified by skv - Mon Jun 16 15:51:21 2003 OCC615 Begin
          else if (Check2dDistance &&
                 RefP2d.Distance(int2d.Value()) < RefP2d.Distance(p2d)) {
            Standard_Real pp1 = int2d.ParamOnFirst();

            if(isper1)
            pp1 = recadre(pp1,pref1,sens1,first1,last1);

            p1  = pp1;
            p2  = pp2;
            p2d = int2d.Value();
          }
//  Modified by skv - Mon Jun 16 15:51:22 2003 OCC615 End
        }
        else if(((int2d.ParamOnFirst() - p1)*sens1 < 0.) &&
              ((int2d.ParamOnSecond() - p2)*sens2 < 0.)) {
          p1 = int2d.ParamOnFirst();
          p2 = int2d.ParamOnSecond();
          p2d = int2d.Value();
        }
        else if((Abs(int2d.ParamOnFirst() - pref1) < Abs(p1 - pref1)) &&
              (Abs(int2d.ParamOnSecond() - pref2) < Abs(p2 - pref2))) {
          p1 = int2d.ParamOnFirst();
          p2 = int2d.ParamOnSecond();
          p2d = int2d.Value();
        }
        else if (Check2dDistance && RefP2d.Distance(int2d.Value()) < RefP2d.Distance(p2d))
          {
            p1 = int2d.ParamOnFirst();
            p2 = int2d.ParamOnSecond();
            p2d = int2d.Value();
          }
      }
      return Standard_True; 
      }
      return Standard_False; 
    }
    else { return Standard_False; }
  }
  else { return Standard_False; }
} 

//=======================================================================
//function : Coefficient
//purpose  : 
//=======================================================================
void ChFi3d_Coefficient(const gp_Vec& V3d,
                  const gp_Vec& D1u,
                  const gp_Vec& D1v,
                  Standard_Real& DU,
                  Standard_Real& DV) 
{
  Standard_Real AA = D1u.SquareMagnitude();
  Standard_Real BB = D1u.Dot(D1v);
  Standard_Real CC = D1v.SquareMagnitude();
  Standard_Real DD = D1u.Dot(V3d);
  Standard_Real EE = D1v.Dot(V3d);
  Standard_Real Delta = AA*CC-BB*BB;
  DU = (DD*CC-EE*BB)/Delta;
  DV = (AA*EE-BB*DD)/Delta;
}

//=======================================================================
//function : ReparamPcurv
//purpose  : Dans le cas ou la pcurve est une BSpline on verifie 
//           ses parametres et on la reparametre eventuellement.
//=======================================================================

void ChFi3d_ReparamPcurv(const Standard_Real Uf, 
                   const Standard_Real Ul, 
                   Handle(Geom2d_Curve)& Pcurv) 
{
  if(Pcurv.IsNull()) return;
  Standard_Real upcf = Pcurv->FirstParameter();
  Standard_Real upcl = Pcurv->LastParameter();
  Handle(Geom2d_Curve) basis = Pcurv;
  Handle(Geom2d_TrimmedCurve) trpc = Handle(Geom2d_TrimmedCurve)::DownCast(Pcurv);
  if(!trpc.IsNull()) basis = trpc->BasisCurve();
  Handle(Geom2d_BSplineCurve) pc = Handle(Geom2d_BSplineCurve)::DownCast(basis);
  if(pc.IsNull()) return;
  if(Abs(upcf - pc->FirstParameter()) > Precision::PConfusion() ||
     Abs(upcl - pc->LastParameter()) > Precision::PConfusion()) {
    pc->Segment(upcf,upcl);
  }
  if(Abs(Uf - pc->FirstParameter()) > Precision::PConfusion() ||
     Abs(Ul - pc->LastParameter()) > Precision::PConfusion()) {
    TColgp_Array1OfPnt2d pol(1,pc->NbPoles());
    pc->Poles(pol);
    TColStd_Array1OfReal kn(1,pc->NbKnots());
    pc->Knots(kn);
    TColStd_Array1OfInteger mu(1,pc->NbKnots());
    pc->Multiplicities(mu);
    Standard_Integer deg = pc->Degree();
    BSplCLib::Reparametrize(Uf,Ul,kn);
    pc = new Geom2d_BSplineCurve(pol,kn,mu,deg);
  }
  Pcurv = pc;
}

//=======================================================================
//function : ProjectPCurv
//purpose  : Calcul la pcurve correspondant a une ligne d intersection
//           3d. Ne doit etre appele que dans les cas analytiques.
//=======================================================================

void ChFi3d_ProjectPCurv(const Handle(Adaptor3d_HCurve)&   HCg, 
                   const Handle(Adaptor3d_HSurface)& HSg, 
                   Handle(Geom2d_Curve)&           Pcurv,
                   const Standard_Real             tol,
                   Standard_Real&                  tolreached) 
{
  if (HSg->GetType() != GeomAbs_BezierSurface &&
      HSg->GetType() != GeomAbs_BSplineSurface) {
    
    ProjLib_ProjectedCurve Projc (HSg,HCg,tol);
    tolreached = Projc.GetTolerance();
    switch (Projc.GetType()) {
    case GeomAbs_Line : 
      {
      Pcurv = new Geom2d_Line(Projc.Line());
      }
      break;
    case GeomAbs_Circle :
      {
      Pcurv = new Geom2d_Circle(Projc.Circle());
      }
      break;
    case GeomAbs_Ellipse :
      {
      Pcurv = new Geom2d_Ellipse(Projc.Ellipse());
      }
      break;
    case GeomAbs_Hyperbola :
      {
      Pcurv = new Geom2d_Hyperbola(Projc.Hyperbola());
      }
      break;
    case GeomAbs_Parabola :
      {
      Pcurv = new Geom2d_Parabola(Projc.Parabola());
      }
      break;
    case GeomAbs_BezierCurve :
      {
      Pcurv = Projc.Bezier(); 
      }
      break;
    case GeomAbs_BSplineCurve :
      {
      Pcurv = Projc.BSpline();
      }
      break;
    default:
      Standard_NotImplemented::Raise("echec approximation de la pcurve ");
    }
  }
}

//=======================================================================
//function : CheckSameParameter
//purpose  : Controle a posteriori que sameparameter a bien fait son boulot
//=======================================================================

Standard_Boolean ChFi3d_CheckSameParameter 
(const Handle(Adaptor3d_HCurve)&   C3d,
 Handle(Geom2d_Curve)&           Pcurv,
 const Handle(Adaptor3d_HSurface)& S,
 const Standard_Real             tol3d,
 Standard_Real&                  tolreached)
{
  tolreached = 0.;
  Standard_Real f = C3d->FirstParameter();
  Standard_Real l = C3d->LastParameter();
  Standard_Integer nbp = 45;
  Standard_Real step = 1./(nbp -1);
  for(Standard_Integer i = 0; i < nbp; i++) {
    Standard_Real t,u,v;
    t = step * i;
    t = (1-t) * f + t * l;
    Pcurv->Value(t).Coord(u,v);
    gp_Pnt pS = S->Value(u,v);
    gp_Pnt pC = C3d->Value(t);
    Standard_Real d2 = pS.SquareDistance(pC);
    tolreached = Max(tolreached,d2);
  }
  tolreached = sqrt(tolreached);
  if(tolreached > tol3d) {
    tolreached *= 2.;
    return Standard_False;
  }
  tolreached *= 2.;
  tolreached = Max(tolreached,Precision::Confusion());
  return Standard_True;
}

//=======================================================================
//function : SameParameter
//purpose  : Encapsulation de Sameparameter
//=======================================================================

Standard_Boolean ChFi3d_SameParameter(const Handle(Adaptor3d_HCurve)&   C3d,
                              Handle(Geom2d_Curve)&           Pcurv,
                              const Handle(Adaptor3d_HSurface)& S,
                              const Standard_Real             tol3d,
                              Standard_Real&                  tolreached)
{
  if(ChFi3d_CheckSameParameter(C3d,Pcurv,S,tol3d,tolreached)) return Standard_True;
  Approx_SameParameter sp(C3d,Pcurv,S,tol3d);
  if(sp.IsDone() && !sp.IsSameParameter()) Pcurv = sp.Curve2d();
  else if(!sp.IsDone() && !sp.IsSameParameter()) {
#ifdef DEB
    cout<<"echec SameParameter"<<endl;
#endif  
    return Standard_False;
  }
  tolreached = sp.TolReached();
  if(tolreached > 1.1*tol3d) {
#ifdef DEB
    cout<<"SameParameter : Tol non atteinte!!!"<<endl;
    cout<<"tol visee : "<<tol3d<<" tol obtenue : "<<tolreached<<endl;
#endif  
  }
  return Standard_True;
}
  
//=======================================================================
//function : SameParameter
//purpose  : Encapsulation de Sameparameter
//=======================================================================

Standard_Boolean ChFi3d_SameParameter(const Handle(Geom_Curve)&   C3d,
                              Handle(Geom2d_Curve)&       Pcurv,
                              const Handle(Geom_Surface)& S,
                              const Standard_Real         Pardeb,
                              const Standard_Real         Parfin,
                              const Standard_Real         tol3d,
                              Standard_Real&              tolreached)
{
  static Handle(GeomAdaptor_HSurface) hs = new GeomAdaptor_HSurface();
  static Handle(GeomAdaptor_HCurve) hc = new GeomAdaptor_HCurve();
  hs->ChangeSurface().Load(S);
  hc->ChangeCurve().Load(C3d,Pardeb,Parfin);
  return ChFi3d_SameParameter(hc,Pcurv,hs,tol3d,tolreached);
}

//=======================================================================
//function : ComputePCurv 
//purpose  : Calcule une droite eventuellement sous forme de BSpline 
//           pour garantir le range et le parametrage identique
//           a une courbe 3d de reference.
//=======================================================================

void ChFi3d_ComputePCurv(const Handle(Adaptor3d_HCurve)&   C3d,
                   const gp_Pnt2d&                 UV1,
                   const gp_Pnt2d&                 UV2,
                   Handle(Geom2d_Curve)&           Pcurv,
                   const Handle(Adaptor3d_HSurface)& S,
                   const Standard_Real             Pardeb,
                   const Standard_Real             Parfin,
                   const Standard_Real             tol3d,
                   Standard_Real&                  tolreached,
                   const Standard_Boolean          reverse)
{
  ChFi3d_ComputePCurv(UV1,UV2,Pcurv,Pardeb,Parfin,reverse);

#ifdef DEB   
  OSD_Chronometer ch;
  ChFi3d_InitChron(ch);//init perf  pour ChFi3d_SameParameter 
#endif

  ChFi3d_SameParameter(C3d,Pcurv,S,tol3d,tolreached);

#ifdef DEB    
  ChFi3d_ResultChron(ch ,t_sameparam); // result perf pour ChFi3d_SameParameter
#endif 
 
}

//=======================================================================
//function : ComputePCurv 
//purpose  : Calcule une droite eventuellement sous forme de BSpline 
//           pour garantir le range et le parametrage identique
//           a une courbe 3d de reference.
//=======================================================================

void ChFi3d_ComputePCurv(const Handle(Geom_Curve)&   C3d,
                   const gp_Pnt2d&             UV1,
                   const gp_Pnt2d&             UV2,
                   Handle(Geom2d_Curve)&       Pcurv,
                   const Handle(Geom_Surface)& S,
                   const Standard_Real         Pardeb,
                   const Standard_Real         Parfin,
                   const Standard_Real         tol3d,
                   Standard_Real&              tolreached,
                   const Standard_Boolean      reverse)
{
  static Handle(GeomAdaptor_HSurface) hs = new GeomAdaptor_HSurface();
  static Handle(GeomAdaptor_HCurve) hc = new GeomAdaptor_HCurve();
  hs->ChangeSurface().Load(S);
  hc->ChangeCurve().Load(C3d,Pardeb,Parfin);
  ChFi3d_ComputePCurv(hc,UV1,UV2,Pcurv,hs,Pardeb,Parfin,tol3d,tolreached,reverse);
}

//=======================================================================
//function : ComputePCurv 
//purpose  : Calcule une droite eventuellement sous forme de BSpline 
//           pour garantir le range.
//=======================================================================

void ChFi3d_ComputePCurv(const gp_Pnt2d& UV1,
                   const gp_Pnt2d& UV2,
                   Handle(Geom2d_Curve)& Pcurv,
                   const Standard_Real Pardeb,
                   const Standard_Real Parfin,
                   const Standard_Boolean reverse)
{
  
  Standard_Real tol = Precision::PConfusion();
  gp_Pnt2d p1,p2;
  if (!reverse) { 
    p1 = UV1;
    p2 = UV2;
  }
  else {
    p1 = UV2;
    p2 = UV1;
  }
  
  if (Abs(p1.X()-p2.X()) <= tol &&
      Abs((p2.Y()-p1.Y())-(Parfin-Pardeb)) <= tol) {
    gp_Pnt2d ppp(p1.X(),p1.Y()-Pardeb);
    Pcurv = new Geom2d_Line(ppp,gp::DY2d());
  }
  else if (Abs(p1.X()-p2.X()) <= tol &&
         Abs((p1.Y()-p2.Y())-(Parfin-Pardeb)) <= tol) {
    gp_Pnt2d ppp(p1.X(),p1.Y()+Pardeb);
    Pcurv = new Geom2d_Line(ppp,gp::DY2d().Reversed());
  }
  else if (Abs(p1.Y()-p2.Y()) <= tol &&
         Abs((p2.X()-p1.X())-(Parfin-Pardeb)) <= tol) {
    gp_Pnt2d ppp(p1.X()-Pardeb,p1.Y());
    Pcurv = new Geom2d_Line(ppp,gp::DX2d());
  }
  else if (Abs(p1.Y()-p2.Y()) <= tol &&
         Abs((p1.X()-p2.X())-(Parfin-Pardeb)) <= tol) {
    gp_Pnt2d ppp(p1.X()+Pardeb,p1.Y());
    Pcurv = new Geom2d_Line(ppp,gp::DX2d().Reversed());
  }
  else{
    TColgp_Array1OfPnt2d p(1,2);
    TColStd_Array1OfReal k(1,2);
    TColStd_Array1OfInteger m(1,2);
    m.Init(2);
    k(1) = Pardeb;
    k(2) = Parfin;
    p(1) = p1;
    p(2) = p2;
    Pcurv = new Geom2d_BSplineCurve(p,k,m,1);
  }
  Pcurv = new Geom2d_TrimmedCurve(Pcurv,Pardeb,Parfin);
}

Handle(GeomFill_Boundary) ChFi3d_mkbound
(const Handle(Adaptor3d_HSurface)& Fac,
 Handle(Geom2d_Curve)& curv, 
 const Standard_Integer sens1,
 const gp_Pnt2d& pfac1,
 const gp_Vec2d& vfac1,
 const Standard_Integer sens2,
 const gp_Pnt2d& pfac2,
 const gp_Vec2d& vfac2,
 const Standard_Real t3d,
 const Standard_Real ta)
{
  gp_Dir2d v1(vfac1);
  if(sens1 == 1) v1.Reverse();
  gp_Dir2d v2(vfac2);
  if(sens2 == 1) v2.Reverse();
  curv = ChFi3d_BuildPCurve(Fac,pfac1,v1,pfac2,v2,Standard_False);
  return ChFi3d_mkbound(Fac,curv,t3d,ta);
}

Handle(GeomFill_Boundary) ChFi3d_mkbound
(const Handle(Adaptor3d_HSurface)& Surf,
 Handle(Geom2d_Curve)& curv,
 const Standard_Integer sens1,
 const gp_Pnt2d& p1,
 gp_Vec&   v1,
 const Standard_Integer sens2,
 const gp_Pnt2d& p2,
 gp_Vec& v2,
 const Standard_Real t3d,
 const Standard_Real ta)
{
  if(sens1 == 1) v1.Reverse();
  if(sens2 == 1) v2.Reverse();
  curv = ChFi3d_BuildPCurve(Surf,p1,v1,p2,v2);
  return ChFi3d_mkbound(Surf,curv,t3d,ta);
}

Handle(GeomFill_Boundary) ChFi3d_mkbound
(const Handle(Geom_Surface)& s,
 const gp_Pnt2d& p1,
 const gp_Pnt2d& p2,
 const Standard_Real t3d,
 const Standard_Real ta,
 const Standard_Boolean isfreeboundary)
{
  Handle(GeomAdaptor_HSurface) HS = new GeomAdaptor_HSurface();
  HS->ChangeSurface().Load(s);
  return ChFi3d_mkbound(HS,p1,p2,t3d,ta,isfreeboundary);
}

Handle(GeomFill_Boundary) ChFi3d_mkbound
(const Handle(Adaptor3d_HSurface)& HS,
 const gp_Pnt2d& p1,
 const gp_Pnt2d& p2,
 const Standard_Real t3d,
 const Standard_Real ta,
 const Standard_Boolean isfreeboundary)
{
  TColgp_Array1OfPnt2d pol(1,2);
  pol(1)=p1;
  pol(2)=p2;
  Handle(Geom2d_Curve) curv = new Geom2d_BezierCurve(pol);
  return ChFi3d_mkbound(HS,curv,t3d,ta,isfreeboundary);
}

Handle(GeomFill_Boundary) ChFi3d_mkbound
(const Handle(Adaptor3d_HSurface)& HS,
 const Handle(Geom2d_Curve)& curv,
 const Standard_Real t3d,
 const Standard_Real ta,
 const Standard_Boolean isfreeboundary)
{
  Handle(Geom2dAdaptor_HCurve) HC = new Geom2dAdaptor_HCurve();
  HC->ChangeCurve2d().Load(curv);  
  Adaptor3d_CurveOnSurface COnS(HC,HS);
 
  if (isfreeboundary) {
    Handle(Adaptor3d_HCurveOnSurface) HCOnS = new Adaptor3d_HCurveOnSurface(COnS);
    return new GeomFill_SimpleBound(HCOnS,t3d,ta); 
  }
  else {
    return new GeomFill_BoundWithSurf(COnS,t3d,ta);
  }
}

Handle(GeomFill_Boundary) ChFi3d_mkbound
(const Handle(Adaptor3d_HSurface)& Fac,
 Handle(Geom2d_Curve)& curv, 
 const gp_Pnt2d& p1,
 const gp_Pnt2d& p2,
 const Standard_Real t3d,
 const Standard_Real ta,
 const Standard_Boolean isfreeboundary)
{
  TColgp_Array1OfPnt2d pol(1,2);
  pol(1)=p1;
  pol(2)=p2;
  curv = new Geom2d_BezierCurve(pol);
  return ChFi3d_mkbound(Fac,curv,t3d,ta,isfreeboundary);
}

Handle(Geom2d_Curve) ChFi3d_BuildPCurve(const gp_Pnt2d& p1,
                              gp_Dir2d& d1,
                              const gp_Pnt2d& p2,
                              gp_Dir2d& d2,
                              const Standard_Boolean redresse)
{
  gp_Vec2d vref(p1,p2);
  gp_Dir2d dref(vref);
  Standard_Real mref = vref.Magnitude();
  if(redresse) {
    if(d1.Dot(dref) < 0.) d1.Reverse();
    if(d2.Dot(dref) > 0.) d2.Reverse();
  }
  //On fait une cubique a la mords moi le noeud
  TColgp_Array1OfPnt2d pol(1,4);
  pol(1)=p1;
  pol(4)=p2;
  Standard_Real Lambda1 = Max(Abs(d2.Dot(d1)),Abs(dref.Dot(d1)));
  Lambda1 = Max(0.5*mref*Lambda1,1.e-5);
  pol(2) = gp_Pnt2d(p1.XY()+Lambda1*d1.XY());
  Standard_Real Lambda2 = Max(Abs(d1.Dot(d2)),Abs(dref.Dot(d2)));
  Lambda2 = Max(0.5*mref*Lambda2,1.e-5);
  pol(3)=gp_Pnt2d(p2.XY()+Lambda2*d2.XY());
  return new Geom2d_BezierCurve(pol);

#if 0
  //essai avec batten
  Standard_Real h =  0.005*mref;
  Standard_Real ang1,ang2;
  ang1=dref.Angle(d1);
  if (dref.Angle(d2) > 0) ang2=PI-dref.Angle(d2);
  else ang2=-PI-dref.Angle(d2);
  Standard_Real angmax = 0.5*PI;
  if((Abs(ang1)>angmax)||(Abs(ang2)>angmax)) {
    //On fait une cubique a la mords moi le noeud
    TColgp_Array1OfPnt2d pol(1,4);
    pol(1)=p1;
    pol(4)=p2;
    Standard_Real Lambda1 = Max(Abs(d2.Dot(d1)),Abs(dref.Dot(d1)));
    Lambda1 = Max(0.5*mref*Lambda1,1.e-5);
    pol(2) = gp_Pnt2d(p1.XY()+Lambda1*d1.XY());
    Standard_Real Lambda2 = Max(Abs(d1.Dot(d2)),Abs(dref.Dot(d2)));
    Lambda2 = Max(0.5*mref*Lambda2,1.e-5);
    pol(3)=gp_Pnt2d(p2.XY()+Lambda2*d2.XY());
    return new Geom2d_BezierCurve(pol);
  }
  else{
    //On lance un Batten
#ifdef DEB
    OSD_Chronometer ch;
    ChFi3d_InitChron(ch);// init performances pour les batten  
#endif  
    FairCurve_Batten Bat(p1,p2,h);
    Bat.SetFreeSliding (Standard_True);
    Bat.SetAngle1(ang1);
    Bat.SetAngle2(ang2);
    FairCurve_AnalysisCode Iana; 
    Standard_Real tol = 1.e-3;
    Standard_Integer nbit;
    Standard_Boolean ok;
//    nbit = 5;
//    ok = Bat.Compute(Iana,nbit,tol);
    nbit = 25;
    ok = Bat.Compute(Iana,nbit,tol);
#ifdef DEB
    ChFi3d_ResultChron( ch,t_batten);  // result performances pour les batten 
#endif
    if(!ok) Standard_Failure::Raise("Echec Batten"); 
    return Bat.Curve();
  } 
#endif
}

Handle(Geom2d_Curve) ChFi3d_BuildPCurve(const Handle(Adaptor3d_HSurface)& Surf,
                              const gp_Pnt2d&                 p1,
                              const gp_Vec2d&                 v1,
                              const gp_Pnt2d&                 p2,
                              const gp_Vec2d&                 v2,
                              const Standard_Boolean          redresse)
{
  gp_Pnt2d pp1 = p1, pp2 = p2;
  gp_Vec2d vv1 = v1, vv2 = v2;
  Standard_Real ures = Surf->UResolution(1.);
  Standard_Real vres = Surf->VResolution(1.);
  Standard_Real invures = 1./ures;
  Standard_Real invvres = 1./vres;
  pp1.SetX(invures*pp1.X()); pp1.SetY(invvres*pp1.Y()); 
  pp2.SetX(invures*pp2.X()); pp2.SetY(invvres*pp2.Y()); 
  vv1.SetX(invures*vv1.X()); vv1.SetY(invvres*vv1.Y()); 
  vv2.SetX(invures*vv2.X()); vv2.SetY(invvres*vv2.Y()); 
  gp_Dir2d d1(vv1), d2(vv2);
  Handle(Geom2d_Curve) g2dc = ChFi3d_BuildPCurve(pp1,d1,pp2,d2,redresse);
  Handle(Geom2d_BezierCurve) pc = 
    Handle(Geom2d_BezierCurve)::DownCast(g2dc);
  Standard_Integer nbp = pc->NbPoles();
  for(Standard_Integer ip = 1; ip <= nbp; ip++) {
    gp_Pnt2d pol = pc->Pole(ip);
    pol.SetX(ures*pol.X()); pol.SetY(vres*pol.Y());
    pc->SetPole(ip,pol);
  }
  return pc;
}

Handle(Geom2d_Curve) ChFi3d_BuildPCurve(const Handle(Adaptor3d_HSurface)& Surf,
                              const gp_Pnt2d&                 p1,
                              const gp_Vec&                   v1,
                              const gp_Pnt2d&                 p2,
                              const gp_Vec&                   v2,
                              const Standard_Boolean          redresse)
{
  gp_Vec D1u,D1v;
  gp_Pnt PP1,PP2;
  Standard_Real DU,DV;
  Surf->D1(p1.X(),p1.Y(),PP1,D1u,D1v);
  ChFi3d_Coefficient(v1,D1u,D1v,DU,DV);
  gp_Vec2d vv1(DU,DV);
  Surf->D1(p2.X(),p2.Y(),PP2,D1u,D1v);
  ChFi3d_Coefficient(v2,D1u,D1v,DU,DV);
  gp_Vec2d vv2(DU,DV);
  gp_Vec Vref(PP1,PP2);
  if(redresse) {
    if(Vref.Dot(v1) < 0.) vv1.Reverse();
    if(Vref.Dot(v2) > 0.) vv2.Reverse();
  }
  return ChFi3d_BuildPCurve(Surf,p1,vv1,p2,vv2,0);
}

//=======================================================================
//function : ComputeArete
//purpose  : 
// pour les remplissages en s.d. un conge avec ces pcurves est construit ainsi
// firstpoint sur S1 -------------edge:courbe3d/pcurves--->lastpoint sur S1
//  |                                                              |
//  |                                                              |
//  |                                                              |
// edge:courbe 3d/pcurves           conge                        edge
//  |         attention il faut tester l orientation du conge avant|
//  |         de determiner les transitions pcurves/conge          |
//  |                                                              |
//  \/                                                             \/
// firstpoint sur S2 -------------edge:courbe3d/pcurves--->lastpoint sur S2
//
//=======================================================================

void  ChFi3d_ComputeArete(const ChFiDS_CommonPoint&   P1,
                    const gp_Pnt2d&             UV1,
                    const ChFiDS_CommonPoint&   P2,
                    const gp_Pnt2d&             UV2,
                    const Handle(Geom_Surface)& Surf,
                    Handle(Geom_Curve)&         C3d,
                    Handle(Geom2d_Curve)&       Pcurv,
                    Standard_Real&              Pardeb,
                    Standard_Real&              Parfin,
                    const Standard_Real         tol3d,
                    const Standard_Real         tol2d,
                    Standard_Real&              tolreached,
                    const Standard_Integer      IFlag)
     // IFlag=0 pcurve et courbe 3d 
     // IFlag>0 pcurve (parametrage impose si IFlag=2)
{
  static Handle(GeomAdaptor_HSurface) hs = new GeomAdaptor_HSurface();
  static Handle(GeomAdaptor_HCurve) hc = new GeomAdaptor_HCurve();

  tolreached = tol3d;

  if (Abs(UV1.X()-UV2.X()) <= tol2d) {
    if (IFlag == 0) {
      Pardeb = UV1.Y();
      Parfin = UV2.Y();
      C3d = Surf->UIso(UV1.X());
      if(Pardeb > Parfin) {
      Pardeb = C3d->ReversedParameter(Pardeb);
      Parfin = C3d->ReversedParameter(Parfin);
      C3d->Reverse();
      }
      Handle(Geom_TrimmedCurve) tc = Handle(Geom_TrimmedCurve)::DownCast(C3d);
      if(!tc.IsNull()) {
      C3d = tc->BasisCurve();
      if (C3d->IsPeriodic()) {
        ElCLib::AdjustPeriodic(C3d->FirstParameter(),C3d->LastParameter(),
                         tol2d,Pardeb,Parfin);
      }
      }
    }
    if(IFlag != 1) {
      hs->ChangeSurface().Load(Surf);
      hc->ChangeCurve().Load(C3d,Pardeb,Parfin);
      ChFi3d_ComputePCurv(hc,UV1,UV2,Pcurv,hs,Pardeb,Parfin,tol3d,tolreached,Standard_False);
    }
    else{
      Pcurv = new Geom2d_Line(UV1,gp_Vec2d(UV1,UV2));
    }
  }
  else if (Abs(UV1.Y()-UV2.Y())<=tol2d) {
    //iso v
    if (IFlag == 0) {
      Pardeb = UV1.X();
      Parfin = UV2.X();
      C3d = Surf->VIso(UV1.Y());
      if(Pardeb > Parfin) {
      Pardeb = C3d->ReversedParameter(Pardeb);
      Parfin = C3d->ReversedParameter(Parfin);
      C3d->Reverse();
      }
      Handle(Geom_TrimmedCurve) tc = Handle(Geom_TrimmedCurve)::DownCast(C3d);
      if(!tc.IsNull()) {
      C3d = tc->BasisCurve();
      if (C3d->IsPeriodic()) {
        ElCLib::AdjustPeriodic(C3d->FirstParameter(),C3d->LastParameter(),
                         tol2d,Pardeb,Parfin);
      }
      }
    }
    if(IFlag != 1) {
      hs->ChangeSurface().Load(Surf);
      hc->ChangeCurve().Load(C3d,Pardeb,Parfin);
      ChFi3d_ComputePCurv(hc,UV1,UV2,Pcurv,hs,Pardeb,Parfin,tol3d,tolreached,Standard_False);
    }
    else{
      Pcurv = new Geom2d_Line(UV1,gp_Vec2d(UV1,UV2));
    }
  }
  else if (IFlag == 0) {

    if (P1.IsVertex() || P2.IsVertex() || !P1.IsOnArc() || !P2.IsOnArc()) {
    // On construit une droite pour ne pas se tromper 
    // d'arc et donc de tangente.
      TColgp_Array1OfPnt2d qoles(1,2);
      qoles(1)=UV1;
      qoles(2)=UV2;
      Pcurv = new Geom2d_BezierCurve(qoles);
    }
    else {
      BRepAdaptor_Curve C1(P1.Arc());
      gp_Pnt Pp;
      gp_Vec Vv1;
      C1.D1(P1.ParameterOnArc(),Pp,Vv1);
      C1.Initialize(P2.Arc());
      gp_Vec Vv2;
      C1.D1(P2.ParameterOnArc(),Pp,Vv2);
      hs->ChangeSurface().Load(Surf);
      Pcurv = ChFi3d_BuildPCurve(hs,UV1,Vv1,UV2,Vv2,Standard_True); 
      // Il y a des cas ou la PCurve ainsi construite sort de la 
      // surface, en particulier lorsque celle-ci provient d un 
      // prolongement. On fait donc un controle a posteriori et
      // si ca sort on la remplace par la droite UV1 UV2, tant
      // pis pour la tangence avec les arcs voisins!!!
      Bnd_Box2d bs;
      Standard_Real umin,umax,vmin,vmax;
      Surf->Bounds(umin,umax,vmin,vmax);
      bs.Update(umin,vmin,umax,vmax);
      Standard_Boolean IN = Standard_True;
      for(Standard_Integer ii = 1; ii <= 4 && IN; ii++) {
      if(bs.IsOut((*((Handle_Geom2d_BezierCurve*) &Pcurv))->Pole(ii))) {
        IN = Standard_False;
        TColgp_Array1OfPnt2d qoles(1,2);
        qoles(1)=UV1;
        qoles(2)=UV2;
        Pcurv = new Geom2d_BezierCurve(qoles);
      }
      }
    }
    Geom2dAdaptor_Curve AC(Pcurv);
    Handle(Geom2dAdaptor_HCurve) AHC = 
      new Geom2dAdaptor_HCurve(AC);
    GeomAdaptor_Surface AS(Surf);
    Handle(GeomAdaptor_HSurface) AHS = 
      new GeomAdaptor_HSurface(AS);
    Adaptor3d_CurveOnSurface Cs(AHC,AHS);
    Pardeb = Cs.FirstParameter();
    Parfin = Cs.LastParameter();
    Standard_Real avtol;
    GeomLib::BuildCurve3d(tol3d,Cs,Pardeb,Parfin,C3d,tolreached,avtol);
  }
  else {
    hs->ChangeSurface().Load(Surf);
    hc->ChangeCurve().Load(C3d,Pardeb,Parfin);
    ChFi3d_ProjectPCurv(hc,hs,Pcurv,tol3d,tolreached);
    gp_Pnt2d p2d = Pcurv->Value(Pardeb);
    if(!UV1.IsEqual(p2d,Precision::PConfusion())) {
      gp_Vec2d v2d(p2d,UV1);
      Pcurv->Translate(v2d);
    }
  }
}

//=======================================================================
//function : FilCurveInDS
//purpose  : 
//=======================================================================


Handle(TopOpeBRepDS_SurfaceCurveInterference)  ChFi3d_FilCurveInDS
(const Standard_Integer Icurv,
 const Standard_Integer Isurf,
 const Handle(Geom2d_Curve)& Pcurv,
 const TopAbs_Orientation Et)
{
  Handle(TopOpeBRepDS_SurfaceCurveInterference) SC1;
  SC1 = new TopOpeBRepDS_SurfaceCurveInterference(TopOpeBRepDS_Transition(Et),
                                      TopOpeBRepDS_SURFACE,
                                      Isurf,TopOpeBRepDS_CURVE,Icurv,
                                      Pcurv);
  return SC1;
}

//=======================================================================
//function : TrsfTrans
//purpose  : 
//           
//=======================================================================
TopAbs_Orientation ChFi3d_TrsfTrans(const IntSurf_TypeTrans T1) 
{
  if (T1 == IntSurf_In) return TopAbs_FORWARD;
  else if(T1 == IntSurf_Out) return TopAbs_REVERSED;
  else return TopAbs_INTERNAL;
}

//=======================================================================
//function : FilCommonPoint
//purpose  : Chargement du common point
//           gestion du fait que l'on est deja sur un vertex existant
//=======================================================================

Standard_EXPORT void ChFi3d_FilCommonPoint(const BRepBlend_Extremity& SP,
                                 const IntSurf_TypeTrans TransLine,
                                 const Standard_Boolean Start,
                                 ChFiDS_CommonPoint& CP,
                                 const Standard_Real Tol)
{
//  BRep_Tool Outil;
  Standard_Real Dist, maxtol = Max(Tol,CP.Tolerance());

  CP.SetPoint(SP.Value()); // On commence par le point, le vecteur
  if (SP.HasTangent()) {
    if (Start) {
      CP.SetVector(SP.Tangent().Reversed()); // On oriente la tangente vers la sortie
    }
    else {
      CP.SetVector(SP.Tangent());
    }
  }
  
  CP.SetParameter(SP.ParameterOnGuide()); //et le parametre de la spine

  if (SP.IsVertex()) { // On charge le Vertex si besoin est
                       // (A l'interieur d'une face)
    TopoDS_Vertex V =  
      Handle(BRepTopAdaptor_HVertex)::DownCast(SP.Vertex())->Vertex();
  
    CP.SetVertex(V);  
    Dist = (SP.Value()).Distance(BRep_Tool::Pnt(V));
    //// modified by jgv, 18.09.02 for OCC571 ////
    //maxtol += Dist;
    maxtol = Max( Dist, maxtol );
    //////////////////////////////////////////////
    CP.SetPoint(BRep_Tool::Pnt(V));
    
    //la sequence d arcs l information est connu par le vertex (ancestor)
    //dans ce cas on ne calculera pas les transitions c est a ce programme
    //de le faire
  }
  
  if (SP.NbPointOnRst() != 0) { //  On charge un arc, et/ou un vertex

    const BRepBlend_PointOnRst& PR = SP.PointOnRst(1);
    Handle(BRepAdaptor_HCurve2d) 
      Harc = Handle(BRepAdaptor_HCurve2d)::DownCast(PR.Arc());
    if(!Harc.IsNull()) {

      Standard_Real DistF, DistL, LeParamAmoi;
      Standard_Integer Index_min;
      TopoDS_Edge E = Harc->ChangeCurve2d().Edge();

      TopoDS_Vertex V[2];
      TopExp::Vertices(E, V[0], V[1]);

      DistF = (SP.Value()).Distance(BRep_Tool::Pnt(V[0]));
      DistL = (SP.Value()).Distance(BRep_Tool::Pnt(V[1]));
      if (DistF<DistL) { Index_min = 0;
                     Dist = DistF; }
      else             { Index_min = 1;
                   Dist = DistL; }

      if (Dist <= maxtol + BRep_Tool::Tolerance(V[Index_min]) ) { 
      // On tombe sur un vertex prexistant
       CP.SetVertex(V[Index_min]); //On Charge l'ancien vertex
         CP.SetPoint( BRep_Tool::Pnt(V[Index_min]) );
         maxtol = Max(BRep_Tool::Tolerance(V[Index_min]),maxtol);
       //// modified by jgv, 18.09.02 for OCC571 ////
       //maxtol += Dist;
       maxtol = Max( Dist, maxtol );
       //////////////////////////////////////////////
         LeParamAmoi = BRep_Tool::Parameter(V[Index_min], E);    
      }
      else {   // Creation d'un arc seul
      maxtol = Max(BRep_Tool::Tolerance(E),maxtol);
      maxtol = Max(SP.Tolerance(),maxtol);
      LeParamAmoi = PR.ParameterOnArc();
      }

      // Definition de l'arc
      TopAbs_Orientation Tr;
      TopAbs_Orientation Or = E.Orientation();
      if (Start) {
      Tr = TopAbs::Reverse(TopAbs::Compose(ChFi3d_TrsfTrans(TransLine),Or));
      }
      else {
      Tr = TopAbs::Compose(ChFi3d_TrsfTrans(TransLine),Or);
      }
      CP.SetArc(maxtol, E, LeParamAmoi, Tr);
    }
  }
  CP.SetTolerance(maxtol); // On finit par la tolerance.
}

//=======================================================================
//function : SolidIndex
//purpose  : 
//=======================================================================

Standard_Integer ChFi3d_SolidIndex(const Handle(ChFiDS_Spine)&  sp,
                           TopOpeBRepDS_DataStructure&  DStr,
                           ChFiDS_Map&                  MapESo,
                           ChFiDS_Map&                  MapESh)
{
  if(sp.IsNull() || sp->NbEdges() == 0) 
    Standard_Failure::Raise("SolidIndex : Spine incomplete");
  TopoDS_Shape edref= sp->Edges(1);
  TopoDS_Shape shellousolid;
  if(!MapESo(edref).IsEmpty()) shellousolid = MapESo(edref).First();
  else shellousolid = MapESh(edref).First();
  Standard_Integer solidindex = DStr.AddShape(shellousolid);
  return solidindex;
}

//=======================================================================
//function : IndexPointInDS
//purpose  : 
//=======================================================================

Standard_Integer  ChFi3d_IndexPointInDS(const ChFiDS_CommonPoint& P1,
                              TopOpeBRepDS_DataStructure& DStr) 
{
  if (P1.IsVertex()) {
    // --------------------------------->  !*!*!* 
    // Attention : On se permet d'ecrabouiller la tolerance
    // il faudrait prevoir un mecanimse propre.
    BRep_Builder B;
    B.UpdateVertex(P1.Vertex(), P1.Point(), P1.Tolerance());
    return DStr.AddShape(P1.Vertex());
  }
  return DStr.AddPoint(TopOpeBRepDS_Point(P1.Point(),P1.Tolerance()));
}

//=======================================================================
//function : FilPointInDS
//purpose  : 
//=======================================================================

Handle(TopOpeBRepDS_CurvePointInterference) ChFi3d_FilPointInDS
(const TopAbs_Orientation Et,
 const Standard_Integer Ic,
 const Standard_Integer Ip,
 const Standard_Real Par,
 const Standard_Boolean IsVertex)
{
  Handle(TopOpeBRepDS_CurvePointInterference) CP1;
  if (IsVertex)    
    CP1 = new TopOpeBRepDS_CurvePointInterference (TopOpeBRepDS_Transition(Et),
                                       TopOpeBRepDS_CURVE,Ic,
                                       TopOpeBRepDS_VERTEX,Ip,Par); 
  else
    CP1 = new TopOpeBRepDS_CurvePointInterference (TopOpeBRepDS_Transition(Et),
                                       TopOpeBRepDS_CURVE,Ic,
                                       TopOpeBRepDS_POINT,Ip,Par);
  return CP1;
}

//=======================================================================
//function : FilVertexInDS
//purpose  : 
//=======================================================================

Handle(TopOpeBRepDS_CurvePointInterference) ChFi3d_FilVertexInDS
(const TopAbs_Orientation Et,
 const Standard_Integer Ic,
 const Standard_Integer Ip,
 const Standard_Real Par)
{
  
  Handle(TopOpeBRepDS_CurvePointInterference) CP1 = new
    TopOpeBRepDS_CurvePointInterference (TopOpeBRepDS_Transition(Et),
                               TopOpeBRepDS_CURVE,Ic,
                               TopOpeBRepDS_VERTEX,Ip,Par);
  return CP1;
}

//=======================================================================
//function : Orientation
//purpose  : retourne l'orientation d'une interference (la premiere trouvee
//           dans la liste).
//=======================================================================

static Standard_Boolean ChFi3d_Orientation
(const TopOpeBRepDS_ListOfInterference& LI,
 const Standard_Integer                 igros,
 const Standard_Integer                 ipetit,
 TopAbs_Orientation&                    Or,
 const Standard_Boolean                 isvertex = Standard_False,
 const Standard_Boolean                 aprendre = Standard_False)
{
  //Dans le cas, ou on veux inserer un point/vertex, on desire savoir si c'est un point ou un vertex, car leur index peuvent etre les memes.
  TopOpeBRepDS_Kind typepetit;
  if (isvertex)
    typepetit =  TopOpeBRepDS_VERTEX;
  else
    typepetit =  TopOpeBRepDS_POINT;
  TopOpeBRepDS_ListIteratorOfListOfInterference itLI(LI);
  for (; itLI.More(); itLI.Next() ) {
    const Handle(TopOpeBRepDS_Interference)& cur = itLI.Value(); 
    TopOpeBRepDS_Kind GK;
    TopOpeBRepDS_Kind SK;
    Standard_Integer S;
    Standard_Integer G;
    cur->GKGSKS(GK,G,SK,S);
    if (aprendre) {
      if ( S == igros && G == ipetit && GK == typepetit) {
      Or = cur->Transition().Orientation(TopAbs_IN);
      return Standard_True;
      }
    }
    else  {
      if ( S == igros && G == ipetit) {
      Or = cur->Transition().Orientation(TopAbs_IN);
      return Standard_True;
      } 
    }
  }
  return Standard_False;
}

//=======================================================================
//function : Contains
//purpose  : Verifie qu une interference n existe pas deja.
//=======================================================================

static Standard_Boolean ChFi3d_Contains
(const TopOpeBRepDS_ListOfInterference& LI,
 const Standard_Integer                 igros,
 const Standard_Integer                 ipetit,
 const Standard_Boolean                 isvertex = Standard_False,
 const Standard_Boolean                 aprendre = Standard_False)                             
{
  TopAbs_Orientation bidOr;
  return ChFi3d_Orientation(LI,igros,ipetit,bidOr,isvertex,aprendre);
}

static void QueryAddVertexInEdge
(      TopOpeBRepDS_ListOfInterference& LI,
 const Standard_Integer                 IC,
 const Standard_Integer                 IV,
 const Standard_Real                    par,
 const TopAbs_Orientation               Or)
{
  TopOpeBRepDS_ListIteratorOfListOfInterference it(LI);
  for (; it.More(); it.Next() ) {
    const Handle(TopOpeBRepDS_Interference)& cur = it.Value();
    const Handle(TopOpeBRepDS_CurvePointInterference)& cpi = 
      Handle(TopOpeBRepDS_CurvePointInterference)::DownCast(cur);
    if(!cpi.IsNull()) {
      Standard_Integer newIV = cpi->Geometry();
      TopOpeBRepDS_Kind kv = cpi->GeometryType();
      TopAbs_Orientation newOr = cpi->Transition().Orientation(TopAbs_IN);
      Standard_Real newpar = cpi->Parameter();
      if(IV == newIV && kv == TopOpeBRepDS_VERTEX && 
       Or == newOr && Abs(par - newpar) < 1.e-10) {
      return;
      }
    }
  }
  Handle(TopOpeBRepDS_CurvePointInterference) interf = 
    ChFi3d_FilVertexInDS(Or,IC,IV,par);
  LI.Append(interf);
}

static void CutEdge(const TopoDS_Vertex&           V,
                const Handle(ChFiDS_SurfData)& SD,
                TopOpeBRepDS_DataStructure&    DStr,
                const Standard_Boolean         /*isfirst*/,
                const Standard_Integer         ons)
{
  if(!SD->IsOnCurve(ons)) return;
  Standard_Integer IC = SD->IndexOfC(ons);
  Standard_Integer IV = DStr.AddShape(V);
  TopOpeBRepDS_ListOfInterference& LI = DStr.ChangeShapeInterferences(IC);
  TopoDS_Edge E = TopoDS::Edge(DStr.Shape(IC));
  E.Orientation(TopAbs_FORWARD);
  TopExp_Explorer ex;

  // les traiter tous en verifiant que ce n est pas 
  // deja fait.
  for(ex.Init(E,TopAbs_VERTEX);ex.More();ex.Next()) {
    const TopoDS_Vertex& vv = TopoDS::Vertex(ex.Current());
    if(vv.IsSame(V)) {
      TopAbs_Orientation Or = TopAbs::Reverse(vv.Orientation());
      Standard_Real par = BRep_Tool::Parameter(vv,E);
      QueryAddVertexInEdge(LI,IC,IV,par,Or);
    }
  }
}

//=======================================================================
//function : findIndexPoint
//purpose  : returns in <ipon> index of point bounding a courve interfering
//           with <Fd> and coinciding with last common point on <OnS> face
//=======================================================================

static Standard_Boolean findIndexPoint(const TopOpeBRepDS_DataStructure& DStr,
                               const Handle(ChFiDS_SurfData)&    Fd,
                               const Standard_Integer            OnS,
                               Standard_Integer&                 ipoin)
{
  ipoin = 0;
  gp_Pnt P = Fd->Vertex(Standard_False,OnS).Point();

  TopOpeBRepDS_ListIteratorOfListOfInterference SCIIt, CPIIt;
  
  SCIIt.Initialize (DStr.SurfaceInterferences(Fd->Surf()));
  for (; SCIIt.More(); SCIIt.Next()) {
    Handle(TopOpeBRepDS_SurfaceCurveInterference) SCI =
      Handle(TopOpeBRepDS_SurfaceCurveInterference)::DownCast(SCIIt.Value());
    if (SCI.IsNull()) continue;
    CPIIt.Initialize (DStr.CurveInterferences(SCI->Geometry()));
    for (; CPIIt.More(); CPIIt.Next()) {
      Handle(TopOpeBRepDS_CurvePointInterference) CPI =
      Handle(TopOpeBRepDS_CurvePointInterference)::DownCast(CPIIt.Value());
      if (CPI.IsNull()) continue;
      Standard_Integer iPoint = CPI->Geometry();
      TopOpeBRepDS_Point tp = DStr.Point(iPoint);
      if (P.IsEqual(tp.Point(), tp.Tolerance())) {
      ipoin = iPoint;
      return Standard_True;
      }
    }
  }
  return Standard_False;
}
// *******************************************************************
//=======================================================================
//function : FilDS
//purpose  : 
//=======================================================================
// *******************************************************************

void  ChFi3d_FilDS(const Standard_Integer       SolidIndex,
               const Handle(ChFiDS_Stripe)& CorDat,
               TopOpeBRepDS_DataStructure&  DStr,
               ChFiDS_Regularities&         reglist,
               const Standard_Real          tol3d,
               const Standard_Real          tol2d) 
{
//  BRep_Tool Outil;
  TopExp_Explorer ex;
  Handle(ChFiDS_Spine) spine = CorDat->Spine();
  Standard_Boolean Closed = Standard_False;
  Standard_Boolean Degene = 0, isVertex1 = 0, isVertex2 = 0, Singulier_en_Bout = 0;
  if(!spine.IsNull()) {
    Closed = spine->IsPeriodic();
  }
  const ChFiDS_SequenceOfSurfData& SeqFil = 
    CorDat->SetOfSurfData()->Sequence();
  Standard_Integer Ipoin1 = CorDat->IndexFirstPointOnS1();
  Standard_Integer Ipoin2 = CorDat->IndexFirstPointOnS2();
  Standard_Integer NumEdge = 1;
  TopoDS_Vertex BoutdeVtx; 
  Standard_Integer Icurv = 0;
  Standard_Integer Iarc1 = 0,Iarc2 = 0;
  TopAbs_Orientation trafil1 = TopAbs_FORWARD, trafil2 = TopAbs_FORWARD;
  Standard_Integer IcFil1,IcFil2,Isurf,Ishape1,Ishape2;
  Standard_Real Pardeb,Parfin;
  TopAbs_Orientation ET1;
  Handle(TopOpeBRepDS_CurvePointInterference) Interfp1,Interfp2;
  Handle(TopOpeBRepDS_SurfaceCurveInterference) Interfc1,Interfc2;
  Handle(TopOpeBRepDS_SurfaceCurveInterference) Interfc3,Interfc4;
  Handle(TopOpeBRepDS_CurvePointInterference) Interfp3,Interfp4;
  Handle(TopOpeBRepDS_CurvePointInterference) Interfp5,Interfp6;
  TopoDS_Face F;
  Handle(Geom2d_Curve) PCurv;
  TopOpeBRepDS_Curve Crv;

  TopOpeBRepDS_ListOfInterference& SolidInterfs = 
    DStr.ChangeShapeInterferences(SolidIndex);

  ChFiDS_Regul regcout; // pour les CD closed and tangent
  ChFiDS_Regul regfilfil; // pour les joints Surf/Surf

  ChFiDS_CommonPoint V3;
  ChFiDS_CommonPoint V4;

  // Nullify degenerated ChFi/Faces interferences, eap occ293
  Standard_Integer j;
  if (SeqFil.Length() > 1) {
    for (j=1; j<=SeqFil.Length(); j++) {
      Handle(ChFiDS_SurfData) Fd = SeqFil(j);
      Standard_Integer onS;
      for (onS=1; onS<=2; onS++) {
      const ChFiDS_FaceInterference& Fi = Fd->Interference(onS);
      IcFil1 = Fi.LineIndex();
      if (!IcFil1) continue;
      Standard_Real FiLen = Abs(Fi.FirstParameter()-Fi.LastParameter());
      if (FiLen > Precision::PConfusion()) continue;
      TopOpeBRepDS_Curve& cc = DStr.ChangeCurve(IcFil1);
      cc.ChangeCurve().Nullify();
#ifdef DEB
      cout << "Nullify degenerated FI of SurfData " << j <<" on S" << onS << endl;
#endif
      // care of CommonPoint, eap occ354
      if (j!=1 && j!=SeqFil.Length()) continue;
      Standard_Boolean isfirst = (j==1);
      Standard_Integer i = isfirst ? j+1 : j-1;
      ChFiDS_CommonPoint& CP1 = SeqFil(i)->ChangeVertex(isfirst,onS);
      if (Fd->Vertex(isfirst,onS).IsOnArc() && CP1.IsOnArc()) {
        ChFiDS_CommonPoint& CP2 = Fd->ChangeVertex(!isfirst,onS);
        CP1.Reset();
        CP1.SetPoint(CP2.Point());
        CP2.Reset();
        CP2.SetPoint(CP1.Point());
        }
      }
    }
  }

  for (j=1; j<=SeqFil.Length(); j++) {

    const Handle(ChFiDS_SurfData)& Fd = SeqFil(j);
    Isurf= Fd->Surf();
    Ishape1 = Fd->IndexOfS1();
    Ishape2 = Fd->IndexOfS2();

    // eap, Apr 29 2002, occ 293
    // now IsInDS() returns nb of surfaces at end being in DS;
    // vars showing which end is in DS
    Standard_Boolean isInDS1 = Standard_False, isInDS2 = Standard_False;
    if (j <= CorDat->IsInDS(Standard_True)) {
      isInDS1 = Standard_True;
      isInDS2 = (j+1 <= CorDat->IsInDS(Standard_True));
    }
    if (SeqFil.Length()-j < CorDat->IsInDS(Standard_False)) {
      isInDS2 = Standard_True;
      isInDS1 = isInDS1 || SeqFil.Length()-j+1 < CorDat->IsInDS(Standard_False);
    }
    
    // creation de la SolidSurfaceInterference
    
    Handle(TopOpeBRepDS_SolidSurfaceInterference) 
      SSI = new TopOpeBRepDS_SolidSurfaceInterference
      (TopOpeBRepDS_Transition(Fd->Orientation()),
       TopOpeBRepDS_SOLID,
       SolidIndex,
       TopOpeBRepDS_SURFACE,
       Isurf);
    
    SolidInterfs.Append(SSI);
    
    const ChFiDS_FaceInterference& Fi1 = Fd->InterferenceOnS1();
    const ChFiDS_FaceInterference& Fi2 = Fd->InterferenceOnS2();     
    const ChFiDS_CommonPoint& V1 = Fd->VertexFirstOnS1();
    const ChFiDS_CommonPoint& V2 = Fd->VertexFirstOnS2();

    // Un petit traitement pour gerer les interference doubles
    if (j>1) {
      if (V1.IsOnArc() && V3.IsOnArc() && V1.Arc().IsSame(V3.Arc())) {
      //on initialise Iarc1
      //Iarc1 = DStr.AddShape(V1.Arc());
      if (ChFi3d_Contains(DStr.ShapeInterferences(Iarc1),Iarc1,Ipoin1) && 
          (V1.TransitionOnArc() != V3.TransitionOnArc()) ) {
        Interfp1= ChFi3d_FilPointInDS(V1.TransitionOnArc(),Iarc1,Ipoin1,
                              V1.ParameterOnArc());
        DStr.ChangeShapeInterferences(V1.Arc()).Append(Interfp1);
      }
      }

      if (V2.IsOnArc() && V4.IsOnArc() && V2.Arc().IsSame(V4.Arc())) {
      //on initialise Iarc2
      //Iarc2 = DStr.AddShape(V2.Arc());
      if ( ChFi3d_Contains(DStr.ShapeInterferences(Iarc2),Iarc2,Ipoin2)  && 
          (V2.TransitionOnArc() != V4.TransitionOnArc()) ) {
        Interfp2= ChFi3d_FilPointInDS(V2.TransitionOnArc(),Iarc2,Ipoin2,
                              V2.ParameterOnArc());
        DStr.ChangeShapeInterferences(V2.Arc()).Append(Interfp2);
      }
      }
    }

    V3 = Fd->VertexLastOnS1();
    V4 = Fd->VertexLastOnS2();
    
    if(Ishape1 != 0) {
      if(Ishape1 > 0) {
      trafil1 = DStr.Shape(Ishape1).Orientation();
      }
      else{
        ChFi3d_Orientation(SolidInterfs,SolidIndex,-Ishape1,trafil1);
      }
      trafil1 = TopAbs::Compose(trafil1,Fd->Orientation());
      trafil1 = TopAbs::Compose(TopAbs::Reverse(Fi1.Transition()),trafil1);
      trafil2 = TopAbs::Reverse(trafil1);
    }
    else{
      if(Ishape2 > 0) {
      trafil2 = DStr.Shape(Ishape2).Orientation();
      }
      else{
        ChFi3d_Orientation(SolidInterfs,SolidIndex,-Ishape2,trafil2);
      }
      trafil2 = TopAbs::Compose(trafil2,Fd->Orientation());
      trafil2 = TopAbs::Compose(TopAbs::Reverse(Fi2.Transition()),trafil2);
      trafil1 = TopAbs::Reverse(trafil2);
    }
    
    ET1 = TopAbs::Reverse(trafil1);
    
    // Un petit paragraphe pour traiter les contacts aretes qui touchent 
    // un vertex de l'obstacle.
    if(V1.IsVertex() && Fd->IsOnCurve1()) {
      const TopoDS_Vertex& vv1 = V1.Vertex();
      CutEdge(vv1,Fd,DStr,1,1);
    }
    if(V2.IsVertex() && Fd->IsOnCurve2()) {
      const TopoDS_Vertex& vv2 = V2.Vertex();
      CutEdge(vv2,Fd,DStr,1,2);
    }
    if(V3.IsVertex() && Fd->IsOnCurve1()) {
      const TopoDS_Vertex& vv3 = V3.Vertex();
      CutEdge(vv3,Fd,DStr,0,1);
    }
    if(V4.IsVertex() && Fd->IsOnCurve2()) {
      const TopoDS_Vertex& vv4 = V4.Vertex();
      CutEdge(vv4,Fd,DStr,0,2);
    }
    
    if (j == 1) {
      isVertex1 = V1.IsVertex();
      isVertex2 = V2.IsVertex();
      Singulier_en_Bout =  (V1.Point().IsEqual(V2.Point(), 0));
      
      if (Singulier_en_Bout) {
         // Queue de Billard
         if ((!V1.IsVertex()) || (!V2.IsVertex())) {
#if DEB
         cout << " Singularite en bout hors vertex " << endl;
         // QUEZAKO ????
#endif
       }
       else {
         isVertex1 = isVertex2 = Standard_True; //precaution...
         // On elimine l'arete de la spine debouchant sur ce vertex.
         TopoDS_Edge Arcspine = spine->Edges(1);
         BoutdeVtx = V1.Vertex();
         Standard_Integer IArcspine = DStr.AddShape(Arcspine);
         Standard_Integer IVtx = CorDat->IndexFirstPointOnS1();
#ifndef DEB
         TopAbs_Orientation OVtx = TopAbs_FORWARD;
#else
         TopAbs_Orientation OVtx;
#endif
         for(ex.Init(Arcspine.Oriented(TopAbs_FORWARD),TopAbs_VERTEX); 
             ex.More(); ex.Next()) {
           if(BoutdeVtx.IsSame(ex.Current())) {
             OVtx = ex.Current().Orientation();
             break;
           }
         }
         OVtx = TopAbs::Reverse(OVtx);
         Standard_Real parVtx = BRep_Tool::Parameter(BoutdeVtx,Arcspine);
         Handle(TopOpeBRepDS_CurvePointInterference) 
           interfv = ChFi3d_FilVertexInDS(OVtx,IArcspine,IVtx,parVtx);
         DStr.ChangeShapeInterferences(IArcspine).Append(interfv);
       }
       }
      else {
      if (V1.IsOnArc()) {
        Iarc1 = DStr.AddShape(V1.Arc()); 
        if ( !ChFi3d_Contains(DStr.ShapeInterferences(Iarc1),Iarc1,Ipoin1) ) {
          Interfp1= ChFi3d_FilPointInDS(V1.TransitionOnArc(),Iarc1,Ipoin1,
                                V1.ParameterOnArc(), isVertex1);
          DStr.ChangeShapeInterferences(V1.Arc()).Append(Interfp1);
        }
      }

      if (V2.IsOnArc()) {
        Iarc2 = DStr.AddShape(V2.Arc());
        if ( !ChFi3d_Contains(DStr.ShapeInterferences(Iarc2),Iarc2,Ipoin2) ) {
          Interfp2= ChFi3d_FilPointInDS(V2.TransitionOnArc(),Iarc2,Ipoin2,
                                V2.ParameterOnArc(),isVertex2);
          DStr.ChangeShapeInterferences(V2.Arc()).Append(Interfp2);
        }
      }
      }

      if (!isInDS1) {
      ET1 = TopAbs::Compose(ET1,CorDat->FirstPCurveOrientation());
      Icurv = CorDat->FirstCurve();
      if(Closed && !Singulier_en_Bout) {
        regcout.SetCurve(Icurv);
        regcout.SetS1(Isurf,Standard_False);
      }
      PCurv = CorDat->FirstPCurve();
      CorDat->FirstParameters(Pardeb,Parfin);
      
      TopOpeBRepDS_ListOfInterference& Li = DStr.ChangeCurveInterferences(Icurv);
      if (Li.IsEmpty()) {
        if(CorDat->FirstPCurveOrientation()==TopAbs_REVERSED) {
          Interfp1=ChFi3d_FilPointInDS
            (TopAbs_REVERSED,Icurv,Ipoin1,Parfin,isVertex1);
          Interfp2=ChFi3d_FilPointInDS
            (TopAbs_FORWARD,Icurv,Ipoin2,Pardeb,isVertex2);
        }
        else{
          Interfp1=ChFi3d_FilPointInDS
            (TopAbs_FORWARD,Icurv,Ipoin1,Pardeb,isVertex1);
          Interfp2=ChFi3d_FilPointInDS
            (TopAbs_REVERSED,Icurv,Ipoin2,Parfin,isVertex2);
        }
        Li.Append(Interfp1);
        Li.Append(Interfp2);
      }
      Interfc1= ChFi3d_FilCurveInDS (Icurv,Isurf,PCurv,ET1);
      DStr.ChangeSurfaceInterferences(Isurf).Append(Interfc1);
        if (Ipoin1 == Ipoin2) {
        TopOpeBRepDS_Curve& TCurv = DStr.ChangeCurve(Icurv);
        TCurv.ChangeCurve().Nullify();
        Handle(TopOpeBRepDS_Interference) bidinterf;
        TCurv.SetSCI(Interfc1,bidinterf);     
      }
      }
    } // Fin du Traitement Initial (j==1)
    else {
      // ---- Interference entre Conges ------
      
      if (!isInDS1) {// eap, Apr 29 2002, occ 293 
      
        if (Degene && isVertex1) {
          // On elimine l'arete de la spine debouchant sur ce vertex.
          NumEdge++; // On a dejas trouve l'arete precedente du vertex
          TopoDS_Edge Arcspine = spine->Edges(NumEdge);
          Standard_Integer IArcspine = DStr.AddShape(Arcspine);
          Standard_Integer IVtx = DStr.AddShape(BoutdeVtx);
#ifndef DEB
          TopAbs_Orientation OVtx = TopAbs_FORWARD;
#else
          TopAbs_Orientation OVtx;
#endif
          for(ex.Init(Arcspine.Oriented(TopAbs_FORWARD),TopAbs_VERTEX); 
              ex.More(); ex.Next()) {
            if(BoutdeVtx.IsSame(ex.Current())) {
              OVtx = ex.Current().Orientation();
              break;
            }
          }
          OVtx = TopAbs::Reverse(OVtx);
          Standard_Real parVtx = BRep_Tool::Parameter(BoutdeVtx,Arcspine);
          Handle(TopOpeBRepDS_CurvePointInterference) 
            interfv = ChFi3d_FilVertexInDS(OVtx,IArcspine,IVtx,parVtx);
          DStr.ChangeShapeInterferences(IArcspine).Append(interfv);
        } // Fin de l'elimination

        gp_Pnt2d UV1 = Fd->InterferenceOnS1().PCurveOnSurf()->
          Value(Fd->InterferenceOnS1().FirstParameter());
        gp_Pnt2d UV2 = Fd->InterferenceOnS2().PCurveOnSurf()-> 
          Value(Fd->InterferenceOnS2().FirstParameter());
        TopOpeBRepDS_Curve& TCurv = DStr.ChangeCurve(Icurv);
        if (Degene) {
          // on associe la pcurve via la SCI a la TopOpeBRepDSCurve.
          ChFi3d_ComputePCurv(UV1,UV2,PCurv,Pardeb,Parfin);       
          Interfc1= ChFi3d_FilCurveInDS (Icurv,Isurf,PCurv,ET1);
          DStr.ChangeSurfaceInterferences(Isurf).Append(Interfc1);
          TCurv.ChangeCurve().Nullify();
          Handle(TopOpeBRepDS_Interference) bidinterf;
          TCurv.SetSCI(Interfc1,bidinterf);           
        }
        else {
          regfilfil.SetS2(Isurf,Standard_False);
          reglist.Append(regfilfil);
          Standard_Real tolreached;
          ChFi3d_ComputePCurv(TCurv.ChangeCurve(),UV1,UV2,PCurv,
                              DStr.Surface(Fd->Surf()).Surface(),
                              Pardeb,Parfin,tol3d,tolreached);
          TCurv.Tolerance(Max(TCurv.Tolerance(),tolreached));
          Interfc1= ChFi3d_FilCurveInDS (Icurv,Isurf,PCurv,ET1);
          DStr.ChangeSurfaceInterferences(Isurf).Append(Interfc1);
        }
      }
    } // Fin Interference entre conges 
   
    // ---- Interference Conges / Faces
    IcFil1 = Fi1.LineIndex();
    
    if (IcFil1!=0 ) {
      Interfc3= ChFi3d_FilCurveInDS (IcFil1,Isurf,
                             Fi1.PCurveOnSurf(),trafil1);
      DStr.ChangeSurfaceInterferences(Isurf).Append(Interfc3);
      Ishape1 = Fd->IndexOfS1();
      // Cas d arete degeneree : on associe la pcurve via la SCI 
      // a la TopOpeBRepDSCurve.
      TopOpeBRepDS_Curve& cc = DStr.ChangeCurve(IcFil1);
      if(cc.Curve().IsNull()) {
      Handle(TopOpeBRepDS_Interference) bidinterf;
      cc.SetSCI(Interfc3,bidinterf);
      }
      else{
      ChFiDS_Regul regon1;
      regon1.SetCurve(IcFil1);
      regon1.SetS1(Isurf,Standard_False);
      if ( Ishape1 < 0 ) {
        Ishape1 = -Ishape1;
        regon1.SetS2(Ishape1,Standard_False);
        Interfc1=ChFi3d_FilCurveInDS(IcFil1,Ishape1,Fi1.PCurveOnFace(),
                               Fi1.Transition()); 
        DStr.ChangeSurfaceInterferences(Ishape1).Append(Interfc1);      
      }
      else if ( Ishape1 > 0 ) {
        regon1.SetS2(Ishape1,Standard_True);
        Interfc1=ChFi3d_FilCurveInDS(IcFil1,Ishape1,Fi1.PCurveOnFace(),
                               Fi1.Transition()); 
        DStr.ChangeShapeInterferences(Ishape1).Append(Interfc1);      
      }
      reglist.Append(regon1);
      }
      // Indice et type du point en Fin
      Standard_Integer ipoin;
      Standard_Boolean isVertex = Fd->VertexLastOnS1().IsVertex();
      if (j == SeqFil.Length()) ipoin = CorDat->IndexLastPointOnS1();
      else if ( j == (SeqFil.Length()-1)  &&  /*Closed &&*/
             (DStr.Curve(SeqFil.Last()->InterferenceOnS1().
                     LineIndex()).Curve().IsNull())) {
      if (Closed) {
        ipoin = CorDat->IndexFirstPointOnS1();
        isVertex = SeqFil(1)->VertexFirstOnS1().IsVertex();
      } else {
        ipoin = CorDat->IndexLastPointOnS1();
        isVertex = SeqFil.Last()->VertexLastOnS1().IsVertex();
      }
      }
      else if(DStr.Curve(IcFil1).Curve().IsNull()) {// Rotation !!
      ipoin = Ipoin1;
      isVertex = isVertex1;
      }
      else if ( ((j==1) || (j== SeqFil.Length()-1)) && 
              ( (Fd->VertexLastOnS1().Point().IsEqual(
              SeqFil(1)->VertexFirstOnS1().Point(), 1.e-7)) ||
              (Fd->VertexLastOnS1().Point().IsEqual(
              SeqFil(SeqFil.Length())->VertexLastOnS1().Point(), 1.e-7))) )
      // Cas des SurfData coupe de facon "Triangulaire"   
      ipoin=CorDat->IndexLastPointOnS1();

      // eap, Apr 29 2002, occ 293
      else if (isInDS2 && findIndexPoint(DStr, Fd, 1, ipoin)) {
#ifdef DEB
      cout << "ChFi3d_FilDS: point " << ipoin <<" is found in DS" << endl;
#endif
      }
      else ipoin = ChFi3d_IndexPointInDS(Fd->VertexLastOnS1(),DStr);
      
      TopOpeBRepDS_ListOfInterference& Li = DStr.ChangeCurveInterferences(IcFil1);
#ifdef DEB
      if (!Li.IsEmpty()) 
      cout<<"FilDS : Ligne de tangence sur plusieurs SurfData"<< endl;
#endif
      if (!ChFi3d_Contains(Li,IcFil1,Ipoin1)) { 
      
      Interfp1 = ChFi3d_FilPointInDS(TopAbs_FORWARD,IcFil1,Ipoin1,
                               Fi1.FirstParameter(),isVertex1);
      DStr.ChangeCurveInterferences(IcFil1).Append(Interfp1);
      }
      if (ipoin == Ipoin1 || !ChFi3d_Contains(Li,IcFil1,ipoin)) { 
      Interfp3 = ChFi3d_FilPointInDS(TopAbs_REVERSED,IcFil1,ipoin,
                               Fi1.LastParameter(), isVertex);
      DStr.ChangeCurveInterferences(IcFil1).Append(Interfp3);
      }
      Ipoin1 = ipoin;
      isVertex1 = isVertex;
    }
    
    IcFil2 = Fi2.LineIndex();
    if (IcFil2!=0) {
      Interfc4=ChFi3d_FilCurveInDS(IcFil2,Isurf,
                           Fi2.PCurveOnSurf(),trafil2);
      DStr.ChangeSurfaceInterferences(Isurf).Append(Interfc4);
      Ishape2 = Fd->IndexOfS2();
      // Cas d arete degeneree : on associe la pcurve via la SCI 
      // a la TopOpeBRepDSCurve.
      TopOpeBRepDS_Curve& cc = DStr.ChangeCurve(IcFil2);
      if(cc.Curve().IsNull()) {
      Handle(TopOpeBRepDS_Interference) bidinterf;
      cc.SetSCI(Interfc4,bidinterf);
      }
      else{
      ChFiDS_Regul regon2;
      regon2.SetCurve(IcFil2);
      regon2.SetS1(Isurf,Standard_False);
      if ( Ishape2 < 0 ) {
        Ishape2 = -Ishape2;
        regon2.SetS2(Ishape2,Standard_False);
        Interfc2=ChFi3d_FilCurveInDS(IcFil2,Ishape2,Fi2.PCurveOnFace(),
                               Fi2.Transition());
        DStr.ChangeSurfaceInterferences(Ishape2).Append(Interfc2);      
      }
      else if ( Ishape2 > 0 ) {
        regon2.SetS2(Ishape2,Standard_True);
        Interfc2=ChFi3d_FilCurveInDS(IcFil2,Ishape2,Fi2.PCurveOnFace(),
                               Fi2.Transition());
        DStr.ChangeShapeInterferences(Ishape2).Append(Interfc2);      
      }
      reglist.Append(regon2);
      }
      // Indice et type du point en Fin
      Standard_Integer ipoin;
      Standard_Boolean isVertex = Fd->VertexLastOnS2().IsVertex();
      if (j == SeqFil.Length() ) ipoin = CorDat->IndexLastPointOnS2();
      else if ( j == (SeqFil.Length()-1)  && /*Closed &&*/
             (DStr.Curve(SeqFil.Last()->InterferenceOnS2().
                     LineIndex()).Curve().IsNull())) {
      if (Closed) {
        ipoin = CorDat->IndexFirstPointOnS2();
        isVertex = SeqFil(1)->VertexFirstOnS2().IsVertex();
      } else {
        ipoin = CorDat->IndexLastPointOnS2();
        isVertex = SeqFil.Last()->VertexLastOnS2().IsVertex();
      }
      }
      else if(DStr.Curve(IcFil2).Curve().IsNull()) { // Rotation !!
      ipoin = Ipoin2;
      isVertex = isVertex2;
      }
      else if(Fd->VertexLastOnS2().Point().IsEqual(
            Fd->VertexLastOnS1().Point(), 0) ) {  //Pincement !!
      ipoin = Ipoin1;
      isVertex = isVertex1;
      }
      else if ( ((j==1) || (j==SeqFil.Length()-1)) && 
              ( (Fd->VertexLastOnS2().Point().IsEqual(
              SeqFil(1)->VertexFirstOnS2().Point(), 1.e-7)) ||
              (Fd->VertexLastOnS2().Point().IsEqual(
              SeqFil(SeqFil.Length())->VertexLastOnS2().Point(), 1.e-7))) )
      // Cas des SurfData coupe de facon "Triangulaire"    
      ipoin=CorDat->IndexLastPointOnS2();

      // eap, Apr 29 2002, occ 293
      else if (isInDS2 && findIndexPoint(DStr, Fd, 2, ipoin)) {
#ifdef DEB
      cout << "ChFi3d_FilDS: point " << ipoin <<" is found in DS" << endl;
#endif
      }
      else ipoin = ChFi3d_IndexPointInDS(Fd->VertexLastOnS2(),DStr);
      
      TopOpeBRepDS_ListOfInterference& Li = DStr.ChangeCurveInterferences(IcFil2);
#ifdef DEB
      if (!Li.IsEmpty()) 
      cout<<"FilDS : Ligne de tangence sur plusieurs SurfData"<< endl;
#endif
      if (!ChFi3d_Contains(Li,IcFil2,Ipoin2)) { 
      Interfp2 = ChFi3d_FilPointInDS(TopAbs_FORWARD,IcFil2,Ipoin2,
                               Fi2.FirstParameter(), isVertex2);
      DStr.ChangeCurveInterferences(IcFil2).Append(Interfp2);
      }
      if (ipoin == Ipoin2 || !ChFi3d_Contains(Li,IcFil2,ipoin)) { 
      Interfp4= ChFi3d_FilPointInDS(TopAbs_REVERSED,IcFil2,ipoin,
                              Fi2.LastParameter(), isVertex );
      DStr.ChangeCurveInterferences(IcFil2).Append(Interfp4);
      }
      Ipoin2 = ipoin;
      isVertex2 = isVertex;      
    }
    
    ET1 = trafil1;
    if (j == SeqFil.Length()) {
      if (!isInDS2) {
      Icurv = CorDat->LastCurve();
      if(Closed && !Singulier_en_Bout && (Ipoin1!=Ipoin2)) {
        regcout.SetS2(Isurf,Standard_False);
        reglist.Append(regcout);
      }
      PCurv = CorDat->LastPCurve();
      ET1 = TopAbs::Compose(ET1,CorDat->LastPCurveOrientation());
      CorDat->LastParameters(Pardeb,Parfin);
      TopOpeBRepDS_ListOfInterference& Li = DStr.ChangeCurveInterferences(Icurv);
      if (Li.IsEmpty()) {
        if(CorDat->LastPCurveOrientation()==TopAbs_REVERSED) {
          Interfp5=ChFi3d_FilPointInDS
            (TopAbs_REVERSED,Icurv,Ipoin1,Parfin, isVertex1);
          Interfp6=ChFi3d_FilPointInDS
            (TopAbs_FORWARD,Icurv,Ipoin2,Pardeb, isVertex2);
        }
        else{
          Interfp5=ChFi3d_FilPointInDS
            (TopAbs_FORWARD,Icurv,Ipoin1,Pardeb, isVertex1);
          Interfp6=ChFi3d_FilPointInDS
            (TopAbs_REVERSED,Icurv,Ipoin2,Parfin, isVertex2);
        }
        Li.Append(Interfp5);
        Li.Append(Interfp6);
      }
      Interfc1= ChFi3d_FilCurveInDS(Icurv,Isurf,PCurv,ET1);
      DStr.ChangeSurfaceInterferences(Isurf).Append(Interfc1);
        if (Ipoin1 == Ipoin2) {
        TopOpeBRepDS_Curve& TCurv = DStr.ChangeCurve(Icurv);
        TCurv.ChangeCurve().Nullify();
        Handle(TopOpeBRepDS_Interference) bidinterf;
        TCurv.SetSCI( Interfc1, bidinterf);
//         bidinterf = TCurv.GetSCI1(); 
//      TCurv.SetSCI(bidinterf, Interfc1);          
      }
      }
    }
    else {
//      Degene = (Fd->VertexLastOnS1().Point().IsEqual(
//                Fd->VertexLastOnS2().Point(), 0) );
      
      // eap, Apr 29 2002, occ 293 
      if (!isInDS2) {
      
        Handle(Geom_Curve) C3d;
        Standard_Real tolreached;
        ChFi3d_ComputeArete(Fd->VertexLastOnS1(),
                            Fd->InterferenceOnS1().PCurveOnSurf()->
                            Value(Fd->InterferenceOnS1().LastParameter()),
                            Fd->VertexLastOnS2(),
                            Fd->InterferenceOnS2().PCurveOnSurf()->
                            Value(Fd->InterferenceOnS2().LastParameter()),
                            DStr.Surface(Fd->Surf()).Surface(),C3d,PCurv,
                            Pardeb,Parfin,tol3d,tol2d,tolreached,0);
        Crv = TopOpeBRepDS_Curve(C3d,tolreached);
        Icurv = DStr.AddCurve(Crv);
        regfilfil.SetCurve(Icurv);
        regfilfil.SetS1(Isurf,Standard_False);
        Interfp5 = ChFi3d_FilPointInDS(TopAbs_FORWARD,Icurv,Ipoin1,Pardeb, isVertex1);
        DStr.ChangeCurveInterferences(Icurv).Append(Interfp5);
        Interfp6= ChFi3d_FilPointInDS(TopAbs_REVERSED,Icurv,Ipoin2,Parfin, isVertex2);
        DStr.ChangeCurveInterferences(Icurv).Append(Interfp6);
        Interfc1= ChFi3d_FilCurveInDS(Icurv,Isurf,PCurv,ET1);
        DStr.ChangeSurfaceInterferences(Isurf).Append(Interfc1);      
      }
    }
   
    Degene = V3.Point().IsEqual(V4.Point(), 0);

    // Traitement des cas degenere     
    if (Degene) {
      // Queue de Billard
      Standard_Boolean Vertex = (V3.IsVertex()) && (V4.IsVertex());
      if (!Vertex) {
#if DEB
      if (j == SeqFil.Length()) {
        cout << " Singularite en bout hors vertex " << endl;
        // QUEZAKO ????
      }
#endif
      }
      else {
      // On elimine l'arete de la spine debouchant sur ce vertex.
      Standard_Boolean Trouve = Standard_False;
      TopoDS_Edge Arcspine;
#ifndef DEB
      TopAbs_Orientation OVtx = TopAbs_FORWARD;
#else
      TopAbs_Orientation OVtx;
#endif
      BoutdeVtx = V3.Vertex();
          
      while (NumEdge<= spine->NbEdges() && !Trouve) { 
        Arcspine = spine->Edges(NumEdge);
        for(ex.Init(Arcspine.Oriented(TopAbs_FORWARD),TopAbs_VERTEX); 
            ex.More() && (!Trouve); ex.Next()) {
          if(BoutdeVtx.IsSame(ex.Current())) {
            OVtx = ex.Current().Orientation();
            if (Closed &&  (NumEdge == 1)) 
               Trouve = (spine->NbEdges() == 1);
              else Trouve = Standard_True;
          }
        }
        if (!Trouve) NumEdge++; // On passe a l'arete suivante
      }
      Standard_Integer IArcspine = DStr.AddShape(Arcspine);
      Standard_Integer IVtx;
      if  (j == SeqFil.Length()) {
        IVtx = CorDat->IndexLastPointOnS1();
      }
      else { IVtx = DStr.AddShape(BoutdeVtx); }
      OVtx = TopAbs::Reverse(OVtx);
      Standard_Real parVtx = BRep_Tool::Parameter(BoutdeVtx,Arcspine);
      Handle(TopOpeBRepDS_CurvePointInterference) 
        interfv = ChFi3d_FilVertexInDS(OVtx,IArcspine,IVtx,parVtx);
      DStr.ChangeShapeInterferences(IArcspine).Append(interfv);
      }
    } // fin du cas Degene
    else if (!(Closed && j == SeqFil.Length())) {
      // Traitement des interference Point / Edges
      if (V3.IsOnArc()) {
      if(!(V3.IsVertex() && Fd->IsOnCurve1())) {
        Iarc1 = DStr.AddShape(V3.Arc());
        if ( !ChFi3d_Contains(DStr.ShapeInterferences(Iarc1),Iarc1,Ipoin1,V3.IsVertex(),Standard_True) ) {
          Handle(TopOpeBRepDS_CurvePointInterference) Interfpp = 
            ChFi3d_FilPointInDS(V3.TransitionOnArc(),
                          Iarc1,Ipoin1,V3.ParameterOnArc(), V3.IsVertex());
          DStr.ChangeShapeInterferences(V3.Arc()).Append(Interfpp);
        }
      }
      }

      if (V4.IsOnArc()) {
      if(!(V4.IsVertex() && Fd->IsOnCurve2())) {
        Iarc2 = DStr.AddShape(V4.Arc());
        if ( !ChFi3d_Contains(DStr.ShapeInterferences(Iarc2),Iarc2,Ipoin2,V4.IsVertex(),Standard_True) ) {
          Handle(TopOpeBRepDS_CurvePointInterference) Intfpp=
            ChFi3d_FilPointInDS(V4.TransitionOnArc(),
                          Iarc2,Ipoin2,V4.ParameterOnArc(), V4.IsVertex());
          DStr.ChangeShapeInterferences(V4.Arc()).Append(Intfpp);
        }
      }
      }
    }
  }
}

//=======================================================================
//function : StripeEdgeInter
//purpose  : This function examines two stripes for an intersection 
//           between curves of interference with faces. If the intersection
//           exists, it will cause bad result, so it's better to quit.
//remark   : If someone somewhen computes the interference between stripes, 
//           this function will become useless.
//author   : akm, 06/02/02. Against bug OCC119.
//=======================================================================

void ChFi3d_StripeEdgeInter (const Handle(ChFiDS_Stripe)& theStripe1,
                       const Handle(ChFiDS_Stripe)& theStripe2,
                       TopOpeBRepDS_DataStructure&  /*DStr*/,
                       const Standard_Real          tol2d)
{
  // Do not check the stripeshaving common corner points
  for (Standard_Integer iSur1=1; iSur1<=2; iSur1++)
    for (Standard_Integer iSur2=1; iSur2<=2; iSur2++)
      if (theStripe1->IndexPoint(0,iSur1)==theStripe2->IndexPoint(0,iSur2) ||
        theStripe1->IndexPoint(0,iSur1)==theStripe2->IndexPoint(1,iSur2) ||
        theStripe1->IndexPoint(1,iSur1)==theStripe2->IndexPoint(0,iSur2) ||
        theStripe1->IndexPoint(1,iSur1)==theStripe2->IndexPoint(1,iSur2))
      return;

  Handle(ChFiDS_HData) aSurDat1 = theStripe1->SetOfSurfData();
  Handle(ChFiDS_HData) aSurDat2 = theStripe2->SetOfSurfData();

  Geom2dInt_GInter anIntersector;
  Standard_Integer iPart1, iPart2;
  Standard_Integer Ishape11, Ishape12, Ishape21, Ishape22;
  // Loop on parts of the first stripe
  for (iPart1=1; iPart1<=aSurDat1->Length(); iPart1++) 
  {
    Handle(ChFiDS_SurfData) aDat1 = aSurDat1->Value(iPart1);
    Ishape11 = aDat1->IndexOfS1();
    Ishape12 = aDat1->IndexOfS2();
    // Loop on parts of the second stripe
    for (iPart2=1; iPart2<=aSurDat2->Length(); iPart2++) 
    {
      Handle(ChFiDS_SurfData) aDat2 = aSurDat2->Value(iPart2);
      Ishape21 = aDat2->IndexOfS1();
      Ishape22 = aDat2->IndexOfS2();

      // Find those FaceInterferences able to intersect
      ChFiDS_FaceInterference aFI1, aFI2;
      if (Ishape11 == Ishape21)
      {
      aFI1 = aDat1->InterferenceOnS1();
      aFI2 = aDat2->InterferenceOnS1();
      } 
      else if (Ishape11 == Ishape22)
      {
      aFI1 = aDat1->InterferenceOnS1();
      aFI2 = aDat2->InterferenceOnS2();
      } 
      else if (Ishape12 == Ishape21)
      {
      aFI1 = aDat1->InterferenceOnS2();
      aFI2 = aDat2->InterferenceOnS1();
      } 
      else if (Ishape12 == Ishape22)
      {
      aFI1 = aDat1->InterferenceOnS2();
      aFI2 = aDat2->InterferenceOnS2();
      }
      else
      {
      // No common faces
      continue;
      }

      if (IsEqual (aFI1.FirstParameter(),aFI1.LastParameter()) ||
        IsEqual (aFI2.FirstParameter(),aFI2.LastParameter()) ||
        aFI1.PCurveOnFace().IsNull() ||
        aFI2.PCurveOnFace().IsNull())
      // Do not waste time on degenerates
      continue;
      // Examine for intersections
      Geom2dAdaptor_Curve aPCurve1 (aFI1.PCurveOnFace(),
                            aFI1.FirstParameter(),
                            aFI1.LastParameter());
      Geom2dAdaptor_Curve aPCurve2 (aFI2.PCurveOnFace(),
                            aFI2.FirstParameter(),
                            aFI2.LastParameter());
      anIntersector.Perform (aPCurve1,
                       aPCurve2,
                       tol2d,
                       Precision::PConfusion());
      if (anIntersector.NbSegments() > 0 ||
        anIntersector.NbPoints() > 0)
      StdFail_NotDone::Raise ("StripeEdgeInter : fillets have too big radiuses");
    }
  }
}

//=======================================================================
//function : IndexOfSurfData
//purpose  : 
//=======================================================================

Standard_Integer ChFi3d_IndexOfSurfData(const TopoDS_Vertex& V1,
                              const Handle(ChFiDS_Stripe)& CD,
                              Standard_Integer& sens)
{
  Handle(ChFiDS_Spine) spine = CD->Spine();
  Standard_Integer Index = 0;
  sens = 1;
  TopoDS_Vertex Vref;
  const TopoDS_Edge& E = spine->Edges(1);
  if (E.Orientation() == TopAbs_REVERSED) Vref = TopExp::LastVertex(E);
  else Vref = TopExp::FirstVertex(E); 
  if (Vref.IsSame(V1)) Index =1;
  else {
    const TopoDS_Edge& E1 = spine->Edges(spine->NbEdges());
    if (E1.Orientation() == TopAbs_REVERSED) Vref = TopExp::FirstVertex(E1);
    else Vref = TopExp::LastVertex(E1); 
    sens = -1;
    if(CD->SetOfSurfData().IsNull()) return 0;
    else if (Vref.IsSame(V1)) Index = CD->SetOfSurfData()->Length();
    else Standard_ConstructionError::Raise("");
  }
  return Index; 
}  

//=======================================================================
//function : EdgeFromV1
//purpose  : 
//=======================================================================

TopoDS_Edge ChFi3d_EdgeFromV1(const TopoDS_Vertex& V1,
                        const Handle(ChFiDS_Stripe)& CD,
                        Standard_Integer& sens)
{
  Handle(ChFiDS_Spine) spine = CD->Spine();
  sens = 1;
  TopoDS_Vertex Vref;
  const TopoDS_Edge& E = spine->Edges(1);
  if (E.Orientation() == TopAbs_REVERSED) Vref = TopExp::LastVertex(E);
  else Vref = TopExp::FirstVertex(E); 
  if (Vref.IsSame(V1)) return E;
  else
    {
      const TopoDS_Edge& E1 = spine->Edges(spine->NbEdges());
      if (E1.Orientation() == TopAbs_REVERSED) Vref = TopExp::FirstVertex(E1);
      else Vref = TopExp::LastVertex(E1); 
      sens = -1;
      if (Vref.IsSame(V1)) return E1;
      else Standard_ConstructionError::Raise("");
    }
  return E;
}

#ifdef DRAW
//=======================================================================
//function : drawline
//purpose  : 
//=======================================================================

static void drawline(const Handle(IntSurf_LineOn2S)& lin)
{
  Handle(Draw_Marker3D) p3d;
  Handle(Draw_Marker2D) p2d;
  Handle(Draw_Marker2D) p2d2;
  
  for(Standard_Integer i = 1; i <= lin->NbPoints(); i++) {
    const IntSurf_PntOn2S& pt = lin->Value(i);
    gp_Pnt point = pt.Value();
    Standard_Real u,v,uu,vv;
    pt.ParametersOnS1(u,v);
    gp_Pnt2d point2d(u,v);
    pt.ParametersOnS2(uu,vv);
    gp_Pnt2d point2d2(uu,vv);
    if(i == 1 || i == lin->NbPoints()) {
      p3d = new Draw_Marker3D(point,Draw_Square,Draw_rouge);
      p2d = new Draw_Marker2D(point2d,Draw_Square,Draw_rouge);
      p2d2 = new Draw_Marker2D(point2d2,Draw_Square,Draw_rouge);
    }
    else {
      p3d = new Draw_Marker3D(point,Draw_Square,Draw_jaune);
      p2d = new Draw_Marker2D(point2d,Draw_Square,Draw_jaune);
      p2d2 = new Draw_Marker2D(point2d2,Draw_Square,Draw_jaune);
    }
    dout<<p3d;
    dout<<p2d;
    dout<<p2d2;
    
    dout.Flush();
  }
}
#endif

//=======================================================================
//function : ConvTol2dToTol3d
//purpose  : Comme son nom l indique.
//=======================================================================

Standard_Real ChFi3d_ConvTol2dToTol3d(const Handle(Adaptor3d_HSurface)& S,
                              const Standard_Real             tol2d)
{
  Standard_Real ures = S->UResolution(1.e-7);
  Standard_Real vres = S->VResolution(1.e-7);
  Standard_Real uresto3d = 1.e-7*tol2d/ures;
  Standard_Real vresto3d = 1.e-7*tol2d/vres;
  return Max(uresto3d,vresto3d);
}

//=======================================================================
//function : EvalTolReached
//purpose  : Comme son nom l indique la fonction ci dessus etant trop
//           dure lorsque le parametrage des surfaces n est pas homogene.
//=======================================================================

Standard_Real ChFi3d_EvalTolReached(const Handle(Adaptor3d_HSurface)& S1,
                            const Handle(Geom2d_Curve)&     pc1,
                            const Handle(Adaptor3d_HSurface)& S2,
                            const Handle(Geom2d_Curve)&     pc2,
                            const Handle(Geom_Curve)&       C)
{
  Standard_Real distmax = 0.;

  Standard_Real f = C->FirstParameter();
  Standard_Real l = C->LastParameter();
  Standard_Integer nbp = 45;
  Standard_Real step = 1./(nbp -1);
  for(Standard_Integer i = 0; i < nbp; i++) {
    Standard_Real t,u,v;
    t = step * i;
    t = (1-t) * f + t * l;
    pc1->Value(t).Coord(u,v);
    gp_Pnt pS1 = S1->Value(u,v);
    pc2->Value(t).Coord(u,v);
    gp_Pnt pS2 = S2->Value(u,v);
    gp_Pnt pC = C->Value(t);
    Standard_Real d = pS1.SquareDistance(pC);
    if(d>distmax) distmax = d;
    d = pS2.SquareDistance(pC);
    if(d>distmax) distmax = d;
    d = pS1.SquareDistance(pS2);
    if(d>distmax) distmax = d;
  }
  distmax = 1.5*sqrt(distmax);
  distmax = Max(distmax, Precision::Confusion());
  return distmax;
}

//=======================================================================
//function : trsfsurf
//purpose  : 
//=======================================================================

Handle(Geom_Surface) trsfsurf(const Handle(Adaptor3d_HSurface)& HS,
                        Handle(Adaptor3d_TopolTool)&      /*dom*/)
{
  //Pour l utilisation des domaines voir avec BUBUCH!!
  Handle(Geom_Surface) res;
  Handle(BRepAdaptor_HSurface) hbs = Handle(BRepAdaptor_HSurface)::DownCast(HS);
  Handle(GeomAdaptor_HSurface) hgs = Handle(GeomAdaptor_HSurface)::DownCast(HS);
  if(!hbs.IsNull()) {
    res = hbs->ChangeSurface().Surface().Surface();
    gp_Trsf trsf = hbs->ChangeSurface().Trsf();
    res = Handle(Geom_Surface)::DownCast(res->Transformed(trsf));
  }
  else if(!hgs.IsNull()) {
    res = hgs->ChangeSurface().Surface();
  }
  Handle(Geom_RectangularTrimmedSurface) 
    tr = Handle(Geom_RectangularTrimmedSurface)::DownCast(res);
  if(!tr.IsNull()) res = tr->BasisSurface();

  Standard_Real U1 = HS->FirstUParameter(), U2 = HS->LastUParameter();
  Standard_Real V1 = HS->FirstVParameter(), V2 = HS->LastVParameter();
  if(!res.IsNull()) {
  // Blindage contre les Construction Error intempestifs
    Standard_Real u1, u2, v1, v2;
    res->Bounds( u1, u2, v1, v2);
    if (!res->IsUPeriodic()) {
      if (U1 < u1) U1 = u1;
      if (U2 > u2) U2 = u2;
    }
    if (!res->IsVPeriodic()) {
      if (V1 < v1) V1 = v1;
      if (V2 > v2) V2 = v2;
    }
    res = new Geom_RectangularTrimmedSurface(res,U1,U2,V1,V2);
  }
//  Handle(GeomAdaptor_HSurface) temp = new GeomAdaptor_HSurface(res,U1,U2,V1,V2);
//  dom = new Adaptor3d_TopolTool(temp);
  return res;
}

//=======================================================================
//function : CurveCleaner
//purpose  : Rend une BSpline le plus continue possible
//           a une tolerance donne
//=======================================================================
static void CurveCleaner(Handle(Geom_BSplineCurve)& BS, 
                   const Standard_Real Tol,
                   const Standard_Integer MultMin)

{
  Standard_Real tol = Tol;
  Standard_Integer Mult, ii, NbK=BS->NbKnots();
  Standard_Boolean Bof;
  

  for (Mult = BS->Degree(); Mult > MultMin; Mult--) {
    tol /= 2; // Reduction progressive
    for (ii=NbK; ii>1; ii--) {
      if (BS->Multiplicity(ii) == Mult) {
      Bof = BS->RemoveKnot(ii, Mult-1, tol);
      }      
    }
  }
}

//=======================================================================
//function : ComputeCurves
//purpose  : Calcule une intersection bornee entre deux HSurfaces.
//           Il faut connaitre les extremites de l intersection et 
//           les surfaces doivent avoir ete retouchees en entree 
//           pour encadrer au mieux (ni trop pres ni trop loin) les
//           points de debut et fin de l intersection.
//           Les intersections analytiques sont traitees a part.
//           <wholeCurv> means that resulting curve is restricted by
//           boundaries of input surfaces (eap 30 May occ354)
//=======================================================================

Standard_Boolean ChFi3d_ComputeCurves(Handle(Adaptor3d_HSurface)&   S1,
                              Handle(Adaptor3d_HSurface)&   S2,
                              const TColStd_Array1OfReal& Pardeb,
                              const TColStd_Array1OfReal& Parfin,
                              Handle(Geom_Curve)&         C3d,
                              Handle(Geom2d_Curve)&       Pc1,
                              Handle(Geom2d_Curve)&       Pc2,
                              const Standard_Real         tol3d,
                              const Standard_Real         tol2d,
                              Standard_Real&              tolreached,
                              const Standard_Boolean      wholeCurv) 
{
  Standard_Real Step = 0.1;

  gp_Pnt pdeb1 = S1->Value(Pardeb(1),Pardeb(2));
  gp_Pnt pfin1 = S1->Value(Parfin(1),Parfin(2));
  gp_Pnt pdeb2 = S2->Value(Pardeb(3),Pardeb(4));
  gp_Pnt pfin2 = S2->Value(Parfin(3),Parfin(4));

  Standard_Real distrefdeb = pdeb1.Distance(pdeb2);//mesure la solidite
  Standard_Real distreffin = pfin1.Distance(pfin2);//des donnees d entree
  if(distrefdeb < tol3d) distrefdeb = tol3d;
  if(distreffin < tol3d) distreffin = tol3d;

  gp_Pnt pdeb,pfin;
  pdeb.SetXYZ(0.5*(pdeb1.XYZ()+pdeb2.XYZ()));
  pfin.SetXYZ(0.5*(pfin1.XYZ()+pfin2.XYZ()));

  Standard_Real distref = 0.005*pdeb.Distance(pfin);
  if(distref < distrefdeb) distref = distrefdeb;
  if(distref < distreffin) distref = distreffin;

  //On traite a part quelques cas analytiques.
  //Pour reorienter eventuellement le resultat de l intersection
  //analytique, on postule que la tangente en debut doit etre dans
  //le sens de la corde deb/fin.
  gp_Vec Vint, Vref(pdeb,pfin);
  gp_Pnt Pbid;
  Standard_Real Udeb,Ufin;
  Standard_Real tolr1,tolr2;
  tolr1 = tolr2 = tolreached = tol3d;
  if((S1->GetType() == GeomAbs_Cylinder && S2->GetType() == GeomAbs_Plane)||
     (S1->GetType() == GeomAbs_Plane && S2->GetType() == GeomAbs_Cylinder)) { 
    gp_Pln pl;
    gp_Cylinder cyl;
    if(S1->GetType() == GeomAbs_Plane) {
      pl = S1->Plane();
      cyl = S2->Cylinder();
    }
    else{
      pl = S2->Plane();
      cyl = S1->Cylinder();
    }
    IntAna_QuadQuadGeo ImpKK(pl,cyl,Precision::Angular(),tol3d);
    if (ImpKK.IsDone()) {
      Standard_Boolean c1line = 0;
      switch  (ImpKK.TypeInter()) {
      case IntAna_Line:
      {
        c1line = 1;
        Standard_Integer nbsol = ImpKK.NbSolutions();
        gp_Lin C1;
        for(Standard_Integer ilin = 1; ilin <= nbsol; ilin++) {
          C1 = ImpKK.Line(ilin);
          Udeb = ElCLib::Parameter(C1,pdeb);
          gp_Pnt ptest = ElCLib::Value(Udeb,C1);
          if(ptest.Distance(pdeb) < tol3d) break;
        }
        Ufin = ElCLib::Parameter(C1,pfin);
        C3d = new Geom_Line(C1);
        ElCLib::D1(Udeb,C1,Pbid,Vint);
      }
      break;
      case IntAna_Circle:
      {
        gp_Circ C1 = ImpKK.Circle(1);
        C3d = new Geom_Circle(C1);
        Udeb = ElCLib::Parameter(C1,pdeb);
        Ufin = ElCLib::Parameter(C1,pfin);
        ElCLib::D1(Udeb,C1,Pbid,Vint);
      }
      break;
      case IntAna_Ellipse:
      {
        gp_Elips C1 = ImpKK.Ellipse(1);
        C3d = new Geom_Ellipse(C1);
        Udeb = ElCLib::Parameter(C1,pdeb);
        Ufin = ElCLib::Parameter(C1,pfin);
        ElCLib::D1(Udeb,C1,Pbid,Vint);
      }
      break;
      default:
      break;
      }
      if (Vint.Dot(Vref)<0) {
      C3d->Reverse();
      if(c1line) {
        Udeb = -Udeb;
        Ufin = -Ufin;
      }
      else{
        Udeb = 2*PI - Udeb;
        Ufin = 2*PI - Ufin;
      }
      }
      if(!c1line) ElCLib::AdjustPeriodic(0.,2*PI,Precision::Angular(),Udeb,Ufin);
      Handle(GeomAdaptor_HCurve) HC = new GeomAdaptor_HCurve();
      HC->ChangeCurve().Load(C3d,Udeb,Ufin);
      ChFi3d_ProjectPCurv(HC,S1,Pc1,tol3d,tolr1);
      if(S1->GetType() == GeomAbs_Cylinder) {
      Standard_Real x,y;
      Pc1->Value(Udeb).Coord(x,y);
      x = Pardeb(1) - x;
      y = Pardeb(2) - y;
      if(Abs(x) >= tol2d || Abs(y) >= tol2d) Pc1->Translate(gp_Vec2d(x,y));
      }
      ChFi3d_ProjectPCurv(HC,S2,Pc2,tol3d,tolr2);
      if(S2->GetType() == GeomAbs_Cylinder) {
      Standard_Real x,y;
      Pc2->Value(Udeb).Coord(x,y);
      x = Pardeb(3) - x;
      y = Pardeb(4) - y;
      if(Abs(x) >= tol2d || Abs(y) >= tol2d) Pc2->Translate(gp_Vec2d(x,y));
      }
      C3d = new Geom_TrimmedCurve(C3d,Udeb,Ufin);
      tolreached = 1.5*Max(tolr1,tolr2);
      tolreached = Min(tolreached,ChFi3d_EvalTolReached(S1,Pc1,S2,Pc2,C3d));
      return Standard_True;
    }
  }    
  else if(S1->GetType() == GeomAbs_Plane && S2->GetType() == GeomAbs_Plane) { 
    IntAna_QuadQuadGeo LInt(S1->Plane(),S2->Plane(),Precision::Angular(),tol3d);
    if (LInt.IsDone()) {
      gp_Lin L = LInt.Line(1);
      C3d = new Geom_Line(L);
      Udeb = ElCLib::Parameter(L,pdeb);
      Ufin = ElCLib::Parameter(L,pfin);
      ElCLib::D1(Udeb,L,Pbid,Vint);
      if (Vint.Dot(Vref)<0) {
      C3d->Reverse();
      Udeb = - Udeb;
      Ufin = - Ufin;
      }
      Handle(GeomAdaptor_HCurve) HC = new GeomAdaptor_HCurve();
      HC->ChangeCurve().Load(C3d,Udeb,Ufin);
      ChFi3d_ProjectPCurv(HC,S1,Pc1,tol3d,tolr1);
      ChFi3d_ProjectPCurv(HC,S2,Pc2,tol3d,tolr2);
      C3d = new Geom_TrimmedCurve(C3d,Udeb,Ufin);
      return Standard_True;
    }
  }
  else {
    // ici on attaque GeomInt.
    //Pour l utilisation des domaines voir avec BUBUCH!!
    Handle(Adaptor3d_TopolTool) dom1,dom2;
    Handle(Geom_Surface) gs1 = trsfsurf(S1,dom1);
    Handle(Geom_Surface) gs2 = trsfsurf(S2,dom2);
    Standard_Integer nbl ;
    if(!gs1.IsNull() && !gs2.IsNull()) {
      GeomInt_IntSS inter;
      //  Modified by skv - Fri Oct 24 14:24:47 2003 OCC4077 Begin
//       Standard_Real tolap = 1.e-7;//car l approx de la wline est faite dans [0,1]
      // Set the lowest tolerance which is used in new boolean operations.
      Standard_Real tolap = 2.e-7;
      //  Modified by skv - Fri Oct 24 14:24:48 2003 OCC4077 End
      inter.Perform(gs1,gs2,tolap,1,1,1);
      if(inter.IsDone()) {
      nbl = inter.NbLines();   
#if defined(IRIX) || defined(__sgi)
        if(nbl==0) {

//  solution de rattrapage pour SGI 
//  si l'intersection de gs1 avec gs2 ne marche pas alors on  tente 
// l'intersection de gs2 avec gs1 

          inter.Perform(gs2,gs1,tolap,1,1,1);
//          inter.Perform(gs2,dom2,gs1,dom1,tolap,1,1,1);
          if(!inter.IsDone()) return Standard_False; 
        nbl = inter.NbLines(); 

//  si GeomInt ne rend pas d'intersection on ne tente pas la solution de 
//  rattrapage 
            if (nbl==0) return Standard_False;
        }
#endif
      GeomAPI_ProjectPointOnCurve proj;
      for(Standard_Integer ilin = 1; ilin <= nbl; ilin++) {
        if(inter.HasLineOnS1(ilin) && inter.HasLineOnS2(ilin)) {
          C3d = inter.Line(ilin);
          Pc1 = inter.LineOnS1(ilin);
          Pc2 = inter.LineOnS2(ilin);
          gp_Pnt ptestdeb, ptestfin;
          Standard_Real Uf=0., Ul=0.;
          if (wholeCurv) {
            Uf = C3d->FirstParameter();
            Ul = C3d->LastParameter();
            ptestdeb = C3d->Value(Uf);
            ptestfin = C3d->Value(Ul);
          }
          else {
              // find end parameters
            Standard_Boolean failedF, failedL;
              failedF = failedL = Standard_False;
              proj.Init( pdeb1, C3d);
              if (proj.NbPoints()==0 && distrefdeb > Precision::Confusion())
                proj.Perform( pdeb2 );
              if (proj.NbPoints()==0)
                failedF = Standard_True;
              else
                Uf = proj.LowerDistanceParameter();
              proj.Perform( pfin1 );
              if (proj.NbPoints()==0 && distreffin > Precision::Confusion())
                proj.Perform( pfin2 );
              if (proj.NbPoints()==0) 
                failedL = Standard_True;
              else
                Ul = proj.LowerDistanceParameter();

              if (failedF && failedL) {
                Uf = C3d->FirstParameter();
                Ul = C3d->LastParameter();
              }
              else if (failedF || failedL) {
                // select right end parameter
                Standard_Real Uok = failedF ? Ul : Uf;
                Standard_Real U1 = C3d->FirstParameter(), U2 = C3d->LastParameter();
                Uok = Abs(Uok-U1) > Abs(Uok-U2) ? U1 : U2;
                if (failedF) Uf = Uok;
                else         Ul = Uok;
              }
              else { // both projected, but where?
                if (Uf == Ul) continue;
              }
              ptestdeb = C3d->Value(Uf);
              ptestfin = C3d->Value(Ul);
              if (C3d->IsPeriodic() && !(failedF && failedL)) {
                // assure the same order of ends, otherwise TrimmedCurve will take
                // the other part of C3d
                gp_Pnt Ptmp;
                gp_Vec DirOld, DirNew(ptestdeb,ptestfin);
                C3d->D1(Uf, Ptmp, DirOld);
                if (DirOld * DirNew < 0) {
                  Standard_Real Utmp = Uf; Uf = Ul; Ul = Utmp;
                  Ptmp = ptestdeb; ptestdeb = ptestfin; ptestfin = Ptmp;
                }
              }
            }
          C3d = new Geom_TrimmedCurve(C3d,Uf,Ul);
          Pc1 = new Geom2d_TrimmedCurve(Pc1,Uf,Ul);
          Pc2 = new Geom2d_TrimmedCurve(Pc2,Uf,Ul);
          //faut il renverser ?
          Standard_Real distdeb = ptestdeb.Distance(pdeb);
          Standard_Real distfin = ptestfin.Distance(pfin);
          if(distdeb > distref || distfin > distref) {
            C3d->Reverse();
            Pc1->Reverse();
            Pc2->Reverse();
            ptestdeb = C3d->Value(C3d->FirstParameter());
            ptestfin = C3d->Value(C3d->LastParameter());
            distdeb = ptestdeb.Distance(pdeb);
            distfin = ptestfin.Distance(pfin);
          }
          if(distdeb < distref && distfin < distref) {
            Uf = C3d->FirstParameter();
            Ul = C3d->LastParameter();
            ChFi3d_ReparamPcurv(Uf,Ul,Pc1);
            ChFi3d_ReparamPcurv(Uf,Ul,Pc2);
            Standard_Real x,y;
            Pc1->Value(Uf).Coord(x,y);
            x = Pardeb(1) - x;
            y = Pardeb(2) - y;
            if(Abs(x) > tol2d || Abs(y) > tol2d) Pc1->Translate(gp_Vec2d(x,y));
            Pc2->Value(Uf).Coord(x,y);
            x = Pardeb(3) - x;
            y = Pardeb(4) - y;
            if(Abs(x) > tol2d || Abs(y) > tol2d) Pc2->Translate(gp_Vec2d(x,y));
            tolreached = ChFi3d_EvalTolReached(S1,Pc1,S2,Pc2,C3d);
            return Standard_True;
          }
        }
      }
      }
    }
  }
#ifdef DEB
  cout<<"echec intersection surf/surf on tente le rattrapage a la main"<<endl;
#endif
  
  // A ce stade : 
  // les intersections classiques ont echouees on attaque le
  // cheminement du desespoir.
//  Standard_Real Step = 0.1;
  while(1) {
    //Attention les parametres de fleche pour le cheminement et 
    //de tolerance pour l'approx ne peuvent etre pris comme ceux 
    //du Builder, on les reestime donc comme on peut.
    Standard_Real fleche = 1.e-3 * pdeb.Distance(pfin);
    Standard_Real tolap = 1.e-7;
    IntPatch_ThePWalkingInterOfThePPIntOfIntersection 
      IntKK(S1,S2,tol3d,tol3d,fleche,Step);
    
    //On connait les extremites de l intersection (Pardeb,Parfin),
    //on essaye de trouver un point de depart franchement au
    //milieu pour ne pas embrouiller le cheminement.
    Standard_Boolean depok = Standard_False;
    IntSurf_PntOn2S pintdep;
    TColStd_Array1OfReal depart(1,4);
    for(Standard_Integer ipdep = 2; ipdep <= 7 && !depok; ipdep++) {
      Standard_Real alpha = 0.1 * ipdep;
      Standard_Real unmoinsalpha = 1. - alpha;
      depart(1) = alpha*Pardeb(1) + unmoinsalpha*Parfin(1);
      depart(2) = alpha*Pardeb(2) + unmoinsalpha*Parfin(2);
      depart(3) = alpha*Pardeb(3) + unmoinsalpha*Parfin(3);
      depart(4) = alpha*Pardeb(4) + unmoinsalpha*Parfin(4);
      depok = IntKK.PerformFirstPoint(depart,pintdep);
    } 
    if(!depok) {
#ifdef DEB 
      cout<<"intersection surf/surf : echec calcul point de depart"<<endl;
#endif
      return Standard_False;
    }
    pintdep.Parameters(depart(1),depart(2),depart(3),depart(4));
    IntKK.Perform(depart);
    if (!IntKK.IsDone()) return Standard_False;
    if (IntKK.NbPoints() <= 30) {
      Step *= 0.5;
      if (Step <= 0.0001) {
#ifdef DEB 
      cout<<"intersection surf/surf : pas assez de points"<<endl;
#endif
      return Standard_False;
      }
    }
    else{
      // A ce stade on a une LineOn2S presentable, on la tronque entre les 
      // points les plus proches des extremites connues on en fait une 
      // WLine et on lance l approx.
      // On retouche ensuite le resultat pour avoir les bons points debut 
      // et fin.
      const Handle(IntSurf_LineOn2S)& L2S = IntKK.Line();
      
      gp_Pnt codeb1 = S1->Value(Pardeb(1),Pardeb(2));
      gp_Pnt codeb2 = S2->Value(Pardeb(3),Pardeb(4));
      Standard_Real tol1 = Max(codeb1.Distance(codeb2),tol3d);
      Standard_Boolean bondeb = (tol1 == tol3d); 
      gp_Pnt pntd(0.5*(codeb1.Coord() + codeb2.Coord()));
      
      gp_Pnt cofin1 = S1->Value(Parfin(1),Parfin(2));
      gp_Pnt cofin2 = S2->Value(Parfin(3),Parfin(4));
      Standard_Real tol2 = Max(cofin1.Distance(cofin2),tol3d);
      Standard_Boolean bonfin = (tol2 == tol3d); 
      gp_Pnt pntf(0.5*(cofin1.Coord() + cofin2.Coord()));
      
      Standard_Integer nbp = L2S->NbPoints(), i;
      Standard_Real ddeb = Precision::Infinite(); 
      Standard_Real dfin = Precision::Infinite();
      Standard_Real dd;
#ifndef DEB
      Standard_Integer indd = 0, indf = 0;
#else
      Standard_Integer indd, indf;
#endif
      for(i = 1; i <= nbp; i++) {
      dd = L2S->Value(i).Value().Distance(pntd);
      if(dd < ddeb) { ddeb = dd; indd = i;}
      dd = L2S->Value(i).Value().Distance(pntf);
      if(dd < dfin) { dfin = dd; indf = i;}
      }
      if(indd > indf) {
      L2S->Reverse();
      indd = nbp - indd + 1;
      indf = nbp - indf + 1;
      }
      for (i = 1; i < indd; i++) { L2S->RemovePoint(1); nbp--; indf--; }
      for (i = indf + 1; i <= nbp; i++) { L2S->RemovePoint(indf + 1); }
      nbp = indf;
      if(nbp==1) return Standard_False;
      //On insere les extremites dans la ligne si les points extremites de
      //celle-ci en sont trop eloignes et si pardeb et parfin sont bons.
      if(ddeb >= tol3d && bondeb) {
      IntSurf_PntOn2S p1 = L2S->Value(1);
      IntSurf_PntOn2S p2 = L2S->Value(2);
         
      gp_Vec v1(pntd,p1.Value());
      gp_Vec v2(p1.Value(),p2.Value());
      gp_Vec v3(pntd,p2.Value());
      p1.SetValue(pntd,Pardeb(1),Pardeb(2),Pardeb(3),Pardeb(4));
      if(v1.Dot(v3) < 0) {
        if(v3.Magnitude() < 0.2*v2.Magnitude()) {
          L2S->RemovePoint(1);
          nbp--;
        }
        L2S->Value(1,p1);
      }
      else if(v1.Magnitude() > 0.2*v2.Magnitude()) {
        L2S->InsertBefore(1,p1);
        nbp++;
      }
      else{
        L2S->Value(1,p1);
      }
      ddeb = 0.;
      }
      if(dfin >= tol3d && bonfin) {
      IntSurf_PntOn2S p1 = L2S->Value(nbp);
      IntSurf_PntOn2S p2 = L2S->Value(nbp - 1);
      gp_Vec v1(pntf,p1.Value());
      gp_Vec v2(p1.Value(),p2.Value());
      gp_Vec v3(pntf,p2.Value());
      p1.SetValue(pntf,Parfin(1),Parfin(2),Parfin(3),Parfin(4));
      if(v1.Dot(v3) < 0) {
        if(v3.Magnitude() < 0.2*v2.Magnitude()) {
          L2S->RemovePoint(nbp);
          nbp--;
        }
        L2S->Value(nbp,p1);
      }
      else if(v1.Magnitude() > 0.2*v2.Magnitude()) {
        L2S->Add(p1);
        nbp++;
      }
      else{
        L2S->Value(nbp,p1);
      }
      dfin = 0.;
      }      
#ifdef DRAW
      if(ChFi3d_GettraceDRAWINT()) drawline(L2S);
#endif
      Handle(IntPatch_TheWLineOfIntersection) 
      WL = new IntPatch_TheWLineOfIntersection(L2S,Standard_False);
      
      GeomInt_WLApprox approx;
      approx.SetParameters(tolap,tol2d,4,8,0,1);
      // gerer ici les approx inutiles sur les plans!!!!!!!!!!!
      approx.Perform(S1,S2,WL,
                 Standard_True,Standard_True,Standard_True,
                 1,nbp);
      if(!approx.IsDone()) return Standard_False;
//      tolreached = approx.TolReached3d();
//      Standard_Real tolr2d = approx.TolReached2d();
//      tolreached = Max(tolreached,ChFi3d_ConvTol2dToTol3d(S1,tolr2d));
//      tolreached = Max(tolreached,ChFi3d_ConvTol2dToTol3d(S2,tolr2d));
      const AppParCurves_MultiBSpCurve& mbs = approx.Value(1);
      Standard_Integer nbpol = mbs.NbPoles();
      TColgp_Array1OfPnt pol3d(1,nbpol);
      mbs.Curve(1,pol3d);
      TColgp_Array1OfPnt2d pol2d1(1,nbpol);
      mbs.Curve(2,pol2d1);
      TColgp_Array1OfPnt2d pol2d2(1,nbpol);
      mbs.Curve(3,pol2d2);
      // On recale les extremites de l intersection sur les points connus.
      if(ddeb >= tol1) {
      pol3d(1) = pntd;
      pol2d1(1).SetCoord(Pardeb(1),Pardeb(2));
      pol2d2(1).SetCoord(Pardeb(3),Pardeb(4));
//    tolreached = Max(tolreached,ddeb);
      }
      
      if(dfin >= tol2) {
      pol3d(nbpol) = pntf;
      pol2d1(nbpol).SetCoord(Parfin(1),Parfin(2));
      pol2d2(nbpol).SetCoord(Parfin(3),Parfin(4));
//    tolreached = Max(tolreached,dfin);
      }
      const TColStd_Array1OfReal& knots = mbs.Knots();
      const TColStd_Array1OfInteger& mults = mbs.Multiplicities();
      Standard_Integer deg = mbs.Degree();
      C3d = new Geom_BSplineCurve(pol3d,knots,mults,deg);
      Pc1 = new Geom2d_BSplineCurve(pol2d1,knots,mults,deg);
      Pc2 = new Geom2d_BSplineCurve(pol2d2,knots,mults,deg);
      tolreached = ChFi3d_EvalTolReached(S1,Pc1,S2,Pc2,C3d);
      return Standard_True;
    }
  }
}

//=======================================================================
//function : IntCS
//purpose  : Calcul rapide de l intersection courbe surface.
//
//=======================================================================

Standard_Boolean ChFi3d_IntCS(Handle(Adaptor3d_HSurface)& S,
                        Handle(Adaptor3d_HCurve)& C,
                        gp_Pnt2d& p2dS,
                        Standard_Real& wc)
{
  IntCurveSurface_HInter Intersection;
  
  Standard_Real uf = C->FirstParameter(), ul = C->LastParameter();
  Standard_Real u1 = S->FirstUParameter(), u2 = S->LastUParameter();
  Standard_Real v1 = S->FirstVParameter(), v2 = S->LastVParameter();
  IntCurveSurface_IntersectionPoint pint;
  Intersection.Perform(C,S);
  Standard_Boolean keepfirst = (wc < -1.e100), keeplast = (wc > 1.e100);
#ifndef DEB
  Standard_Real temp = 0.;
#else
  Standard_Real temp;
#endif
  if(keepfirst) temp = 1.e100;
  if(keeplast) temp = -1.e100;
  Standard_Real dist = 2.e100;
  if(Intersection.IsDone()) {
    Standard_Integer nbp = Intersection.NbPoints(),i,isol = 0;
    for (i = 1; i <= nbp; i++) {
      pint = Intersection.Point(i);
      Standard_Real up = pint.U();
      Standard_Real vp = pint.V();
      if(S->IsUPeriodic()) up = ChFi3d_InPeriod(up,u1,u1+S->UPeriod(),1.e-8);
      if(S->IsVPeriodic()) vp = ChFi3d_InPeriod(vp,v1,v1+S->VPeriod(),1.e-8);
      if(uf <= pint.W() && ul >= pint.W() &&
       u1 <= up && u2 >= up &&
       v1 <= vp && v2 >= vp) {
      if(keepfirst && pint.W() < temp) {
        temp = pint.W();
        isol = i;
      }
      else if(keeplast && pint.W() > temp) {
        temp = pint.W();
        isol = i;
      }
      else if(Abs(pint.W() - wc) < dist) {
        dist = Abs(pint.W() - wc);
        isol = i;
      }
      }
    }
    if(isol == 0) return Standard_False;
    pint = Intersection.Point(isol);
    Standard_Real up = pint.U();
    Standard_Real vp = pint.V();
    if(S->IsUPeriodic()) up = ChFi3d_InPeriod(up,u1,u1+S->UPeriod(),1.e-8);
    if(S->IsVPeriodic()) vp = ChFi3d_InPeriod(vp,v1,v1+S->VPeriod(),1.e-8);
    p2dS.SetCoord(up,vp);
    wc = pint.W();
    return Standard_True;
  }
  return Standard_False;
}

//=======================================================================
//function : ComputesIntPC
//purpose  : Intersection de deux PCurves de type FaceInterference
//           les parametres sur les pcurves du point solution sont 
//           UInt1,UInt2
//=======================================================================

void ChFi3d_ComputesIntPC (const ChFiDS_FaceInterference&      Fi1,
                     const ChFiDS_FaceInterference&      Fi2,
                     const Handle(GeomAdaptor_HSurface)& HS1,
                     const Handle(GeomAdaptor_HSurface)& HS2,
                     Standard_Real&                      UInt1, 
                     Standard_Real&                      UInt2)
{
  gp_Pnt bid;
  ChFi3d_ComputesIntPC(Fi1,Fi2,HS1,HS2,UInt1,UInt2,bid);
}

void ChFi3d_ComputesIntPC (const ChFiDS_FaceInterference&      Fi1,
                     const ChFiDS_FaceInterference&      Fi2,
                     const Handle(GeomAdaptor_HSurface)& HS1,
                     const Handle(GeomAdaptor_HSurface)& HS2,
                     Standard_Real&                      UInt1, 
                     Standard_Real&                      UInt2,
                     gp_Pnt&                             P)
{    
  // Une seule intersection a realiser, on prend tout de meme
  // le soin de valider les extremites par un extrema c3d/c3d
  // realise sur les pcurveonsurf des conges.
  
  Standard_Real x,y,distref;
  Fi1.PCurveOnSurf()->Value(UInt1).Coord(x,y);
  gp_Pnt p3d1 = HS1->Value(x,y);
  Fi2.PCurveOnSurf()->Value(UInt2).Coord(x,y);
  gp_Pnt p3d2 = HS2->Value(x,y);
  distref = p3d1.Distance(p3d2);
  P.SetXYZ(0.5*(p3d1.XYZ() + p3d2.XYZ()));
  // recalcul de l'extrema
  Standard_Real delt1 = 
    Min(0.1,0.05*(Fi1.LastParameter() - Fi1.FirstParameter()));
  Handle(Geom2dAdaptor_HCurve) hc2d1 = 
    new Geom2dAdaptor_HCurve(Fi1.PCurveOnSurf(),UInt1-delt1,UInt1+delt1);
  Adaptor3d_CurveOnSurface cons1(hc2d1,HS1);
  Standard_Real delt2 = 
    Min(0.1,0.05*(Fi2.LastParameter() - Fi2.FirstParameter()));
  Handle(Geom2dAdaptor_HCurve) hc2d2 = 
    new Geom2dAdaptor_HCurve(Fi2.PCurveOnSurf(),UInt2-delt2,UInt2+delt2);
  Adaptor3d_CurveOnSurface cons2(hc2d2,HS2);
  Extrema_LocateExtCC ext(cons1,cons2,UInt1,UInt2);
  if(ext.IsDone()) {
    Standard_Real dist = ext.Value();
    if(dist<distref) {
      Extrema_POnCurv ponc1,ponc2;
      ext.Point(ponc1,ponc2);
      UInt1 = ponc1.Parameter();
      UInt2 = ponc2.Parameter();
      gp_Pnt Pnt1 = ponc1.Value();
      gp_Pnt Pnt2 = ponc2.Value();
      P.SetXYZ(0.5*(Pnt1.XYZ() + Pnt2.XYZ()));
    }
  }
} 

//=======================================================================
//function : BoundSurf
//purpose  : computes a GeomAdaptor_Surface from the surface of the 
//           SurfData Fd1 and trims it to allow the intersection computation

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

Handle(GeomAdaptor_HSurface) ChFi3d_BoundSurf(TopOpeBRepDS_DataStructure&    DStr,
                                    const Handle(ChFiDS_SurfData)& Fd1,
                                    const Standard_Integer&        IFaCo1,
                                    const Standard_Integer&        IFaArc1)
{
  //rmq : comme en fait les 2 interferences de Fd1 ne servent qu'a donner les 
  //      bornes, les indices  IFaCo1 et IFaArc1 sont inutiles.
  //      On les garde ici en option au cas ou il faudrait borner de facon plus
  //      restrictive (avec des points d'intersection en argument en plus).
  
  Handle(GeomAdaptor_HSurface) HS1 = new GeomAdaptor_HSurface();
  GeomAdaptor_Surface& S1 = HS1->ChangeSurface();
  S1.Load(DStr.Surface(Fd1->Surf()).Surface());
  
  if ((IFaCo1 == 0)||(IFaArc1 == 0)) 
    return HS1;
  
  const ChFiDS_FaceInterference& FiCo1 = Fd1->Interference(IFaCo1);
  const ChFiDS_FaceInterference& FiArc1 = Fd1->Interference(IFaArc1);
  
  Standard_Real Du,Dv,mu,Mu,mv,Mv;
  gp_Pnt2d UVf1,UVf2,UVl1,UVl2;
  
  UVf1 = FiCo1.PCurveOnSurf()->Value(FiCo1.FirstParameter());
  UVl1 = FiCo1.PCurveOnSurf()->Value(FiCo1.LastParameter());
  UVf2 = FiArc1.PCurveOnSurf()->Value(FiArc1.FirstParameter());
  UVl2 = FiArc1.PCurveOnSurf()->Value(FiArc1.LastParameter());
  ChFi3d_Boite(UVf1,UVf2,UVl1,UVl2,Du,Dv,mu,Mu,mv,Mv);
  GeomAbs_SurfaceType styp = S1.GetType();
  if (styp == GeomAbs_Cylinder) { 
    Dv = Max(0.5*Dv,4.*S1.Cylinder().Radius());
    Du = 0.;
    S1.Load(DStr.Surface(Fd1->Surf()).Surface(),
          mu,Mu,mv-Dv,Mv+Dv);
  }
  //Dans le cas d'un tore ou cone, il ne faut pas que l'agrandissement des bounds engendrent une surface avec une periode plus grande que 2PI. lvt
  else if (styp == GeomAbs_Torus ||
         styp == GeomAbs_Cone) {
    Du = Min(PI-0.5*Du,0.1*Du);
    Dv = 0.;
    S1.Load(DStr.Surface(Fd1->Surf()).Surface(),
          mu-Du,Mu+Du,mv,Mv);
  }
  else if (styp == GeomAbs_Plane) {
    Du = Max(0.5*Du,4.*Dv);
    Dv = 0.;
    S1.Load(DStr.Surface(Fd1->Surf()).Surface(),
          mu-Du,Mu+Du,mv,Mv);
  }
  return HS1;
}



//=======================================================================
//function : SearchPivot
//purpose  : 
//=======================================================================
Standard_Integer ChFi3d_SearchPivot(Standard_Integer* s,
                            Standard_Real u[3][3],
                            const Standard_Real t)
{
  //Cette fonction recherche comme pivot une cd dont les sections 
  //ne se croisent pas sur la face opposee. 
  //         - il y aura peut etre des cas suffisamment asymetriques
  //           pour qu aucun des trois conges ne convienne!! A VOIR.
  //         - dans le cas ou plusieurs conviennent on prend le 
  //           premier qui n est pas forcement le meilleur, prevoir 
  //           d affiner cela en comparant les parametres sur les
  //           lignes guide et (/ou) les rayons.
  
  Standard_Boolean bondeb,bonfin;
  for(Standard_Integer i = 0; i <= 2; i++) {
    if(s[(i+1)%3] == 1) {bondeb = (u[(i+1)%3][i]-u[(i+1)%3][(i+2)%3] >= -t);}
    else {bondeb = (u[(i+1)%3][i]-u[(i+1)%3][(i+2)%3] <= t);}
    if(s[(i+2)%3] == 1) {bonfin = (u[(i+2)%3][i]-u[(i+2)%3][(i+1)%3] >= -t);}
    else {bonfin = (u[(i+2)%3][i]-u[(i+2)%3][(i+1)%3] <= t);}
    if (bondeb && bonfin) { return i; }
  }
  return -1;
}



//=======================================================================
//function : SearchFD
//purpose  : 
//=======================================================================
Standard_Boolean ChFi3d_SearchFD(TopOpeBRepDS_DataStructure& DStr,
                         const Handle(ChFiDS_Stripe)& cd1, 
                         const Handle(ChFiDS_Stripe)& cd2,
                         const Standard_Integer sens1,
                         const Standard_Integer sens2,
                         Standard_Integer& i1,
                         Standard_Integer& i2,
                         Standard_Real& p1,
                         Standard_Real& p2,
                         const Standard_Integer ind1,
                         const Standard_Integer ind2,
                         TopoDS_Face& face,
                         Standard_Boolean& sameside,
                         Standard_Integer& jf1,
                         Standard_Integer& jf2)
{
  Standard_Boolean found = Standard_False;
  Standard_Integer id1 = ind1, id2 = ind2;
  Standard_Integer if1 = ind1, if2 = ind2;
  Standard_Integer l1 = cd1->SetOfSurfData()->Length();
  Standard_Integer l2 = cd2->SetOfSurfData()->Length();
  Standard_Integer i;
  Standard_Boolean fini1 = Standard_False, fini2 = Standard_False;
  Standard_Boolean visavis,visavisok = Standard_False;
  TopoDS_Vertex Vtx;
  while( !found ) {
    for(i = id1; (i*sens1) <= (if1*sens1) && !found && !fini2; i = i+sens1 ) { 
      if(ChFi3d_IsInFront(DStr,cd1,cd2,i,if2,sens1,sens2,p1,p2,face,sameside,jf1,jf2,visavis,Vtx,Standard_False,0)) {
      i1 = i;
      i2 = if2;
      found = Standard_True;
      }
      else if (visavis && !visavisok) {
      visavisok = Standard_True;
      i1 = i;
      i2 = if2;
      }
    }
    if(!fini1) {
      if1 = if1 + sens1;
      if(if1 < 1 || if1 > l1) { if1 = if1 - sens1; fini1 = Standard_True; }
    }
    
    for(i = id2; (i*sens2) <= (if2*sens2) && !found && !fini1; i = i+sens2 ) { 
      if(ChFi3d_IsInFront(DStr,cd1,cd2,if1,i,sens1,sens2,p1,p2,face,sameside,jf1,jf2,visavis,Vtx,Standard_False,0)) {
      i1 = if1;
      i2 = i;
      found = Standard_True;
      }
      else if (visavis && !visavisok) {
      visavisok = Standard_True;
      i1 = if1;
      i2 = i;
      }
    }
    if(!fini2) {
      if2 = if2 + sens2;
      if(if2 < 1 || if2 > l2) { if2 = if2 - sens2; fini2 = Standard_True; }
    }
    if(fini1 && fini2) break;
  }
  return found;
}

//=======================================================================
//function : Parameters
//purpose  : compute the parameters <u> and <v> of the 3d point <p3d> 
//           on the surface <S> if it's an analytic surface
//=======================================================================

void ChFi3d_Parameters(const Handle(Geom_Surface)& S,
                   const gp_Pnt& p3d,
                   Standard_Real& u,
                   Standard_Real& v)
{
  GeomAdaptor_Surface gas(S);
  switch ( gas.GetType() ) {
  case GeomAbs_Plane :
    ElSLib::Parameters(gas.Plane(),p3d,u,v);
    break;
  case GeomAbs_Cylinder :
    ElSLib::Parameters(gas.Cylinder(),p3d,u,v);
    break;
  case GeomAbs_Cone :
    ElSLib::Parameters(gas.Cone(),p3d,u,v);
    break;
  case GeomAbs_Sphere :
    ElSLib::Parameters(gas.Sphere(),p3d,u,v);
    break;
  case GeomAbs_Torus :
    ElSLib::Parameters(gas.Torus(),p3d,u,v);
    break;
  case GeomAbs_BezierSurface :
  case GeomAbs_BSplineSurface :  
  default :
    {
      GeomAPI_ProjectPointOnSurf tool(p3d,S);
      if ( tool.NbPoints() != 1 )
      StdFail_NotDone::Raise("");
      else
      tool.Parameters(1,u,v);
    }
  }
}

//=======================================================================
//function : TrimCurve
//purpose  : trims the curve <gc> between the points <FirstP> and
//           <LastP>. The trimmed curve is <gtc>
//=======================================================================

void ChFi3d_TrimCurve(const Handle(Geom_Curve)& gc,
                  const gp_Pnt& FirstP,
                  const gp_Pnt& LastP,
                  Handle(Geom_TrimmedCurve)& gtc)
{
#ifndef DEB
  Standard_Real uf = 0.,ul = 0.;
#else
  Standard_Real uf,ul;
#endif
  GeomAdaptor_Curve gac(gc);
  switch ( gac.GetType() ) {
  case GeomAbs_Line :
    {
      uf = ElCLib::Parameter(gac.Line(),FirstP);
      ul = ElCLib::Parameter(gac.Line(),LastP);
    }
    break;
  case GeomAbs_Circle :
    {
      uf = ElCLib::Parameter(gac.Circle(),FirstP);
      ul = ElCLib::Parameter(gac.Circle(),LastP);
    }
    break;
  case GeomAbs_Ellipse :
    {
      uf = ElCLib::Parameter(gac.Ellipse(),FirstP);
      ul = ElCLib::Parameter(gac.Ellipse(),LastP);
    }
    break;
  case GeomAbs_Hyperbola :
    {
      uf = ElCLib::Parameter(gac.Hyperbola(),FirstP);
      ul = ElCLib::Parameter(gac.Hyperbola(),LastP);
    }
    break;
  case GeomAbs_Parabola :
    {
      uf = ElCLib::Parameter(gac.Parabola(),FirstP);
      ul = ElCLib::Parameter(gac.Parabola(),LastP);
    }
    break;
  default :
    {
    GeomAPI_ProjectPointOnCurve tool(FirstP,gc);
    if ( tool.NbPoints() != 1 )
      StdFail_NotDone::Raise("");
    else
      uf = tool.Parameter(1);
    tool.Init(LastP,gc);
    if ( tool.NbPoints() != 1 )
      StdFail_NotDone::Raise("");
    else
      ul = tool.Parameter(1);
  }
  }
  gtc = new Geom_TrimmedCurve(gc,uf,ul);
}

//=======================================================================
//function : PerformElSpine
//purpose  : 
//=======================================================================

static Standard_Boolean GoodExt(const Handle(Geom_Curve)& C,
                        const gp_Vec&             V,
                        const Standard_Real       f,
                        const Standard_Real       l,
                        const Standard_Real       a)
{
  for(Standard_Integer i = 0; i < 6; i++) {
    gp_Pnt d0; gp_Vec d1;
    Standard_Real t = i * 0.2;
    C->D1(((1-t)*f+t*l),d0,d1);
    Standard_Real ang = d1.Angle(V);
    Standard_Real angref = a*t + 0.002;
    if(ang > angref) return 0;
  }
  return 1;
}

Standard_EXPORT void ChFi3d_PerformElSpine(Handle(ChFiDS_HElSpine)& HES,
                                 Handle(ChFiDS_Spine)&    Spine,
                                 const GeomAbs_Shape      continuity,
                                 const Standard_Real      tol) 
{
  ChFiDS_ElSpine& ES = HES->ChangeCurve();
  
  Standard_Real WF = ES.FirstParameter();
  Standard_Real WL = ES.LastParameter();
  Standard_Real Wrefdeb = WF;
  Standard_Real Wreffin = WL;
  Standard_Integer IF,IL;
  Standard_Boolean periodic = Spine->IsPeriodic();
  Standard_Real nwf = WF, nwl = WL,period;
  Standard_Integer nbed = Spine->NbEdges();
  Standard_Real pared,tolpared;
  Handle(Geom_BSplineCurve) BSpline;
  Handle(Geom_BezierCurve) Bezier;
  
  if(periodic) {
    period = Spine->Period();
    nwf = ElCLib::InPeriod(WF,-tol,period-tol);
    IF = Spine->Index(nwf,1);
    nwl = ElCLib::InPeriod(WL,tol,period+tol);
    IL = Spine->Index(nwl,0);
    if(nwl<nwf+tol) IL += nbed;
  }
  else{
    IF = Spine->Index(WF,1);
    IL = Spine->Index(WL,0);
    Wrefdeb = Max(Spine->FirstParameter(IF),WF);
    Wreffin = Min(Spine->LastParameter(IL),WL);
  }
  
  gp_Pnt PDeb, PFin, Bout;
  gp_Vec VrefDeb, VrefFin;
  Spine->D1(WF,PDeb,VrefDeb);
  Spine->D1(WL,PFin,VrefFin);
  VrefDeb.Normalize();
  VrefFin.Normalize();
  
  Standard_Boolean   Bof;
  Standard_Integer IEdge;
  
  Standard_Real  First, Last, epsV;
  TColgp_Array1OfPnt ExtrapPole(1, 5);
  TColgp_Array1OfPnt ExtraCoeffs(1, 5);
  TColgp_Array1OfXYZ Cont(1,5);
  
  TopoDS_Edge E, Eold;
  TopoDS_Vertex V;
  
  // Attention on segmente eventuellement la premiere et la
  // derniere arete.
  
  // Traitment de la premiere arete
  E = Spine->Edges(IF);
  Bof = BRepLib::BuildCurve3d(E);
  const BRepAdaptor_Curve& edc = Spine->CurrentElementarySpine(IF);
  tolpared = edc.Resolution(tol);
  Handle(Geom_Curve) Cv = BRep_Tool::Curve(E, First, Last);
  Standard_Real urefdeb = Spine->FirstParameter(IF);
  Standard_Boolean checkdeb = (nwf > urefdeb),cepadur = 0;
  if(checkdeb) Spine->Parameter(IF,nwf,pared,0);
  if(E.Orientation() == TopAbs_REVERSED) {
    Standard_Real sov = First;
    First = Cv->ReversedParameter(Last);
    Last = Cv->ReversedParameter(sov);
    if(checkdeb) pared = Cv->ReversedParameter(pared);
    else pared = First;
    if(First < pared) First = pared; 
    if(IL == IF) {
      Standard_Real ureffin = Spine->LastParameter(IL);
      Standard_Boolean checkfin = (nwl < ureffin);
      if(checkfin) {
      Spine->Parameter(IL,nwl,pared,0);
      pared = Cv->ReversedParameter(pared);
      }
      else pared = Last;
      if(pared < Last) Last = pared; 
    }
    Cv = Cv->Reversed();
  }
  else{
    if(!checkdeb) pared = First;
    if(First < pared) First = pared; 
    if(IL == IF) {
      Standard_Real ureffin = Spine->LastParameter(IL);
      Standard_Boolean checkfin = (nwl < ureffin);
      if(checkfin) Spine->Parameter(IL,nwl,pared,0);
      else pared = Last;
      if(pared < Last) Last = pared; 
    }
  }

  if(Abs(Last-First) < tolpared) cepadur = 1;

  //Petite veru pour les cas ou un KPart a bouffe l arete
  //sans parvenir a terminer. On tire une droite.
  if(cepadur) {
    Handle(Geom_Line) L;
    gp_Pnt ptemp; gp_Vec vtemp;
    if(WL < Spine->FirstParameter(1) + tol) {
      ES.LastPointAndTgt(ptemp,vtemp);
      gp_Dir d(vtemp);
      gp_Pnt olin;
      olin.ChangeCoord().SetLinearForm(-WL,d.XYZ(),PFin.XYZ());
      L = new Geom_Line(olin,d);
      ES.SetCurve(L);
    }
    else if(WF > Spine->LastParameter(nbed) - tol) {
      ES.FirstPointAndTgt(ptemp,vtemp);
      gp_Dir d(vtemp);
      gp_Pnt olin;
      olin.ChangeCoord().SetLinearForm(-WF,d.XYZ(),PDeb.XYZ());
      L = new Geom_Line(olin,d);
      ES.SetCurve(L);
    }
#ifdef DRAW
    if (ChFi3d_GettraceDRAWSPINE()) {
//      POP pour NT
//      char name[100];
      char* name = new char[100];
      sprintf(name,"elspine_%d",NbSp++);
      Handle(Geom_Curve) bid = new Geom_TrimmedCurve(L,WF,WL);
      DrawTrSurf::Set(name,bid);
    }
#endif
    return;
  }
  Handle(Geom_BoundedCurve) TC = new (Geom_TrimmedCurve)(Cv, First, Last);
  Handle(Geom_BSplineCurve) BS;
  BS = GeomConvert::CurveToBSplineCurve(TC);
  CurveCleaner(BS, Abs(WL-WF)*1.e-4, 0);
  //Smoothing of the curve
  if (TC->Continuity() < GeomAbs_C2 && !ChFi3d_IsSmooth( TC ))
    {
      BS = ChFi3d_ApproxByC2( TC );
      TC = BS;
    }

  //  Concatenation des aretes suivantes
  GeomConvert_CompCurveToBSplineCurve Concat( TC, Convert_QuasiAngular );
  Eold = E;
  
  for (IEdge=IF+1; IEdge<=IL; IEdge++) {
    Standard_Integer iloc = IEdge;
    if(periodic) iloc = (IEdge - 1)%nbed + 1;
    E = Spine->Edges(iloc);
    if (! BRep_Tool::Degenerated(E)) {
      Bof = TopExp::CommonVertex(Eold, E, V);
      if (Bof) {epsV = BRep_Tool::Tolerance(V);}
      else {epsV = tol;}
      Bof = BRepLib::BuildCurve3d(E);
      if (!Bof) {Standard_ConstructionError::
               Raise("PerformElSpine : BuildCurve3d error");}
      Cv = BRep_Tool::Curve(E, First, Last);
      if(IEdge == IL) {
      Standard_Real ureffin = Spine->LastParameter(iloc);
      Standard_Boolean checkfin = (nwl < ureffin);
      if(checkfin) Spine->Parameter(iloc,nwl,pared,0);
      else pared = Last;
      if(E.Orientation() == TopAbs_REVERSED) {
        Standard_Real sov = First;
        First = Cv->ReversedParameter(Last);
        Last = Cv->ReversedParameter(sov);
        if(checkfin) pared = Cv->ReversedParameter(pared);
        else pared = Last;
        Cv = Cv->Reversed();
      }
      if(pared < Last) Last = pared; 
      }
      TC = new (Geom_TrimmedCurve)(Cv, First, Last);
      BS = GeomConvert::CurveToBSplineCurve(TC);
      CurveCleaner(BS, Abs(WL-WF)*1.e-4, 0);
      //Smoothing of the curve
      if (TC->Continuity() < GeomAbs_C2 && !ChFi3d_IsSmooth( TC ))
      {
        BS = ChFi3d_ApproxByC2( TC );
        TC = BS;
      }
      Standard_Real tolrac = Min(tol,epsV);
      Bof = Concat.Add( TC, 2*tolrac, Standard_True );
      // si l'ajout ne s'est pas bien passe on essai d'augmenter la tolerance
      if (!Bof) {
      Bof = Concat.Add( TC, 2*epsV, Standard_True );
      }
      if (!Bof) {
      Bof = Concat.Add( TC, 200*epsV, Standard_True );
#ifdef DEB
      cout << "Tolerance de chainage non respecte" << endl;
#endif
      if (!Bof) {Standard_ConstructionError::
                 Raise("PerformElSpine: spine merged error");}
      }
    }
    Eold = E;
  }
  
  // On a la portion d elspine calculee sans prolongements sur la partie 
  // valide des aretes du chemin.
  BSpline = Concat.BSplineCurve();
  // On reparametre ici pour coller au mieux a l abscisse des aretes.
  TColStd_Array1OfReal BSNoeuds (1, BSpline->NbKnots());
  BSpline->Knots(BSNoeuds);
  BSplCLib::Reparametrize (Wrefdeb, Wreffin, BSNoeuds);
  BSpline->SetKnots(BSNoeuds);

  // Traitement des Extremites
  Standard_Real LocalWL = WL;
  Standard_Real LocalWF = WF; 
  GeomAdaptor_Curve gacurve;
  Standard_Integer caredeb = 0, carefin = 0;
  Handle(Geom_BSplineCurve) newc;
  Standard_Real Angle = PI*0.75;
  if (!ES.IsPeriodic() && !PDeb.IsEqual(BSpline->Pole(1), tol) ) {
    // Prolongement C3 au debut
    // afin d'eviter des pts d'inflexions dans la partie utile de la
    // spine le prolongement se fait jusqu'a un point eloigne.
    if(BSpline->IsRational()) caredeb = 1;
    Standard_Real rabdist = Wrefdeb - WF;
    Bout = PDeb.Translated(-20*rabdist * VrefDeb);
    Standard_Boolean goodext = 0;
    for(Standard_Integer icont = 3; icont>=1 && !goodext; icont--) {
      newc = BSpline;
      GeomLib::ExtendCurveToPoint( newc, Bout, icont, Standard_False);
      gacurve.Load(newc);
      GCPnts_AbscissaPoint GCP(gacurve,-rabdist,Wrefdeb,WF);
      if(GCP.IsDone()) {
      WF = GCP.Parameter();
      goodext = GoodExt(newc,VrefDeb,Wrefdeb,WF,Angle);
      }
    }
    if(caredeb) caredeb = newc->NbKnots() - BSpline->NbKnots();
    BSpline = newc;
    LocalWF = BSpline->FirstParameter();
  }
  
  if (!ES.IsPeriodic() && !PFin.IsEqual(BSpline->Pole(BSpline->NbPoles()), tol) ) {
    // Prolongement C3 en fin
    if(BSpline->IsRational()) carefin = 1;
    Standard_Real rabdist = WL - Wreffin;
    Bout = PFin.Translated(20*rabdist * VrefFin);
    Standard_Boolean goodext = 0;
    for(Standard_Integer icont = 3; icont>=1 && !goodext; icont--) {
      newc = BSpline;
      GeomLib::ExtendCurveToPoint( newc, Bout, icont, Standard_True);
      gacurve.Load(newc);
      GCPnts_AbscissaPoint GCP(gacurve,rabdist,Wreffin,WL);
      if(GCP.IsDone()) {
      WL = GCP.Parameter();
      goodext = GoodExt(newc, VrefFin, Wreffin,WL,Angle);
      }
    }
    if(carefin) carefin = newc->NbKnots() - BSpline->NbKnots();
    BSpline = newc;
    LocalWL = BSpline->LastParameter(); 
  }
  
  //Reparametrisation et segmentation sur le domaine de la Spine.
  if(Abs(BSpline->FirstParameter() - WF)<tol) WF = BSpline->FirstParameter();
  if(Abs(BSpline->LastParameter() - WL)<tol) WL = BSpline->LastParameter();
  if ( (LocalWF<WF) ||  (LocalWL>WL)) {   // pour eviter des pb avec segment!
    BSpline->Segment(WF, WL);
    ES.FirstParameter(WF);
    ES.LastParameter(WL);
  }
  if (BSpline->IsRational()) {
    Handle(Geom_BSplineCurve) C1;
    C1 =  Handle(Geom_BSplineCurve)::DownCast(BSpline->Copy());
    GeomConvert::C0BSplineToC1BSplineCurve(C1, tol, 0.1);
    // Il faut s'assurer que l'origine n'a pas bouge (cts21158)
    if (C1->FirstParameter() == BSpline->FirstParameter()) BSpline = C1;
    else {
#if DEB
      cout << "Attention : Echec de C0BSplineToC1 !" << endl;
#endif
    }
  }
  
  // Deformation eventuelle pour rendre la spine C2.
  // ou C3 pour des approx C2
  if((caredeb || carefin) && BSpline->Degree() < 8) BSpline->IncreaseDegree(8);
  Standard_Integer fk = 2, lk = BSpline->NbKnots()-1;
  if(BSpline->IsPeriodic()) fk = 1;
  if(caredeb) fk += caredeb;
  if(carefin) lk -= carefin;
  Standard_Integer   MultMax, ii; 
  if (continuity == GeomAbs_C3) {
    if (BSpline->Degree() < 7) BSpline->IncreaseDegree(7);
    MultMax = BSpline->Degree() - 3;
  }
  else {
    if (BSpline->Degree() < 5) BSpline->IncreaseDegree(5);
    MultMax = BSpline->Degree() - 2;
  }
  // correction C2 ou C3 (si possible)
  CurveCleaner(BSpline, Abs(WL-WF)*1.e-4, 1);
  CurveCleaner(BSpline, Abs(WL-WF)*1.e-2, MultMax);
  Standard_Integer MultMin = Max(BSpline->Degree() - 4, 1);
  for (ii = fk; ii <= lk; ii++) {
    if( BSpline->Multiplicity(ii) > MultMax ) {
      Bof = BSpline->RemoveKnot(ii, MultMax, Abs(WL-WF)/10);
    }
    // Voir C4
    if( BSpline->Multiplicity(ii) > MultMin ) {
      Bof = BSpline->RemoveKnot(ii, MultMin, Abs(WL-WF)*1.e-4);
    }       
  }
  // elspine periodique => BSpline Periodique
  if(ES.IsPeriodic()) {
    if(!BSpline->IsPeriodic()) {
      BSpline->SetPeriodic();
      Bof = BSpline->RemoveKnot(1, MultMax, Abs(WL-WF)/10);
    }
  }
  else { 
    // Sinon faut il bouger les poles pour les adapter 
    // aux nouvelles tangentes ?
    Standard_Boolean adjust = Standard_False; 
    gp_Pnt P1, P2;
    gp_Vec V1, V2;
    BSpline->D1(WF, P1, V1);
    V1.Normalize();
    ES.FirstPointAndTgt(PDeb,VrefDeb);
    Standard_Real scaldeb = VrefDeb.Dot(V1);
    Standard_Real disdeb = PDeb.Distance(P1);
    if((Abs(WF-LocalWF) < 1.e-12) &&
       ((scaldeb <= 0.9999999) ||
      disdeb >= tol)) {   
      // Oui s'il n'y as pas eu de prolongement et que la tangente n'est pas
      // la bonne.
      adjust = Standard_True;
    } 
    BSpline->D1(WL, P2, V2);
    V2.Normalize();
    ES.LastPointAndTgt(PFin,VrefFin);
    Standard_Real scalfin = VrefFin.Dot(V2); 
    Standard_Real disfin = PFin.Distance(P2);
    if((Abs(WL-LocalWL) < 1.e-12) && 
       ((scalfin <= 0.9999999)||
      disfin >= tol)) {
      // de meme a la fin
      adjust = Standard_True;
    }
    if(adjust) GeomLib::AdjustExtremity(BSpline, PDeb, PFin, VrefDeb, VrefFin);
  }

  // Le Resultat       
#if DEB
  if (! BSpline->IsCN(2)) cout << "Attention Spine non C2 !" << endl;
#endif

  ES.SetCurve(BSpline);
#ifdef DRAW
  if (ChFi3d_GettraceDRAWSPINE()) {
    char name[100];
    sprintf(name,"elspine_%d",NbSp++);
    char* Temp = name;
    DrawTrSurf::Set(Temp,BSpline);
//    DrawTrSurf::Set(name,BSpline);
  }
#endif
}

//=======================================================================
//function : cherche_face1
//purpose  : cherche la face F differente de F1 dans la map.
// La map  contient  les deux faces adjacentes a une edge 
//=======================================================================
void ChFi3d_cherche_face1 (const TopTools_ListOfShape & map,
                     const TopoDS_Face & F1,
                     TopoDS_Face &  F)   
{ 
  TopoDS_Face Fcur;
  Standard_Boolean trouve=Standard_False;
  TopTools_ListIteratorOfListOfShape It;
  for (It.Initialize(map);It.More()&&!trouve;It.Next()) { 
    Fcur=TopoDS::Face (It.Value());
    if (!Fcur.IsSame(F1)) {
      F=Fcur;trouve=Standard_True;}
  }
} 
//=======================================================================
//function : cherche_element
//purpose  : cherche l'edge E  de F1 differente de E1 et contenant le vertex V
// Vtx est l'autre vertex de E 
//=======================================================================

void ChFi3d_cherche_element(const TopoDS_Vertex & V,
                      const TopoDS_Edge & E1,
                      const TopoDS_Face & F1,
                      TopoDS_Edge & E , 
                      TopoDS_Vertex  & Vtx )
{ 
  Standard_Integer ie; 
  TopoDS_Vertex V1,V2;
  Standard_Boolean trouve=Standard_False;
  TopoDS_Edge Ecur;
  TopTools_IndexedMapOfShape  MapE;
  TopExp::MapShapes( F1,TopAbs_EDGE,MapE);
  for ( ie=1; ie<= MapE.Extent()&& !trouve; ie++) {  
    Ecur = TopoDS::Edge (MapE(ie));
    if (!Ecur.IsSame(E1)) { 
      TopTools_IndexedMapOfShape  MapV;
      TopExp::MapShapes(Ecur, TopAbs_VERTEX, MapV);
      if (MapV.Extent()==2) {
      V1 = TopoDS::Vertex (MapV(1));
      V2 = TopoDS::Vertex (MapV(2));
      if (V1.IsSame(V)) { 
        Vtx=V2;
        E=Ecur;
        trouve=Standard_True;
      }
      else if (V2.IsSame(V)) {
        Vtx=V1;
        E=Ecur;
        trouve=Standard_True;
      }
      }
    }
  }
} 
//=======================================================================
//function : cherche_edge
//purpose  : cherche l'edge E  de F1 differente de la liste d'edges E1 et
//            contenant le vertex V  Vtx est l'autre vertex de E 
//=======================================================================
 
void ChFi3d_cherche_edge(const TopoDS_Vertex & V,
                      const  TopTools_Array1OfShape & E1,
                      const TopoDS_Face & F1,
                      TopoDS_Edge & E , 
                      TopoDS_Vertex  & Vtx )
{ 
  Standard_Integer ie,i; 
  TopoDS_Vertex V1,V2;
  Standard_Boolean trouve=Standard_False;
  TopoDS_Edge Ecur;
  Standard_Boolean same;
  TopTools_IndexedMapOfShape  MapE;
  TopExp::MapShapes( F1,TopAbs_EDGE,MapE);
  for ( ie=1; ie<= MapE.Extent()&& !trouve; ie++) {  
    Ecur=TopoDS::Edge (MapE(ie));
    same=Standard_False;
    for (i=E1.Lower();i<=E1.Upper() ;i++) {
       if (Ecur.IsSame(E1.Value(i))) same=Standard_True;
    }
      if (!same) {
      TopTools_IndexedMapOfShape  MapV;
      TopExp::MapShapes(Ecur, TopAbs_VERTEX, MapV);
      if (MapV.Extent()==2) {
      V1 = TopoDS::Vertex (MapV(1));
      V2 = TopoDS::Vertex (MapV(2));
      if (V1.IsSame(V)) { 
        Vtx=V2;
        E=Ecur;
        trouve=Standard_True;
      }
      else if (V2.IsSame(V)) {
        Vtx=V1;
        E=Ecur;
        trouve=Standard_True;
      }
      }
    }
  }
}

//=======================================================================
//function : nbface
//purpose  : calcule le nombre de faces communes a un vertex
//           
//=======================================================================
Standard_Integer  ChFi3d_nbface (const TopTools_ListOfShape & mapVF )
{ Standard_Integer nface=0;     
  TopTools_ListIteratorOfListOfShape ItF,JtF;
  Standard_Integer  fj = 0;
  for (ItF.Initialize(mapVF); ItF.More(); ItF.Next()) {
    fj++;
    Standard_Integer kf = 1;
    const TopoDS_Shape& cur = ItF.Value();
    for (JtF.Initialize(mapVF); JtF.More( )&&(kf<fj); JtF.Next(), kf++) {
      if(cur.IsSame(JtF.Value())) break;
    }
    if(kf == fj) nface++;
  }
  return nface;  
}

//=======================================================================
//function : edge_common_faces 
//purpose  :  determine les deux faces partageant une edge.
//            F1 =F2 si on a une arete de couure 
//=======================================================================

void ChFi3d_edge_common_faces (const TopTools_ListOfShape & mapEF,
                               TopoDS_Face & F1,
                               TopoDS_Face &  F2)   
{ TopTools_ListIteratorOfListOfShape It;
  TopoDS_Face F;
  Standard_Boolean trouve;
  It.Initialize(mapEF);  
  F1=TopoDS::Face(It.Value());
  trouve=Standard_False;
  for(It.Initialize(mapEF);It.More()&&!trouve;It.Next()) { 
     F=TopoDS::Face (It.Value());
     if (!F.IsSame(F1)) {
       F2=F;trouve=Standard_True;
     }
  }
  if (!trouve) F2=F1;
}

/***********************************************************/
// donne l'angle entre les edges E1 et E2 . Vtx est le vertex
// commun aux edges
/************************************************************/
Standard_Real ChFi3d_AngleEdge (const TopoDS_Vertex & Vtx,
                                  const TopoDS_Edge&  E1,
                                  const TopoDS_Edge &  E2)
{   Standard_Real angle;
    BRepAdaptor_Curve BCurv1(E1);
    BRepAdaptor_Curve BCurv2(E2);
    Standard_Real parE1,parE2;
    gp_Vec dir1,dir2 ;
    gp_Pnt P1,P2 ;
    parE1=BRep_Tool::Parameter(Vtx,E1);
    parE2=BRep_Tool::Parameter(Vtx,E2);
    BCurv1.D1(parE1,P1,dir1);
    BCurv2.D1(parE2,P2,dir2);
    if (!Vtx.IsSame(TopExp::FirstVertex(E1))) dir1.Reverse();
    if (!Vtx.IsSame(TopExp::FirstVertex(E2)))  dir2.Reverse();
    angle=Abs(dir1.Angle(dir2));
    return angle;
}

//==================================================================
// ChercheBordsLibres
// determine si le vertex V1 a des aretes de bords libres
// edgelibre1 et edgelibre2 .
// On suppose qu'un sommet ne peut avoir que 2 aretes de bords libres
//===================================================================
void ChFi3d_ChercheBordsLibres(const  ChFiDS_Map & myVEMap,
                         const TopoDS_Vertex & V1,
                         Standard_Boolean & bordlibre,
                         TopoDS_Edge & edgelibre1,
                         TopoDS_Edge & edgelibre2) 
{ 
  bordlibre=Standard_False;
  TopTools_ListIteratorOfListOfShape ItE,ItE1;
  Standard_Integer nboccur;
  for (ItE.Initialize(myVEMap(V1)); ItE.More()&&!bordlibre; ItE.Next()) {
    nboccur=0;
    const TopoDS_Edge& cur = TopoDS::Edge(ItE.Value());
    if (!BRep_Tool::Degenerated(cur)) {
      for (ItE1.Initialize(myVEMap(V1)); ItE1.More(); ItE1.Next()) {
        const TopoDS_Edge& cur1 = TopoDS::Edge(ItE1.Value());
        if (cur1.IsSame(cur)) nboccur++;
      }
    }
    if (nboccur==1) {
      edgelibre1=cur;
      bordlibre=Standard_True;
    }
  }
  if (bordlibre) {
    bordlibre=Standard_False;
    for (ItE.Initialize(myVEMap(V1)); ItE.More()&&!bordlibre; ItE.Next()) {
      nboccur=0;
      const TopoDS_Edge& cur = TopoDS::Edge(ItE.Value());
      if (!BRep_Tool::Degenerated(cur)&&!cur.IsSame(edgelibre1)) {
      for (ItE1.Initialize(myVEMap(V1)); ItE1.More(); ItE1.Next()) {
        const TopoDS_Edge& cur1 = TopoDS::Edge(ItE1.Value());
        if (cur1.IsSame(cur)) nboccur++;
      }
      }
      if (nboccur==1) {
      edgelibre2=cur;
      bordlibre=Standard_True;
      }
    }
  }
}

//=======================================================================
//function : NbNotDegeneratedEdges
//purpose  : calcule le nb d'aretes non degenerees de la Map VEMap(Vtx)
// Attention les aretes de jointures sont comptees deux fois
//=======================================================================
Standard_Integer ChFi3d_NbNotDegeneratedEdges (const TopoDS_Vertex& Vtx,
                              const ChFiDS_Map& VEMap)
{
  TopTools_ListIteratorOfListOfShape ItE;
  Standard_Integer nba=VEMap(Vtx).Extent();
  for (ItE.Initialize(VEMap(Vtx)); ItE.More(); ItE.Next()) {
    const TopoDS_Edge& cur = TopoDS::Edge(ItE.Value());
    if (BRep_Tool::Degenerated(cur)) nba--;
  }
  return nba;
}

//=======================================================================
//function : NumberOfEdges
//purpose  : calcule le nombre d'aretes arrivant au sommet Vtx
// les aretes degenerees ne sont pas comptees. 
//=======================================================================
Standard_Integer ChFi3d_NumberOfEdges(const TopoDS_Vertex& Vtx,
                              const ChFiDS_Map& VEMap)
{
  Standard_Integer nba;
  Standard_Boolean bordlibre;
  TopoDS_Edge edgelibre1,edgelibre2;
  nba=ChFi3d_NbNotDegeneratedEdges(Vtx, VEMap);
  ChFi3d_ChercheBordsLibres(VEMap,Vtx,bordlibre,edgelibre1,edgelibre2);
  if (bordlibre) nba=(nba-2)/2 +2;
  else  nba=nba/2;
  return nba;
}
//=====================================================
// function cherche_vertex
// cherche le vertex commun entre deux edges 
//=====================================================

void ChFi3d_cherche_vertex (const TopoDS_Edge & E1,
                      const TopoDS_Edge & E2,
                      TopoDS_Vertex & vertex,
                      Standard_Boolean & trouve)
{ Standard_Integer i,j; 
  TopoDS_Vertex Vcur1,Vcur2;
  trouve=Standard_False;
  TopTools_IndexedMapOfShape  MapV1,MapV2;
  TopExp::MapShapes( E1,TopAbs_VERTEX,MapV1);
  TopExp::MapShapes( E2,TopAbs_VERTEX,MapV2);
  for ( i=1; i<= MapV1.Extent()&&!trouve; i++) {
    TopoDS_Shape alocalshape = TopoDS_Shape (MapV1(i));
    Vcur1=TopoDS::Vertex(alocalshape);
//    Vcur1=TopoDS::Vertex(TopoDS_Shape (MapV1(i)));
    for ( j=1; j<= MapV2.Extent()&&!trouve; j++) {
      TopoDS_Shape aLocalShape = TopoDS_Shape (MapV2(j));
      Vcur2=TopoDS::Vertex(aLocalShape);
//      Vcur2=TopoDS::Vertex(TopoDS_Shape (MapV2(j)));
      if (Vcur2.IsSame(Vcur1)) {
      vertex=Vcur1;trouve=Standard_True;
      }
    }
  }
}     
//============================================================
// determine si F a une arete de couture
//=============================================================
void ChFi3d_Couture( const TopoDS_Face & F,
                     Standard_Boolean & couture,
                     TopoDS_Edge & edgecouture)
{   TopoDS_Edge Ecur;
    couture=Standard_False;
    TopTools_IndexedMapOfShape  MapE1;
    TopExp::MapShapes( F,TopAbs_EDGE,MapE1);
    TopLoc_Location Loc;
    Handle(Geom_Surface) Surf =BRep_Tool::Surface(F,Loc);
    for ( Standard_Integer i=1; i<= MapE1.Extent()&&!couture; i++) {
      TopoDS_Shape aLocalShape = TopoDS_Shape (MapE1(i));
      Ecur=TopoDS::Edge(aLocalShape);
//      Ecur=TopoDS::Edge(TopoDS_Shape (MapE1(i)));
      if (BRep_Tool::IsClosed(Ecur,Surf,Loc)) {
      couture=Standard_True;
      edgecouture=Ecur;
      } 
    }
}

void ChFi3d_CoutureOnVertex( const TopoDS_Face & F,
                       const TopoDS_Vertex & V,
                       Standard_Boolean & couture,
                       TopoDS_Edge & edgecouture)
{   TopoDS_Edge Ecur;
    couture = Standard_False;
    TopTools_IndexedMapOfShape  MapE1;
    TopExp::MapShapes( F,TopAbs_EDGE,MapE1);
    TopLoc_Location Loc;
    Handle(Geom_Surface) Surf = BRep_Tool::Surface(F,Loc);
    for ( Standard_Integer i=1; i <= MapE1.Extent(); i++) {
      TopoDS_Shape aLocalShape = TopoDS_Shape (MapE1(i));
      Ecur=TopoDS::Edge(aLocalShape);
//      Ecur=TopoDS::Edge(TopoDS_Shape (MapE1(i)));
      if (BRep_Tool::IsClosed(Ecur,Surf,Loc)) {
      TopoDS_Vertex Vf, Vl;
      TopExp::Vertices( Ecur, Vf, Vl );
      if (Vf.IsSame(V) || Vl.IsSame(V))
        {
          couture = Standard_True;
          edgecouture = Ecur;
          break;
        }
      } 
    }
}

Standard_Boolean ChFi3d_IsPseudoSeam( const TopoDS_Edge& E,
                              const TopoDS_Face& F )
{
  if (! BRep_Tool::IsClosed( E, F ))
    return Standard_False;

  Standard_Boolean NeighborSeamFound = Standard_False;
  TopoDS_Vertex Vf, Vl, V1, V2;
  TopExp::Vertices( E, Vf, Vl );
  TopExp_Explorer Explo( F, TopAbs_EDGE );
  for (; Explo.More(); Explo.Next())
    {
      TopoDS_Edge Ecur = TopoDS::Edge( Explo.Current() );
      if (! Ecur.IsSame(E))
      {
        TopExp::Vertices( Ecur, V1, V2 );
        if ((V1.IsSame(Vf) || V1.IsSame(Vl) || V2.IsSame(Vf) || V2.IsSame(Vl)) &&
            BRepTools::IsReallyClosed( Ecur, F ))
          {
            NeighborSeamFound = Standard_True;
            break;
          }
      }
    }
  return NeighborSeamFound;
}

Handle(Geom_BSplineCurve) ChFi3d_ApproxByC2( const Handle(Geom_Curve)& C )
{
  Standard_Real First = C->FirstParameter(), Last = C->LastParameter();
  Standard_Integer NbPoints = 101;

  TColgp_Array1OfPnt Points( 1, NbPoints );
  Standard_Real delta = (Last - First) / (NbPoints-1);
  for (Standard_Integer i = 1; i <= NbPoints-1; i++)
    Points(i) = C->Value(First + (i-1)*delta);
  Points(NbPoints) = C->Value(Last);

  GeomAPI_PointsToBSpline Approx( Points );
  Handle(Geom_BSplineCurve) BS = Approx.Curve();
  return BS;
}

Standard_Boolean ChFi3d_IsSmooth( const Handle(Geom_Curve)& C )
{
  GeomAdaptor_Curve GAC( C );

  Standard_Integer ii;
  Standard_Integer intrv, nbintv = GAC.NbIntervals(GeomAbs_CN);
  TColStd_Array1OfReal TI(1,nbintv+1);
  GAC.Intervals(TI,GeomAbs_CN);
  Standard_Real Resolution = gp::Resolution(), Curvature;
  GeomLProp_CLProps LProp(C, 2, Resolution);
  gp_Pnt P1, P2;
  Standard_Integer Discretisation = 30;
  
//   Standard_Real t = TI(1);
//   LProp.SetParameter(t);
//   if (!LProp.IsTangentDefined())
//     return Standard_False;
//   C->D0(t, P1);
//   LProp.CentreOfCurvature(P2);
//   gp_Vec PrevVec(P1, P2);

  gp_Vec PrevVec;
  Standard_Boolean prevVecFound = Standard_False;
  Standard_Integer intrvFound = 0;
  for( intrv = 1; intrv <= nbintv; intrv++) {
    Standard_Real t = TI(intrv);
    Standard_Real step = (TI(intrv+1) - t) / Discretisation;
    for (ii = 1; ii <= Discretisation; ii++) {
      LProp.SetParameter(t);
      if (!LProp.IsTangentDefined())
      return Standard_False;
      Curvature = Abs(LProp.Curvature());
      if (Curvature > Resolution) {
      C->D0(t, P1);
      LProp.CentreOfCurvature(P2);
      PrevVec = gp_Vec(P1, P2);
      prevVecFound = Standard_True;
      break;
      }
      t += step;
    }
    if( prevVecFound ) {
      intrvFound = intrv;
      break;
    }
  }

  if( !prevVecFound )
    return Standard_True;

  //for (intrv = 1; intrv <= nbintv; intrv++) {
  for (intrv = intrvFound; intrv <= nbintv; intrv++) {
    Standard_Real t = TI(intrv);
    Standard_Real step = (TI(intrv+1) - t) / Discretisation;
    for (ii = 1; ii <= Discretisation; ii++)
      {      
      LProp.SetParameter(t);
      if (!LProp.IsTangentDefined())
        return Standard_False;
      Curvature = Abs(LProp.Curvature());
      if (Curvature > Resolution)
        {
          C->D0(t, P1);
          LProp.CentreOfCurvature(P2);
          gp_Vec Vec(P1, P2);
          Standard_Real Angle = PrevVec.Angle( Vec );
          if (Angle > PI/3.)
            return Standard_False;
          Standard_Real Ratio = Vec.Magnitude() / PrevVec.Magnitude();
          if (Ratio < 1.)
            Ratio = 1. / Ratio;
          if (Ratio > 2. && (intrv != nbintv || ii != Discretisation))
            return Standard_False;
          PrevVec = Vec;
        }
      t += step;
      }
  }

  return Standard_True;
}

Generated by  Doxygen 1.6.0   Back to index