]> OCCT Git - occt-copy.git/commitdiff
Implement the new class that allows merging of planar faces constituting a cylindrica...
authormsv <msv@opencascade.com>
Fri, 5 Aug 2016 08:13:44 +0000 (11:13 +0300)
committerakz <akz@opencascade.com>
Fri, 12 Aug 2016 10:17:55 +0000 (13:17 +0300)
The algorithm has front end in draw as a command combinetocylinder.

src/SWDRAW/SWDRAW_ShapeUpgrade.cxx
src/ShapeUpgrade/FILES
src/ShapeUpgrade/ShapeUpgrade_CombineToCylinder.cxx [new file with mode: 0644]
src/ShapeUpgrade/ShapeUpgrade_CombineToCylinder.hxx [new file with mode: 0644]

index 4735f8fc16d5f3235f4d8c26a5309985c2f1b77d..ff34103de9980ca4bc41eec563dad97a0012ceff 100644 (file)
@@ -41,6 +41,7 @@
 #include <ShapeFix.hxx>
 #include <ShapeFix_ComposeShell.hxx>
 #include <ShapeUpgrade.hxx>
+#include <ShapeUpgrade_CombineToCylinder.hxx>
 #include <ShapeUpgrade_RemoveInternalWires.hxx>
 #include <ShapeUpgrade_RemoveLocations.hxx>
 #include <ShapeUpgrade_ShapeConvertToBezier.hxx>
@@ -1476,6 +1477,69 @@ Standard_Integer reshape(Draw_Interpretor& di,
   return 0;
 }
 
+//=======================================================================
+// CombineToCylinder
+//=======================================================================
+static Standard_Integer CombineToCylinder(Draw_Interpretor& di, Standard_Integer n, const char** a)
+{
+  if (n < 3)
+  {
+    di << "Use combinetocyl result shape [-t val] [-a val]\n";
+    di << "options:\n";
+    di << "-t val to set linear tolerance (default is 1e-7)\n";
+    di << "-a val to set angular tolerance (default is 1e-12)\n";
+    return 1;
+  }
+
+  TopoDS_Shape aShape = DBRep::Get(a[2]);
+  if (aShape.IsNull())
+    return 1;
+
+  // default values
+  Standard_Real aLinTol = Precision::Confusion();
+  Standard_Real aAngTol = Precision::Angular();
+
+  if (n > 3)
+  {
+    for (int i = 3; i < n; i++)
+    {
+      if (!strcmp(a[i], "-t") || !strcmp(a[i], "-a"))
+      {
+        if (++i < n)
+        {
+          if (a[i - 1][1] == 't')
+            aLinTol = Draw::Atof(a[i]);
+          else
+            aAngTol = Draw::Atof(a[i]) * M_PI / 180.;
+        }
+        else
+        {
+          di << "value expected after " << a[i - 1];
+          return 1;
+        }
+      }
+    }
+  }
+
+  ShapeUpgrade_CombineToCylinder aTool;
+  aTool.SetShape(aShape);
+  aTool.SetLinearTolerance(aLinTol);
+  aTool.SetAngularTolerance(aAngTol);
+  aTool.Build();
+  if (!aTool.IsDone())
+  {
+    di << "Not done";
+  }
+  else
+  {
+    di << "Number of new merged faces: " << aTool.NbNewFaces();
+    di << "\nNumber of replaced faces: " << aTool.NbReplacedFaces();
+    const TopoDS_Shape& Result = aTool.Shape();
+    DBRep::Set(a[1], Result);
+  }
+  return 0;
+}
+
 //=======================================================================
 //function : InitCommands
 //purpose  : 
@@ -1589,4 +1653,7 @@ Standard_Integer reshape(Draw_Interpretor& di,
     "\n      '-remove what'         Removes 'what' sub-shape"
     "\n    Requests '-replace' and '-remove' can be repeated many times.",
     __FILE__, reshape, g);
+
+  theCommands.Add("combinetocylinder",
+                  "run w/o arguments to get help", __FILE__, CombineToCylinder, g);
 }
index 74b07ba1cc3710397a4f69a701358e12a0acb6aa..3ea7ef3e092972e76cd54f967764e00f12599965 100644 (file)
@@ -4,6 +4,8 @@ ShapeUpgrade_ClosedEdgeDivide.cxx
 ShapeUpgrade_ClosedEdgeDivide.hxx
 ShapeUpgrade_ClosedFaceDivide.cxx
 ShapeUpgrade_ClosedFaceDivide.hxx
+ShapeUpgrade_CombineToCylinder.cxx
+ShapeUpgrade_CombineToCylinder.hxx
 ShapeUpgrade_ConvertCurve2dToBezier.cxx
 ShapeUpgrade_ConvertCurve2dToBezier.hxx
 ShapeUpgrade_ConvertCurve3dToBezier.cxx
