[occt.git] / src / BRepFill / BRepFill_TrimSurfaceTool.cxx

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


#include <GeomAdaptor_SurfaceOfRevolution.hxx>
#include <AppParCurves_MultiCurve.hxx>
#include <BRep_Tool.hxx>
#include <BRepFill_ApproxSeewing.hxx>
#include <BRepFill_ComputeCLine.hxx>
#include <BRepFill_MultiLine.hxx>
#include <BRepFill_TrimSurfaceTool.hxx>
#include <BRepIntCurveSurface_Inter.hxx>
#include <BSplCLib.hxx>
#include <ElCLib.hxx>
#include <Geom2d_BSplineCurve.hxx>
#include <Geom2d_Curve.hxx>
#include <Geom2d_TrimmedCurve.hxx>
#include <Geom2dAdaptor_Curve.hxx>
#include <Geom2dAPI_ProjectPointOnCurve.hxx>
#include <Geom2dInt_GInter.hxx>
#include <Geom_BSplineCurve.hxx>
#include <Geom_Curve.hxx>
#include <Geom_Plane.hxx>
#include <Geom_SurfaceOfRevolution.hxx>
#include <Geom_TrimmedCurve.hxx>
#include <GeomAdaptor_HCurve.hxx>
#include <GeomAdaptor_Surface.hxx>
#include <GeomProjLib.hxx>
#include <gp.hxx>
#include <gp_Cone.hxx>
#include <gp_Cylinder.hxx>
#include <gp_Pnt2d.hxx>
#include <gp_Sphere.hxx>
#include <gp_Torus.hxx>
#include <IntRes2d_IntersectionPoint.hxx>
#include <IntRes2d_IntersectionSegment.hxx>
#include <PLib.hxx>
#include <Precision.hxx>
#include <Standard_NoSuchObject.hxx>
#include <Standard_NotImplemented.hxx>
#include <StdFail_NotDone.hxx>
#include <TColgp_Array1OfPnt.hxx>
#include <TColStd_Array1OfInteger.hxx>
#include <TColStd_Array1OfReal.hxx>
#include <TopAbs_State.hxx>
#include <TopExp.hxx>
#include <TopExp_Explorer.hxx>
#include <TopoDS.hxx>
#include <TopoDS_Edge.hxx>
#include <TopoDS_Face.hxx>
#include <TopoDS_Vertex.hxx>

//#define DRAW
#include <stdio.h>
#ifdef DRAW
#include <DrawTrSurf.hxx>
#include <DBRep.hxx>
#endif
#ifdef OCCT_DEBUG
static Standard_Boolean Affich       = Standard_False;
static Standard_Integer NBCALL  = 0;
#endif

//=======================================================================
//function : BRepFill_TrimSurfaceTool
//purpose  : Initialisation with two neighbor faces
//          Edge1 and Edge2 are parallel edges corresponding
//          to minimum iso on F1 and F2 respectively.
//          ie Edge1 is Umin or VMin on F1.
//          Inv1 and Inv2 show if Edge1 and Edge2 are
//          returned parallel.   
//=======================================================================

BRepFill_TrimSurfaceTool::BRepFill_TrimSurfaceTool
  (const Handle(Geom2d_Curve)& Bis, 
  const TopoDS_Face&          Face1, 
  const TopoDS_Face&          Face2,
  const TopoDS_Edge&          Edge1,
  const TopoDS_Edge&          Edge2,
  const Standard_Boolean      Inv1,
  const Standard_Boolean      Inv2 ) :
myFace1(Face1),
  myFace2(Face2),
  myEdge1(Edge1),
  myEdge2(Edge2),
  myInv1(Inv1),
  myInv2(Inv2),
  myBis  (Bis)
{
#ifdef OCCT_DEBUG
  if ( Affich) {
    NBCALL++;
    cout << " ---------->TrimSurfaceTool : NBCALL = " << NBCALL << endl;
#ifdef DRAW
    char name[256];

    sprintf(name,"FACE1_%d",NBCALL);
    DBRep::Set(name,myFace1);

    sprintf(name,"FACE2_%d",NBCALL);
    DBRep::Set(name,myFace2);

    sprintf(name,"EDGE1_%d",NBCALL);
    DBRep::Set(name,myEdge1);

    sprintf(name,"EDGE2_%d",NBCALL);
    DBRep::Set(name,myEdge2);

    sprintf(name,"BISSEC_%d",NBCALL);
    DrawTrSurf::Set(name,myBis);
#endif    
  }
#endif
}


//=======================================================================
//function : Bubble
//purpose  : Order the sequence of points by increasing x. 
//=======================================================================

static void Bubble(TColgp_SequenceOfPnt& Seq) 
{
  Standard_Boolean Invert = Standard_True;
  Standard_Integer NbPoints = Seq.Length();
  while (Invert) {
    Invert = Standard_False;
    for ( Standard_Integer i = 1; i < NbPoints; i++) {
      gp_Pnt P1 = Seq.Value(i);
      gp_Pnt P2 = Seq.Value(i+1);
      if (P2.X()<P1.X())  {
        Seq.Exchange(i,i+1);
        Invert = Standard_True;
      }
    }
  }
}


//=======================================================================
//function : EvalPhase
//purpose  : 
//=======================================================================

static Standard_Real EvalPhase(const TopoDS_Edge& Edge,
  const TopoDS_Face& Face,
  const GeomAdaptor_Surface& GAS,
  const gp_Ax3& Axis)
{
  gp_Pnt2d PE1,PE2,PF1,PF2;
  Standard_Real VDeg;
  Standard_Real V = 0.;
  BRep_Tool::UVPoints(Edge,Face,PE1,PE2);
  VDeg = PE1.Y();
  TopExp_Explorer Exp(Face,TopAbs_EDGE);
  for (; Exp.More(); Exp.Next()) {
    if ( !TopoDS::Edge(Exp.Current()).IsSame(Edge)) {
      BRep_Tool::UVPoints(TopoDS::Edge(Exp.Current()),Face,PF1,PF2);
      V = ( Abs(PF1.Y() - VDeg) > Abs(PF2.Y() - VDeg)) ? PF1.Y() : PF2.Y();
      break;
    }
  }
  gp_Pnt P = GAS.Value(0., V);

  if ( gp_Vec(Axis.Location(), P).Dot(Axis.XDirection()) < 0.) 
    return M_PI;
  else
    return 0.;
}


