0023367: New functionality restoring the middle path of pipe-like shape
[occt.git] / src / BRepOffsetAPI / BRepOffsetAPI_MiddlePath.cxx
1 // File:      BRepOffsetAPI_MiddlePath.cxx
2 // Created:   06.08.12 16:53:16
3 // Author:    jgv@ROLEX
4 // Copyright: Open CASCADE 2012
5
6 #include <BRepOffsetAPI_MiddlePath.ixx>
7 #include <BRepOffsetAPI_MiddlePath.hxx>
8
9 #include <ShapeUpgrade_UnifySameDomain.hxx>
10
11 #include <gp_Lin.hxx>
12 #include <Geom_Curve.hxx>
13 #include <Geom_TrimmedCurve.hxx>
14 #include <Geom_Line.hxx>
15 #include <Geom_BezierCurve.hxx>
16 #include <Geom_BSplineCurve.hxx>
17 #include <BRep_Tool.hxx>
18 #include <gce_MakeLin.hxx>
19
20 #include <BRepLib_MakeWire.hxx>
21
22 #include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
23 #include <TopTools_MapIteratorOfMapOfShape.hxx>
24 #include <TopExp.hxx>
25 #include <TopTools_ListIteratorOfListOfShape.hxx>
26 #include <TopoDS.hxx>
27 #include <BRepTools.hxx>
28 #include <TopTools_SequenceOfShape.hxx>
29 #include <BRepTools_WireExplorer.hxx>
30 #include <TopoDS_Iterator.hxx>
31 #include <BRep_Builder.hxx>
32 #include <Precision.hxx>
33 #include <TopExp_Explorer.hxx>
34 #include <BRepExtrema_DistShapeShape.hxx>
35 #include <Geom2d_Curve.hxx>
36 #include <Geom2d_Line.hxx>
37 #include <GCE2d_MakeLine.hxx>
38 #include <BRepLib_MakeEdge.hxx>
39 #include <BRepLib.hxx>
40 #include <GeomAbs_CurveType.hxx>
41 #include <BRepAdaptor_Curve.hxx>
42 #include <TopTools_Array1OfShape.hxx>
43 #include <BRepLib_MakeFace.hxx>
44 #include <TColgp_Array1OfPnt.hxx>
45 #include <TColgp_HArray1OfPnt.hxx>
46 #include <TColgp_Array1OfVec.hxx>
47 #include <TColStd_HArray1OfBoolean.hxx>
48 #include <GProp_GProps.hxx>
49 #include <BRepGProp.hxx>
50 #include <Geom_Circle.hxx>
51 #include <gp_Circ.hxx>
52 #include <GC_MakeCircle.hxx>
53 #include <TColgp_SequenceOfPnt.hxx>
54 #include <GeomLib.hxx>
55 #include <GeomAPI_Interpolate.hxx>
56
57 static Standard_Boolean IsClosed(const TopoDS_Wire& aWire)
58 {
59   TopoDS_Vertex V1, V2;
60   TopExp::Vertices(aWire, V1, V2);
61   return (V1.IsSame(V2));
62 }
63
64 static Standard_Boolean IsLinear(const TopoDS_Edge& anEdge,
65                                  gp_Lin& aLine)
66 {
67   Standard_Real fpar, lpar;
68   Handle(Geom_Curve) aCurve = BRep_Tool::Curve(anEdge, fpar, lpar);
69   if (aCurve->IsInstance(STANDARD_TYPE(Geom_TrimmedCurve)))
70     aCurve = ((Handle(Geom_TrimmedCurve)&) aCurve)->BasisCurve();
71
72   gp_Pnt Pnt1, Pnt2;
73   if (aCurve->IsKind(STANDARD_TYPE(Geom_Line)))
74   {
75     aLine = ((Handle(Geom_Line)&) aCurve)->Lin();
76     return Standard_True;
77   }
78   else if (aCurve->IsKind(STANDARD_TYPE(Geom_BezierCurve)))
79   {
80     Handle(Geom_BezierCurve) theBezier = (Handle(Geom_BezierCurve)&) aCurve;
81     if (theBezier->NbPoles() == 2)
82     {
83       Pnt1 = theBezier->Pole(1);
84       Pnt2 = theBezier->Pole(2);
85       aLine = gce_MakeLin(Pnt1, Pnt2);
86       return Standard_True;
87     }
88   }
89   else if (aCurve->IsKind(STANDARD_TYPE(Geom_BSplineCurve)))
90   {
91     Handle(Geom_BSplineCurve) theBSpline = (Handle(Geom_BSplineCurve)&) aCurve;
92     if (theBSpline->NbPoles() == 2)
93     {
94       Pnt1 = theBSpline->Pole(1);
95       Pnt2 = theBSpline->Pole(2);
96       aLine = gce_MakeLin(Pnt1, Pnt2);
97       return Standard_True;
98     }
99   }
100
101   return Standard_False;
102 }
103
104 static GeomAbs_CurveType TypeOfEdge(const TopoDS_Edge& anEdge)
105 {
106   gp_Lin aLin;
107   if (IsLinear(anEdge, aLin))
108     return GeomAbs_Line;
109
110   BRepAdaptor_Curve BAcurve(anEdge);
111   return BAcurve.GetType();
112 }
113
114 static gp_Vec TangentOfEdge(const TopoDS_Shape& aShape,
115                             const Standard_Boolean OnFirst)
116 {
117   TopoDS_Edge anEdge = TopoDS::Edge(aShape);
118   TopAbs_Orientation anOr = anEdge.Orientation();
119   
120   Standard_Real fpar, lpar;
121   Handle(Geom_Curve) aCurve = BRep_Tool::Curve(anEdge, fpar, lpar);
122   Standard_Real thePar;
123   if (OnFirst)
124     thePar = (anOr == TopAbs_FORWARD)? fpar : lpar;
125   else
126     thePar = (anOr == TopAbs_FORWARD)? lpar : fpar;
127
128   gp_Pnt thePoint;
129   gp_Vec theTangent;
130   aCurve->D1(thePar, thePoint, theTangent);
131   if (anOr == TopAbs_REVERSED)
132     theTangent.Reverse();
133
134   return theTangent;
135 }
136
137
138 static Standard_Boolean IsValidEdge(const TopoDS_Edge& theEdge,
139                                     const TopoDS_Face& theFace)
140 {
141   TopoDS_Vertex V1, V2;
142   TopExp::Vertices(theEdge, V1, V2);
143
144   Standard_Real Tol = Precision::Confusion();
145   Standard_Integer i;
146   
147   TopExp_Explorer Explo(theFace, TopAbs_EDGE);
148   for (; Explo.More(); Explo.Next())
149   {
150     const TopoDS_Shape& anEdge = Explo.Current();
151     BRepExtrema_DistShapeShape DistMini(theEdge, anEdge);
152     if (DistMini.Value() <= Tol)
153     {
154       for (i = 1; i <= DistMini.NbSolution(); i++)
155       {
156         BRepExtrema_SupportType theType = DistMini.SupportTypeShape2(i);
157         if (theType == BRepExtrema_IsOnEdge)
158           return Standard_False;
159         //theType is "IsVertex"
160         TopoDS_Shape aVertex = DistMini.SupportOnShape2(i);
161         if (!(aVertex.IsSame(V1) || aVertex.IsSame(V2)))
162           return Standard_False;
163       }
164     }
165   }
166
167   return Standard_True;
168 }
169
170 /*
171 //=======================================================================
172 //function : BRepOffsetAPI_MiddlePath
173 //purpose  : Constructor
174 //=======================================================================
175
176 BRepOffsetAPI_MiddlePath::BRepOffsetAPI_MiddlePath(const TopoDS_Shape& aShape,
177                                                    const TopoDS_Wire&  StartWire)
178 {
179   myInitialShape = aShape;
180   myStartWire    = StartWire;
181   myClosedSection = myStartWire.Closed();
182 }
183
184 //=======================================================================
185 //function : BRepOffsetAPI_MiddlePath
186 //purpose  : Constructor
187 //=======================================================================
188
189 BRepOffsetAPI_MiddlePath::BRepOffsetAPI_MiddlePath(const TopoDS_Shape& aShape,
190                                                    const TopoDS_Edge&  StartEdge)
191 {
192   myInitialShape = aShape;
193
194   BRepLib_MakeWire MW(StartEdge);
195   
196   //BB.Add(myStartWire, StartEdge);
197
198   TopTools_IndexedDataMapOfShapeListOfShape EFmap;
199   TopTools_IndexedDataMapOfShapeListOfShape VEmap;
200   TopExp::MapShapesAndAncestors(myInitialShape, TopAbs_EDGE,   TopAbs_FACE, EFmap);
201   TopExp::MapShapesAndAncestors(myInitialShape, TopAbs_VERTEX, TopAbs_EDGE, VEmap);
202   
203   //Standard_Boolean Start = Standard_True;
204   //if (Start)
205   //{
206   //TopExp::Vertices(CurEdge, V1, V2);
207   //  StartVertex = V1;
208   //  CurVertex   = V2;
209   //  if (VEmap(CurVertex).Extent() == 2) //end: two free edges
210   //  {
211   //    StartVertex = V2;
212   //    CurVertex   = V1;
213   //    if (VEmap(CurVertex).Extent() == 2) //end: two free edges
214   //      break;
215   //  }
216   //  Start = Standard_False;
217   //  continue;
218   //}
219
220   TopoDS_Vertex StartVertex, CurVertex, V1, V2;
221   TopExp::Vertices(StartEdge, StartVertex, CurVertex);
222   TopoDS_Edge CurEdge = StartEdge;
223   Standard_Integer i;
224   for (i = 1; i <= 2; i++)
225   {
226     for (;;)
227     {
228       const TopTools_ListOfShape& LE = VEmap.FindFromKey(CurVertex);
229       if (LE.Extent() == 2) //end: two free edges or one closed free edge
230         break;
231       TopTools_ListIteratorOfListOfShape itl(LE);
232       TopoDS_Edge anEdge;
233       for (; itl.More(); itl.Next())
234       {
235         anEdge = TopoDS::Edge(itl.Value());
236         if (anEdge.IsSame(CurEdge))
237           continue;
238         if (EFmap.FindFromKey(anEdge).Extent() == 1) //another free edge found
239           break;
240       }
241       //BB.Add(myStartWire, anEdge);
242       MW.Add(anEdge);
243       TopExp::Vertices(anEdge, V1, V2);
244       CurVertex = (V1.IsSame(CurVertex))? V2 : V1;
245       CurEdge = anEdge;
246       if (CurVertex.IsSame(StartVertex))
247         break;
248     }
249     if (CurVertex.IsSame(StartVertex))
250       break;
251     CurVertex = StartVertex;
252     CurEdge = StartEdge;
253   }
254   
255   myStartWire = MW.Wire();
256   myClosedSection = myStartWire.Closed();
257 }
258 */
259
260 //=======================================================================
261 //function : BRepOffsetAPI_MiddlePath
262 //purpose  : Constructor
263 //=======================================================================
264
265 BRepOffsetAPI_MiddlePath::BRepOffsetAPI_MiddlePath(const TopoDS_Shape& aShape,
266                                                    const TopoDS_Shape& StartShape,
267                                                    const TopoDS_Shape& EndShape)
268 {
269   ShapeUpgrade_UnifySameDomain Unifier(aShape);
270   Unifier.Build();
271   myInitialShape = Unifier.Shape();
272
273   TopoDS_Wire aStartWire, anEndWire;
274   if (StartShape.ShapeType() == TopAbs_FACE)
275   {
276     const TopoDS_Face& StartFace = TopoDS::Face(StartShape);
277     aStartWire = BRepTools::OuterWire(StartFace);
278   }
279   else
280     aStartWire = TopoDS::Wire(StartShape);
281   
282   if (EndShape.ShapeType() == TopAbs_FACE)
283   {
284     const TopoDS_Face& EndFace = TopoDS::Face(EndShape);
285     anEndWire = BRepTools::OuterWire(EndFace);
286   }
287   else
288     anEndWire = TopoDS::Wire(EndShape);
289
290   BRepLib_MakeWire MWstart;
291   //TopTools_MapOfShape MapEdges;
292   BRepTools_WireExplorer wexp(aStartWire);
293   for (; wexp.More(); wexp.Next())
294   {
295     TopoDS_Shape anEdge = wexp.Current();
296     TopoDS_Shape NewEdge = Unifier.Generated(anEdge);
297     if (!NewEdge.IsNull())
298       MWstart.Add(TopoDS::Edge(NewEdge));
299   }
300   myStartWire = MWstart.Wire();
301
302   BRepLib_MakeWire MWend;
303   //MapEdges.Clear();
304   for (wexp.Init(anEndWire); wexp.More(); wexp.Next())
305   {
306     TopoDS_Shape anEdge = wexp.Current();
307     TopoDS_Shape NewEdge = Unifier.Generated(anEdge);
308     if (!NewEdge.IsNull())
309       MWend.Add(TopoDS::Edge(NewEdge));
310   }
311   myEndWire = MWend.Wire();
312
313   myClosedSection = myStartWire.Closed();
314   myClosedRing    = myStartWire.IsSame(myEndWire);
315 }
316
317 //=======================================================================
318 //function : Build
319 //purpose  : 
320 //=======================================================================
321
322 void BRepOffsetAPI_MiddlePath::Build()
323 {
324   TopTools_ListIteratorOfListOfShape itl;
325   
326   TopTools_SequenceOfShape StartVertices;
327   TopTools_MapOfShape EndVertices;
328   TopTools_MapOfShape EndEdges;
329   BRepOffsetAPI_SequenceOfSequenceOfShape SectionsEdges;
330   
331   BRepTools_WireExplorer wexp(myStartWire);
332   TopTools_SequenceOfShape EdgeSeq;
333   for (; wexp.More(); wexp.Next())
334   {
335     StartVertices.Append(wexp.CurrentVertex());
336     EdgeSeq.Append(wexp.Current());
337   }
338   if (!myClosedSection)
339     StartVertices.Append(wexp.CurrentVertex());
340   SectionsEdges.Append(EdgeSeq);
341   
342   for (wexp.Init(myEndWire); wexp.More(); wexp.Next())
343   {
344     EndVertices.Add(wexp.CurrentVertex());
345     EndEdges.Add(wexp.Current());
346   }
347   if (!myClosedSection)
348     EndVertices.Add(wexp.CurrentVertex());
349   
350
351   TopoDS_Iterator itw(myStartWire);
352   for (; itw.More(); itw.Next())
353     myStartWireEdges.Add(itw.Value());
354   for (itw.Initialize(myEndWire); itw.More(); itw.Next())
355     myEndWireEdges.Add(itw.Value());
356
357   TopTools_IndexedDataMapOfShapeListOfShape VEmap;
358   TopTools_IndexedDataMapOfShapeListOfShape EFmap;
359   TopExp::MapShapesAndAncestors(myInitialShape, TopAbs_VERTEX, TopAbs_EDGE, VEmap);
360   TopExp::MapShapesAndAncestors(myInitialShape, TopAbs_EDGE,   TopAbs_FACE, EFmap);
361
362   TopTools_MapOfShape CurVertices;
363   
364   Standard_Integer i, j, k;
365   TopoDS_Edge anEdge;
366   TopoDS_Vertex V1, V2, NextVertex;
367   //Initialization of <myPaths>
368   for (i = 1; i <= StartVertices.Length(); i++)
369   {
370     TopTools_SequenceOfShape Edges;
371     const TopTools_ListOfShape& LE = VEmap.FindFromKey(StartVertices(i));
372     for (itl.Initialize(LE); itl.More(); itl.Next())
373     {
374       anEdge = TopoDS::Edge(itl.Value());
375       if (!myStartWireEdges.Contains(anEdge))
376       {
377         TopExp::Vertices(anEdge, V1, V2, Standard_True);
378         if (V1.IsSame(StartVertices(i)))
379           CurVertices.Add(V2);
380         else
381         {
382           anEdge.Reverse();
383           CurVertices.Add(V1);
384         }
385         Edges.Append(anEdge);
386         break;
387       }
388     }
389     myPaths.Append(Edges);
390   }
391
392   //Filling of "myPaths"
393   TopTools_ListOfShape NextVertices;
394   for (;;)
395   {
396     for (i = 1; i <= myPaths.Length(); i++)
397     {
398       const TopoDS_Shape& theShape = myPaths(i).Last();
399       TopoDS_Edge theEdge;
400       TopoDS_Vertex theVertex;
401       if (theShape.ShapeType() == TopAbs_EDGE)
402       {
403         theEdge = TopoDS::Edge(theShape);
404         theVertex = TopExp::LastVertex(theEdge, Standard_True);
405       }
406       else //last segment of path was punctual
407       {
408         theEdge = TopoDS::Edge(myPaths(i)(myPaths(i).Length()-1));
409         theVertex = TopoDS::Vertex(theShape);
410       }
411       
412       if (EndVertices.Contains(theVertex))
413         continue;
414       const TopTools_ListOfShape& LE = VEmap.FindFromKey(theVertex);
415       TopTools_MapOfShape NextEdgeCandidates;
416       for (itl.Initialize(LE); itl.More(); itl.Next())
417       {
418         anEdge = TopoDS::Edge(itl.Value());
419         if (anEdge.IsSame(theEdge))
420           continue;
421         TopExp::Vertices(anEdge, V1, V2, Standard_True);
422         if (V1.IsSame(theVertex))
423           NextVertex = V2;
424         else
425         {
426           anEdge.Reverse();
427           NextVertex = V1;
428         }
429         if (!CurVertices.Contains(NextVertex))
430           NextEdgeCandidates.Add(anEdge);
431       }
432       if (!NextEdgeCandidates.IsEmpty())
433       {
434         if (NextEdgeCandidates.Extent() > 1)
435           myPaths(i).Append(theVertex); //punctual segment of path
436         else
437         {
438           TopTools_MapIteratorOfMapOfShape mapit(NextEdgeCandidates);
439           anEdge = TopoDS::Edge(mapit.Key());
440           myPaths(i).Append(anEdge);
441           NextVertex = TopExp::LastVertex(anEdge, Standard_True);
442           NextVertices.Append(NextVertex);
443         }
444       }
445     }
446     if (NextVertices.IsEmpty())
447       break;
448     for (itl.Initialize(NextVertices); itl.More(); itl.Next())
449       CurVertices.Add(itl.Value());
450     NextVertices.Clear();
451   }
452
453   //Temporary
454   /*
455   TopoDS_Compound aCompound, aCmp1;
456   BRep_Builder BB;
457   BB.MakeCompound(aCompound);
458   BB.MakeCompound(aCmp1);
459   for (i = 1; i <= myPaths.Length(); i++)
460   {
461     TopoDS_Compound aCmp;
462     BB.MakeCompound(aCmp);
463     for (j = 1; j <= myPaths(i).Length(); j++)
464       BB.Add(aCmp, myPaths(i)(j));
465     BB.Add(aCmp1, aCmp);
466   }
467   BB.Add(aCompound, aCmp1);
468          
469   myShape = aCompound;
470
471   Done();
472   return;
473   */
474   ////////////
475   
476   //Building of set of sections
477   Standard_Integer NbE = EdgeSeq.Length();
478   Standard_Integer NbPaths = myPaths.Length();
479   Standard_Integer NbVer = myPaths.Length();
480   if (myClosedSection)
481     NbVer++;
482   i = 1;
483   for (;;)
484   {
485     for (j = 1; j <= EdgeSeq.Length(); j++)
486       EdgeSeq(j).Nullify();
487
488     Standard_Boolean ToInsertVertex = Standard_False;
489     
490     for (j = 2; j <= NbVer; j++)
491     {
492       if (!EdgeSeq(j-1).IsNull())
493         continue;
494
495       //for the end of initial shape
496       if (myPaths(j-1).Length() < i)
497       {
498         TopoDS_Edge aE1 = TopoDS::Edge(myPaths(j-1)(i-1));
499         TopoDS_Shape LastVer = TopExp::LastVertex(aE1, Standard_True);
500         myPaths(j-1).Append(LastVer);
501       }
502       if (myPaths((j<=NbPaths)? j : 1).Length() < i)
503       {
504         TopoDS_Edge aE2 = TopoDS::Edge(myPaths((j<=NbPaths)? j : 1)(i-1));
505         TopoDS_Shape LastVer = TopExp::LastVertex(aE2, Standard_True);
506         myPaths((j<=NbPaths)? j : 1).Append(LastVer);
507       }
508       //////////////////////////////
509       
510       if (ToInsertVertex)
511       {
512         if (myPaths(j-1)(i).ShapeType() == TopAbs_EDGE)
513         {
514           TopoDS_Edge aE1 = TopoDS::Edge(myPaths(j-1)(i));
515           TopoDS_Shape fver = TopExp::FirstVertex(aE1, Standard_True);
516           myPaths(j-1).InsertBefore(i, fver);
517         }
518         if (myPaths((j<=NbPaths)? j : 1)(i).ShapeType() == TopAbs_EDGE)
519         {
520           TopoDS_Edge aE2 = TopoDS::Edge(myPaths((j<=NbPaths)? j : 1)(i));
521           TopoDS_Shape fver = TopExp::FirstVertex(aE2, Standard_True);
522           myPaths((j<=NbPaths)? j : 1).InsertBefore(i, fver);
523         }
524         ToInsertVertex = Standard_False;
525       }
526       
527       TopoDS_Edge E1, E2;
528       if (myPaths(j-1)(i).ShapeType() == TopAbs_EDGE)
529         E1 = TopoDS::Edge(myPaths(j-1)(i));
530       if (myPaths((j<=NbPaths)? j : 1)(i).ShapeType() == TopAbs_EDGE)
531         E2 = TopoDS::Edge(myPaths((j<=NbPaths)? j : 1)(i));
532       TopoDS_Edge E12 = TopoDS::Edge(SectionsEdges(i)(j-1));
533       
534       //TopoDS_Vertex PrevVertex = TopoDS::Vertex(VerSeq(j-1));
535       //TopoDS_Vertex CurVertex  = TopoDS::Vertex(VerSeq(j));
536       TopoDS_Vertex PrevVertex = (E1.IsNull())? TopoDS::Vertex(myPaths(j-1)(i))
537         : TopExp::LastVertex(E1, Standard_True);
538       TopoDS_Vertex CurVertex = (E2.IsNull())? TopoDS::Vertex(myPaths((j<=NbPaths)? j : 1)(i))
539         : TopExp::LastVertex(E2, Standard_True);
540       
541       TopoDS_Edge ProperEdge;
542       const TopTools_ListOfShape& LE = VEmap.FindFromKey(PrevVertex);
543       //Temporary
544       Standard_Integer LenList = LE.Extent();
545       ///////////
546       for (itl.Initialize(LE); itl.More(); itl.Next())
547       {
548         anEdge = TopoDS::Edge(itl.Value());
549         TopExp::Vertices(anEdge, V1, V2);
550         if ((V1.IsSame(PrevVertex) && V2.IsSame(CurVertex) ||
551              V1.IsSame(CurVertex) && V2.IsSame(PrevVertex)) &&
552             !anEdge.IsSame(E1))
553         {
554           ProperEdge = anEdge;
555           break;
556         }
557       }
558       
559       if ((myPaths(j-1)(i)).ShapeType() == TopAbs_VERTEX &&
560           (myPaths((j<=NbPaths)? j : 1)(i)).ShapeType() == TopAbs_VERTEX)
561       {
562         EdgeSeq(j-1) = ProperEdge;
563         continue;
564       }
565
566       TopoDS_Vertex PrevPrevVer = (E1.IsNull())? PrevVertex
567         : TopExp::FirstVertex(E1, Standard_True);
568       TopoDS_Vertex PrevCurVer  = (E2.IsNull())? CurVertex
569         : TopExp::FirstVertex(E2, Standard_True);
570       if (ProperEdge.IsNull()) //no connection between these two vertices
571       {
572         //Find the face on which E1, E2 and E12 lie
573         //ToInsertVertex = Standard_False;
574         const TopoDS_Shape& EE1 = (E1.IsNull())?
575           myPaths(j-1)(i-1) : E1;
576         const TopoDS_Shape& EE2 = (E2.IsNull())?
577           myPaths((j<=NbPaths)? j : 1)(i-1) : E2;
578         const TopTools_ListOfShape& LF = EFmap.FindFromKey(EE1);
579         TopoDS_Face theFace;
580         for (itl.Initialize(LF); itl.More(); itl.Next())
581         {
582           const TopoDS_Shape& aFace = itl.Value();
583           TopExp_Explorer Explo(aFace, TopAbs_EDGE);
584           for (; Explo.More(); Explo.Next())
585           {
586             if (EE2.IsSame(Explo.Current()))
587             {
588               const TopTools_ListOfShape& LFsec = EFmap.FindFromKey(E12);
589               TopTools_ListIteratorOfListOfShape itlsec(LFsec);
590               for (; itlsec.More(); itlsec.Next())
591                 if (aFace.IsSame(itlsec.Value()))
592                 {
593                   theFace = TopoDS::Face(aFace);
594                   break;
595                 }
596               if (!theFace.IsNull())
597                 break;
598             }
599           }
600           if (!theFace.IsNull())
601             break;
602         }
603         TopTools_ListOfShape ListOneFace;
604         ListOneFace.Append(theFace);
605
606         if (E1.IsNull() || E2.IsNull())
607         {
608           if (E1.IsNull())
609             E1 = TopoDS::Edge(myPaths(j-1)(i-1));
610           if (E2.IsNull())
611             E2 = TopoDS::Edge(myPaths((j<=NbPaths)? j : 1)(i-1));
612           Standard_Real fpar1, lpar1, fpar2, lpar2;
613           Standard_Real FirstPar1, LastPar1, FirstPar2, LastPar2;
614           Handle(Geom2d_Curve) PCurve1 = BRep_Tool::CurveOnSurface(E1, theFace, fpar1, lpar1);
615           Handle(Geom2d_Curve) PCurve2 = BRep_Tool::CurveOnSurface(E2, theFace, fpar2, lpar2);
616           if (E1.Orientation() == TopAbs_FORWARD)
617           { FirstPar1 = fpar1; LastPar1 = lpar1; }
618           else
619           { FirstPar1 = lpar1; LastPar1 = fpar1; }
620           if (E2.Orientation() == TopAbs_FORWARD)
621           { FirstPar2 = fpar2; LastPar2 = lpar2; }
622           else
623           { FirstPar2 = lpar2; LastPar2 = fpar2; }
624           gp_Pnt2d FirstPnt2d = PCurve1->Value(LastPar1);
625           gp_Pnt2d LastPnt2d  = PCurve2->Value(LastPar2);
626           Handle(Geom_Surface) theSurf = BRep_Tool::Surface(theFace);
627           Handle(Geom2d_Line) theLine = GCE2d_MakeLine(FirstPnt2d, LastPnt2d);
628           Standard_Real len_ne = FirstPnt2d.Distance(LastPnt2d);
629           TopoDS_Edge NewEdge = BRepLib_MakeEdge(theLine, theSurf,
630                                                  PrevVertex, CurVertex,
631                                                  0., len_ne);
632           BRepLib::BuildCurve3d(NewEdge);
633           EdgeSeq(j-1) = NewEdge;
634           EFmap.Add(NewEdge, ListOneFace);
635         }
636         else //E1 is edge
637         {
638           //Extract points 2d
639           Standard_Real fpar1, lpar1, fpar2, lpar2;
640           Standard_Real FirstPar1, LastPar1, FirstPar2, LastPar2;
641           Handle(Geom2d_Curve) PCurve1 = BRep_Tool::CurveOnSurface(E1, theFace, fpar1, lpar1);
642           Handle(Geom2d_Curve) PCurve2 = BRep_Tool::CurveOnSurface(E2, theFace, fpar2, lpar2);
643           if (E1.Orientation() == TopAbs_FORWARD)
644           { FirstPar1 = fpar1; LastPar1 = lpar1; }
645           else
646           { FirstPar1 = lpar1; LastPar1 = fpar1; }
647           if (E2.Orientation() == TopAbs_FORWARD)
648           { FirstPar2 = fpar2; LastPar2 = lpar2; }
649           else
650           { FirstPar2 = lpar2; LastPar2 = fpar2; }
651           gp_Pnt2d FirstPnt2d = PCurve1->Value(LastPar1);
652           gp_Pnt2d LastPnt2d  = PCurve2->Value(LastPar2);
653           Handle(Geom_Surface) theSurf = BRep_Tool::Surface(theFace);
654           Handle(Geom2d_Line) theLine = GCE2d_MakeLine(FirstPnt2d, LastPnt2d);
655           Standard_Real len_ne = FirstPnt2d.Distance(LastPnt2d);
656           TopoDS_Edge NewEdge = BRepLib_MakeEdge(theLine, theSurf,
657                                                  PrevVertex, CurVertex,
658                                                  0., len_ne);
659           BRepLib::BuildCurve3d(NewEdge);
660           gp_Pnt2d PrevFirstPnt2d = PCurve1->Value(FirstPar1);
661           gp_Pnt2d PrevLastPnt2d  = PCurve2->Value(FirstPar2);
662           Handle(Geom2d_Line) Line1 = GCE2d_MakeLine(PrevFirstPnt2d, LastPnt2d);
663           Handle(Geom2d_Line) Line2 = GCE2d_MakeLine(FirstPnt2d, PrevLastPnt2d);
664           Standard_Real len_ne1 = PrevFirstPnt2d.Distance(LastPnt2d);
665           TopoDS_Edge NewEdge1 = BRepLib_MakeEdge(Line1, theSurf,
666                                                   PrevPrevVer, CurVertex,
667                                                   0., len_ne1);
668           BRepLib::BuildCurve3d(NewEdge1);
669           Standard_Real len_ne2 = FirstPnt2d.Distance(PrevLastPnt2d);
670           TopoDS_Edge NewEdge2 = BRepLib_MakeEdge(Line2, theSurf,
671                                                   PrevVertex, PrevCurVer,
672                                                   0., len_ne2);
673           BRepLib::BuildCurve3d(NewEdge2);
674           Standard_Boolean good_ne  = IsValidEdge(NewEdge, theFace);
675           Standard_Boolean good_ne1 = IsValidEdge(NewEdge1, theFace);
676           Standard_Boolean good_ne2 = IsValidEdge(NewEdge2, theFace);
677
678           GeomAbs_CurveType type_E1 = TypeOfEdge(E1);
679           GeomAbs_CurveType type_E2 = TypeOfEdge(E2);
680
681           Standard_Integer ChooseEdge = 0;
682           
683           if (!good_ne || type_E1 != type_E2)
684           {
685             if (type_E1 == type_E2) //!good_ne
686             {
687               if (good_ne1)
688                 ChooseEdge = 1;
689               else
690                 ChooseEdge = 2;
691             }
692             else //types are different
693             {
694               if (type_E1 == GeomAbs_Line)
695                 ChooseEdge = 1;
696               else if (type_E2 == GeomAbs_Line)
697                 ChooseEdge = 2;
698               else //to be developed later...
699               {}
700             }
701           }
702           
703           if (ChooseEdge == 0)
704           {
705             EdgeSeq(j-1) = NewEdge;
706             EFmap.Add(NewEdge, ListOneFace);
707           }
708           else if (ChooseEdge == 1)
709           {
710             EdgeSeq(j-1) = NewEdge1;
711             EFmap.Add(NewEdge1, ListOneFace);
712             for (k = 1; k < j-1; k++)
713               EdgeSeq(k).Nullify();
714             for (k = 1; k <= j-1; k++)
715             {
716               TopoDS_Edge aLastEdge = TopoDS::Edge(myPaths(k)(i));
717               TopoDS_Shape VertexAsEdge = TopExp::FirstVertex(aLastEdge, Standard_True);
718               myPaths(k).InsertBefore(i, VertexAsEdge);
719             }
720             j = 1; //start from beginning
721           }
722           else if (ChooseEdge == 2)
723           {
724             EdgeSeq(j-1) = NewEdge2;
725             EFmap.Add(NewEdge2, ListOneFace);
726             ToInsertVertex = Standard_True;
727           }
728         } //else //E1 is edge
729       } //if (ProperEdge.IsNull())
730       else //connecting edge exists
731       {
732         /*
733         if (ToInsertVertex)
734         {
735           myPaths(j-1).InsertBefore(i, PrevPrevVer);
736           myPaths((j<=NbPaths)? j : 1).InsertBefore(i, PrevCurVer);
737           EdgeSeq(j-1) = E12;
738         }
739         else
740         */
741           EdgeSeq(j-1) = ProperEdge;
742       }
743     } //for (j = 2; j <= NbVer; j++)
744     SectionsEdges.Append(EdgeSeq);
745
746     //check for exit from for(;;)
747     Standard_Integer NbEndEdges = 0;
748     for (j = 1; j <= EdgeSeq.Length(); j++)
749       if (EndEdges.Contains(EdgeSeq(j)))
750         NbEndEdges++;
751     if (NbEndEdges == NbE)
752       break;
753     
754     i++;
755   } //for (;;)
756   
757
758   //final phase: building of middle path
759   Standard_Integer NbSecFaces = SectionsEdges.Length();
760   TopTools_Array1OfShape SecFaces(1, NbSecFaces);
761   for (i = 1; i <= NbSecFaces; i++)
762   {
763     BRepLib_MakeWire MW;
764     for (j = 1; j <= NbE; j++)
765     {
766       anEdge = TopoDS::Edge(SectionsEdges(i)(j));
767       MW.Add(anEdge);
768     }
769     if (!myClosedSection)
770     {
771       TopExp::Vertices(MW.Wire(), V1, V2);
772       anEdge = BRepLib_MakeEdge(V2, V1);
773       MW.Add(anEdge);
774     }
775     TopoDS_Wire aWire = MW.Wire();
776     BRepLib_MakeFace MF(aWire, Standard_True); //Only plane
777     if (MF.IsDone())
778       SecFaces(i) = MF.Face();
779     else
780       SecFaces(i) = aWire;
781   }
782
783   TColgp_Array1OfPnt Centers(1, NbSecFaces);
784   for (i = 1; i <= NbSecFaces; i++)
785   {
786     GProp_GProps Properties;
787     if (SecFaces(i).ShapeType() == TopAbs_FACE)
788       BRepGProp::SurfaceProperties(SecFaces(i), Properties);
789     else //wire
790       BRepGProp::LinearProperties(SecFaces(i), Properties);
791       
792     Centers(i) = Properties.CentreOfMass();
793   }
794
795   TopTools_Array1OfShape MidEdges(1, NbSecFaces-1);
796   Standard_Real LinTol = 1.e-5;
797   Standard_Real AngTol = 1.e-7;
798   gp_Pnt Pnt1, Pnt2;
799   for (i = 1; i < NbSecFaces; i++)
800   {
801     GeomAbs_CurveType TypeOfMidEdge = GeomAbs_OtherCurve;
802     for (j = 1; j <= myPaths.Length(); j++)
803     {
804       const TopoDS_Shape& aShape = myPaths(j)(i);
805       if (aShape.ShapeType() == TopAbs_VERTEX)
806       {
807         TypeOfMidEdge = GeomAbs_OtherCurve;
808         break;
809       }
810       anEdge = TopoDS::Edge(aShape);
811       GeomAbs_CurveType aType = TypeOfEdge(anEdge);
812       if (j == 1)
813         TypeOfMidEdge = aType;
814       else
815       {
816         if (aType != TypeOfMidEdge)
817         {
818           TypeOfMidEdge = GeomAbs_OtherCurve;
819           break;
820         }
821       }
822     }
823     if (TypeOfMidEdge == GeomAbs_Line)
824       MidEdges(i) = BRepLib_MakeEdge(Centers(i), Centers(i+1));
825     else if (TypeOfMidEdge == GeomAbs_Circle)
826     {
827       gp_Ax1 theAxis;
828       gp_Dir theDir1, theDir2;
829       Standard_Real theAngle;
830       gp_Vec theTangent;
831       Standard_Boolean SimilarArcs = Standard_True;
832       for (j = 1; j <= myPaths.Length(); j++)
833       {
834         anEdge = TopoDS::Edge(myPaths(j)(i));
835         Standard_Real fpar, lpar;
836         Handle(Geom_Curve) aCurve = BRep_Tool::Curve(anEdge, fpar, lpar);
837         if (aCurve->IsInstance(STANDARD_TYPE(Geom_TrimmedCurve)))
838           aCurve = ((Handle(Geom_TrimmedCurve)&) aCurve)->BasisCurve();
839         Pnt1 = aCurve->Value(fpar);
840         Pnt2 = aCurve->Value(lpar);
841         Handle(Geom_Circle) aCircle = Handle(Geom_Circle)::DownCast(aCurve);
842         gp_Circ aCirc = aCircle->Circ();
843         if (j == 1)
844         {
845           theAxis = aCirc.Axis();
846           theDir1 = gp_Vec(aCirc.Location(), Pnt1);
847           theDir2 = gp_Vec(aCirc.Location(), Pnt2);
848           theAngle = lpar - fpar;
849           Standard_Real theParam = (anEdge.Orientation() == TopAbs_FORWARD)?
850             fpar : lpar;
851           aCurve->D1(theParam, Pnt1, theTangent);
852           if (anEdge.Orientation() == TopAbs_REVERSED)
853             theTangent.Reverse();
854         }
855         else
856         {
857           gp_Ax1 anAxis = aCirc.Axis();
858           gp_Lin aLin(anAxis);
859           if (!aLin.Contains(theAxis.Location(), LinTol) ||
860               !anAxis.IsParallel(theAxis, AngTol))
861           {
862             SimilarArcs = Standard_False;
863             break;
864           }
865           gp_Dir aDir1 = gp_Vec(aCirc.Location(), Pnt1);
866           gp_Dir aDir2 = gp_Vec(aCirc.Location(), Pnt2);
867           if (!(aDir1.IsEqual(theDir1, AngTol) && aDir2.IsEqual(theDir2, AngTol) ||
868                 aDir1.IsEqual(theDir2, AngTol) && aDir2.IsEqual(theDir1, AngTol)))
869           {
870             SimilarArcs = Standard_False;
871             break;
872           }
873         }
874       }
875       if (SimilarArcs)
876       {
877         gp_XYZ AxisLoc = theAxis.Location().XYZ();
878         gp_XYZ AxisDir = theAxis.Direction().XYZ();
879         Standard_Real Parameter = (Centers(i).XYZ() - AxisLoc) * AxisDir;
880         gp_Pnt theCenterOfCirc(AxisLoc + Parameter*AxisDir);
881         
882         gp_Vec Vec1(theCenterOfCirc, Centers(i));
883         gp_Vec Vec2(theCenterOfCirc, Centers(i+1));
884         /*
885         gp_Dir VecProd = Vec1 ^ Vec2;
886         if (theAxis.Direction() * VecProd < 0.)
887           theAxis.Reverse();
888         */
889         
890         Standard_Real anAngle = Vec1.AngleWithRef(Vec2, theAxis.Direction());
891         if (anAngle < 0.)
892           anAngle += 2.*M_PI;
893         if (Abs(anAngle - theAngle) > AngTol)
894           theAxis.Reverse();
895         gp_Ax2 theAx2(theCenterOfCirc, theAxis.Direction(), Vec1);
896         Handle(Geom_Circle) theCircle = GC_MakeCircle(theAx2, Vec1.Magnitude());
897         gp_Vec aTangent;
898         theCircle->D1( 0., Pnt1, aTangent );
899         if (aTangent * theTangent < 0.)
900         {
901           theAxis.Reverse();
902           theAx2 = gp_Ax2(theCenterOfCirc, theAxis.Direction(), Vec1);
903           theCircle = GC_MakeCircle(theAx2, Vec1.Magnitude());
904         }
905         MidEdges(i) = BRepLib_MakeEdge(theCircle, 0., theAngle);
906       }
907     }
908   }
909
910   //Build missed edges
911   for (i = 1; i < NbSecFaces; i++)
912   {
913     if (MidEdges(i).IsNull())
914     {
915       for (j = i+1; j < NbSecFaces; j++)
916       {
917         if (!MidEdges(j).IsNull())
918           break;
919       }
920       //from i to j-1 all edges are null
921       Handle(TColgp_HArray1OfPnt) thePoints = new TColgp_HArray1OfPnt(1, j-i+1);
922       TColgp_Array1OfVec theTangents(1, j-i+1);
923       Handle(TColStd_HArray1OfBoolean) theFlags = new TColStd_HArray1OfBoolean(1, j-i+1);
924       for (k = i; k <= j; k++)
925         thePoints->SetValue(k-i+1, Centers(k));
926       for (k = i; k <= j; k++)
927       {
928         TColgp_SequenceOfPnt PntSeq;
929         for (Standard_Integer indp = 1; indp <= myPaths.Length(); indp++)
930         {
931           gp_Vec aTangent;
932           if (k == i)
933           {
934             if (myPaths(indp)(k).ShapeType() == TopAbs_VERTEX)
935               continue;
936             aTangent = TangentOfEdge(myPaths(indp)(k), Standard_True); //at begin
937           }
938           else if (k == j)
939           {
940             if (myPaths(indp)(k-1).ShapeType() == TopAbs_VERTEX)
941               continue;
942             aTangent = TangentOfEdge(myPaths(indp)(k-1), Standard_False); //at end
943           }
944           else
945           {
946             if (myPaths(indp)(k-1).ShapeType() == TopAbs_VERTEX ||
947                 myPaths(indp)(k).ShapeType() == TopAbs_VERTEX)
948               continue;
949             gp_Vec Tangent1 = TangentOfEdge(myPaths(indp)(k-1), Standard_False); //at end
950             gp_Vec Tangent2 = TangentOfEdge(myPaths(indp)(k), Standard_True); //at begin
951             aTangent = Tangent1 + Tangent2;
952           }
953           aTangent.Normalize();
954           gp_Pnt aPnt(aTangent.XYZ());
955           PntSeq.Append(aPnt);
956         }
957         TColgp_Array1OfPnt PntArray(1, PntSeq.Length());
958         for (Standard_Integer ip = 1; ip <= PntSeq.Length(); ip++)
959           PntArray(ip) = PntSeq(ip);
960         gp_Pnt theBary;
961         gp_Dir xdir, ydir;
962         Standard_Real xgap, ygap, zgap;
963         GeomLib::Inertia(PntArray, theBary, xdir, ydir, xgap, ygap, zgap);
964         gp_Vec theTangent(theBary.XYZ());
965         theTangents(k-i+1) = theTangent;
966       }
967       theFlags->Init(Standard_True);
968
969       GeomAPI_Interpolate Interpol(thePoints, Standard_False, LinTol);
970       Interpol.Load(theTangents, theFlags);
971       Interpol.Perform();
972       if (!Interpol.IsDone())
973       {
974         cout<<endl<<"Interpolation failed"<<endl;
975       }
976       Handle(Geom_Curve) InterCurve = Interpol.Curve();
977       MidEdges(i) = BRepLib_MakeEdge(InterCurve);
978       i = j;
979     }
980   }  
981
982   BRepLib_MakeWire MakeFinalWire;
983   for (i = 1; i < NbSecFaces; i++)
984     if (!MidEdges(i).IsNull())
985       MakeFinalWire.Add(TopoDS::Edge(MidEdges(i)));
986
987   TopoDS_Wire FinalWire = MakeFinalWire.Wire();
988   myShape = MakeFinalWire.Wire();
989   
990   //Temporary
991   /*
992   TopoDS_Compound aCompound, aCmp1, aCmp2;
993   BRep_Builder BB;
994   BB.MakeCompound(aCompound);
995   BB.MakeCompound(aCmp1);
996   BB.MakeCompound(aCmp2);
997   for (i = 1; i <= myPaths.Length(); i++)
998   {
999     TopoDS_Compound aCmp;
1000     BB.MakeCompound(aCmp);
1001     for (j = 1; j <= myPaths(i).Length(); j++)
1002       BB.Add(aCmp, myPaths(i)(j));
1003     BB.Add(aCmp1, aCmp);
1004   }
1005   for (i = 1; i <= SectionsEdges.Length(); i++)
1006   {
1007     TopoDS_Wire aWire;
1008     BB.MakeWire(aWire);
1009     for (j = 1; j <= SectionsEdges(i).Length(); j++)
1010       BB.Add(aWire, SectionsEdges(i)(j));
1011     BB.Add(aCmp2, aWire);
1012   }
1013   BB.Add(aCompound, aCmp1);
1014   BB.Add(aCompound, aCmp2);
1015   BB.Add(aCompound, FinalWire);
1016          
1017   myShape = aCompound;
1018   */
1019   ////////////
1020
1021   Done();
1022 }