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