//=======================================================================
//function : EvalParameters
//purpose  : 
//=======================================================================

static void EvalParameters(const TopoDS_Edge&          Edge,
  const TopoDS_Face&          Face,
  const Handle(Geom2d_Curve)& Bis ,
  TColgp_SequenceOfPnt& Seq  ) 
{
  Standard_Boolean Degener = BRep_Tool::Degenerated(Edge);
  // return curves 3d associated to edges.
  TopLoc_Location L;
  Standard_Real   f,l;

  Handle(Geom_TrimmedCurve) CT;
  Handle(Geom_Plane) Plane = new Geom_Plane(0,0,1,0);

  Geom2dInt_GInter Intersector;

  Standard_Integer NbPoints, NbSegments;
  Standard_Real U1, U2;
  gp_Pnt  P;//,PSeq;

  //  Standard_Real Tol = Precision::Intersection();
  //  modified by NIZHNY-EAP Wed Dec 22 15:00:51 1999 ___BEGIN___
  Standard_Real Tol = 1.e-6;  // BRepFill_Precision();
  Standard_Real TolC = 0.;

  if ( !Degener) {
    Handle(Geom_Curve) C = BRep_Tool::Curve(Edge,L,f,l);
    CT = new Geom_TrimmedCurve(C,f,l);
    CT->Transform(L.Transformation());
    // projection of 3d curves in the plane xOy
    Handle(Geom2d_Curve) C2d = GeomProjLib::Curve2d(CT,Plane);

    Geom2dAdaptor_Curve AC(C2d);
    Geom2dAdaptor_Curve ABis(Bis);

    Intersector = Geom2dInt_GInter(ABis, AC, TolC, Tol);

    if ( !Intersector.IsDone()) {
      StdFail_NotDone::Raise("BRepFill_TrimSurfaceTool::IntersectWith");
    }

    NbPoints = Intersector.NbPoints();

    if (NbPoints < 1) {
      // try to elongate curves and enlarge tolerance

      // don't do it rightaway from the beginning in order not to get
      // extra solutions those would cause *Exception*: incoherent intersection

      GeomAbs_CurveType CType = AC.GetType(), BisType = ABis.GetType();
      Standard_Boolean  canElongateC = !(CType == GeomAbs_BezierCurve  ||
                                         CType == GeomAbs_BSplineCurve ||
                                         CType == GeomAbs_OffsetCurve  ||
                                         CType == GeomAbs_OtherCurve);
      Standard_Boolean canElongateBis = !(BisType == GeomAbs_BezierCurve  ||
                                          BisType == GeomAbs_BSplineCurve ||
                                          BisType == GeomAbs_OffsetCurve  ||
                                          BisType == GeomAbs_OtherCurve);
      
      Handle(Geom2d_TrimmedCurve) TBis = Handle(Geom2d_TrimmedCurve)::DownCast(Bis);
      Handle(Geom2d_TrimmedCurve) TC2d = Handle(Geom2d_TrimmedCurve)::DownCast(C2d);

      if (canElongateC) {
        TC2d->SetTrim(TC2d->FirstParameter() - Tol, TC2d->LastParameter() + Tol);
        AC.Load(TC2d);
      }
      if (canElongateBis) {
        TBis->SetTrim(TBis->FirstParameter() - Tol, TBis->LastParameter() + Tol);
        ABis.Load(TBis);
      }
      Intersector = Geom2dInt_GInter(ABis, AC, TolC, Tol*10);

      if ( !Intersector.IsDone()) {
        StdFail_NotDone::Raise("BRepFill_TrimSurfaceTool::IntersectWith");
      }

      NbPoints = Intersector.NbPoints();
    }
    //  modified by NIZHNY-EAP Wed Dec 22 15:00:56 1999 ___END___
    if (NbPoints > 0) {

      for ( Standard_Integer i = 1; i <= NbPoints; i++) {
        U1 = Intersector.Point(i).ParamOnFirst();
        U2 = Intersector.Point(i).ParamOnSecond();
        P = gp_Pnt(U1,U2,0.);
        Seq.Append(P);
      }
    }

    NbSegments = Intersector.NbSegments();

    if (NbSegments > 0) {
#ifdef OCCT_DEBUG
      cout << " IntersectWith : " << NbSegments  
        << " Segments of intersection" << endl;
#endif
      IntRes2d_IntersectionSegment Seg;
      for ( Standard_Integer i = 1; i <= NbSegments; i++) {
        Seg = Intersector.Segment(i);
        U1  = Seg.FirstPoint().ParamOnFirst();
        U1 += Seg.LastPoint().ParamOnFirst();
        U1 /= 2.;
        U2  = Seg.FirstPoint().ParamOnSecond();
        U2 += Seg.LastPoint().ParamOnSecond();
        U2 /= 2.;
        P = gp_Pnt(U1,U2,0.);
        Seq.Append(P);
      }
    }
    // Order the sequence by increasing parameter on the bissectrice.
    Bubble( Seq);

    //  modified by NIZHNY-EAP Fri Dec 24 18:47:24 1999 ___BEGIN___
    // Remove double points
    gp_Pnt P1, P2;
    for ( Standard_Integer i = 1; i < NbPoints; i++) {
      P1 = Seq.Value(i);
      P2 = Seq.Value(i+1);
      if ( P2.X()-P1.X() < Tol ) {
        //	cout<<"REMOVE "<<P1.X()<<endl;
        Seq.Remove(i--);
        NbPoints--;
      }
    }
    //  modified by NIZHNY-EAP Fri Dec 24 18:47:28 1999 ___END___
  }
  else {
    // the edge is degenerated : the point and it is found if it is
    // on the bissectrice.

    gp_Pnt P3d = BRep_Tool::Pnt( TopExp::FirstVertex(Edge));
    gp_Pnt2d P2d( P3d.X(), P3d.Y());

    Standard_Real UBis = Bis->FirstParameter();
    gp_Pnt2d PBis = Bis->Value( UBis);

    //  modified by NIZHNY-EAP Wed Jan 12 11:41:30 2000 ___BEGIN___
    // inside gp_Pnt2d::Distance
    // Infinite * Infinite => Exception: DefaultNumericError
    // Case encounered: UBis < Precision::Infinite()
    // but PBis.X() > Precision::Infinite()
    if (Precision::IsPositiveInfinite(Abs(PBis.X())) ||
      Precision::IsPositiveInfinite(Abs(PBis.Y())) ||
      PBis.Distance(P2d) > Tol) {
        //  modified by NIZHNY-EAP Wed Jan 12 11:41:40 2000 ___END___
        UBis = Bis->LastParameter();
        if (UBis >= Precision::Infinite()) return;
        PBis = Bis->Value( UBis);
        if ( PBis.Distance(P2d) > Tol) return;
    }

    // evaluate parameter intersection.
    Handle(Geom_Surface) GS = BRep_Tool::Surface(Face);
    GeomAdaptor_Surface GAS(GS);

    gp_Ax3 Axis;
    Standard_Real Phase = 0.;

    switch ( GAS.GetType()) {

    case GeomAbs_Sphere: 
      Axis = GAS.Sphere().Position(); break;
    case GeomAbs_Cone: {
      //----------------------------------------------------------
      // if myFace1 is not at the same side of the apex as the point
      // of parameter 0 0 on the cone => phase = M_PI.
      //----------------------------------------------------------
      Axis = GAS.Cone().Position(); 
      Phase = EvalPhase(Edge,Face,GAS,Axis);
      break;
                       }
    case GeomAbs_Torus:
      Axis = GAS.Torus().Position(); break;
    case GeomAbs_Cylinder:
      Axis = GAS.Cylinder().Position(); break;
    case GeomAbs_SurfaceOfRevolution: {      
      //----------------------------------------------------------
      // if myFace1 is not at the same side of the apex as the point
      // of parameter 0 0 on the cone => phase = M_PI.
      //----------------------------------------------------------
      Handle(Geom_SurfaceOfRevolution) GSRev = 
        Handle(Geom_SurfaceOfRevolution)::DownCast(GS);
      Handle(GeomAdaptor_HCurve) HC = 
        new GeomAdaptor_HCurve(GSRev->BasisCurve());
      GeomAdaptor_SurfaceOfRevolution ASRev(HC,GAS.AxeOfRevolution());
      Axis  = ASRev.Axis();       
      Phase = EvalPhase(Edge,Face,GAS,Axis);
      break;
                                      }
    default:
      Standard_NotImplemented::Raise(" BRepFill_TrimSurfaceTool");
    }

    gp_Vec2d D12d = Bis->DN(UBis,1);
    gp_Vec   D1( D12d.X(), D12d.Y(), 0.);

    Standard_Real U = Axis.XDirection().
      AngleWithRef(D1,Axis.XDirection()^Axis.YDirection());
    U += Phase;
    if ( U < 0.) U += 2*M_PI;

    P = gp_Pnt(Bis->FirstParameter(), U, 0.);
    Seq.Append(P);
  }
}


