dd62446d6b807d983524d177aaf076bdf74daac0
[occt.git] / 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 }