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