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