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