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