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