//=======================================================================
//function : IntersectWith
//purpose  : 
//=======================================================================

void BRepFill_TrimSurfaceTool::IntersectWith
  (const TopoDS_Edge&          EdgeOnF1,
  const TopoDS_Edge&          EdgeOnF2,
  TColgp_SequenceOfPnt& Points   ) 
  const 
{
#ifdef DRAW
  if ( Affich) {
    char name[256];  
    Standard_Integer i1 = 0, i2 = 2;
    sprintf(name,"EdgeOnF1_%d_%d",i1, NBCALL);
    DBRep::Set(name,EdgeOnF1);

    sprintf(name,"EdgeOnF2_%d_%d",i2, NBCALL);
    DBRep::Set(name,EdgeOnF2);
  }

#endif 
  Points.Clear();
  TColgp_SequenceOfPnt Points2;

  EvalParameters(EdgeOnF1, myFace1, myBis, Points);
  EvalParameters(EdgeOnF2, myFace2, myBis, Points2);

  StdFail_NotDone_Raise_if
    ( Points.Length() != Points2.Length(),
    "BRepFill_TrimSurfaceTool::IntersectWith: incoherent intersection");

  gp_Pnt PSeq;
  Standard_Integer NbPoints = Points.Length();
  for ( Standard_Integer i = 1; i <= NbPoints; i++) {
    PSeq = Points(i);
    PSeq.SetZ((Points2.Value(i)).Y());
    Points.SetValue(i,PSeq);
    //    cout<<"BisPar "<<PSeq.X()<<endl;
  }
}


//=======================================================================
//function : IsOnFace
//purpose  : 
//=======================================================================

Standard_Boolean BRepFill_TrimSurfaceTool::IsOnFace
  (const gp_Pnt2d& Point) const 
{
  gp_Pnt P( Point.X(), Point.Y(), 0.);
  gp_Lin Line( P, gp::DZ());

  BRepIntCurveSurface_Inter Inter;

  // eval if is on face 1
  //  modified by NIZHNY-EAP Fri Jan 21 09:49:09 2000 ___BEGIN___
  Inter.Init(myFace1, Line,1e-6);//Precision::PConfusion());
  if (Inter.More()) return Standard_True;

  // eval if is on face 2
  Inter.Init(myFace2, Line, 1e-6);//Precision::PConfusion());
  return Inter.More();
  //  modified by NIZHNY-EAP Fri Jan 21 09:49:14 2000 ___END___
}


