src/BRepFill/BRepFill_TrimSurfaceTool.cxx

   1 // Created on: 1994-10-21
   2 // Created by: Bruno DUMORTIER
   3 // Copyright (c) 1994-1999 Matra Datavision
   4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   5 //
   6 // This file is part of Open CASCADE Technology software library.
   7 //
   8 // This library is free software; you can redistribute it and/or modify it under
   9 // the terms of the GNU Lesser General Public License version 2.1 as published
  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.
  13 //
  14 // Alternatively, this file may be used under the terms of Open CASCADE
  15 // commercial license or contractual agreement.
  16
  17
  18 #include <Adaptor3d_SurfaceOfRevolution.hxx>
  19 #include <AppParCurves_MultiCurve.hxx>
  20 #include <BRep_Tool.hxx>
  21 #include <BRepFill_ApproxSeewing.hxx>
  22 #include <BRepFill_ComputeCLine.hxx>
  23 #include <BRepFill_MultiLine.hxx>
  24 #include <BRepFill_TrimSurfaceTool.hxx>
  25 #include <BRepIntCurveSurface_Inter.hxx>
  26 #include <BSplCLib.hxx>
  27 #include <ElCLib.hxx>
  28 #include <Geom2d_BSplineCurve.hxx>
  29 #include <Geom2d_Curve.hxx>
  30 #include <Geom2d_TrimmedCurve.hxx>
  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>
  36 #include <Geom_Plane.hxx>
  37 #include <Geom_SurfaceOfRevolution.hxx>
  38 #include <Geom_TrimmedCurve.hxx>
  39 #include <GeomAdaptor_HCurve.hxx>
  40 #include <GeomAdaptor_Surface.hxx>
  41 #include <GeomProjLib.hxx>
  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>
  48 #include <IntRes2d_IntersectionPoint.hxx>
  49 #include <IntRes2d_IntersectionSegment.hxx>
  50 #include <PLib.hxx>
  51 #include <Precision.hxx>
  52 #include <Standard_NoSuchObject.hxx>
  53 #include <Standard_NotImplemented.hxx>
  54 #include <StdFail_NotDone.hxx>
  55 #include <TColgp_Array1OfPnt.hxx>
  56 #include <TColStd_Array1OfInteger.hxx>
  57 #include <TColStd_Array1OfReal.hxx>
  58 #include <TopAbs_State.hxx>
  59 #include <TopExp.hxx>
  60 #include <TopExp_Explorer.hxx>
  61 #include <TopoDS.hxx>
  62 #include <TopoDS_Edge.hxx>
  63 #include <TopoDS_Face.hxx>
  64 #include <TopoDS_Vertex.hxx>
  65
  66 //#define DRAW
  67 #include <stdio.h>
  68 #ifdef DRAW
  69 #include <DrawTrSurf.hxx>
  70 #include <DBRep.hxx>
  71 #endif
  72 #ifdef OCCT_DEBUG
  73 static Standard_Boolean Affich       = Standard_False;
  74 static Standard_Integer NBCALL  = 0;
  75 #endif
  76
  77 //=======================================================================
  78 //function : BRepFill_TrimSurfaceTool
  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.
  85 //=======================================================================
  86
  87 BRepFill_TrimSurfaceTool::BRepFill_TrimSurfaceTool
  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 ) :
  95 myFace1(Face1),
  96   myFace2(Face2),
  97   myEdge1(Edge1),
  98   myEdge2(Edge2),
  99   myInv1(Inv1),
 100   myInv2(Inv2),
 101   myBis  (Bis)
 102 {
 103 #ifdef OCCT_DEBUG
 104   if ( Affich) {
 105     NBCALL++;
 106     cout << " ---------->TrimSurfaceTool : NBCALL = " << NBCALL << endl;
 107 #ifdef DRAW
 108     char name[256];
 109
 110     sprintf(name,"FACE1_%d",NBCALL);
 111     DBRep::Set(name,myFace1);
 112
 113     sprintf(name,"FACE2_%d",NBCALL);
 114     DBRep::Set(name,myFace2);
 115
 116     sprintf(name,"EDGE1_%d",NBCALL);
 117     DBRep::Set(name,myEdge1);
 118
 119     sprintf(name,"EDGE2_%d",NBCALL);
 120     DBRep::Set(name,myEdge2);
 121
 122     sprintf(name,"BISSEC_%d",NBCALL);
 123     DrawTrSurf::Set(name,myBis);
 124 #endif
 125   }
 126 #endif
 127 }
 128
 129
 130 //=======================================================================
 131 //function : Bubble
 132 //purpose  : Order the sequence of points by increasing x.
 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())  {
 145         Seq.Exchange(i,i+1);
 146         Invert = Standard_True;
 147       }
 148     }
 149   }
 150 }
 151
 152
 153 //=======================================================================
 154 //function : EvalPhase
 155 //purpose  :
 156 //=======================================================================
 157
 158 static Standard_Real EvalPhase(const TopoDS_Edge& Edge,
 159   const TopoDS_Face& Face,
 160   const GeomAdaptor_Surface& GAS,
 161   const gp_Ax3& Axis)
 162 {
 163   gp_Pnt2d PE1,PE2,PF1,PF2;
 164   Standard_Real VDeg;
 165   Standard_Real V = 0.;
 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);
 177
 178   if ( gp_Vec(Axis.Location(), P).Dot(Axis.XDirection()) < 0.)
 179     return M_PI;
 180   else
 181     return 0.;
 182 }
 183
 184
 185 //=======================================================================
 186 //function : EvalParameters
 187 //purpose  :
 188 //=======================================================================
 189
 190 static void EvalParameters(const TopoDS_Edge&          Edge,
 191   const TopoDS_Face&          Face,
 192   const Handle(Geom2d_Curve)& Bis ,
 193   TColgp_SequenceOfPnt& Seq  )
 194 {
 195   Standard_Boolean Degener = BRep_Tool::Degenerated(Edge);
 196   // return curves 3d associated to edges.
 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
 209   //  Standard_Real Tol = Precision::Intersection();
 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.;
 213
 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());
 218     // projection of 3d curves in the plane xOy
 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);
 225
 226     if ( !Intersector.IsDone()) {
 227       StdFail_NotDone::Raise("BRepFill_TrimSurfaceTool::IntersectWith");
 228     }
 229
 230     NbPoints = Intersector.NbPoints();
 231
 232     if (NbPoints < 1) {
 233       // try to elongate curves and enlarge tolerance
 234
 235       // don't do it rightaway from the beginning in order not to get
 236       // extra solutions those would cause *Exception*: incoherent intersection
 237
 238       GeomAbs_CurveType CType = AC.GetType(), BisType = ABis.GetType();
 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);
 247
 248       Handle(Geom2d_TrimmedCurve) TBis = Handle(Geom2d_TrimmedCurve)::DownCast(Bis);
 249       Handle(Geom2d_TrimmedCurve) TC2d = Handle(Geom2d_TrimmedCurve)::DownCast(C2d);
 250
 251       if (canElongateC) {
 252         TC2d->SetTrim(TC2d->FirstParameter() - Tol, TC2d->LastParameter() + Tol);
 253         AC.Load(TC2d);
 254       }
 255       if (canElongateBis) {
 256         TBis->SetTrim(TBis->FirstParameter() - Tol, TBis->LastParameter() + Tol);
 257         ABis.Load(TBis);
 258       }
 259       Intersector = Geom2dInt_GInter(ABis, AC, TolC, Tol*10);
 260
 261       if ( !Intersector.IsDone()) {
 262         StdFail_NotDone::Raise("BRepFill_TrimSurfaceTool::IntersectWith");
 263       }
 264
 265       NbPoints = Intersector.NbPoints();
 266     }
 267     //  modified by NIZHNY-EAP Wed Dec 22 15:00:56 1999 ___END___
 268     if (NbPoints > 0) {
 269
 270       for ( Standard_Integer i = 1; i <= NbPoints; i++) {
 271         U1 = Intersector.Point(i).ParamOnFirst();
 272         U2 = Intersector.Point(i).ParamOnSecond();
 273         P = gp_Pnt(U1,U2,0.);
 274         Seq.Append(P);
 275       }
 276     }
 277
 278     NbSegments = Intersector.NbSegments();
 279
 280     if (NbSegments > 0) {
 281 #ifdef OCCT_DEBUG
 282       cout << " IntersectWith : " << NbSegments
 283         << " Segments of intersection" << endl;
 284 #endif
 285       IntRes2d_IntersectionSegment Seg;
 286       for ( Standard_Integer i = 1; i <= NbSegments; i++) {
 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);
 296       }
 297     }
 298     // Order the sequence by increasing parameter on the bissectrice.
 299     Bubble( Seq);
 300
 301     //  modified by NIZHNY-EAP Fri Dec 24 18:47:24 1999 ___BEGIN___
 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 ) {
 308         //      cout<<"REMOVE "<<P1.X()<<endl;
 309         Seq.Remove(i--);
 310         NbPoints--;
 311       }
 312     }
 313     //  modified by NIZHNY-EAP Fri Dec 24 18:47:28 1999 ___END___
 314   }
 315   else {
 316     // the edge is degenerated : the point and it is found if it is
 317     // on the bissectrice.
 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);
 324
 325     //  modified by NIZHNY-EAP Wed Jan 12 11:41:30 2000 ___BEGIN___
 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())) ||
 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;
 338     }
 339
 340     // evaluate parameter intersection.
 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()) {
 348
 349     case GeomAbs_Sphere:
 350       Axis = GAS.Sphere().Position(); break;
 351     case GeomAbs_Cone: {
 352       //----------------------------------------------------------
 353       // if myFace1 is not at the same side of the apex as the point
 354       // of parameter 0 0 on the cone => phase = M_PI.
 355       //----------------------------------------------------------
 356       Axis = GAS.Cone().Position();
 357       Phase = EvalPhase(Edge,Face,GAS,Axis);
 358       break;
 359                        }
 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       //----------------------------------------------------------
 366       // if myFace1 is not at the same side of the apex as the point
 367       // of parameter 0 0 on the cone => phase = M_PI.
 368       //----------------------------------------------------------
 369       Handle(Geom_SurfaceOfRevolution) GSRev =
 370         Handle(Geom_SurfaceOfRevolution)::DownCast(GS);
 371       Handle(GeomAdaptor_HCurve) HC =
 372         new GeomAdaptor_HCurve(GSRev->BasisCurve());
 373       Adaptor3d_SurfaceOfRevolution ASRev(HC,GAS.AxeOfRevolution());
 374       Axis  = ASRev.Axis();
 375       Phase = EvalPhase(Edge,Face,GAS,Axis);
 376       break;
 377                                       }
 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.);
 384
 385     Standard_Real U = Axis.XDirection().
 386       AngleWithRef(D1,Axis.XDirection()^Axis.YDirection());
 387     U += Phase;
 388     if ( U < 0.) U += 2*M_PI;
 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
 402   (const TopoDS_Edge&          EdgeOnF1,
 403   const TopoDS_Edge&          EdgeOnF2,
 404   TColgp_SequenceOfPnt& Points   )
 405   const
 406 {
 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
 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(),
 427     "BRepFill_TrimSurfaceTool::IntersectWith: incoherent intersection");
 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);
 435     //    cout<<"BisPar "<<PSeq.X()<<endl;
 436   }
 437 }
 438
 439
 440 //=======================================================================
 441 //function : IsOnFace
 442 //purpose  :
 443 //=======================================================================
 444
 445 Standard_Boolean BRepFill_TrimSurfaceTool::IsOnFace
 446   (const gp_Pnt2d& Point) const
 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
 454   //  modified by NIZHNY-EAP Fri Jan 21 09:49:09 2000 ___BEGIN___
 455   Inter.Init(myFace1, Line,1e-6);//Precision::PConfusion());
 456   if (Inter.More()) return Standard_True;
 457
 458   // eval if is on face 2
 459   Inter.Init(myFace2, Line, 1e-6);//Precision::PConfusion());
 460   return Inter.More();
 461   //  modified by NIZHNY-EAP Fri Jan 21 09:49:14 2000 ___END___
 462 }
 463
 464
 465 //=======================================================================
 466 //function : ProjOn
 467 //purpose  :
 468 //=======================================================================
 469
 470 Standard_Real BRepFill_TrimSurfaceTool::ProjOn(const gp_Pnt2d& Point,
 471   const TopoDS_Edge& Edge) const
 472 {
 473   TopLoc_Location L;
 474   Standard_Real   f,l;
 475
 476
 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());
 481
 482   // projection of curves 3d in the plane xOy
 483   Handle(Geom_Plane) Plane = new Geom_Plane(0,0,1,0);
 484   Handle(Geom2d_Curve) C2d = GeomProjLib::Curve2d(CT,Plane);
 485
 486   // evaluate the projection of the point on the curve.
 487   Geom2dAPI_ProjectPointOnCurve Projector(Point, C2d);
 488 #ifdef OCCT_DEBUG
 489   Standard_Real Dist = Projector.LowerDistance();
 490   if ( Dist > Precision::Confusion() ) {
 491     cout << " *** WARNING  TrimSurfaceTool:  *** " << endl;
 492     cout << "      --> the point is not on the edge" <<endl;
 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
 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
 514 {
 515   Handle(Geom2d_TrimmedCurve) CT =
 516     new Geom2d_TrimmedCurve(myBis,U1,U2);
 517   BRepFill_MultiLine ML(myFace1,myFace2,
 518     myEdge1,myEdge2,myInv1,myInv2,CT);
 519
 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 }