diff --git a/src/ShapeUpgrade/ShapeUpgrade_CombineToCylinder.cxx b/src/ShapeUpgrade/ShapeUpgrade_CombineToCylinder.cxx
new file mode 100644 (file)
index 0000000..3f993c2
--- /dev/null
@@ -0,0 +1,885 @@
+//-Copyright: Open CASCADE 2016
+// Created on: 2016-08-01
+// Created by: Mikhail Sazonov
+// Copyright (c) 2016 OPEN CASCADE SAS
+//
+// This file is part of Open CASCADE Technology software library.
+//
+// This library is free software; you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License version 2.1 as published
+// by the Free Software Foundation, with special exception defined in the file
+// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
+// distribution for complete text of the license and disclaimer of any warranty.
+//
+// Alternatively, this file may be used under the terms of Open CASCADE
+// commercial license or contractual agreement.
+
+#include <ShapeUpgrade_CombineToCylinder.hxx>
+#include <BRep_Builder.hxx>
+#include <BRepAdaptor_Curve.hxx>
+#include <BRepAdaptor_Surface.hxx>
+#include <BRepBuilderAPI_MakeEdge.hxx>
+#include <BRepBuilderAPI_MakeFace.hxx>
+#include <BRepTools.hxx>
+#include <BRepTools_ReShape.hxx>
+#include <ElCLib.hxx>
+#include <GCPnts_AbscissaPoint.hxx>
+#include <Geom_Circle.hxx>
+#include <Geom_BSplineCurve.hxx>
+#include <GeomAPI_PointsToBSpline.hxx>
+#include <gce_MakeCirc.hxx>
+#include <gp_Ax1.hxx>
+#include <gp_Ax3.hxx>
+#include <gp_Circ.hxx>
+#include <gp_Cylinder.hxx>
+#include <gp_Pln.hxx>
+#include <IntAna_IntConicQuad.hxx>
+#include <ShapeAnalysis_Curve.hxx>
+#include <ShapeFix_Shape.hxx>
+#include <TColStd_MapOfInteger.hxx>
+#include <TopoDS.hxx>
+#include <TopoDS_Compound.hxx>
+#include <TopoDS_Edge.hxx>
+#include <TopoDS_Face.hxx>
+#include <TopExp.hxx>
+#include <TopExp_Explorer.hxx>
+#include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
+#include <TopTools_IndexedMapOfShape.hxx>
+#include <TopTools_MapOfShape.hxx>
+#include <Precision.hxx>
+
+//=======================================================================
+//function : ShapeUpgrade_CombineToCylinder
+//purpose  : 
+//=======================================================================
+ShapeUpgrade_CombineToCylinder::ShapeUpgrade_CombineToCylinder()
+  : myLinTol(Precision::Confusion()),
+    myAngTol(Precision::Angular()),
+    myNbNewFaces(0),
+    myNbReplacedFaces(0),
+    myDone(Standard_False)
+{
+  myContext = new BRepTools_ReShape;
+}
+
+//=======================================================================
+//function : sortEdges
+//purpose  : 
+//=======================================================================
+static void sortEdges(const TopoDS_Shape& theShape,
+                            const TopTools_IndexedDataMapOfShapeListOfShape& theMapEF,
+                            const Standard_Real theLinTol,
+                            const Standard_Real theAngTol,
+                            TopTools_MapOfShape& theStopEdges,
+                            NCollection_DataMap<TopoDS_Shape,gp_Cylinder>& theMapECyl)
+{
+  TopExp_Explorer anEx(theShape, TopAbs_EDGE);
+  for (; anEx.More(); anEx.Next())
+  {
+    const TopoDS_Edge& anEdge = TopoDS::Edge(anEx.Current());
+    const TopTools_ListOfShape* pFaces = theMapEF.Seek(anEdge);
+    if (!pFaces)
+      continue;
+
+    // check for manifold edge
+    if (pFaces->Extent() != 2)
+    {
+      theStopEdges.Add(anEdge);
+      continue;
+    }
+
+    // Get curve of the edge and check if it is a line
+    BRepAdaptor_Curve aCurve(anEdge);
+    if (aCurve.GetType() != GeomAbs_Line)
+    {
+      theStopEdges.Add(anEdge);
+      continue;
+    }
+    // Check if the edge has enough length
+    Standard_Real aLength = GCPnts_AbscissaPoint::Length(aCurve);
+    if (aLength < 2.*theLinTol)
+    {
+      theStopEdges.Add(anEdge);
+      continue;
+    }
+
+    // check that both faces are planar, if no then add the edge to stop list
+    TopoDS_Face aFaces[2] = { TopoDS::Face(pFaces->First()), TopoDS::Face(pFaces->Last()) };
+    BRepAdaptor_Surface aSurf[2] = {
+      BRepAdaptor_Surface(aFaces[0], Standard_False),
+      BRepAdaptor_Surface(aFaces[1], Standard_False) };
+    if (aSurf[0].GetType() != GeomAbs_Plane || aSurf[1].GetType() != GeomAbs_Plane)
+    {
+      theStopEdges.Add(anEdge);
+      continue;
+    }
+
+    // check if edge is sharp, if yes then add it to stop list
+    gp_Pln aPl1 = aSurf[0].Plane();
+    gp_Pln aPl2 = aSurf[1].Plane();
+    if (aPl1.Axis().Direction().Angle(aPl2.Axis().Direction()) > theAngTol)
+    {
+      theStopEdges.Add(anEdge);
+      continue;
+    }
+
+    // Create plane perpendicular to the edge
+    gp_Lin aLin = aCurve.Line();
+    gp_Pln aPln(aLin.Location(), aLin.Direction());
+
+    // Get opposite edges on each face
+    TopoDS_Edge anOpEdges[2];
+    gp_Lin anOpLines[2];
+    Standard_Integer i;
+    for (i = 0; i < 2; i++)
+    {
+      TopExp_Explorer anEx1(aFaces[i], TopAbs_EDGE);
+      Standard_Integer nbE = 0;
+      for (; anEx1.More(); anEx1.Next(), nbE++)
+      {
+        if (anOpEdges[i].IsNull())
+        {
+          const TopoDS_Edge& aE = TopoDS::Edge(anEx1.Current());
+          if (!aE.IsSame(anEdge))
+          {
+            BRepAdaptor_Curve aC(aE);
+            if (aC.GetType() == GeomAbs_Line)
+            {
+              gp_Lin aL = aC.Line();
+              if (aL.Direction().IsParallel(aLin.Direction(), theAngTol))
+              {
+                anOpEdges[i] = aE;
+                anOpLines[i] = aL;
+              }
+            }
+          }
+        }
+      }
+      if (nbE != 4)
+        break; // face has non-rectangular form
+    }
+    if (anOpEdges[0].IsNull() || anOpEdges[1].IsNull())
+    {
+      theStopEdges.Add(anEdge);
+      continue;
+    }
+
+    // intersect all 3 edges with the plane
+    gp_Lin aLines[3] = { anOpLines[0], aLin, anOpLines[1] };
+    gp_Pnt anIntPnt[3];
+    for (i = 0; i < 3; i++)
+    {
+      IntAna_IntConicQuad anInter(aLines[i], aPln, theAngTol);
+      if (!anInter.IsDone() || anInter.NbPoints() != 1)
+        break;
+      anIntPnt[i] = anInter.Point(1);
+    }
+    if (i < 3)
+    {
+      theStopEdges.Add(anEdge);
+      continue;
+    }
+
+    // check that distances to points on opposite edges are equal
+    Standard_Real aDist1 = anIntPnt[0].Distance(anIntPnt[1]);
+    Standard_Real aDist2 = anIntPnt[2].Distance(anIntPnt[1]);
+    if (fabs(aDist1 - aDist2) > theLinTol)
+    {
+      theStopEdges.Add(anEdge);
+      continue;
+    }
+
+    // compute the circle radius
+    gce_MakeCirc aCircMaker(anIntPnt[0], anIntPnt[1], anIntPnt[2]);
+    if (!aCircMaker.IsDone())
+    {
+      theStopEdges.Add(anEdge);
+      continue;
+    }
+    Standard_Real aRadius = aCircMaker.Value().Radius();
+    const gp_Ax1& aCylAxis = aCircMaker.Value().Axis();
+    gp_Cylinder aCyl(gp_Ax3(aCylAxis.Location(), aCylAxis.Direction()), aRadius);
+    theMapECyl.Bind(anEdge, aCyl);
+  }
+}
+
+//=======================================================================
+//function : updateBoundary
+//purpose  : 
+//=======================================================================
+static void updateBoundary(TopTools_ListOfShape& theBoundary,
+                           const TopoDS_Shape& theFace,
+                           const TopoDS_Shape& theConnectingEdge,
+                           const TopTools_MapOfShape& theStopEdges)
+{
+  // iterate on all edges of the face
+  TopExp_Explorer anEx(theFace, TopAbs_EDGE);
+  for (; anEx.More(); anEx.Next())
+  {
+    const TopoDS_Edge& anEdge = TopoDS::Edge(anEx.Current());
+    if (anEdge.IsSame(theConnectingEdge) || theStopEdges.Contains(anEdge))
+      continue; // skip stop edges
+    // find this edge in the boundary
+    TopTools_ListIteratorOfListOfShape anIt(theBoundary);
+    while (anIt.More())
+    {
+      if (anIt.Value().IsSame(anEdge))
+        break;
+      anIt.Next();
+    }
+    if (anIt.More())
+    {
+      // found, remove it from boundary
+      theBoundary.Remove(anIt);
+    }
+    else
+    {
+      // not found, add to boundary
+      theBoundary.Append(anEdge);
+    }
+  }
+}
+
+//=======================================================================
+//function : Auxiliary structures
+//purpose  : 
+//=======================================================================
+namespace CombineToCylinder
+{
+  struct Area
+  {
+    TopoDS_Compound    Faces;
+    const gp_Cylinder* pCylinder;
+    Standard_Integer   nbFaces;
+
+    Area() : pCylinder(0L), nbFaces(0)
+    {}
+    void Init()
+    {
+      Faces.Nullify();
+      nbFaces = 0;
+      BRep_Builder().MakeCompound(Faces);
+    }
+    void AddFace(const TopoDS_Shape& theFace)
+    {
+      BRep_Builder().Add(Faces, theFace);
+      nbFaces++;
+    }
+  };
+}
+
+//=======================================================================
+//function : EqualCylinders
+//purpose  : 
+//=======================================================================
+inline Standard_Boolean EqualCylinders(const gp_Cylinder& theCyl1,
+                                       const gp_Cylinder& theCyl2,
+                                       const Standard_Real theLinTol,
+                                       const Standard_Real theAngTol)
+{
+  // check radii
+  if (fabs(theCyl1.Radius() - theCyl2.Radius()) > theLinTol)
+    return Standard_False;
+  // check if parallel
+  if (!theCyl1.Axis().Direction().IsParallel(theCyl2.Axis().Direction(), theAngTol))
+    return Standard_False;
+  // check distance between axes
+  gp_XYZ aVecO1O2 = theCyl2.Location().XYZ() - theCyl1.Location().XYZ();
+  Standard_Real aDist1 = aVecO1O2.CrossSquareMagnitude(theCyl2.Axis().Direction().XYZ());
+  Standard_Real aDist2 = aVecO1O2.CrossSquareMagnitude(theCyl1.Axis().Direction().XYZ());
+  Standard_Real aTol = theLinTol * theLinTol;
+  return aDist1 <= aTol && aDist2 < aTol;
+}
+
+//=======================================================================
+//function : collectAreas
+//purpose  : 
+//=======================================================================
+static void collectAreas(const TopoDS_Shape& theShape,
+                         const TopTools_IndexedDataMapOfShapeListOfShape& theMapEF,
+                         const TopTools_MapOfShape& theStopEdges,
+                         const NCollection_DataMap<TopoDS_Shape, gp_Cylinder>& theMapECyl,
+                         const Standard_Real theLinTol,
+                         const Standard_Real theAngTol,
+                         NCollection_List<CombineToCylinder::Area>& theAreas)
+{
+  TopTools_MapOfShape aPassedFaces;
+
+  // Repeat while there are unprocessed faces
+  for (;;)
+  {
+    // get first unprocessed face
+    TopoDS_Shape aCurFace;
+    TopExp_Explorer anEx(theShape, TopAbs_FACE);
+    for (; anEx.More(); anEx.Next())
+    {
+      aCurFace = anEx.Current();
+      if (!aPassedFaces.Contains(aCurFace))
+        break;
+    }
+    if (!anEx.More())
+      break; // all are processed, stop
+
+    // Start area with it
+    CombineToCylinder::Area anArea;
+    anArea.Init();
+    anArea.AddFace(aCurFace);
+    aPassedFaces.Add(aCurFace);
+
+    // Define edges of boundary of area
+    TopTools_ListOfShape aBoundary;
+    updateBoundary(aBoundary, aCurFace, TopoDS_Shape(), theStopEdges);
+
+    // Add all neighbors to the area while there are unprocessed edges of boundary
+    while (!aBoundary.IsEmpty())
+    {
+      TopTools_ListIteratorOfListOfShape anIt(aBoundary);
+      while (anIt.More())
+      {
+        const TopoDS_Shape& anEdge = anIt.Value();
+        const gp_Cylinder* pCyl = theMapECyl.Seek(anEdge);
+        const TopTools_ListOfShape* pFaces = theMapEF.Seek(anEdge);
+        Standard_ASSERT_RAISE(pCyl && pFaces && pFaces->Extent() == 2, "collectAreas");
+        if (!anArea.pCylinder || EqualCylinders(*anArea.pCylinder, *pCyl, theLinTol, theAngTol))
+        {
+          anArea.pCylinder = pCyl;
+          aCurFace = (aPassedFaces.Contains(pFaces->First()) ? pFaces->Last() : pFaces->First());
+          if (!aPassedFaces.Contains(aCurFace))
+          {
+            // Update boundary adding the next face
+            updateBoundary(aBoundary, aCurFace, anEdge, theStopEdges);
+            anArea.AddFace(aCurFace);
+            aPassedFaces.Add(aCurFace);
+          }
+        }
+        aBoundary.Remove(anIt);
+      }
+    }
+
+    // check if area is well formed
+    if (anArea.pCylinder && anArea.nbFaces > 2)
+    {
+      // Add area to the result list
+      theAreas.Append(anArea);
+    }
+  }
+}
+
+//=======================================================================
+//function : getNextEdge
+//purpose  : 
+//=======================================================================
+static Standard_Boolean getNextEdge(const TopTools_ListOfShape& theChain,
+                                    const Standard_Boolean theFwd,
+                                    const TopTools_IndexedDataMapOfShapeListOfShape& theMapVE,
+                                    const Standard_Real theSmoothAng,
+                                    TopoDS_Edge& theNextEdge)
+{
+  // Get end edge and vertex
+  TopoDS_Edge anEndE, anOppoEdge;
+  TopoDS_Vertex anEndV;
+  if (theFwd)
+  {
+    anEndE = TopoDS::Edge(theChain.Last());
+    anOppoEdge = TopoDS::Edge(theChain.First());
+    anEndV = TopExp::LastVertex(anEndE, Standard_True);
+  }
+  else
+  {
+    anEndE = TopoDS::Edge(theChain.First());
+    anOppoEdge = TopoDS::Edge(theChain.Last());
+    anEndV = TopExp::FirstVertex(anEndE, Standard_True);
+  }
+
+  // Get next edge using topology map; reject the case of multi-connected vertex,
+  // when there are more than one candidate edge
+  const TopTools_ListOfShape& aConnEdges = theMapVE.FindFromKey(anEndV);
+  // Note! aConnEdges can contain duplicates
+  theNextEdge.Nullify();
+  TopTools_ListIteratorOfListOfShape itL(aConnEdges);
+  for (; itL.More(); itL.Next())
+  {
+    const TopoDS_Shape& aE = itL.Value();
+    if (aE.IsSame(anEndE))
+      continue;
+    if (theNextEdge.IsNull())
+      theNextEdge = TopoDS::Edge(aE);
+    else if(!aE.IsSame(theNextEdge))
+      return Standard_False; // multi-connected vertex
+  }
+
+  // check for closed chain
+  if (theNextEdge.IsSame(anOppoEdge))
+    return Standard_False;
+
+  // Orient the next edge
+  TopoDS_Vertex aNextEndV = (theFwd ? TopExp::FirstVertex(theNextEdge, Standard_True) :
+    TopExp::LastVertex(theNextEdge, Standard_True));
+  if (!aNextEndV.IsSame(anEndV))
+    theNextEdge.Reverse();
+
+  if (theSmoothAng > 0.)
+  {
+    // Check that the edge is linear
+    BRepAdaptor_Curve aNextC(theNextEdge);
+    if (aNextC.GetType() != GeomAbs_Line)
+      return Standard_False;
+
+    // Verify smoothness of connection
+    BRepAdaptor_Curve anEndC(anEndE);
+    gp_Lin anEndLin = anEndC.Line();
+    gp_Lin aNextLin = aNextC.Line();
+    const gp_Dir& anEndDir = anEndLin.Direction();
+    const gp_Dir& aNextDir = aNextLin.Direction();
+    if (!anEndDir.IsParallel(aNextDir, theSmoothAng))
+      return Standard_False;
+  }
+  return Standard_True;
+}
+
+//=======================================================================
+//function : makeEdgesChain
+//purpose  : 
+//=======================================================================
+static void makeEdgesChain(TopTools_ListOfShape& theChain,
+                           const TopTools_IndexedDataMapOfShapeListOfShape& theMapVE,
+                           const Standard_Real theSmoothAng)
+{
+  // Cycle for both ends of the chain
+  Standard_Boolean isFwd[2] = { Standard_True, Standard_False };
+  for (Standard_Integer i = 0; i < 2; i++)
+  {
+    TopoDS_Edge aNextEdge;
+    while (getNextEdge(theChain, isFwd[i], theMapVE, theSmoothAng, aNextEdge))
+    {
+      // add next edge to the chain
+      if (isFwd[i])
+        theChain.Append(aNextEdge);
+      else
+        theChain.Prepend(aNextEdge);
+    }
+  }
+}
+
+//=======================================================================
+//function : createCylFaceFromArea
+//purpose  : 
+//=======================================================================
+static TopoDS_Face createCylFaceFromArea(const TopoDS_Compound& theFaces,
+                                         const gp_Cylinder& theCylinder)
+{
+  // Find boundary edges of the area
+  TopTools_IndexedDataMapOfShapeListOfShape aMapEF, aMapVE;
+  TopExp::MapShapesAndAncestors(theFaces, TopAbs_EDGE, TopAbs_FACE, aMapEF);
+  TopTools_ListOfShape aBoundary;
+  Standard_Integer i;
+  for (i = 1; i <= aMapEF.Extent(); i++)
+  {
+    if (aMapEF(i).Extent() == 1)
+    {
+      const TopoDS_Shape& aE = aMapEF.FindKey(i);
+      aBoundary.Append(aE);
+      TopExp::MapShapesAndAncestors(aE, TopAbs_VERTEX, TopAbs_EDGE, aMapVE);
+    }
+  }
+
+  // Check the case of closed cylinder. In this case there will be two 
+  // distinct chains of boundary edges, and we need to create a seam edge
+  TopoDS_Edge aSeamEdge;
+  TopTools_ListOfShape aChain;
+  aChain.Append(aBoundary.First());
+  makeEdgesChain(aChain, aMapVE, 0.);
+  TopTools_MapOfShape aMapChain;
+  TopTools_ListIteratorOfListOfShape itC(aChain);
+  for (; itC.More(); itC.Next())
+    aMapChain.Add(itC.Value());
+  TopTools_ListIteratorOfListOfShape itB(aBoundary);
+  for (; itB.More(); itB.Next())
+  {
+    if (!aMapChain.Contains(itB.Value()))
+      break;
+  }
+  if (itB.More())
+  {
+    // It is the case.
+    // Get any multi-connected edge and use it as a seam
+    for (i = 1; i <= aMapEF.Extent(); i++)
+      if (aMapEF(i).Extent() == 2)
+        aSeamEdge = TopoDS::Edge(aMapEF.FindKey(i));
+  }
+
+  // Build new cylindrical face
+  TopoDS_Wire aWire;
+  BRep_Builder().MakeWire(aWire);
+  for (itB.Init(aBoundary); itB.More(); itB.Next())
+  {
+    BRep_Builder().Add(aWire, itB.Value());
+  }
+  if (!aSeamEdge.IsNull())
+  {
+    BRep_Builder().Add(aWire, aSeamEdge);
+    BRep_Builder().Add(aWire, aSeamEdge.Reversed());
+  }
+  TopoDS_Face aCylFace = BRepBuilderAPI_MakeFace(theCylinder, aWire);
+
+  // Check orientation of an initial face, so that possibly to reverse the new face.
+  // The new cylindrical surface always has normal directed outside of cylinder.
+  // If an initial face has normal (taking face orientation into account) directed 
+  // inside then we must reverse the new face.
+  const TopoDS_Edge& aE = TopoDS::Edge(aBoundary.First());
+  const TopoDS_Face& aF = TopoDS::Face(aMapEF.FindFromKey(aE).First());
+  BRepAdaptor_Surface aSurf(aF, Standard_False);
+  Standard_Real aParF, aParL;
+  Handle(Geom2d_Curve) aPCurve = BRep_Tool::CurveOnSurface(aE, aF, aParF, aParL);
+  Standard_ASSERT_RAISE(!aPCurve.IsNull(), "createCylFaceFromArea");
+  gp_Pnt2d aP2d = aPCurve->Value(aParF);
+  gp_Pnt aP;
+  gp_Vec aDU, aDV;
+  aSurf.D1(aP2d.X(), aP2d.Y(), aP, aDU, aDV);
+  gp_XYZ aN = aDU.XYZ() ^ aDV.XYZ();
+  if (aF.Orientation() == TopAbs_REVERSED)
+    aN.Reverse();
+  gp_XYZ aVecFromAxis = aP.XYZ() - theCylinder.Location().XYZ();
+  if (aVecFromAxis * aN < 0.)
+  {
+    // need reverse
+    aCylFace.Reverse();
+  }
+  return aCylFace;
+}
+
+//=======================================================================
+//function : updateVertexOnCirc
+//purpose  : 
+//=======================================================================
+static void updateVertexOnCirc(const TopoDS_Vertex& theV,
+                               const gp_Circ& theCirc,
+                               Standard_Real& thePar)
+{
+  Standard_Real aTol = BRep_Tool::Tolerance(theV);
+  gp_Pnt aP = BRep_Tool::Pnt(theV);
+  thePar = ElCLib::Parameter(theCirc, aP);
+  gp_Pnt aProj = ElCLib::Value(thePar, theCirc);
+  Standard_Real aDist = aProj.SquareDistance(aP);
+  if (aDist > aTol*aTol)
+    BRep_Builder().UpdateVertex(theV, sqrt(aDist)*1.001);
+}
+
+//=======================================================================
+//function : createCurveEdge
+//purpose  : 
+//=======================================================================
+static TopoDS_Edge createCurveEdge(const TopTools_ListOfShape& theChain,
+                                   const TopTools_IndexedDataMapOfShapeListOfShape& theMapEF,
+                                   const Standard_Real theLinTol)
+{
+  // Find a cylindrical face common for all edges
+  TopTools_ListOfShape aCylFaces;
+  TopTools_ListIteratorOfListOfShape itL(theChain);
+  for (; itL.More(); itL.Next())
+  {
+    const TopoDS_Shape& aE = itL.Value();
+    const TopTools_ListOfShape& aFaces = theMapEF.FindFromKey(aE);
+    // get cyl faces of the current edge
+    TopTools_ListOfShape aCurCylFaces;
+    TopTools_ListIteratorOfListOfShape itF(aFaces);
+    for (; itF.More(); itF.Next())
+    {
+      const TopoDS_Face& aF = TopoDS::Face(itF.Value());
+      BRepAdaptor_Surface aSurf(aF, Standard_False);
+      if (aSurf.GetType() == GeomAbs_Cylinder)
+        aCurCylFaces.Append(aF);
+    }
+
+    if (aCylFaces.IsEmpty())
+      // initialize cyl faces with the first edge
+      aCylFaces.Append(aCurCylFaces);
+    else
+    {
+      // filter cyl faces leaving only common ones
+      TopTools_ListIteratorOfListOfShape itCF(aCylFaces);
+      while (itCF.More())
+      {
+        const TopoDS_Shape& aF = itCF.Value();
+        TopTools_ListIteratorOfListOfShape itCF1(aCurCylFaces);
+        for (; itCF1.More(); itCF1.Next())
+          if (itCF1.Value().IsSame(aF))
+            break;
+        if (!itCF1.More())
+          aCylFaces.Remove(itCF);
+        else
+          itCF.Next();
+      }
+    }
+    if (aCylFaces.IsEmpty())
+      return TopoDS_Edge(); // no common cyl faces
+  }
+  if (aCylFaces.IsEmpty())
+    return TopoDS_Edge();
+  const TopoDS_Face& aCylFace = TopoDS::Face(aCylFaces.First());
+  BRepAdaptor_Surface aCylSurf(aCylFace, Standard_False);
+  gp_Cylinder aCyl = aCylSurf.Cylinder();
+  const gp_Dir& aCylDir = aCyl.Axis().Direction();
+  const gp_Pnt& aCylLoc = aCyl.Location();
+  const Standard_Real aRadius = aCyl.Radius();
+
+  // Check if the chain is closed
+  TopoDS_Vertex aVF, aVL;
+  aVF = TopoDS::Vertex(TopExp::FirstVertex(TopoDS::Edge(theChain.First()), Standard_True));
+  aVL = TopoDS::Vertex(TopExp::LastVertex(TopoDS::Edge(theChain.Last()), Standard_True));
+  Standard_Boolean isClosed = aVF.IsSame(aVL);
+
+  // Check if all points are in plane perpendicular to the cylinder axis;
+  // in such case the curve is a circle
+  TColgp_Array1OfPnt aPnts(1, theChain.Extent()+1);
+  Standard_Integer i = 1;
+  for (itL.Init(theChain); itL.More(); itL.Next(), i++)
+  {
+    const TopoDS_Edge& aE = TopoDS::Edge(itL.Value());
+    TopoDS_Shape aV = TopExp::FirstVertex(aE, Standard_True);
+    aPnts(i) = BRep_Tool::Pnt(TopoDS::Vertex(aV));
+  }
+  aPnts(i) = BRep_Tool::Pnt(aVL);
+  gp_XYZ aMidPnt(0, 0, 0);
+  for (i = 1; i <= aPnts.Length(); i++)
+    aMidPnt += aPnts(i).XYZ();
+  aMidPnt /= aPnts.Length();
+  for (i = 1; i <= aPnts.Length(); i++)
+  {
+    Standard_Real aDist = fabs((aPnts(i).XYZ() - aMidPnt) * aCylDir.XYZ());
+    if (aDist > theLinTol)
+      break;
+  }
+  Standard_Boolean isCircle = (i > aPnts.Length());
+
+  // Create new edge
+  TopoDS_Edge aNewEdge;
+  if (isCircle)
+  {
+    // Create circular edge.
+    // Find center by projecting aMidPnt to the cylinder axis
+    gp_XYZ aVecOM = aMidPnt - aCylLoc.XYZ();
+    gp_Pnt aCenter = aCylLoc.Translated(aCylDir.XYZ() * (aVecOM * aCylDir.XYZ()));
+    gp_XYZ aXDir = aPnts(1).XYZ() - aCenter.XYZ();
+    gp_XYZ aCircDir = aCylDir.XYZ();
+    gp_XYZ aTanDir = aPnts(2).XYZ() - aPnts(1).XYZ();
+    if ((aXDir ^ aTanDir)* aCylDir.XYZ() < 0.)
+      aCircDir.Reverse();
+    gp_Circ aCirc(gp_Ax2(aCenter, aCircDir, aXDir), aRadius);
+    Handle(Geom_Circle) aCurve = new Geom_Circle(aCirc);
+    // update vertices tolerance
+    Standard_Real aPar;
+    updateVertexOnCirc(aVF, aCirc, aPar);
+    if (isClosed)
+    {
+      // full circle
+      aNewEdge = BRepBuilderAPI_MakeEdge(aCurve, aVF, aVL, 0, M_PI*2.);
+    }
+    else
+    {
+      // an arc
+      updateVertexOnCirc(aVL, aCirc, aPar);
+      aNewEdge = BRepBuilderAPI_MakeEdge(aCurve, aVF, aVL, 0., aPar);
+    }
+  }
+  else
+  {
+    gp_Pnt aPnts1[3];
+    const gp_Pnt* pPnts = &aPnts(1);
+    Standard_Integer np = aPnts.Length();
+    if (np == 2)
+    {
+      gp_XYZ anEdgeDir = aPnts(2).XYZ() - aPnts(1).XYZ();
+      const Standard_Real anEpsAng = 1e-4;
+      if (aCylDir.IsParallel(anEdgeDir, anEpsAng))
+        return aNewEdge;
+      // we have only two points; create additional point in the middle of the segment
+      aPnts1[0] = aPnts(1);
+      aPnts1[2] = aPnts(2);
+      aPnts1[1] = (aPnts1[0].XYZ() + aPnts1[2].XYZ()) * 0.5;
+      pPnts = aPnts1;
+      np = 3;
+    }
+    TColgp_Array1OfPnt aPntA(*pPnts, 1, np);
+    // Create b-spline curve
+    Standard_Integer aDegMin = 3;
+    Standard_Integer aDegMax = 3;
+    GeomAbs_Shape aCont = GeomAbs_C2;
+    GeomAPI_PointsToBSpline aCreator(aPntA, aDegMin, aDegMax, aCont, theLinTol);
+    if (aCreator.IsDone())
+    {
+      const Handle(Geom_BSplineCurve)& aCurve = aCreator.Curve();
+      Standard_Real aTol = BRep_Tool::Tolerance(aVF);
+      Standard_Real aParF = aCurve->FirstParameter();
+      Standard_Real aParL = aCurve->LastParameter();
+      Standard_Real aDist = aCurve->Value(aParF).SquareDistance(aPntA(1));
+      if (aDist > aTol*aTol)
+        BRep_Builder().UpdateVertex(aVF, sqrt(aDist));
+      if (!isClosed)
+      {
+        aTol = BRep_Tool::Tolerance(aVL);
+        aDist = aCurve->Value(aParL).SquareDistance(aPntA(np));
+        if (aDist > aTol*aTol)
+          BRep_Builder().UpdateVertex(aVL, sqrt(aDist));
+      }
+      aNewEdge = BRepBuilderAPI_MakeEdge(aCurve, aVF, aVL, aParF, aParL);
+    }
+  }
+  return aNewEdge;
+}
+
+//=======================================================================
+//function : mergeEdges
+//purpose  : 
+//=======================================================================
+static void mergeEdges(const TopoDS_Shape& theShape,
+                       const Standard_Real theLinTol,
+                       const Standard_Real theAngTol,
+                       const Handle(BRepTools_ReShape)& theContext)
+{
+  // Create topological data maps
+  TopTools_IndexedDataMapOfShapeListOfShape aMapVE, aMapEF;
+  TopExp::MapShapesAndAncestors(theShape, TopAbs_VERTEX, TopAbs_EDGE, aMapVE);
+  TopExp::MapShapesAndAncestors(theShape, TopAbs_EDGE, TopAbs_FACE, aMapEF);
+
+  // The idea is to gather chains of edges in which all intermediate vertices
+  // have only two neighboring edges, and do not have connection to any other edges
+  // from the shape.
+
+  // Prepare set of all edges
+  TopTools_IndexedMapOfShape anEdges;
+  TopExp::MapShapes(theShape, TopAbs_EDGE, anEdges);
+  TColStd_MapOfInteger anEdgesInd;  // set of indices of left edges
+  Standard_Integer i;
+  for (i = 1; i <= anEdges.Extent(); i++)
+    anEdgesInd.Add(i);
+
+  // Process while there are unprocessed edges
+  while (!anEdgesInd.IsEmpty())
+  {
+    // Get an unprocessed edge
+    Standard_Integer aEInd = TColStd_MapIteratorOfMapOfInteger(anEdgesInd).Value();
+    anEdgesInd.Remove(aEInd);
+    TopoDS_Edge aFirstEdge = TopoDS::Edge(anEdges(aEInd));
+    BRepAdaptor_Curve aCurve(aFirstEdge);
+    if (aCurve.GetType() != GeomAbs_Line)
+      continue; // only linear edges are to be merged
+
+    // Start a new chain
+    TopTools_ListOfShape aChain;
+    aChain.Append(aFirstEdge.Oriented(TopAbs_FORWARD));
+    makeEdgesChain(aChain, aMapVE, theAngTol);
+
+    // Remove it from the left edges
+    TopTools_ListIteratorOfListOfShape itL(aChain);
+    for (; itL.More(); itL.Next())
+    {
+      const TopoDS_Shape& aE = itL.Value();
+      Standard_Integer aEI = anEdges.FindIndex(aE);
+      anEdgesInd.Remove(aEI);
+    }
+
+    if (aChain.Extent() > 0)
+    {
+      // Create a new curvilinear edge from the polyline
+      TopoDS_Edge aNewEdge = createCurveEdge(aChain, aMapEF, theLinTol);
+      if (!aNewEdge.IsNull())
+      {
+        aNewEdge.Orientation(aFirstEdge.Orientation());
+        // Replace edges from the chain with the new edge.
+        Standard_Integer nbMerged = 0;
+        for (itL.Init(aChain); itL.More(); itL.Next(), nbMerged++)
+        {
+          const TopoDS_Shape& aE = itL.Value();
+          if (nbMerged == 0)
+            theContext->Replace(aE, aNewEdge);
+          else
+            theContext->Remove(aE);
+        }
+      }
+    }
+  }
+}
+
+//=======================================================================
+//function : Build
+//purpose  : 
+//=======================================================================
+void ShapeUpgrade_CombineToCylinder::Build()
+{
+  myDone = Standard_False;
+  myNbNewFaces = myNbReplacedFaces = 0;
+  myShape.Nullify();
+
+  TopTools_IndexedDataMapOfShapeListOfShape aMapEF;
+  TopExp::MapShapesAndAncestors(myInitShape, TopAbs_EDGE, TopAbs_FACE, aMapEF);
+
+  // Sort edges in the shape on two groups - stop edges and smooth edges that are 
+  // to be located on a cylinder wall parallel to its axis.
+  TopTools_MapOfShape aStopEdges;
+  NCollection_DataMap<TopoDS_Shape, gp_Cylinder> aMapECyl;
+  sortEdges(myInitShape, aMapEF, myLinTol, myAngTol, aStopEdges, aMapECyl);
+
+#ifdef DEBUG_CombineToCylinder_StopEdges
+  TopoDS_Compound aCompE;
+  BRep_Builder().MakeCompound(aCompE);
+  TopTools_MapIteratorOfMapOfShape aMapIt(aStopEdges);
+  for (; aMapIt.More(); aMapIt.Next())
+    BRep_Builder().Add(aCompE, aMapIt.Value());
+  BRepTools::Write(aCompE, "stope.brep");
+#endif
+
+  // Collect faces forming same-cylinder areas
+  NCollection_List<CombineToCylinder::Area> anAreas;
+  collectAreas(myInitShape, aMapEF, aStopEdges, aMapECyl, myLinTol, myAngTol, anAreas);
+
+#ifdef DEBUG_CombineToCylinder_Areas
+  TopoDS_Compound aComp;
+  BRep_Builder().MakeCompound(aComp);
+  TopTools_ListIteratorOfListOfShape anItA(anAreas);
+  for (; anItA.More(); anItA.Next())
+    BRep_Builder().Add(aComp, anItA.Value());
+  BRepTools::Write(aComp, "areas.brep");
+#endif
+
+  // Process each area
+  NCollection_List<CombineToCylinder::Area>::Iterator anItAreas(anAreas);
+  for (; anItAreas.More(); anItAreas.Next())
+  {
+    const CombineToCylinder::Area& anArea = anItAreas.Value();
+    if (!anArea.pCylinder)
+      continue;
+
+    TopoDS_Shape aCylFace = createCylFaceFromArea(anArea.Faces, *anArea.pCylinder);
+
+    // Replace faces from the list with the new face.
+    TopExp_Explorer anEx(anArea.Faces, TopAbs_FACE);
+    Standard_Integer nbMerged = 0;
+    for (; anEx.More(); anEx.Next(), nbMerged++)
+    {
+      const TopoDS_Shape& aFace = anEx.Current();
+      if (nbMerged == 0)
+        myContext->Replace(aFace, aCylFace);
+      else
+        myContext->Remove(aFace);
+    }
+    myNbNewFaces++;
+    myNbReplacedFaces += nbMerged;
+  }
+
+  // Build intermediate result shape using replacement context.
+  // This shape will contain chains of linear edges along cylinders
+  // instead of circles or whatever curves.
+  // The next step is to merge them.
+  TopoDS_Shape aTempShape = myContext->Apply(myInitShape);
+  mergeEdges(aTempShape, myLinTol, myAngTol, myContext);
+
+  // Rebuild the result shape using replacement context.
+  myShape = myContext->Apply(myInitShape);
+
+  // build pcurves
+  ShapeFix_Shape aFix(myShape);
+  aFix.Perform();
+  myShape = aFix.Shape();
+  myDone = Standard_True;
+}
diff --git a/src/ShapeUpgrade/ShapeUpgrade_CombineToCylinder.hxx b/src/ShapeUpgrade/ShapeUpgrade_CombineToCylinder.hxx
new file mode 100644 (file)
index 0000000..4b37b86
--- /dev/null
@@ -0,0 +1,93 @@
+//-Copyright: Open CASCADE 2016
+// Created on: 2016-08-01
+// Created by: Mikhail Sazonov
+// Copyright (c) 2016 OPEN CASCADE SAS
+//
+// This file is part of Open CASCADE Technology software library.
+//
+// This library is free software; you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License version 2.1 as published
+// by the Free Software Foundation, with special exception defined in the file
+// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
+// distribution for complete text of the license and disclaimer of any warranty.
+//
+// Alternatively, this file may be used under the terms of Open CASCADE
+// commercial license or contractual agreement.
+
+#ifndef _ShapeUpgrade_CombineToCylinder_HeaderFile
+#define _ShapeUpgrade_CombineToCylinder_HeaderFile
+
+#include <Standard_TypeDef.hxx>
+#include <TopoDS_Shape.hxx>
+
+class BRepTools_ReShape;
+
+//! This tool tries to merge planar faces constituting a cylindrical surface into one face.
+//! The algorithm checks that all merged faces have smooth change of normal within the given
+//! angular tolerance.
+class ShapeUpgrade_CombineToCylinder
+{
+public:
+  
+  //! empty constructor
+  Standard_EXPORT ShapeUpgrade_CombineToCylinder();
+  
+  //! Sets input shape
+  void SetShape(const TopoDS_Shape& theShape)
+  {
+    myInitShape = theShape;
+  }
+  
+  //! Sets the linear tolerance. Default value is Precision::Confusion().
+  void SetLinearTolerance(const Standard_Real theValue)
+  {
+    myLinTol = theValue;
+  }
+
+  //! Sets the angular tolerance. Default value is Precision::Angular().
+  void SetAngularTolerance(const Standard_Real theValue)
+  {
+    myAngTol = theValue;
+  }
+
+  //! Builds the resulting shape
+  Standard_EXPORT void Build();
+  
+  //! Returns true if the result shape was successfully built.
+  Standard_Boolean IsDone() const
+  {
+    return myDone;
+  }
+
+  //! Returns the result shape.
+  const TopoDS_Shape& Shape() const
+  {
+    return myShape;
+  }
+
+  //! Returns the number of new faces built by merging
+  Standard_Integer NbNewFaces() const
+  {
+    return myNbNewFaces;
+  }
+
+  //! Returns the number of replaced faces
+  Standard_Integer NbReplacedFaces() const
+  {
+    return myNbReplacedFaces;
+  }
+
+private:
+
+  TopoDS_Shape myInitShape;
+  Standard_Real myLinTol;
+  Standard_Real myAngTol;
+  TopoDS_Shape myShape;
+  Handle(BRepTools_ReShape) myContext;
+  Standard_Integer myNbNewFaces;
+  Standard_Integer myNbReplacedFaces;
+  Standard_Boolean myDone;
+
+};
+
+#endif // _ShapeUpgrade_CombineToCylinder_HeaderFile