//=======================================================================
//function : ProjOn
//purpose  : 
//=======================================================================

Standard_Real BRepFill_TrimSurfaceTool::ProjOn(const gp_Pnt2d& Point,
  const TopoDS_Edge& Edge) const 
{
  TopLoc_Location L;
  Standard_Real   f,l;


  Handle(Geom_Curve) C1 = BRep_Tool::Curve(Edge,L,f,l);
  Handle(Geom_TrimmedCurve) CT = new Geom_TrimmedCurve(C1,f,l);
  CT->Transform(L.Transformation());

  // projection of curves 3d in the plane xOy
  Handle(Geom_Plane) Plane = new Geom_Plane(0,0,1,0);
  Handle(Geom2d_Curve) C2d = GeomProjLib::Curve2d(CT,Plane);

  // evaluate the projection of the point on the curve.
  Geom2dAPI_ProjectPointOnCurve Projector(Point, C2d);
#ifdef OCCT_DEBUG
  Standard_Real Dist = Projector.LowerDistance();
  if ( Dist > Precision::Confusion() ) {
    cout << " *** WARNING  TrimSurfaceTool:  *** " << endl;
    cout << "      --> the point is not on the edge" <<endl;
    cout << "          distance  = " << Dist << endl;
  }
#endif

  Standard_Real U = Projector.LowerDistanceParameter();
  return U;
}


//=======================================================================
//function : Project
//purpose  : 
//=======================================================================

void BRepFill_TrimSurfaceTool::Project
  (const Standard_Real U1,
  const Standard_Real U2,
  Handle(Geom_Curve)&   Curve,
  Handle(Geom2d_Curve)& PCurve1,
  Handle(Geom2d_Curve)& PCurve2,
  GeomAbs_Shape&        Cont) const 
{
  Handle(Geom2d_TrimmedCurve) CT = 
    new Geom2d_TrimmedCurve(myBis,U1,U2);
  BRepFill_MultiLine ML(myFace1,myFace2,
    myEdge1,myEdge2,myInv1,myInv2,CT);

  Cont = ML.Continuity();

  if ( ML.IsParticularCase()) {
    ML.Curves(Curve,PCurve1,PCurve2);
  }
  else {
    BRepFill_ApproxSeewing AppSeew(ML);

    Curve   = AppSeew.Curve();
    PCurve1 = AppSeew.CurveOnF1();
    PCurve2 = AppSeew.CurveOnF2();
  }
}
Commit	Line	Data
b311480e	1	// Created on: 1994-10-21
	2	// Created by: Bruno DUMORTIER
	3	// Copyright (c) 1994-1999 Matra Datavision
973c2be1	4	// Copyright (c) 1999-2014 OPEN CASCADE SAS
b311480e	5	//
973c2be1	6	// This file is part of Open CASCADE Technology software library.
b311480e	7	//
d5f74e42	8	// This library is free software; you can redistribute it and/or modify it under
d5f74e42	9	// the terms of the GNU Lesser General Public License version 2.1 as published
973c2be1	10	// by the Free Software Foundation, with special exception defined in the file
	11	// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
	12	// distribution for complete text of the license and disclaimer of any warranty.
b311480e	13	//
973c2be1	14	// Alternatively, this file may be used under the terms of Open CASCADE
973c2be1	15	// commercial license or contractual agreement.
7fd59977	16
7fd59977	17
6b84c3f7	18	#include <GeomAdaptor_SurfaceOfRevolution.hxx>
7fd59977	19	#include <AppParCurves_MultiCurve.hxx>
42cf5bc1	20	#include <BRep_Tool.hxx>
42cf5bc1	21	#include <BRepFill_ApproxSeewing.hxx>
7fd59977	22	#include <BRepFill_ComputeCLine.hxx>
7fd59977	23	#include <BRepFill_MultiLine.hxx>
42cf5bc1	24	#include <BRepFill_TrimSurfaceTool.hxx>
7fd59977	25	#include <BRepIntCurveSurface_Inter.hxx>
7fd59977	26	#include <BSplCLib.hxx>
42cf5bc1	27	#include <ElCLib.hxx>
7fd59977	28	#include <Geom2d_BSplineCurve.hxx>
42cf5bc1	29	#include <Geom2d_Curve.hxx>
7fd59977	30	#include <Geom2d_TrimmedCurve.hxx>
42cf5bc1	31	#include <Geom2dAdaptor_Curve.hxx>
	32	#include <Geom2dAPI_ProjectPointOnCurve.hxx>
	33	#include <Geom2dInt_GInter.hxx>
	34	#include <Geom_BSplineCurve.hxx>
	35	#include <Geom_Curve.hxx>
7fd59977	36	#include <Geom_Plane.hxx>
7fd59977	37	#include <Geom_SurfaceOfRevolution.hxx>
7fd59977	38	#include <Geom_TrimmedCurve.hxx>
42cf5bc1	39	#include <GeomAdaptor_HCurve.hxx>
42cf5bc1	40	#include <GeomAdaptor_Surface.hxx>
7fd59977	41	#include <GeomProjLib.hxx>
42cf5bc1	42	#include <gp.hxx>
	43	#include <gp_Cone.hxx>
	44	#include <gp_Cylinder.hxx>
	45	#include <gp_Pnt2d.hxx>
	46	#include <gp_Sphere.hxx>
	47	#include <gp_Torus.hxx>
