0024899: Time of computation of intersection points with help of class BRepIntCurveSu...
[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 #include <TopExp_Explorer.hxx>
63
64 #ifdef DRAW
65 #include <DrawTrSurf.hxx>
66 #include <DBRep.hxx>
67 #endif
68 #ifdef DEB
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 DEB
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 DEB
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 DEB  
472   Standard_Real Dist = 
473 #endif
474     Projector.LowerDistance();
475 #ifdef DEB
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 }