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