7fd59977	48	#include <IntRes2d_IntersectionPoint.hxx>
7fd59977	49	#include <IntRes2d_IntersectionSegment.hxx>
42cf5bc1	50	#include <PLib.hxx>
7fd59977	51	#include <Precision.hxx>
42cf5bc1	52	#include <Standard_NoSuchObject.hxx>
7fd59977	53	#include <Standard_NotImplemented.hxx>
42cf5bc1	54	#include <StdFail_NotDone.hxx>
7fd59977	55	#include <TColgp_Array1OfPnt.hxx>
7fd59977	56	#include <TColStd_Array1OfInteger.hxx>
42cf5bc1	57	#include <TColStd_Array1OfReal.hxx>
7fd59977	58	#include <TopAbs_State.hxx>
7fd59977	59	#include <TopExp.hxx>
42cf5bc1	60	#include <TopExp_Explorer.hxx>
7fd59977	61	#include <TopoDS.hxx>
42cf5bc1	62	#include <TopoDS_Edge.hxx>
42cf5bc1	63	#include <TopoDS_Face.hxx>
7fd59977	64	#include <TopoDS_Vertex.hxx>
7fd59977	65
df573a26	66	//#define DRAW
42cf5bc1	67	#include <stdio.h>
7fd59977	68	#ifdef DRAW
	69	#include <DrawTrSurf.hxx>
	70	#include <DBRep.hxx>
	71	#endif
0797d9d3	72	#ifdef OCCT_DEBUG
7fd59977	73	static Standard_Boolean Affich = Standard_False;
df573a26	74	static Standard_Integer NBCALL = 0;
7fd59977	75	#endif
	76
	77	//=======================================================================
	78	//function : BRepFill_TrimSurfaceTool
0d969553 Y	79	//purpose : Initialisation with two neighbor faces
	80	// Edge1 and Edge2 are parallel edges corresponding
	81	// to minimum iso on F1 and F2 respectively.
	82	// ie Edge1 is Umin or VMin on F1.
	83	// Inv1 and Inv2 show if Edge1 and Edge2 are
	84	// returned parallel.
7fd59977	85	//=======================================================================
	86
	87	BRepFill_TrimSurfaceTool::BRepFill_TrimSurfaceTool
df573a26	88	(const Handle(Geom2d_Curve)& Bis,
	89	const TopoDS_Face& Face1,
	90	const TopoDS_Face& Face2,
	91	const TopoDS_Edge& Edge1,
	92	const TopoDS_Edge& Edge2,
	93	const Standard_Boolean Inv1,
	94	const Standard_Boolean Inv2 ) :
7fd59977	95	myFace1(Face1),
df573a26	96	myFace2(Face2),
	97	myEdge1(Edge1),
	98	myEdge2(Edge2),
	99	myInv1(Inv1),
	100	myInv2(Inv2),
	101	myBis (Bis)
7fd59977	102	{
0797d9d3	103	#ifdef OCCT_DEBUG
7fd59977	104	if ( Affich) {
df573a26	105	NBCALL++;
7fd59977	106	cout << " ---------->TrimSurfaceTool : NBCALL = " << NBCALL << endl;
7fd59977	107	#ifdef DRAW
7fd59977	108	char name[256];
df573a26	109
7fd59977	110	sprintf(name,"FACE1_%d",NBCALL);
7fd59977	111	DBRep::Set(name,myFace1);
df573a26	112
7fd59977	113	sprintf(name,"FACE2_%d",NBCALL);
7fd59977	114	DBRep::Set(name,myFace2);
df573a26	115
7fd59977	116	sprintf(name,"EDGE1_%d",NBCALL);
7fd59977	117	DBRep::Set(name,myEdge1);
df573a26	118
7fd59977	119	sprintf(name,"EDGE2_%d",NBCALL);
7fd59977	120	DBRep::Set(name,myEdge2);
df573a26	121
7fd59977	122	sprintf(name,"BISSEC_%d",NBCALL);
	123	DrawTrSurf::Set(name,myBis);
	124	#endif
7fd59977	125	}
	126	#endif
	127	}
	128
	129
	130	//=======================================================================
	131	//function : Bubble
0d969553	132	//purpose : Order the sequence of points by increasing x.
7fd59977	133	//=======================================================================
	134
	135	static void Bubble(TColgp_SequenceOfPnt& Seq)
	136	{
	137	Standard_Boolean Invert = Standard_True;
	138	Standard_Integer NbPoints = Seq.Length();
	139	while (Invert) {
	140	Invert = Standard_False;
	141	for ( Standard_Integer i = 1; i < NbPoints; i++) {
	142	gp_Pnt P1 = Seq.Value(i);
	143	gp_Pnt P2 = Seq.Value(i+1);
	144	if (P2.X()<P1.X()) {
df573a26	145	Seq.Exchange(i,i+1);
df573a26	146	Invert = Standard_True;
7fd59977	147	}
	148	}
	149	}
	150	}
	151
	152
	153	//=======================================================================
	154	//function : EvalPhase
	155	//purpose :
	156	//=======================================================================
	157
	158	static Standard_Real EvalPhase(const TopoDS_Edge& Edge,
df573a26	159	const TopoDS_Face& Face,
	160	const GeomAdaptor_Surface& GAS,
	161	const gp_Ax3& Axis)
7fd59977	162	{
	163	gp_Pnt2d PE1,PE2,PF1,PF2;
	164	Standard_Real VDeg;
7fd59977	165	Standard_Real V = 0.;
7fd59977	166	BRep_Tool::UVPoints(Edge,Face,PE1,PE2);
	167	VDeg = PE1.Y();
	168	TopExp_Explorer Exp(Face,TopAbs_EDGE);
	169	for (; Exp.More(); Exp.Next()) {
	170	if ( !TopoDS::Edge(Exp.Current()).IsSame(Edge)) {
	171	BRep_Tool::UVPoints(TopoDS::Edge(Exp.Current()),Face,PF1,PF2);
	172	V = ( Abs(PF1.Y() - VDeg) > Abs(PF2.Y() - VDeg)) ? PF1.Y() : PF2.Y();
	173	break;
	174	}
	175	}
	176	gp_Pnt P = GAS.Value(0., V);
df573a26	177
7fd59977	178	if ( gp_Vec(Axis.Location(), P).Dot(Axis.XDirection()) < 0.)
c6541a0c	179	return M_PI;
7fd59977	180	else
	181	return 0.;
	182	}
	183
	184
	185	//=======================================================================
	186	//function : EvalParameters
	187	//purpose :
	188	//=======================================================================
	189
	190	static void EvalParameters(const TopoDS_Edge& Edge,
df573a26	191	const TopoDS_Face& Face,
	192	const Handle(Geom2d_Curve)& Bis ,
	193	TColgp_SequenceOfPnt& Seq )
