0024624: Lost word in license statement in source files
[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 #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 }