0026252: GeomAdaptor_Surface should use inner adaptor to calculate values of complex...
[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 <GeomAdaptor_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 //#define DRAW
67 #include <stdio.h>
68 #ifdef DRAW
69 #include <DrawTrSurf.hxx>
70 #include <DBRep.hxx>
71 #endif
72 #ifdef OCCT_DEBUG
73 static Standard_Boolean Affich       = Standard_False;
74 static Standard_Integer NBCALL  = 0;
75 #endif
76
77 //=======================================================================
78 //function : BRepFill_TrimSurfaceTool
79 //purpose  : Initialisation with two neighbor faces
80 //          Edge1 and Edge2 are parallel edges corresponding
81 //          to minimum iso on F1 and F2 respectively.
82 //          ie Edge1 is Umin or VMin on F1.
83 //          Inv1 and Inv2 show if Edge1 and Edge2 are
84 //          returned parallel.   
85 //=======================================================================
86
87 BRepFill_TrimSurfaceTool::BRepFill_TrimSurfaceTool
88   (const Handle(Geom2d_Curve)& Bis, 
89   const TopoDS_Face&          Face1, 
90   const TopoDS_Face&          Face2,
91   const TopoDS_Edge&          Edge1,
92   const TopoDS_Edge&          Edge2,
93   const Standard_Boolean      Inv1,
94   const Standard_Boolean      Inv2 ) :
95 myFace1(Face1),
96   myFace2(Face2),
97   myEdge1(Edge1),
98   myEdge2(Edge2),
99   myInv1(Inv1),
100   myInv2(Inv2),
101   myBis  (Bis)
102 {
103 #ifdef OCCT_DEBUG
104   if ( Affich) {
105     NBCALL++;
106     cout << " ---------->TrimSurfaceTool : NBCALL = " << NBCALL << endl;
107 #ifdef DRAW
108     char name[256];
109
110     sprintf(name,"FACE1_%d",NBCALL);
111     DBRep::Set(name,myFace1);
112
113     sprintf(name,"FACE2_%d",NBCALL);
114     DBRep::Set(name,myFace2);
115
116     sprintf(name,"EDGE1_%d",NBCALL);
117     DBRep::Set(name,myEdge1);
118
119     sprintf(name,"EDGE2_%d",NBCALL);
120     DBRep::Set(name,myEdge2);
121
122     sprintf(name,"BISSEC_%d",NBCALL);
123     DrawTrSurf::Set(name,myBis);
124 #endif    
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_OffsetCurve  ||
242                                          CType == GeomAbs_OtherCurve);
243       Standard_Boolean canElongateBis = !(BisType == GeomAbs_BezierCurve  ||
244                                           BisType == GeomAbs_BSplineCurve ||
245                                           BisType == GeomAbs_OffsetCurve  ||
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 OCCT_DEBUG
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       GeomAdaptor_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 #ifdef DRAW
408   if ( Affich) {
409     char name[256];  
410     Standard_Integer i1 = 0, i2 = 2;
411     sprintf(name,"EdgeOnF1_%d_%d",i1, NBCALL);
412     DBRep::Set(name,EdgeOnF1);
413
414     sprintf(name,"EdgeOnF2_%d_%d",i2, NBCALL);
415     DBRep::Set(name,EdgeOnF2);
416   }
417
418 #endif 
419   Points.Clear();
420   TColgp_SequenceOfPnt Points2;
421
422   EvalParameters(EdgeOnF1, myFace1, myBis, Points);
423   EvalParameters(EdgeOnF2, myFace2, myBis, Points2);
424
425   StdFail_NotDone_Raise_if
426     ( Points.Length() != Points2.Length(),
427     "BRepFill_TrimSurfaceTool::IntersectWith: incoherent intersection");
428
429   gp_Pnt PSeq;
430   Standard_Integer NbPoints = Points.Length();
431   for ( Standard_Integer i = 1; i <= NbPoints; i++) {
432     PSeq = Points(i);
433     PSeq.SetZ((Points2.Value(i)).Y());
434     Points.SetValue(i,PSeq);
435     //    cout<<"BisPar "<<PSeq.X()<<endl;
436   }
437 }
438
439
440 //=======================================================================
441 //function : IsOnFace
442 //purpose  : 
443 //=======================================================================
444
445 Standard_Boolean BRepFill_TrimSurfaceTool::IsOnFace
446   (const gp_Pnt2d& Point) const 
447 {
448   gp_Pnt P( Point.X(), Point.Y(), 0.);
449   gp_Lin Line( P, gp::DZ());
450
451   BRepIntCurveSurface_Inter Inter;
452
453   // eval if is on face 1
454   //  modified by NIZHNY-EAP Fri Jan 21 09:49:09 2000 ___BEGIN___
455   Inter.Init(myFace1, Line,1e-6);//Precision::PConfusion());
456   if (Inter.More()) return Standard_True;
457
458   // eval if is on face 2
459   Inter.Init(myFace2, Line, 1e-6);//Precision::PConfusion());
460   return Inter.More();
461   //  modified by NIZHNY-EAP Fri Jan 21 09:49:14 2000 ___END___
462 }
463
464
465 //=======================================================================
466 //function : ProjOn
467 //purpose  : 
468 //=======================================================================
469
470 Standard_Real BRepFill_TrimSurfaceTool::ProjOn(const gp_Pnt2d& Point,
471   const TopoDS_Edge& Edge) const 
472 {
473   TopLoc_Location L;
474   Standard_Real   f,l;
475
476
477
478   Handle(Geom_Curve) C1 = BRep_Tool::Curve(Edge,L,f,l);
479   Handle(Geom_TrimmedCurve) CT = new Geom_TrimmedCurve(C1,f,l);
480   CT->Transform(L.Transformation());
481
482   // projection of curves 3d in the plane xOy
483   Handle(Geom_Plane) Plane = new Geom_Plane(0,0,1,0);
484   Handle(Geom2d_Curve) C2d = GeomProjLib::Curve2d(CT,Plane);
485
486   // evaluate the projection of the point on the curve.
487   Geom2dAPI_ProjectPointOnCurve Projector(Point, C2d);
488 #ifdef OCCT_DEBUG
489   Standard_Real Dist = Projector.LowerDistance();
490   if ( Dist > Precision::Confusion() ) {
491     cout << " *** WARNING  TrimSurfaceTool:  *** " << endl;
492     cout << "      --> the point is not on the edge" <<endl;
493     cout << "          distance  = " << Dist << endl;
494   }
495 #endif
496
497   Standard_Real U = Projector.LowerDistanceParameter();
498   return U;
499 }
500
501
502 //=======================================================================
503 //function : Project
504 //purpose  : 
505 //=======================================================================
506
507 void BRepFill_TrimSurfaceTool::Project
508   (const Standard_Real U1,
509   const Standard_Real U2,
510   Handle(Geom_Curve)&   Curve,
511   Handle(Geom2d_Curve)& PCurve1,
512   Handle(Geom2d_Curve)& PCurve2,
513   GeomAbs_Shape&        Cont) const 
514 {
515   Handle(Geom2d_TrimmedCurve) CT = 
516     new Geom2d_TrimmedCurve(myBis,U1,U2);
517   BRepFill_MultiLine ML(myFace1,myFace2,
518     myEdge1,myEdge2,myInv1,myInv2,CT);
519
520   Cont = ML.Continuity();
521
522   if ( ML.IsParticularCase()) {
523     ML.Curves(Curve,PCurve1,PCurve2);
524   }
525   else {
526     BRepFill_ApproxSeewing AppSeew(ML);
527
528     Curve   = AppSeew.Curve();
529     PCurve1 = AppSeew.CurveOnF1();
530     PCurve2 = AppSeew.CurveOnF2();
531   }
532 }