7fd59977	194	{
7fd59977	195	Standard_Boolean Degener = BRep_Tool::Degenerated(Edge);
0d969553	196	// return curves 3d associated to edges.
7fd59977	197	TopLoc_Location L;
	198	Standard_Real f,l;
	199
	200	Handle(Geom_TrimmedCurve) CT;
	201	Handle(Geom_Plane) Plane = new Geom_Plane(0,0,1,0);
	202
	203	Geom2dInt_GInter Intersector;
	204
	205	Standard_Integer NbPoints, NbSegments;
	206	Standard_Real U1, U2;
	207	gp_Pnt P;//,PSeq;
	208
df573a26	209	// Standard_Real Tol = Precision::Intersection();
7fd59977	210	// modified by NIZHNY-EAP Wed Dec 22 15:00:51 1999 ___BEGIN___
	211	Standard_Real Tol = 1.e-6; // BRepFill_Precision();
	212	Standard_Real TolC = 0.;
df573a26	213
7fd59977	214	if ( !Degener) {
	215	Handle(Geom_Curve) C = BRep_Tool::Curve(Edge,L,f,l);
	216	CT = new Geom_TrimmedCurve(C,f,l);
	217	CT->Transform(L.Transformation());
0d969553	218	// projection of 3d curves in the plane xOy
7fd59977	219	Handle(Geom2d_Curve) C2d = GeomProjLib::Curve2d(CT,Plane);
	220
	221	Geom2dAdaptor_Curve AC(C2d);
	222	Geom2dAdaptor_Curve ABis(Bis);
	223
	224	Intersector = Geom2dInt_GInter(ABis, AC, TolC, Tol);
df573a26	225
7fd59977	226	if ( !Intersector.IsDone()) {
	227	StdFail_NotDone::Raise("BRepFill_TrimSurfaceTool::IntersectWith");
	228	}
df573a26	229
7fd59977	230	NbPoints = Intersector.NbPoints();
	231
	232	if (NbPoints < 1) {
	233	// try to elongate curves and enlarge tolerance
df573a26	234
7fd59977	235	// don't do it rightaway from the beginning in order not to get
7fd59977	236	// extra solutions those would cause Exception: incoherent intersection
df573a26	237
7fd59977	238	GeomAbs_CurveType CType = AC.GetType(), BisType = ABis.GetType();
1aec3320	239	Standard_Boolean canElongateC = !(CType == GeomAbs_BezierCurve \|\|
	240	CType == GeomAbs_BSplineCurve \|\|
	241	CType == GeomAbs_OffsetCurve \|\|
	242	CType == GeomAbs_OtherCurve);
	243	Standard_Boolean canElongateBis = !(BisType == GeomAbs_BezierCurve \|\|
	244	BisType == GeomAbs_BSplineCurve \|\|
	245	BisType == GeomAbs_OffsetCurve \|\|
	246	BisType == GeomAbs_OtherCurve);
7fd59977	247
	248	Handle(Geom2d_TrimmedCurve) TBis = Handle(Geom2d_TrimmedCurve)::DownCast(Bis);
	249	Handle(Geom2d_TrimmedCurve) TC2d = Handle(Geom2d_TrimmedCurve)::DownCast(C2d);
df573a26	250
7fd59977	251	if (canElongateC) {
df573a26	252	TC2d->SetTrim(TC2d->FirstParameter() - Tol, TC2d->LastParameter() + Tol);
df573a26	253	AC.Load(TC2d);
7fd59977	254	}
7fd59977	255	if (canElongateBis) {
df573a26	256	TBis->SetTrim(TBis->FirstParameter() - Tol, TBis->LastParameter() + Tol);
df573a26	257	ABis.Load(TBis);
7fd59977	258	}
7fd59977	259	Intersector = Geom2dInt_GInter(ABis, AC, TolC, Tol*10);
df573a26	260
7fd59977	261	if ( !Intersector.IsDone()) {
df573a26	262	StdFail_NotDone::Raise("BRepFill_TrimSurfaceTool::IntersectWith");
	263	}
	264
7fd59977	265	NbPoints = Intersector.NbPoints();
	266	}
	267	// modified by NIZHNY-EAP Wed Dec 22 15:00:56 1999 ___END___
	268	if (NbPoints > 0) {
df573a26	269
7fd59977	270	for ( Standard_Integer i = 1; i <= NbPoints; i++) {
df573a26	271	U1 = Intersector.Point(i).ParamOnFirst();
	272	U2 = Intersector.Point(i).ParamOnSecond();
	273	P = gp_Pnt(U1,U2,0.);
	274	Seq.Append(P);
7fd59977	275	}
	276	}
	277
	278	NbSegments = Intersector.NbSegments();
df573a26	279
7fd59977	280	if (NbSegments > 0) {
0797d9d3	281	#ifdef OCCT_DEBUG
7fd59977	282	cout << " IntersectWith : " << NbSegments
df573a26	283	<< " Segments of intersection" << endl;
7fd59977	284	#endif
	285	IntRes2d_IntersectionSegment Seg;
	286	for ( Standard_Integer i = 1; i <= NbSegments; i++) {
df573a26	287	Seg = Intersector.Segment(i);
	288	U1 = Seg.FirstPoint().ParamOnFirst();
	289	U1 += Seg.LastPoint().ParamOnFirst();
	290	U1 /= 2.;
	291	U2 = Seg.FirstPoint().ParamOnSecond();
	292	U2 += Seg.LastPoint().ParamOnSecond();
	293	U2 /= 2.;
	294	P = gp_Pnt(U1,U2,0.);
	295	Seq.Append(P);
7fd59977	296	}
7fd59977	297	}
0d969553	298	// Order the sequence by increasing parameter on the bissectrice.
7fd59977	299	Bubble( Seq);
df573a26	300
df573a26	301	// modified by NIZHNY-EAP Fri Dec 24 18:47:24 1999 ___BEGIN___
7fd59977	302	// Remove double points
	303	gp_Pnt P1, P2;
	304	for ( Standard_Integer i = 1; i < NbPoints; i++) {
	305	P1 = Seq.Value(i);
	306	P2 = Seq.Value(i+1);
	307	if ( P2.X()-P1.X() < Tol ) {
df573a26	308	// cout<<"REMOVE "<<P1.X()<<endl;
	309	Seq.Remove(i--);
	310	NbPoints--;
7fd59977	311	}
7fd59977	312	}
df573a26	313	// modified by NIZHNY-EAP Fri Dec 24 18:47:28 1999 ___END___
7fd59977	314	}
7fd59977	315	else {
0d969553 Y	316	// the edge is degenerated : the point and it is found if it is
0d969553 Y	317	// on the bissectrice.
7fd59977	318
	319	gp_Pnt P3d = BRep_Tool::Pnt( TopExp::FirstVertex(Edge));
	320	gp_Pnt2d P2d( P3d.X(), P3d.Y());
	321
	322	Standard_Real UBis = Bis->FirstParameter();
	323	gp_Pnt2d PBis = Bis->Value( UBis);
df573a26	324
df573a26	325	// modified by NIZHNY-EAP Wed Jan 12 11:41:30 2000 ___BEGIN___
7fd59977	326	// inside gp_Pnt2d::Distance
	327	// Infinite * Infinite => Exception: DefaultNumericError
	328	// Case encounered: UBis < Precision::Infinite()
	329	// but PBis.X() > Precision::Infinite()
	330	if (Precision::IsPositiveInfinite(Abs(PBis.X())) \|\|
df573a26	331	Precision::IsPositiveInfinite(Abs(PBis.Y())) \|\|
	332	PBis.Distance(P2d) > Tol) {
	333	// modified by NIZHNY-EAP Wed Jan 12 11:41:40 2000 ___END___
	334	UBis = Bis->LastParameter();
	335	if (UBis >= Precision::Infinite()) return;
	336	PBis = Bis->Value( UBis);
	337	if ( PBis.Distance(P2d) > Tol) return;
7fd59977	338	}
7fd59977	339
0d969553	340	// evaluate parameter intersection.
7fd59977	341	Handle(Geom_Surface) GS = BRep_Tool::Surface(Face);
	342	GeomAdaptor_Surface GAS(GS);
	343
	344	gp_Ax3 Axis;
	345	Standard_Real Phase = 0.;
	346
	347	switch ( GAS.GetType()) {
df573a26	348
7fd59977	349	case GeomAbs_Sphere:
	350	Axis = GAS.Sphere().Position(); break;
	351	case GeomAbs_Cone: {
	352	//----------------------------------------------------------
0d969553	353	// if myFace1 is not at the same side of the apex as the point
c6541a0c	354	// of parameter 0 0 on the cone => phase = M_PI.
7fd59977	355	//----------------------------------------------------------
	356	Axis = GAS.Cone().Position();
	357	Phase = EvalPhase(Edge,Face,GAS,Axis);
	358	break;
df573a26	359	}
7fd59977	360	case GeomAbs_Torus:
	361	Axis = GAS.Torus().Position(); break;
	362	case GeomAbs_Cylinder:
	363	Axis = GAS.Cylinder().Position(); break;
	364	case GeomAbs_SurfaceOfRevolution: {
	365	//----------------------------------------------------------
0d969553	366	// if myFace1 is not at the same side of the apex as the point
c6541a0c	367	// of parameter 0 0 on the cone => phase = M_PI.
7fd59977	368	//----------------------------------------------------------
7fd59977	369	Handle(Geom_SurfaceOfRevolution) GSRev =
df573a26	370	Handle(Geom_SurfaceOfRevolution)::DownCast(GS);
7fd59977	371	Handle(GeomAdaptor_HCurve) HC =
df573a26	372	new GeomAdaptor_HCurve(GSRev->BasisCurve());
6b84c3f7	373	GeomAdaptor_SurfaceOfRevolution ASRev(HC,GAS.AxeOfRevolution());
7fd59977	374	Axis = ASRev.Axis();
	375	Phase = EvalPhase(Edge,Face,GAS,Axis);
	376	break;
df573a26	377	}
7fd59977	378	default:
	379	Standard_NotImplemented::Raise(" BRepFill_TrimSurfaceTool");
	380	}
	381
	382	gp_Vec2d D12d = Bis->DN(UBis,1);
	383	gp_Vec D1( D12d.X(), D12d.Y(), 0.);
df573a26	384
7fd59977	385	Standard_Real U = Axis.XDirection().
	386	AngleWithRef(D1,Axis.XDirection()^Axis.YDirection());
	387	U += Phase;
c6541a0c	388	if ( U < 0.) U += 2*M_PI;
7fd59977	389
	390	P = gp_Pnt(Bis->FirstParameter(), U, 0.);
	391	Seq.Append(P);
	392	}
	393	}
	394
	395
	396	//=======================================================================
	397	//function : IntersectWith
	398	//purpose :
	399	//=======================================================================
	400
	401	void BRepFill_TrimSurfaceTool::IntersectWith
