Salome HOME
Merge branch 'Dev_0.6.1' of newgeom:newgeom into Dev_0.6.1
[modules/shaper.git] / src / GeomAlgoAPI / GeomAlgoAPI_SketchBuilder.cpp
1 // Copyright (C) 2014-20xx CEA/DEN, EDF R&D
2
3 // File:        GeomAlgoAPI_SketchBuilder.cpp
4 // Created:     02 Jun 2014
5 // Author:      Artem ZHIDKOV
6
7 #include <GeomAlgoAPI_SketchBuilder.h>
8 #include <GeomAPI_PlanarEdges.h>
9
10 #include <set>
11
12 #include <gp_Dir.hxx>
13 #include <gp_Pln.hxx>
14
15 #include <BOPAlgo_Builder.hxx>
16 #include <BOPAlgo_Operation.hxx>
17 #include <BOPAlgo_PaveFiller.hxx>
18 #include <BOPTools.hxx>
19 #include <BRep_Builder.hxx>
20 #include <BRep_Curve3D.hxx>
21 #include <BRep_Tool.hxx>
22 #include <BRepAlgoAPI_Cut.hxx>
23 #include <BRepBuilderAPI_MakeFace.hxx>
24 #include <BRepClass_FaceClassifier.hxx>
25 #include <Geom_Curve.hxx>
26 #include <Geom_Plane.hxx>
27 #include <TopExp.hxx>
28 #include <TopExp_Explorer.hxx>
29 #include <TopoDS_Edge.hxx>
30 #include <TopoDS_Face.hxx>
31 #include <TopoDS_Vertex.hxx>
32 #include <TopoDS_Wire.hxx>
33
34 #include <Precision.hxx>
35
36 #include <assert.h>
37
38 #ifndef DBL_MAX
39 #define DBL_MAX 1.7976931348623158e+308 
40 #endif
41
42 const double tolerance = Precision::Confusion();
43 // This value helps to find direction on the boundaries of curve.
44 // It is not significant for lines, but is used for circles to avoid 
45 // wrong directions of movement (when two edges are tangent on the certain vertex)
46 const double shift = acos(1.0 - 3.0 * tolerance);
47
48 /// \brief Search first vertex - the vertex with lowest x coordinate, which is used in 2 edges at least
49 static const TopoDS_Vertex& findStartVertex(const BOPCol_IndexedDataMapOfShapeListOfShape& theMapVE,
50                                            const gp_Dir& theDirX, const gp_Dir& theDirY);
51
52 /// \brief Search the vertex on the sketch candidate to be the next one in the loop
53 static void findNextVertex(const TopoDS_Vertex& theStartVertex,
54                            const BOPCol_IndexedDataMapOfShapeListOfShape& theVertexEdgeMap,
55                            const gp_Dir& theStartDir, const gp_Dir& theNormal,
56                            TopoDS_Vertex& theNextVertex, TopoDS_Edge& theNextEdge,
57                            gp_Dir& theNextDir);
58
59 /// \brief Create planar face using the edges surrounding it
60 static void createFace(const TopoDS_Vertex& theStartVertex,
61                        const std::list<TopoDS_Edge>::iterator& theStartEdge,
62                        const std::list<TopoDS_Edge>::iterator& theEndOfEdges,
63                        const gp_Pln& thePlane, TopoDS_Face& theResFace);
64
65 /// \bief Create planar wire
66 static void createWireList(const TopoDS_Vertex& theStartVertex,
67                            const std::list<TopoDS_Edge>::iterator& theStartEdge,
68                            const std::list<TopoDS_Edge>::iterator& theEndOfEdges,
69                            const std::set<TopoDS_Edge*>& theEdgesInLoops,
70                            std::list<TopoDS_Wire>& theResWires);
71
72 /// \brief Calculate outer tengency on the edge in specified vertex
73 static gp_Dir getOuterEdgeDirection(const TopoDS_Edge& theEdge, const TopoDS_Vertex& theVertex);
74
75 /// \brief Unnecessary edges will be removed from the map.
76 ///        Positions of iterator will be updated
77 static void removeWasteEdges(std::list<TopoDS_Vertex>::iterator& theStartVertex,
78                              std::list<TopoDS_Edge>::iterator& theStartEdge,
79                              const std::list<TopoDS_Vertex>::iterator& theEndOfVertexes,
80                              const std::list<TopoDS_Edge>::iterator& theEndOfEdges,
81                              BOPCol_IndexedDataMapOfShapeListOfShape& theMapVE);
82
83
84 void GeomAlgoAPI_SketchBuilder::createFaces(
85     const std::shared_ptr<GeomAPI_Pnt>& theOrigin, const std::shared_ptr<GeomAPI_Dir>& theDirX,
86     const std::shared_ptr<GeomAPI_Dir>& theDirY, const std::shared_ptr<GeomAPI_Dir>& theNorm,
87     const std::list<std::shared_ptr<GeomAPI_Shape> >& theFeatures,
88     std::list<std::shared_ptr<GeomAPI_Shape> >& theResultFaces,
89     std::list<std::shared_ptr<GeomAPI_Shape> >& theResultWires)
90 {
91   if (theFeatures.empty())
92     return;
93
94   // Create the list of edges with shared vertexes
95   BOPAlgo_Builder aBuilder;
96   BOPAlgo_PaveFiller aPF;
97   TopoDS_Shape aFeaturesCompound;
98
99   // Obtain only edges from the features list
100   std::list<std::shared_ptr<GeomAPI_Shape> > anEdges;
101   std::list<std::shared_ptr<GeomAPI_Shape> >::const_iterator aFeatIt = theFeatures.begin();
102   for (; aFeatIt != theFeatures.end(); aFeatIt++) {
103     std::shared_ptr<GeomAPI_Shape> aShape(*aFeatIt);
104     const TopoDS_Edge& anEdge = aShape->impl<TopoDS_Edge>();
105     if (anEdge.ShapeType() == TopAbs_EDGE)
106       anEdges.push_back(aShape);
107   }
108
109   if (anEdges.size() == 1) {  // If there is only one feature, BOPAlgo_Builder will decline to work. Need to process it anyway
110     aFeaturesCompound = anEdges.front()->impl<TopoDS_Shape>();
111   } else {
112     std::list<std::shared_ptr<GeomAPI_Shape> >::const_iterator anIt = anEdges.begin();
113     for (; anIt != anEdges.end(); anIt++) {
114       std::shared_ptr<GeomAPI_Shape> aPreview(*anIt);
115       aBuilder.AddArgument(aPreview->impl<TopoDS_Edge>());
116     }
117     aPF.SetArguments(aBuilder.Arguments());
118     aPF.Perform();
119     int aErr = aPF.ErrorStatus();
120     if (aErr)
121       return;
122     aBuilder.PerformWithFiller(aPF);
123     aErr = aBuilder.ErrorStatus();
124     if (aErr)
125       return;
126     aFeaturesCompound = aBuilder.Shape();
127   }
128
129   BOPCol_IndexedDataMapOfShapeListOfShape aMapVE;  // map between vertexes and edges
130   BOPTools::MapShapesAndAncestors(aFeaturesCompound, TopAbs_VERTEX, TopAbs_EDGE, aMapVE);
131   if (aMapVE.IsEmpty())  // in case of not-initialized circle
132     return;
133
134   gp_Dir aDirX = theDirX->impl<gp_Dir>();
135   gp_Dir aDirY = theDirY->impl<gp_Dir>();
136   gp_Dir aNorm = theNorm->impl<gp_Dir>();
137
138   gp_Pln aPlane(theOrigin->impl<gp_Pnt>(), aNorm);
139
140   // Set of edges used in loops
141   std::set<TopoDS_Edge*> anEdgesInLoops;
142   // Lists for processed vertexes and edges
143   std::list<TopoDS_Vertex> aProcVertexes;
144   std::list<TopoDS_Edge> aProcEdges;
145
146   // Search the start vertex
147   TopoDS_Vertex aStartVertex = findStartVertex(aMapVE, aDirX, aDirY);
148   aProcVertexes.push_back(aStartVertex);
149
150   TopoDS_Vertex aCurVertex = aStartVertex;
151   gp_Dir aCurDir = aDirY.Reversed();
152   gp_Dir aCurNorm = aNorm.Reversed();
153
154   // Go through the edges and find loops
155   TopoDS_Vertex aNextVertex;
156   TopoDS_Edge aBindingEdge;
157   gp_Dir aNextDir;
158   while (aMapVE.Extent() > 0) {
159     if (aCurVertex.IsNull())
160       return;
161     if (!aProcEdges.empty())
162       aBindingEdge = aProcEdges.back();
163     findNextVertex(aCurVertex, aMapVE, aCurDir, aCurNorm, aNextVertex, aBindingEdge, aNextDir);
164     aCurNorm = aNorm;
165
166     // Try to find next vertex in the list of already processed
167     bool isLoopFound = false;
168     std::list<TopoDS_Vertex>::iterator aVertIter = aProcVertexes.begin();
169     std::list<TopoDS_Edge>::iterator anEdgeIter = aProcEdges.begin();
170     for (; aVertIter != aProcVertexes.end(); aVertIter++) {
171       if (aVertIter->IsSame(aNextVertex)) {
172         isLoopFound = true;
173         break;
174       }
175       if (anEdgeIter != aProcEdges.end())
176         anEdgeIter++;
177     }
178
179     bool isCircleFound = (isLoopFound && anEdgeIter == aProcEdges.end());
180     aProcVertexes.push_back(aNextVertex);
181     aProcEdges.push_back(aBindingEdge);
182
183     if (isLoopFound) {
184       // If the binding edge is a full circle, then it may be added recently. Need to update edge iterator
185       if (isCircleFound) {
186         anEdgeIter = aProcEdges.end();
187         anEdgeIter--;
188       }
189       else if (aVertIter != aProcVertexes.begin()) {
190         // Check the orientation of the loop
191         gp_Dir aCN = getOuterEdgeDirection(*anEdgeIter, *aVertIter);
192         gp_Dir aCP = getOuterEdgeDirection(aProcEdges.back(), *aVertIter);
193         aCN.Reverse();
194         aCP.Reverse();
195         if (aCN.DotCross(aCP, aNorm) < -tolerance) {
196           // The found loop has wrong orientation and may contain sub-loops.
197           // Need to check it once again with another initial direction.
198           aCurVertex = *aVertIter;
199           do {
200             aProcVertexes.pop_back();
201             aProcEdges.pop_back();
202           } while (aCurVertex != aProcVertexes.back());
203           aCurDir = aCN.Reversed();
204           aCurNorm = aNorm.Reversed();
205           continue;
206         }
207       }
208
209       if (!isCircleFound && anEdgeIter != aProcEdges.end() && 
210           anEdgeIter->IsSame(aProcEdges.back())) { // The loop on the same edge taken twice
211         aProcVertexes.pop_back();
212         aProcEdges.pop_back();
213         aCurVertex = aProcVertexes.back();
214         aCurDir = getOuterEdgeDirection(aProcEdges.back(), aCurVertex);
215         aCurNorm = aNorm.Reversed();
216         continue;
217       }
218
219       // When the orientation is correct or the edges looped through
220       // the first element, create new face and remove unnecessary edges.
221       TopoDS_Face aPatch;
222       createFace(*aVertIter, anEdgeIter, aProcEdges.end(), aPlane, aPatch);
223       if (!aPatch.IsNull()) {
224         std::shared_ptr<GeomAPI_Shape> aFace(new GeomAPI_Shape);
225         aFace->setImpl(new TopoDS_Face(aPatch));
226         theResultFaces.push_back(aFace);
227       }
228       // push the edges used in the loop to the map
229       std::list<TopoDS_Edge>::iterator anIter;
230       for (anIter = anEdgeIter; anIter != aProcEdges.end(); anIter++)
231         anEdgesInLoops.insert(&(*anIter));
232       // remove unnecessary edges
233       std::list<TopoDS_Vertex>::iterator aCopyVLoop = aVertIter;
234       std::list<TopoDS_Edge>::iterator aCopyELoop = anEdgeIter;
235       removeWasteEdges(aVertIter, anEdgeIter, aProcVertexes.end(), aProcEdges.end(), aMapVE);
236
237       // revert the list of remaining edges
238       std::list<TopoDS_Vertex> aRemainVertexes;
239       for (; aVertIter != aProcVertexes.end(); aVertIter++)
240         aRemainVertexes.push_front(*aVertIter);
241       std::list<TopoDS_Edge> aRemainEdges;
242       for (; anEdgeIter != aProcEdges.end(); anEdgeIter++)
243         aRemainEdges.push_front(*anEdgeIter);
244       // remove edges and vertexes used in the loop and add remaining ones
245       if (aCopyVLoop != aProcVertexes.begin()) {
246         aVertIter = aCopyVLoop;
247         aVertIter--;
248       } else
249         aVertIter = aProcVertexes.end();
250       aProcVertexes.erase(aCopyVLoop, aProcVertexes.end());
251       aProcVertexes.insert(aProcVertexes.end(), aRemainVertexes.begin(), aRemainVertexes.end());
252       if (aCopyELoop != aProcEdges.begin()) {
253         anEdgeIter = aCopyELoop;
254         anEdgeIter--;
255       } else
256         anEdgeIter = aProcEdges.end();
257       aProcEdges.erase(aCopyELoop, aProcEdges.end());
258       aProcEdges.insert(aProcEdges.end(), aRemainEdges.begin(), aRemainEdges.end());
259
260       if (aVertIter == aProcVertexes.end())
261         aVertIter = aProcVertexes.begin();
262       else
263         aVertIter++;
264       if (anEdgeIter == aProcEdges.end())
265         anEdgeIter = aProcEdges.begin();
266       else
267         anEdgeIter++;
268       aCopyVLoop = aVertIter;
269       aCopyELoop = anEdgeIter;
270
271       if (aVertIter != aProcVertexes.end() && 
272           aMapVE.Contains(*aVertIter) && aMapVE.FindFromKey(*aVertIter).Extent() <= 1)
273         removeWasteEdges(aVertIter, anEdgeIter, aProcVertexes.end(), aProcEdges.end(), aMapVE);
274       if (aCopyVLoop != aVertIter)
275         aProcVertexes.erase(aCopyVLoop, aVertIter);
276       if (aCopyELoop != anEdgeIter)
277         aProcEdges.erase(aCopyELoop, anEdgeIter);
278
279       // Check whether the next vertex already exists
280       if (aVertIter != aProcVertexes.end())
281         aVertIter++;
282       if (aVertIter != aProcVertexes.end() && anEdgeIter != aProcEdges.end()) {
283         aNextVertex = *aVertIter;
284         aNextDir = getOuterEdgeDirection(*anEdgeIter, aNextVertex);
285         aProcVertexes.erase(++aVertIter, aProcVertexes.end());
286         aProcEdges.erase(++anEdgeIter, aProcEdges.end());
287       } else {
288         // Recalculate current vertex and current direction
289         aProcEdges.clear();
290         aProcVertexes.clear();
291         if (aMapVE.Extent() > 0) {
292           aNextVertex = findStartVertex(aMapVE, aDirX, aDirY);
293           aProcVertexes.push_back(aNextVertex);
294         }
295         aNextDir = aDirY.Reversed();
296         aCurNorm = aNorm.Reversed();
297       }
298     }
299
300     // if next vertex connected only to alone edge, this is a part of wire (not a closed loop),
301     // we need to go back through the list of already checked edges to find a branching vertex
302     if (!aMapVE.IsEmpty() && aMapVE.Contains(aNextVertex)
303         && aMapVE.FindFromKey(aNextVertex).Size() == 1) {
304       std::list<TopoDS_Vertex>::reverse_iterator aVRIter = aProcVertexes.rbegin();
305       std::list<TopoDS_Edge>::reverse_iterator aERIter = aProcEdges.rbegin();
306       if (aVRIter != aProcVertexes.rend())
307         aVRIter++;
308       if (aERIter != aProcEdges.rend())
309         aERIter++;
310
311       for (; aERIter != aProcEdges.rend(); aERIter++, aVRIter++)
312         if (aMapVE.FindFromKey(*aVRIter).Size() > 2)
313           break;
314       if (aERIter != aProcEdges.rend()
315           || (aVRIter != aProcVertexes.rend() && aMapVE.FindFromKey(*aVRIter).Size() == 1)) {  // the branching vertex was found or current list of edges is a wire without branches
316         std::list<TopoDS_Edge>::iterator aEIter;
317         TopoDS_Edge aCurEdge;
318         if (aERIter != aProcEdges.rend()) {
319           aEIter = aERIter.base();
320           aCurEdge = *aERIter;
321         } else
322           aEIter = aProcEdges.begin();
323         std::list<TopoDS_Wire> aTail;
324         createWireList(*aVRIter, aEIter, aProcEdges.end(), anEdgesInLoops, aTail);
325         std::list<TopoDS_Wire>::const_iterator aTailIter = aTail.begin();
326         for (; aTailIter != aTail.end(); aTailIter++)
327           if (!aTailIter->IsNull()) {
328             std::shared_ptr<GeomAPI_Shape> aWire(new GeomAPI_Shape);
329             aWire->setImpl(new TopoDS_Shape(*aTailIter));
330             theResultWires.push_back(aWire);
331           }
332         std::list<TopoDS_Vertex>::iterator aVIter = aVRIter.base();
333         std::list<TopoDS_Edge>::iterator aEItCopy = aEIter;
334         removeWasteEdges(--aVIter, aEItCopy, aProcVertexes.end(), aProcEdges.end(), aMapVE);
335
336         aProcEdges.erase(aEIter, aProcEdges.end());
337         aVIter = aVRIter.base();
338         aProcVertexes.erase(aVIter, aProcVertexes.end());
339
340         if (!aProcVertexes.empty()) {
341           aNextVertex = aProcVertexes.back();
342           if (!aCurEdge.IsNull())
343             aNextDir = getOuterEdgeDirection(aCurEdge, aNextVertex);
344         }
345       } else {  // there is no branching vertex in the list of proceeded,
346                 // so we should revert the list and go the opposite way
347         aProcVertexes.reverse();
348         aProcEdges.reverse();
349         aNextVertex = aProcVertexes.back();
350         aNextDir =
351             aProcEdges.empty() ? aDirY : getOuterEdgeDirection(aProcEdges.back(), aNextVertex);
352       }
353     }
354
355     // When all edges of current component of connectivity are processed,
356     // we should repeat procedure for finding initial vertex in the remaining
357     if (!aMapVE.IsEmpty() && !aMapVE.Contains(aNextVertex)) {
358       aProcVertexes.clear();
359       aProcEdges.clear();
360
361       aStartVertex = findStartVertex(aMapVE, aDirX, aDirY);
362       aProcVertexes.push_back(aStartVertex);
363       aNextVertex = aStartVertex;
364       aNextDir = aDirY.Reversed();
365       aCurNorm = aNorm.Reversed();
366     }
367
368     aCurVertex = aNextVertex;
369     aCurDir = aNextDir;
370   }
371
372   if (theResultFaces.size() > 1)
373     fixIntersections(theResultFaces);
374 }
375
376 void GeomAlgoAPI_SketchBuilder::createFaces(const std::shared_ptr<GeomAPI_Pnt>& theOrigin,
377                                             const std::shared_ptr<GeomAPI_Dir>& theDirX,
378                                             const std::shared_ptr<GeomAPI_Dir>& theDirY,
379                                             const std::shared_ptr<GeomAPI_Dir>& theNorm,
380                                             const std::shared_ptr<GeomAPI_Shape>& theWire,
381                                             std::list<std::shared_ptr<GeomAPI_Shape> >& theResultFaces)
382 {
383   std::shared_ptr<GeomAPI_PlanarEdges> aWire = std::dynamic_pointer_cast<GeomAPI_PlanarEdges>(theWire);
384   if(!aWire)
385     return;
386   // Filter wires, return only faces.
387   std::list<std::shared_ptr<GeomAPI_Shape> > aFilteredWires;
388   createFaces(theOrigin, theDirX, theDirY, theNorm,
389               aWire->getEdges(), theResultFaces, aFilteredWires);
390 }
391
392
393 void GeomAlgoAPI_SketchBuilder::fixIntersections(
394     std::list<std::shared_ptr<GeomAPI_Shape> >& theFaces)
395 {
396   BRepClass_FaceClassifier aClassifier;
397
398   std::list<std::shared_ptr<GeomAPI_Shape> >::iterator anIter1 = theFaces.begin();
399   std::list<std::shared_ptr<GeomAPI_Shape> >::iterator anIter2;
400   for (; anIter1 != theFaces.end(); anIter1++) {
401     anIter2 = anIter1;
402     for (++anIter2; anIter2 != theFaces.end(); anIter2++) {
403       const TopoDS_Face& aF1 = (*anIter1)->impl<TopoDS_Face>();
404       assert(aF1.ShapeType() == TopAbs_FACE);  // all items in result list should be faces
405       TopExp_Explorer aVert2((*anIter2)->impl<TopoDS_Shape>(), TopAbs_VERTEX);
406       for (; aVert2.More(); aVert2.Next()) {
407         const TopoDS_Vertex& aV = (const TopoDS_Vertex&)aVert2.Current();
408         aClassifier.Perform(aF1, BRep_Tool::Pnt(aV), tolerance);
409         TopAbs_State aState = aClassifier.State();
410         if (aState != TopAbs_IN && aState != TopAbs_ON)
411           break;
412       }
413       if (aVert2.More()) {  // second shape is not inside first, change the shapes order and repeat comparision
414         const TopoDS_Face& aF2 = (*anIter2)->impl<TopoDS_Face>();
415         assert(aF2.ShapeType() == TopAbs_FACE);  // all items in result list should be faces
416         TopExp_Explorer aVert1((*anIter1)->impl<TopoDS_Shape>(), TopAbs_VERTEX);
417         for (; aVert1.More(); aVert1.Next()) {
418           const TopoDS_Vertex& aV = (const TopoDS_Vertex&)aVert2.Current();
419           aClassifier.Perform(aF2, BRep_Tool::Pnt(aV), tolerance);
420           TopAbs_State aState = aClassifier.State();
421           if (aState != TopAbs_IN && aState != TopAbs_ON)
422             break;
423         }
424         if (!aVert1.More()) {  // first shape should be cut from the second
425           BRepAlgoAPI_Cut aCut((*anIter2)->impl<TopoDS_Shape>(), (*anIter1)->impl<TopoDS_Shape>());
426           aCut.Build();
427           TopExp_Explorer anExp(aCut.Shape(), TopAbs_FACE);
428           bool isFirstFace = true;
429           for (; anExp.More(); anExp.Next()) {
430             if (anExp.Current().ShapeType() != TopAbs_FACE) continue;
431             if (isFirstFace) {
432               (*anIter2)->setImpl(new TopoDS_Shape(anExp.Current()));
433               isFirstFace = false;
434             } else {
435               std::shared_ptr<GeomAPI_Shape> aShape(new GeomAPI_Shape);
436               aShape->setImpl(new TopoDS_Shape(anExp.Current()));
437               theFaces.push_back(aShape);
438             }
439           }
440         }
441       } else {  // second shape should be cut from the first
442         BRepAlgoAPI_Cut aCut((*anIter1)->impl<TopoDS_Shape>(), (*anIter2)->impl<TopoDS_Shape>());
443         aCut.Build();
444         TopExp_Explorer anExp(aCut.Shape(), TopAbs_FACE);
445         bool isFirstFace = true;
446         for (; anExp.More(); anExp.Next()) {
447           if (anExp.Current().ShapeType() != TopAbs_FACE) continue;
448           if (isFirstFace) {
449             (*anIter1)->setImpl(new TopoDS_Shape(anExp.Current()));
450             isFirstFace = false;
451           } else {
452             std::shared_ptr<GeomAPI_Shape> aShape(new GeomAPI_Shape);
453             aShape->setImpl(new TopoDS_Shape(anExp.Current()));
454             theFaces.push_back(aShape);
455           }
456         }
457       }
458     }
459   }
460 }
461
462 // =================== Auxiliary functions ====================================
463 const TopoDS_Vertex& findStartVertex(const BOPCol_IndexedDataMapOfShapeListOfShape& theMapVE,
464                                     const gp_Dir& theDirX, const gp_Dir& theDirY)
465 {
466   int aStartVertexInd = 1;
467   double aMaxX = -DBL_MAX;
468   double aMaxY = -DBL_MAX;
469   int aNbVert = theMapVE.Extent();
470   for (int i = 1; i <= aNbVert; i++) {
471     const TopoDS_Vertex& aV = (const TopoDS_Vertex&) theMapVE.FindKey(i);
472     const gp_Pnt& aVertPnt = BRep_Tool::Pnt(aV);
473
474     double aX = aVertPnt.XYZ().Dot(theDirX.XYZ());
475     double aY = aVertPnt.XYZ().Dot(theDirY.XYZ());
476     if ((aX > aMaxX || (fabs(aX - aMaxX) < tolerance && aY > aMaxY))
477         && theMapVE.FindFromIndex(i).Extent() > 1) {
478       aMaxX = aX;
479       aMaxY = aY;
480       aStartVertexInd = i;
481     }
482   }
483   return static_cast<const TopoDS_Vertex&>(theMapVE.FindKey(aStartVertexInd));
484 }
485
486 void findNextVertex(const TopoDS_Vertex& theStartVertex,
487                     const BOPCol_IndexedDataMapOfShapeListOfShape& theVertexEdgeMap,
488                     const gp_Dir& theStartDir, const gp_Dir& theNormal, TopoDS_Vertex& theNextVertex,
489                     TopoDS_Edge& theNextEdge, gp_Dir& theNextDir)
490 {
491   theNextVertex = TopoDS_Vertex();
492   const BOPCol_ListOfShape& anEdgesList = theVertexEdgeMap.FindFromKey(theStartVertex);
493   int anEdgesNum = anEdgesList.Extent();
494   BOPCol_ListOfShape::Iterator aEdIter(anEdgesList);
495   double aBestEdgeProj = DBL_MAX;
496   for (; aEdIter.More(); aEdIter.Next()) {
497     const TopoDS_Edge& anEdge = static_cast<const TopoDS_Edge&>(aEdIter.Value());
498     gp_Dir aTang = getOuterEdgeDirection(anEdge, theStartVertex);
499     aTang.Reverse();
500
501     // The projection is normalized in segment (-1, 1),
502     // where (-1, 0] corresponds to the angles (pi/2, 0] between theStartDir and aTang
503     // and [0, 1) corresponds to the angles [0, -pi/2)
504     double aProj = (aTang.Dot(theStartDir) - 1.0) * 0.5;
505     if (anEdgesNum > 1 && fabs(fabs(aProj) - 1) < tolerance)
506       continue;
507     if (theStartDir.DotCross(aTang, theNormal) < tolerance)
508       aProj *= -1.0;
509
510     if (aProj < aBestEdgeProj) {
511       aBestEdgeProj = aProj;
512       theNextEdge = anEdge;
513       TopExp_Explorer aVertExp(theNextEdge, TopAbs_VERTEX);
514       for (; aVertExp.More(); aVertExp.Next())
515         if (!aVertExp.Current().IsSame(theStartVertex)) {
516           theNextVertex = static_cast<const TopoDS_Vertex&>(aVertExp.Current());
517           theNextDir = getOuterEdgeDirection(anEdge, theNextVertex);
518           break;
519         }
520       if (!aVertExp.More()) {  // This edge is a full circle
521         TopoDS_Vertex aV1, aV2;
522         TopExp::Vertices(theNextEdge, aV1, aV2);
523         if (aV1.Orientation() == theStartVertex.Orientation())
524           theNextVertex = aV2;
525         else
526           theNextVertex = aV1;
527         theNextDir = getOuterEdgeDirection(anEdge, theNextVertex);
528       }
529     }
530   }
531
532   // Probably there are two tangent edges. We will take the edge differs from current one
533   if (theNextVertex.IsNull() && anEdgesNum == 2) {
534     BOPCol_ListOfShape::Iterator aEdIter(anEdgesList);
535     if (aEdIter.Value() == theNextEdge)
536       aEdIter.Next();
537     theNextEdge = static_cast<const TopoDS_Edge&>(aEdIter.Value());
538     TopoDS_Vertex aV1, aV2;
539     TopExp::Vertices(theNextEdge, aV1, aV2);
540     theNextVertex = theStartVertex.IsSame(aV1) ? aV2 : aV1;
541     theNextDir = getOuterEdgeDirection(theNextEdge, theNextVertex);
542   }
543 }
544
545 static void addEdgeToWire(const TopoDS_Edge& theEdge, const BRep_Builder& theBuilder,
546                           TopoDS_Shape& theSpliceVertex, TopoDS_Wire& theWire)
547 {
548   TopoDS_Edge anEdge = theEdge;
549   bool isCurVertChanged = false;
550   TopoDS_Shape aCurVertChanged;
551
552   TopExp_Explorer aVertExp(theEdge, TopAbs_VERTEX);
553   for (; aVertExp.More(); aVertExp.Next()) {
554     const TopoDS_Shape& aVertex = aVertExp.Current();
555     if (aVertex.IsSame(theSpliceVertex)
556         && aVertex.Orientation() != theEdge.Orientation()) {  // Current vertex is the last for the edge, so its orientation is wrong, need to revert the edge
557       anEdge.Reverse();
558       break;
559     }
560     if (!aVertex.IsSame(theSpliceVertex)) {
561       aCurVertChanged = aVertex;
562       isCurVertChanged = true;
563     }
564   }
565   theSpliceVertex = isCurVertChanged ? aCurVertChanged : aVertExp.Current();
566
567   theBuilder.Add(theWire, anEdge);
568 }
569
570 void createFace(const TopoDS_Vertex& theStartVertex,
571                 const std::list<TopoDS_Edge>::iterator& theStartEdge,
572                 const std::list<TopoDS_Edge>::iterator& theEndOfEdges,
573                 const gp_Pln& thePlane,
574                 TopoDS_Face& theResFace)
575 {
576   TopoDS_Wire aResWire;
577   BRep_Builder aBuilder;
578   aBuilder.MakeWire(aResWire);
579
580   TopoDS_Vertex aCurVertex = theStartVertex;
581   std::list<TopoDS_Edge>::const_iterator anEdgeIter = theStartEdge;
582   for (; anEdgeIter != theEndOfEdges; anEdgeIter++) {
583     if (!anEdgeIter->IsNull())
584       addEdgeToWire(*anEdgeIter, aBuilder, aCurVertex, aResWire);
585   }
586
587   BRepBuilderAPI_MakeFace aFaceBuilder(thePlane, aResWire);
588   if (aFaceBuilder.Error() == BRepBuilderAPI_FaceDone)
589     theResFace = aFaceBuilder.Face();
590 }
591
592 void createWireList(const TopoDS_Vertex& theStartVertex,
593                     const std::list<TopoDS_Edge>::iterator& theStartEdge,
594                     const std::list<TopoDS_Edge>::iterator& theEndOfEdges,
595                     const std::set<TopoDS_Edge*>& theEdgesInLoops,
596                     std::list<TopoDS_Wire>& theResWires)
597 {
598   BRep_Builder aBuilder;
599   bool needNewWire = true;
600   TopoDS_Vertex aCurVertex = theStartVertex;
601
602   std::list<TopoDS_Edge>::iterator anIter = theStartEdge;
603   while (anIter != theEndOfEdges) {
604     while (anIter != theEndOfEdges && needNewWire && theEdgesInLoops.count(&(*anIter)) != 0) {
605       TopExp_Explorer aVertExp(*anIter, TopAbs_VERTEX);
606       for (; aVertExp.More(); aVertExp.Next())
607         if (!aVertExp.Current().IsSame(aCurVertex)) {
608           aCurVertex = static_cast<const TopoDS_Vertex&>(aVertExp.Current());
609           break;
610         }
611       anIter++;
612     }
613     if (anIter == theEndOfEdges)
614       break;
615
616     if (needNewWire) {  // The new wire should be created
617       TopoDS_Wire aWire;
618       aBuilder.MakeWire(aWire);
619       theResWires.push_back(aWire);
620       needNewWire = false;
621     } else if (theEdgesInLoops.count(&(*anIter)) != 0) {  // There was found the edge already used in loop.
622                                                           // Current wire should be released and new one should started
623       needNewWire = true;
624       continue;
625     }
626
627     addEdgeToWire(*anIter, aBuilder, aCurVertex, theResWires.back());
628     anIter++;
629   }
630 }
631
632 gp_Dir getOuterEdgeDirection(const TopoDS_Edge& theEdge, const TopoDS_Vertex& theVertex)
633 {
634   gp_Pnt aVertexPnt = BRep_Tool::Pnt(theVertex);
635
636   // Convert the edge to the curve to calculate the tangency.
637   // There should be only one curve in the edge.
638   double aFirst, aLast;
639   Handle(Geom_Curve) aCurve = BRep_Tool::Curve(theEdge, aFirst, aLast);
640
641   gp_Pnt aPnt;
642   gp_Vec aTang;
643   // A direction is determined not in the boundary points but in the points with small shift.
644   // It was done to avoid tangency between circle and other edge in the shared vertex.
645   aCurve->D1(aFirst + shift > aLast ? aFirst : aFirst + shift, aPnt, aTang);
646   aCurve->D0(aFirst, aPnt);
647   if (aVertexPnt.IsEqual(aPnt, tolerance))
648     return gp_Dir(aTang.Reversed());
649
650   aCurve->D1(aLast - shift < aFirst ? aLast : aLast - shift, aPnt, aTang);
651   return gp_Dir(aTang);
652 }
653
654 void removeWasteEdges(std::list<TopoDS_Vertex>::iterator& theStartVertex,
655                       std::list<TopoDS_Edge>::iterator& theStartEdge,
656                       const std::list<TopoDS_Vertex>::iterator& theEndOfVertexes,
657                       const std::list<TopoDS_Edge>::iterator& theEndOfEdges,
658                       BOPCol_IndexedDataMapOfShapeListOfShape& theMapVE)
659 {
660   bool isVertStep = true;
661   while (theStartVertex != theEndOfVertexes && theStartEdge != theEndOfEdges) {
662     BOPCol_ListOfShape& aBunch = theMapVE.ChangeFromKey(*theStartVertex);
663     BOPCol_ListOfShape::Iterator anApprEdge(aBunch);
664     for (; anApprEdge.More(); anApprEdge.Next())
665       if (anApprEdge.Value().IsSame(*theStartEdge))
666         break;
667     if (anApprEdge.More())
668       aBunch.Remove(anApprEdge);
669
670     if (isVertStep)
671       theStartVertex++;
672     else {
673       theStartEdge++;
674       // check current vertex to be a branching point
675       // if so, it will be a new starting point to find a loop
676       if (aBunch.Size() > 1)
677         break;
678     }
679     isVertStep = !isVertStep;
680   }
681
682   // The map of vertex-edges may be changed
683   BOPCol_IndexedDataMapOfShapeListOfShape aMapVECopy;
684   BOPCol_IndexedDataMapOfShapeListOfShape::Iterator aMapIter(theMapVE);
685   for (int ind = 1; aMapIter.More(); aMapIter.Next(), ind++)
686     if (!aMapIter.Value().IsEmpty())
687       aMapVECopy.Add(theMapVE.FindKey(ind), aMapIter.Value());
688   theMapVE.Clear();
689   theMapVE.Exchange(aMapVECopy);
690 }
691