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   Standard_Real V = 0.;
 168   BRep_Tool::UVPoints(Edge,Face,PE1,PE2);
 169   VDeg = PE1.Y();
 170   TopExp_Explorer Exp(Face,TopAbs_EDGE);
 171   for (; Exp.More(); Exp.Next()) {
 172     if ( !TopoDS::Edge(Exp.Current()).IsSame(Edge)) {
 173       BRep_Tool::UVPoints(TopoDS::Edge(Exp.Current()),Face,PF1,PF2);
 174       V = ( Abs(PF1.Y() - VDeg) > Abs(PF2.Y() - VDeg)) ? PF1.Y() : PF2.Y();
 175       break;
 176     }
 177   }
 178   gp_Pnt P = GAS.Value(0., V);
 179
 180   if ( gp_Vec(Axis.Location(), P).Dot(Axis.XDirection()) < 0.)
 181     return M_PI;
 182   else
 183     return 0.;
 184 }
 185
 186
 187 //=======================================================================
 188 //function : EvalParameters
 189 //purpose  :
 190 //=======================================================================
 191
 192 static void EvalParameters(const TopoDS_Edge&          Edge,
 193                            const TopoDS_Face&          Face,
 194                            const Handle(Geom2d_Curve)& Bis ,
 195                                  TColgp_SequenceOfPnt& Seq  )
 196 {
 197   Standard_Boolean Degener = BRep_Tool::Degenerated(Edge);
 198   // return curves 3d associated to edges.
 199   TopLoc_Location L;
 200   Standard_Real   f,l;
 201
 202   Handle(Geom_TrimmedCurve) CT;
 203   Handle(Geom_Plane) Plane = new Geom_Plane(0,0,1,0);
 204
 205   Geom2dInt_GInter Intersector;
 206
 207   Standard_Integer NbPoints, NbSegments;
 208   Standard_Real U1, U2;
 209   gp_Pnt  P;//,PSeq;
 210
 211 //  Standard_Real Tol = Precision::Intersection();
 212   //  modified by NIZHNY-EAP Wed Dec 22 15:00:51 1999 ___BEGIN___
 213   Standard_Real Tol = 1.e-6;  // BRepFill_Precision();
 214   Standard_Real TolC = 0.;
 215
 216   if ( !Degener) {
 217     Handle(Geom_Curve) C = BRep_Tool::Curve(Edge,L,f,l);
 218     CT = new Geom_TrimmedCurve(C,f,l);
 219     CT->Transform(L.Transformation());
 220     // projection of 3d curves in the plane xOy
 221     Handle(Geom2d_Curve) C2d = GeomProjLib::Curve2d(CT,Plane);
 222
 223     Geom2dAdaptor_Curve AC(C2d);
 224     Geom2dAdaptor_Curve ABis(Bis);
 225
 226     Intersector = Geom2dInt_GInter(ABis, AC, TolC, Tol);
 227
 228     if ( !Intersector.IsDone()) {
 229       StdFail_NotDone::Raise("BRepFill_TrimSurfaceTool::IntersectWith");
 230     }
 231
 232     NbPoints = Intersector.NbPoints();
 233
 234     if (NbPoints < 1) {
 235       // try to elongate curves and enlarge tolerance
 236
 237       // don't do it rightaway from the beginning in order not to get
 238       // extra solutions those would cause *Exception*: incoherent intersection
 239
 240       GeomAbs_CurveType CType = AC.GetType(), BisType = ABis.GetType();
 241       Standard_Boolean  canElongateC = !(CType == GeomAbs_BezierCurve ||
 242                                          CType == GeomAbs_BSplineCurve ||
 243                                          CType == GeomAbs_OtherCurve);
 244       Standard_Boolean canElongateBis = !(BisType == GeomAbs_BezierCurve ||
 245                                           BisType == GeomAbs_BSplineCurve ||
 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 DEB
 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   Points.Clear();
 408   TColgp_SequenceOfPnt Points2;
 409
 410   EvalParameters(EdgeOnF1, myFace1, myBis, Points);
 411   EvalParameters(EdgeOnF2, myFace2, myBis, Points2);
 412
 413   StdFail_NotDone_Raise_if
 414     ( Points.Length() != Points2.Length(),
 415      "BRepFill_TrimSurfaceTool::IntersectWith: incoherent intersection");
 416
 417   gp_Pnt PSeq;
 418   Standard_Integer NbPoints = Points.Length();
 419   for ( Standard_Integer i = 1; i <= NbPoints; i++) {
 420     PSeq = Points(i);
 421     PSeq.SetZ((Points2.Value(i)).Y());
 422     Points.SetValue(i,PSeq);
 423 //    cout<<"BisPar "<<PSeq.X()<<endl;
 424   }
 425 }
 426
 427
 428 //=======================================================================
 429 //function : IsOnFace
 430 //purpose  :
 431 //=======================================================================
 432
 433 Standard_Boolean BRepFill_TrimSurfaceTool::IsOnFace
 434 (const gp_Pnt2d& Point) const
 435 {
 436   gp_Pnt P( Point.X(), Point.Y(), 0.);
 437   gp_Lin Line( P, gp::DZ());
 438
 439   BRepIntCurveSurface_Inter Inter;
 440
 441   // eval if is on face 1
 442 //  modified by NIZHNY-EAP Fri Jan 21 09:49:09 2000 ___BEGIN___
 443   Inter.Init(myFace1, Line, 1e-6);//Precision::PConfusion());
 444   if (Inter.More()) return Standard_True;
 445
 446   // eval if is on face 2
 447   Inter.Init(myFace2, Line, 1e-6);//Precision::PConfusion());
 448   return Inter.More();
 449 //  modified by NIZHNY-EAP Fri Jan 21 09:49:14 2000 ___END___
 450 }
 451
 452
 453 //=======================================================================
 454 //function : ProjOn
 455 //purpose  :
 456 //=======================================================================
 457
 458 Standard_Real BRepFill_TrimSurfaceTool::ProjOn(const gp_Pnt2d& Point,
 459                                                const TopoDS_Edge& Edge) const
 460 {
 461   TopLoc_Location L;
 462   Standard_Real   f,l;
 463
 464
 465
 466   Handle(Geom_Curve) C1 = BRep_Tool::Curve(Edge,L,f,l);
 467   Handle(Geom_TrimmedCurve) CT = new Geom_TrimmedCurve(C1,f,l);
 468   CT->Transform(L.Transformation());
 469
 470   // projection of curves 3d in the plane xOy
 471   Handle(Geom_Plane) Plane = new Geom_Plane(0,0,1,0);
 472   Handle(Geom2d_Curve) C2d = GeomProjLib::Curve2d(CT,Plane);
 473
 474   // evaluate the projection of the point on the curve.
 475   Geom2dAPI_ProjectPointOnCurve Projector(Point, C2d);
 476 #ifdef DEB
 477   Standard_Real Dist =
 478 #endif
 479     Projector.LowerDistance();
 480 #ifdef DEB
 481   if ( Dist > Precision::Confusion() ) {
 482     cout << " *** WARNING  TrimSurfaceTool:  *** " << endl;
 483     cout << "      --> the point is not on the edge" <<endl;
 484     cout << "          distance  = " << Dist << endl;
 485   }
 486 #endif
 487
 488   Standard_Real U = Projector.LowerDistanceParameter();
 489   return U;
 490 }
 491
 492
 493 //=======================================================================
 494 //function : Project
 495 //purpose  :
 496 //=======================================================================
 497
 498 void BRepFill_TrimSurfaceTool::Project
 499 (const Standard_Real U1,
 500  const Standard_Real U2,
 501        Handle(Geom_Curve)&   Curve,
 502        Handle(Geom2d_Curve)& PCurve1,
 503        Handle(Geom2d_Curve)& PCurve2,
 504        GeomAbs_Shape&        Cont) const
 505 {
 506   Standard_Integer Deg1, Deg2;
 507   Deg1 = Deg2 = 8;
 508
 509   Handle(Geom2d_TrimmedCurve) CT =
 510     new Geom2d_TrimmedCurve(myBis,U1,U2);
 511   BRepFill_MultiLine ML(myFace1,myFace2,
 512                         myEdge1,myEdge2,myInv1,myInv2,CT);
 513
 514   Cont = ML.Continuity();
 515
 516   if ( ML.IsParticularCase()) {
 517     ML.Curves(Curve,PCurve1,PCurve2);
 518   }
 519   else {
 520     BRepFill_ApproxSeewing AppSeew(ML);
 521
 522     Curve   = AppSeew.Curve();
 523     PCurve1 = AppSeew.CurveOnF1();
 524     PCurve2 = AppSeew.CurveOnF2();
 525   }
 526 }