df573a26	402	(const TopoDS_Edge& EdgeOnF1,
	403	const TopoDS_Edge& EdgeOnF2,
	404	TColgp_SequenceOfPnt& Points )
	405	const
7fd59977	406	{
df573a26	407	#ifdef DRAW
	408	if ( Affich) {
	409	char name[256];
	410	Standard_Integer i1 = 0, i2 = 2;
	411	sprintf(name,"EdgeOnF1_%d_%d",i1, NBCALL);
	412	DBRep::Set(name,EdgeOnF1);
	413
	414	sprintf(name,"EdgeOnF2_%d_%d",i2, NBCALL);
	415	DBRep::Set(name,EdgeOnF2);
	416	}
	417
	418	#endif
7fd59977	419	Points.Clear();
	420	TColgp_SequenceOfPnt Points2;
	421
	422	EvalParameters(EdgeOnF1, myFace1, myBis, Points);
	423	EvalParameters(EdgeOnF2, myFace2, myBis, Points2);
	424
	425	StdFail_NotDone_Raise_if
	426	( Points.Length() != Points2.Length(),
df573a26	427	"BRepFill_TrimSurfaceTool::IntersectWith: incoherent intersection");
7fd59977	428
	429	gp_Pnt PSeq;
	430	Standard_Integer NbPoints = Points.Length();
	431	for ( Standard_Integer i = 1; i <= NbPoints; i++) {
	432	PSeq = Points(i);
	433	PSeq.SetZ((Points2.Value(i)).Y());
	434	Points.SetValue(i,PSeq);
df573a26	435	// cout<<"BisPar "<<PSeq.X()<<endl;
7fd59977	436	}
	437	}
	438
	439
	440	//=======================================================================
	441	//function : IsOnFace
	442	//purpose :
	443	//=======================================================================
	444
	445	Standard_Boolean BRepFill_TrimSurfaceTool::IsOnFace
