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