df573a26	446	(const gp_Pnt2d& Point) const
7fd59977	447	{
	448	gp_Pnt P( Point.X(), Point.Y(), 0.);
	449	gp_Lin Line( P, gp::DZ());
	450
	451	BRepIntCurveSurface_Inter Inter;
	452
	453	// eval if is on face 1
df573a26	454	// modified by NIZHNY-EAP Fri Jan 21 09:49:09 2000 ___BEGIN___
004e8466	455	Inter.Init(myFace1, Line,1e-6);//Precision::PConfusion());
7fd59977	456	if (Inter.More()) return Standard_True;
df573a26	457
7fd59977	458	// eval if is on face 2
	459	Inter.Init(myFace2, Line, 1e-6);//Precision::PConfusion());
	460	return Inter.More();
df573a26	461	// modified by NIZHNY-EAP Fri Jan 21 09:49:14 2000 ___END___
7fd59977	462	}
	463
	464
	465	//=======================================================================
	466	//function : ProjOn
	467	//purpose :
	468	//=======================================================================
	469
	470	Standard_Real BRepFill_TrimSurfaceTool::ProjOn(const gp_Pnt2d& Point,
df573a26	471	const TopoDS_Edge& Edge) const
7fd59977	472	{
	473	TopLoc_Location L;
	474	Standard_Real f,l;
	475
df573a26	476
7fd59977	477
	478	Handle(Geom_Curve) C1 = BRep_Tool::Curve(Edge,L,f,l);
	479	Handle(Geom_TrimmedCurve) CT = new Geom_TrimmedCurve(C1,f,l);
	480	CT->Transform(L.Transformation());
df573a26	481
0d969553	482	// projection of curves 3d in the plane xOy
7fd59977	483	Handle(Geom_Plane) Plane = new Geom_Plane(0,0,1,0);
	484	Handle(Geom2d_Curve) C2d = GeomProjLib::Curve2d(CT,Plane);
	485
0d969553	486	// evaluate the projection of the point on the curve.
7fd59977	487	Geom2dAPI_ProjectPointOnCurve Projector(Point, C2d);
0797d9d3	488	#ifdef OCCT_DEBUG
63c629aa	489	Standard_Real Dist = Projector.LowerDistance();
7fd59977	490	if ( Dist > Precision::Confusion() ) {
7fd59977	491	cout << " * WARNING TrimSurfaceTool: * " << endl;
0d969553	492	cout << " --> the point is not on the edge" <<endl;
7fd59977	493	cout << " distance = " << Dist << endl;
	494	}
	495	#endif
	496
	497	Standard_Real U = Projector.LowerDistanceParameter();
	498	return U;
	499	}
	500
	501
	502	//=======================================================================
	503	//function : Project
	504	//purpose :
	505	//=======================================================================
	506
	507	void BRepFill_TrimSurfaceTool::Project
df573a26	508	(const Standard_Real U1,
	509	const Standard_Real U2,
	510	Handle(Geom_Curve)& Curve,
	511	Handle(Geom2d_Curve)& PCurve1,
	512	Handle(Geom2d_Curve)& PCurve2,
	513	GeomAbs_Shape& Cont) const
7fd59977	514	{
7fd59977	515	Handle(Geom2d_TrimmedCurve) CT =
	516	new Geom2d_TrimmedCurve(myBis,U1,U2);
	517	BRepFill_MultiLine ML(myFace1,myFace2,
df573a26	518	myEdge1,myEdge2,myInv1,myInv2,CT);
df573a26	519
7fd59977	520	Cont = ML.Continuity();
	521
	522	if ( ML.IsParticularCase()) {
	523	ML.Curves(Curve,PCurve1,PCurve2);
	524	}
	525	else {
	526	BRepFill_ApproxSeewing AppSeew(ML);
	527
	528	Curve = AppSeew.Curve();
	529	PCurve1 = AppSeew.CurveOnF1();
	530	PCurve2 = AppSeew.CurveOnF2();
	531	}
	532	}