Move to Eigen sparse matrix for coherency coefficients.

Progress indication.

src/BRepMesh/BRepMesh_Block.hxx

@@ -0,0 +1,270 @@
+// Created on: 2015-11-27
+// Created by: Danila ULYANOV
+// Copyright (c) 2015 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 _BRepMesh_Block_HeaderFile
+#define _BRepMesh_Block_HeaderFile
+
+#include <stdlib.h>
+#include <stdexcept>
+
+template <class Type> class Storage
+{
+public:
+  Storage (int theSize)
+  {
+    m_first = m_last = NULL;
+    m_BlockSize = theSize;
+  }
+
+  ~Storage()
+  {
+    while (m_first)
+    {
+      Block* next = m_first->Next;
+      delete[] ( (char*) m_first);
+      m_first = next;
+    }
+  }
+
+  Type* New (int aNum = 1)
+  {
+    Type* aResult;
+
+    if (!m_last || m_last->Current + aNum > m_last->Last)
+    {
+      if (m_last && m_last->Next)
+      {
+        m_last = m_last->Next;
+      }
+
+      else
+      {
+        Block* aNext = (Block*) new char [sizeof (Block) + (m_BlockSize - 1) *sizeof (Type)];
+
+        if (!aNext)
+        {
+          throw std::runtime_error ("Not enough memory!");
+        }
+
+        if (m_last)
+        {
+          m_last->Next = aNext;
+        }
+
+        else
+        {
+          m_first = aNext;
+        }
+
+        m_last = aNext;
+        m_last->Current = & (m_last->Data[0]);
+        m_last->Last = m_last->Current + m_BlockSize;
+        m_last->Next = NULL;
+      }
+    }
+
+    aResult = m_last->Current;
+    m_last->Current += aNum;
+    return aResult;
+  }
+
+  Type* ScanFirst()
+  {
+    for (m_scanCurrentBlock = m_first; m_scanCurrentBlock; m_scanCurrentBlock = m_scanCurrentBlock->Next)
+    {
+      m_scanCurrentData = & (m_scanCurrentBlock->Data[0]);
+
+      if (m_scanCurrentData < m_scanCurrentBlock->Current)
+      {
+        return m_scanCurrentData ++;
+      }
+    }
+
+    return NULL;
+  }
+
+  Type* ScanNext()
+  {
+    while (m_scanCurrentData >= m_scanCurrentBlock->Current)
+    {
+      m_scanCurrentBlock = m_scanCurrentBlock->Next;
+
+      if (!m_scanCurrentBlock)
+      {
+        return NULL;
+      }
+
+      m_scanCurrentData = & (m_scanCurrentBlock->Data[0]);
+    }
+
+    return m_scanCurrentData++;
+  }
+
+  struct BlockIterator;
+  Type* ScanFirst (BlockIterator& anIter)
+  {
+    for (anIter.ScanCurrentBlock = m_first; anIter.ScanCurrentBlock; anIter.ScanCurrentBlock = anIter.ScanCurrentBlock->Next)
+    {
+      anIter.ScanCurrentData = & (anIter.ScanCurrentBlock->Data[0]);
+
+      if (anIter.ScanCurrentData < anIter.ScanCurrentBlock->Current)
+      {
+        return anIter.ScanCurrentData ++;
+      }
+    }
+
+    return NULL;
+  }
+  Type* ScanNext (BlockIterator& anIter)
+  {
+    while (anIter.ScanCurrentData >= anIter.ScanCurrentBlock->Current)
+    {
+      anIter.ScanCurrentBlock = anIter.ScanCurrentBlock->Next;
+
+      if (!anIter.ScanCurrentBlock)
+      {
+        return NULL;
+      }
+
+      anIter.ScanCurrentData = & (anIter.ScanCurrentBlock->Data[0]);
+    }
+
+    return anIter.ScanCurrentData ++;
+  }
+
+  void Reset()
+  {
+    Block* aBlock;
+
+    if (!m_first)
+    {
+      return;
+    }
+
+    for (aBlock = m_first; ; aBlock = aBlock->Next)
+    {
+      aBlock->Current = & (aBlock->Data[0]);
+
+      if (aBlock == m_last)
+      {
+        break;
+      }
+    }
+
+    m_last = m_first;
+  }
+
+private:
+
+  typedef struct BlockStruct
+  {
+    Type* Current, *Last;
+    struct BlockStruct* Next;
+    Type Data[1];
+  } Block;
+
+  int m_BlockSize;
+  Block* m_first;
+  Block* m_last;
+public:
+
+  struct BlockIterator
+  {
+    Block* ScanCurrentBlock;
+    Type* ScanCurrentData;
+  };
+
+private:
+  Block* m_scanCurrentBlock;
+  Type* m_scanCurrentData;
+};
+
+template <class Type> class DataBlock
+{
+public:
+  DataBlock (int size)
+  {
+    m_first = NULL;
+    m_firstFree = NULL;
+    m_blockSize = size;
+  }
+
+  ~DataBlock()
+  {
+    while (m_first)
+    {
+      Block* next = m_first->Next;
+      delete[] ( (char*) m_first);
+      m_first = next;
+    }
+  }
+
+  Type* New()
+  {
+    BlockItem* anItem;
+
+    if (!m_firstFree)
+    {
+      Block* next = m_first;
+      m_first = (Block*) new char [sizeof (Block) + (m_blockSize - 1) *sizeof (BlockItem)];
+
+      if (!m_first)
+      {
+        throw std::runtime_error ("Not enough memory!");
+      }
+
+      m_firstFree = & (m_first->Data[0]);
+
+      for (anItem = m_firstFree; anItem < m_firstFree + m_blockSize - 1; anItem++)
+      {
+        anItem->NextFree = anItem + 1;
+      }
+
+      anItem->NextFree = NULL;
+      m_first->Next = next;
+    }
+
+    anItem = m_firstFree;
+    m_firstFree = anItem->NextFree;
+    return (Type*) anItem;
+  }
+
+  void Delete (Type* t)
+  {
+    ( (BlockItem*) t)->NextFree = m_firstFree;
+    m_firstFree = (BlockItem*) t;
+  }
+
+private:
+
+  typedef union BlockItemStruct
+  {
+    Type Val;
+    BlockItemStruct* NextFree;
+  } BlockItem;
+
+  typedef struct BlockStruct
+  {
+    struct BlockStruct* Next;
+    BlockItem Data[1];
+  } Block;
+
+  int m_blockSize;
+  Block* m_first;
+  BlockItem* m_firstFree;
+};
+
+#endif // _BRepMesh_Block_HeaderFile
+

src/BRepMesh/BRepMesh_EdgeSet.hxx

@@ -0,0 +1,98 @@
+// Created on: 2015-08-28
+// Created by: Danila ULYANOV
+// Copyright (c) 2015 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 _BRepMesh_EdgeSet_Header
+#define _BRepMesh_EdgeSet_Header
+
+#include <BVH_PrimitiveSet.hxx>
+#include <BVH_Types.hxx>
+
+//! Class representing set of boundary edges.
+template<class T, int N>
+class BRepMesh_EdgeSet : public BVH_PrimitiveSet<T, N>
+{
+public:
+
+  typedef typename BVH::VectorType<T, N>::Type BVH_VecNt;
+  typedef typename BVH::ArrayType<T, N>::Type BVH_ArrayNt;
+
+public:
+
+  //! Creates empty edge set.
+  BRepMesh_EdgeSet (const BVH_ArrayNt& theTriangulationVertices);
+
+  //! Releases resources of edge set.
+  virtual ~BRepMesh_EdgeSet();
+
+private:
+
+  BRepMesh_EdgeSet (const BRepMesh_EdgeSet& theOther);
+
+  BRepMesh_EdgeSet& operator= (const BRepMesh_EdgeSet& theOther);
+
+public:
+
+  //! Array of indices of edge vertices.
+  BVH_Array4i Elements;
+
+private:
+
+  const BVH_ArrayNt& Vertices;
+
+public:
+
+  //! Returns index of edge which is closest to the given one.
+  Standard_Integer ClosestEdge (const BVH_VecNt& theEdgeV1,
+                                const BVH_VecNt& theEdgeV2,
+                                T* theDist = NULL);
+
+  //! Returns index of edge which is closest to the given one (accounting for edge orientation).
+  Standard_Integer CoherentEdge (const BVH_VecNt& theEdgeV1,
+                                 const BVH_VecNt& theEdgeV2,
+                                 T* theDist = NULL);
+
+  //! Returns the distance between to edges.
+  static T EdgeToEdgeDistance (const BVH_VecNt& theE1V1,
+                               const BVH_VecNt& theE1V2,
+                               const BVH_VecNt& theE2V1,
+                               const BVH_VecNt& theE2V2);
+
+  //! Returns the distance between to edges.
+  static T EdgeToBoxDistance (const BVH_VecNt& theEV1,
+                              const BVH_VecNt& theEV2,
+                              const BVH_VecNt& theMinPoint,
+                              const BVH_VecNt& theMaxPoint);
+
+  //! Returns total number of edges.
+  virtual Standard_Integer Size() const;
+
+  //! Returns AABB of entire set of objects.
+  using BVH_PrimitiveSet<T, N>::Box;
+
+  //! Returns AABB of the given edge.
+  virtual BVH_Box<T, N> Box (const Standard_Integer theIndex) const;
+
+  //! Returns centroid position along the given axis.
+  virtual T Center (const Standard_Integer theIndex,
+                    const Standard_Integer theAxis) const;
+
+  //! Performs transposing the two given edges in the set.
+  virtual void Swap (const Standard_Integer theIndex1,
+                     const Standard_Integer theIndex2);
+
+};
+
+#include <BRepMesh_EdgeSet.lxx>
+#endif // _BRepMesh_EdgeSet_Header

src/BRepMesh/BRepMesh_EdgeSet.lxx

@@ -0,0 +1,452 @@
+// Created on: 2015-08-28
+// Created by: Danila ULYANOV
+// Copyright (c) 2015 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 <BRepMesh_LineBoxDistance.hxx>
+
+// =======================================================================
+// function : BRepMesh_EdgeSet
+// purpose  :
+// =======================================================================
+template<class T, int N>
+BRepMesh_EdgeSet<T, N>::BRepMesh_EdgeSet (const BVH_ArrayNt& theTriangulationVertices)
+  : Vertices (theTriangulationVertices)
+{
+  //
+}
+
+// =======================================================================
+// function : ~BRepMesh_EdgeSet
+// purpose  :
+// =======================================================================
+template<class T, int N>
+BRepMesh_EdgeSet<T, N>::~BRepMesh_EdgeSet()
+{
+  //
+}
+
+// =======================================================================
+// function : Size
+// purpose  :
+// =======================================================================
+template<class T, int N>
+Standard_Integer BRepMesh_EdgeSet<T, N>::Size() const
+{
+  return BVH::Array<Standard_Integer, 4>::Size (Elements);
+}
+
+// =======================================================================
+// function : Box
+// purpose  :
+// =======================================================================
+template<class T, int N>
+BVH_Box<T, N> BRepMesh_EdgeSet<T, N>::Box (const Standard_Integer theIndex) const
+{
+  const BVH_Vec4i& anIndex = BVH::Array<Standard_Integer, 4>::Value (Elements, theIndex);
+
+  const BVH_VecNt& aPoint0 = BVH::Array<T, N>::Value (Vertices, anIndex.x());
+  const BVH_VecNt& aPoint1 = BVH::Array<T, N>::Value (Vertices, anIndex.y());
+  
+  BVH_VecNt aMinPoint = aPoint0;
+  BVH_VecNt aMaxPoint = aPoint0;
+
+  BVH::BoxMinMax<T, N>::CwiseMin (aMinPoint, aPoint1);
+  BVH::BoxMinMax<T, N>::CwiseMax (aMaxPoint, aPoint1);
+
+  return BVH_Box<T, N> (aMinPoint, aMaxPoint);
+}
+
+// =======================================================================
+// function : Center
+// purpose  :
+// =======================================================================
+template<class T, int N>
+T BRepMesh_EdgeSet<T, N>::Center (const Standard_Integer theIndex,
+                                  const Standard_Integer theAxis) const
+{
+  const BVH_Vec4i& anIndex = BVH::Array<Standard_Integer, 4>::Value (Elements, theIndex);
+
+  const BVH_VecNt& aPoint0 = BVH::Array<T, N>::Value (Vertices, anIndex.x());
+  const BVH_VecNt& aPoint1 = BVH::Array<T, N>::Value (Vertices, anIndex.y());
+
+  return ( BVH::VecComp<T, N>::Get (aPoint0, theAxis) +
+           BVH::VecComp<T, N>::Get (aPoint1, theAxis) ) * static_cast<T> (1.0 / 2.0);
+}
+
+// =======================================================================
+// function : Swap
+// purpose  :
+// =======================================================================
+template<class T, int N>
+void BRepMesh_EdgeSet<T, N>::Swap (const Standard_Integer theIndex1,
+                                    const Standard_Integer theIndex2)
+{
+  BVH_Vec4i& anIndices1 = BVH::Array<Standard_Integer, 4>::ChangeValue (Elements, theIndex1);
+  BVH_Vec4i& anIndices2 = BVH::Array<Standard_Integer, 4>::ChangeValue (Elements, theIndex2);
+
+  std::swap (anIndices1, anIndices2);
+}
+
+// =======================================================================
+// function : ClosestEdge
+// purpose  :
+// =======================================================================
+template<class T, int N>
+Standard_Integer BRepMesh_EdgeSet<T, N>::ClosestEdge (const BVH_VecNt& theEdgeV1,
+                                                      const BVH_VecNt& theEdgeV2,
+                                                      T* theDist)
+{
+  Standard_Integer aStack[32];
+
+  const NCollection_Handle<BVH_Tree<T, N> >& aBVH = BVH();
+
+  if (aBVH.IsNull())
+  {
+    return -1;
+  }
+
+  Standard_Integer aNode = 0;
+
+  T aMinDist = std::numeric_limits<T>::max();
+  Standard_Integer aClosestEdge = 0;
+
+  Standard_Integer aHead = -1; // stack head position
+
+  for (;;)
+  {
+    if (aBVH->IsOuter (aNode))
+    {
+      for (Standard_Integer anElemNb = aBVH->BegPrimitive (aNode); anElemNb <= aBVH->EndPrimitive (aNode); ++anElemNb)
+      {
+        const BVH_Vec4i& anEdge = Elements[anElemNb];
+        const BVH_VecNt& aVrt0 = Vertices[anEdge.x()];
+        const BVH_VecNt& aVrt1 = Vertices[anEdge.y()];
+        T aDist = EdgeToEdgeDistance (theEdgeV1, theEdgeV2, aVrt0, aVrt1);
+
+        if (aDist < aMinDist)
+        {
+          aMinDist = aDist;
+          aClosestEdge = anElemNb;
+        }
+      }
+
+      if (aHead < 0)
+      {
+        break;
+      }
+
+      aNode = aStack[aHead--];
+    }
+    else
+    {
+      T aLftDist = EdgeToBoxDistance (theEdgeV1,
+                                      theEdgeV2,
+                                      aBVH->MinPoint (aBVH->LeftChild (aNode)),
+                                      aBVH->MaxPoint (aBVH->LeftChild (aNode)));
+
+      T aRghDist = EdgeToBoxDistance (theEdgeV1,
+                                      theEdgeV2,
+                                      aBVH->MinPoint (aBVH->RightChild (aNode)),
+                                      aBVH->MaxPoint (aBVH->RightChild (aNode)));
+
+      const bool aProceedLft = aLftDist < aMinDist;
+      const bool aProceedRgh = aRghDist < aMinDist;
+
+      if (aProceedLft && aProceedRgh)
+      {
+        int aLeft = aBVH->LeftChild (aNode);
+        int aRight = aBVH->RightChild (aNode);
+
+        aNode = (aLftDist < aRghDist) ? aLeft : aRight;
+        aStack[++aHead] = (aRghDist <= aLftDist) ? aLeft : aRight;
+      }
+      else
+      {
+        if (aProceedLft || aProceedRgh)
+        {
+          aNode = aProceedLft ? aBVH->LeftChild (aNode) : aBVH->RightChild (aNode);
+        }
+        else
+        {
+          if (aHead < 0)
+          {
+            break;
+          }
+
+          aNode = aStack[aHead--];
+        }
+      }
+    }
+  }
+
+  if (theDist != NULL)
+  {
+    *theDist = aMinDist;
+  }
+
+  return aClosestEdge;
+}
+
+// =======================================================================
+// function : CoherentEdge
+// purpose  :
+// =======================================================================
+template<class T, int N>
+Standard_Integer BRepMesh_EdgeSet<T, N>::CoherentEdge (const BVH_VecNt& theEdgeV1,
+                                                       const BVH_VecNt& theEdgeV2,
+                                                       T* theDist)
+{
+  Standard_Integer aStack[32];
+
+  const NCollection_Handle<BVH_Tree<T, N> >& aBVH = BVH();
+
+  if (aBVH.IsNull())
+  {
+    return -1;
+  }
+
+  Standard_Integer aNode = 0;
+
+  T aMinDist = std::numeric_limits<T>::max();
+  Standard_Integer aClosestEdge = 0;
+
+  Standard_Integer aHead = -1; // stack head position
+
+  for (;;)
+  {
+    if (aBVH->IsOuter (aNode))
+    {
+      for (Standard_Integer anElemNb = aBVH->BegPrimitive (aNode); anElemNb <= aBVH->EndPrimitive (aNode); ++anElemNb)
+      {
+        const BVH_Vec4i& anEdge = Elements[anElemNb];
+        const BVH_VecNt& aVrt0 = Vertices[anEdge.x()];
+        const BVH_VecNt& aVrt1 = Vertices[anEdge.y()];
+        T aDist = EdgeToEdgeDistance (theEdgeV1, theEdgeV2, aVrt0, aVrt1);
+
+        // add scaled vertex distances
+        aDist += (1e-6 * (theEdgeV1 - aVrt0).Modulus());
+        aDist += (1e-6 * (theEdgeV1 - aVrt1).Modulus());
+        aDist += (1e-6 * (theEdgeV2 - aVrt0).Modulus());
+        aDist += (1e-6 * (theEdgeV2 - aVrt1).Modulus());
+
+        if (aDist < aMinDist)
+        {
+          aMinDist = aDist;
+          aClosestEdge = anElemNb;
+        }
+      }
+
+      if (aHead < 0)
+      {
+        break;
+      }
+
+      aNode = aStack[aHead--];
+    }
+    else
+    {
+      T aLftDist = EdgeToBoxDistance (theEdgeV1,
+                                      theEdgeV2,
+                                      aBVH->MinPoint (aBVH->LeftChild (aNode)),
+                                      aBVH->MaxPoint (aBVH->LeftChild (aNode)));
+
+      T aRghDist = EdgeToBoxDistance (theEdgeV1,
+                                      theEdgeV2,
+                                      aBVH->MinPoint (aBVH->RightChild (aNode)),
+                                      aBVH->MaxPoint (aBVH->RightChild (aNode)));
+
+      const bool aProceedLft = aLftDist < aMinDist;
+      const bool aProceedRgh = aRghDist < aMinDist;
+
+      if (aProceedLft && aProceedRgh)
+      {
+        int aLeft = aBVH->LeftChild (aNode);
+        int aRight = aBVH->RightChild (aNode);
+
+        aNode = (aLftDist < aRghDist) ? aLeft : aRight;
+        aStack[++aHead] = (aRghDist <= aLftDist) ? aLeft : aRight;
+      }
+      else
+      {
+        if (aProceedLft || aProceedRgh)
+        {
+          aNode = aProceedLft ? aBVH->LeftChild (aNode) : aBVH->RightChild (aNode);
+        }
+        else
+        {
+          if (aHead < 0)
+          {
+            break;
+          }
+
+          aNode = aStack[aHead--];
+        }
+      }
+    }
+  }
+
+  if (theDist != NULL)
+  {
+    *theDist = aMinDist;
+  }
+
+  return aClosestEdge;
+}
+
+#define PARALLEL_THRESOLD 1e-8
+
+// =======================================================================
+// function : EdgeToEdgeDistance
+// purpose  :
+// =======================================================================
+template<class T, int N>
+T BRepMesh_EdgeSet<T, N>::EdgeToEdgeDistance (const BVH_VecNt& theE1V1,
+                                              const BVH_VecNt& theE1V2,
+                                              const BVH_VecNt& theE2V1,
+                                              const BVH_VecNt& theE2V2)
+{
+  BVH_VecNt u = theE1V2 - theE1V1;
+  BVH_VecNt v = theE2V2 - theE2V1;
+  BVH_VecNt w = theE1V1 - theE2V1;
+  Standard_Real a = u.Dot (u);       // always >= 0
+  Standard_Real b = u.Dot (v);
+  Standard_Real c = v.Dot (v);       // always >= 0
+  Standard_Real d = u.Dot (w);
+  Standard_Real e = v.Dot (w);
+  Standard_Real D = a * c - b * b;    // always >= 0
+  Standard_Real sc, sN, sD = D;       // sc = sN / sD, default sD = D >= 0
+  Standard_Real tc, tN, tD = D;       // tc = tN / tD, default tD = D >= 0
+
+  // compute the line parameters of the two closest points
+  if (D < PARALLEL_THRESOLD)   // the lines are almost parallel
+  {
+    sN = 0.0;         // force using point P0 on segment S1
+    sD = 1.0;         // to prevent possible division by 0.0 later
+    tN = e;
+    tD = c;
+  }
+  else                   // get the closest points on the infinite lines
+  {
+    sN = (b * e - c * d);
+    tN = (a * e - b * d);
+
+    if (sN < 0.0)          // sc < 0 => the s=0 edge is visible
+    {
+      sN = 0.0;
+      tN = e;
+      tD = c;
+    }
+    else if (sN > sD)    // sc > 1  => the s=1 edge is visible
+    {
+      sN = sD;
+      tN = e + b;
+      tD = c;
+    }
+  }
+
+  if (tN < 0.0)              // tc < 0 => the t=0 edge is visible
+  {
+    tN = 0.0;
+
+    // recompute sc for this edge
+    if (-d < 0.0)
+    {
+      sN = 0.0;
+    }
+    else if (-d > a)
+    {
+      sN = sD;
+    }
+    else
+    {
+      sN = -d;
+      sD = a;
+    }
+  }
+  else if (tN > tD)        // tc > 1  => the t=1 edge is visible
+  {
+    tN = tD;
+
+    // recompute sc for this edge
+    if ((-d + b) < 0.0)
+    {
+      sN = 0;
+    }
+    else if ((-d + b) > a)
+    {
+      sN = sD;
+    }
+    else
+    {
+      sN = (-d +  b);
+      sD = a;
+    }
+  }
+
+  // finally do the division to get sc and tc
+  sc = (abs (sN) < PARALLEL_THRESOLD ? 0.0 : sN / sD);
+  tc = (abs (tN) < PARALLEL_THRESOLD ? 0.0 : tN / tD);
+
+  // get the difference of the two closest points
+  BVH_VecNt dP = w + (u * sc) - (v * tc);  // =  S1(sc) - S2(tc)
+
+  return dP.Modulus(); // return the closest distance
+}
+
+#undef PARALLEL_THRESOLD
+
+// =======================================================================
+// function : EdgeToBoxDistance
+// purpose  :
+// =======================================================================
+template<class T, int N>
+T BRepMesh_EdgeSet<T, N>::EdgeToBoxDistance (const BVH_VecNt& theEV1,
+                                             const BVH_VecNt& theEV2,
+                                             const BVH_VecNt& theMinPoint,
+                                             const BVH_VecNt& theMaxPoint)
+{
+  T anSqrDist = std::numeric_limits<T>::max();
+
+  BVH_VecNt aVClosest = theEV1;
+  aVClosest = aVClosest.cwiseMax (theMinPoint);
+  aVClosest = aVClosest.cwiseMin (theMaxPoint);
+
+  anSqrDist = (theEV1 - aVClosest).SquareModulus();
+
+  if (anSqrDist == (T) 0.0)
+  {
+    return anSqrDist;
+  }
+
+  aVClosest = theEV2;
+  aVClosest = aVClosest.cwiseMax (theMinPoint);
+  aVClosest = aVClosest.cwiseMin (theMaxPoint);
+
+  anSqrDist = std::min (anSqrDist, (theEV2 - aVClosest).SquareModulus());
+
+  if (anSqrDist == (T) 0.0)
+  {
+    return anSqrDist;
+  }
+
+  LineToBoxDistance<T, N> aDist;
+  LineToBoxDistance<Standard_Real, 3>::Result aResult = aDist (theEV1, theEV2 - theEV1, theMinPoint, theMaxPoint);
+
+  if (aResult.lineParameter > (T) 0.0 && aResult.lineParameter < (T) 1.0)
+  {
+    anSqrDist = std::min (anSqrDist, aResult.sqrDistance);
+  }
+
+  return sqrt (anSqrDist);
+}

src/BRepMesh/BRepMesh_LineBoxDistance.hxx

@@ -0,0 +1,605 @@
+// Created on: 2015-11-27
+// Created by: Danila ULYANOV
+// Copyright (c) 2015 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 _BRepMesh_LineBoxDistance_Header
+#define _BRepMesh_LineBoxDistance_Header
+
+#include <BVH_Types.hxx>
+
+template <class T, int N>
+class LineToBoxDistance
+{
+public:
+
+  typedef typename BVH::VectorType<T, N>::Type BVH_VecNt;
+
+public:
+  struct Result
+  {
+    T distance, sqrDistance;
+    T lineParameter;
+    BVH_VecNt closestPoint[2];
+  };
+
+  Result operator() (const BVH_VecNt& theOrigin,
+                     const BVH_VecNt& theDirect,
+                     const BVH_VecNt& theMinPoint,
+                     const BVH_VecNt& theMaxPoint);
+
+protected:
+  // Compute the distance and closest point between a line and an
+  // axis-aligned box whose center is the origin.  On input, 'point' is the
+  // line origin and 'direction' is the line direction.  On output, 'point'
+  // is the point on the box closest to the line.  The 'direction' is
+  // non-const to allow transforming the problem into the first octant.
+  void DoQuery (BVH_VecNt& point, BVH_VecNt& direction,
+                BVH_VecNt const& boxExtent, Result& result);
+
+private:
+  void Face (int i0, int i1, int i2, BVH_VecNt& pnt,
+             BVH_VecNt const& dir, BVH_VecNt const& PmE,
+             BVH_VecNt const& boxExtent, Result& result);
+
+  void CaseNoZeros (BVH_VecNt& pnt, BVH_VecNt const& dir,
+                    BVH_VecNt const& boxExtent, Result& result);
+
+  void Case0 (int i0, int i1, int i2, BVH_VecNt& pnt,
+              BVH_VecNt const& dir, BVH_VecNt const& boxExtent,
+              Result& result);
+
+  void Case00 (int i0, int i1, int i2, BVH_VecNt& pnt,
+               BVH_VecNt const& dir, BVH_VecNt const& boxExtent,
+               Result& result);
+
+  void Case000 (BVH_VecNt& pnt, BVH_VecNt const& boxExtent,
+                Result& result);
+};
+
+// =======================================================================
+// function :
+// purpose  :
+// =======================================================================
+template<class T, int N>
+typename LineToBoxDistance<T, N>::Result
+LineToBoxDistance<T, N>::operator() (const BVH_VecNt& theOrigin,
+                                     const BVH_VecNt& theDirect,
+                                     const BVH_VecNt& theMinPoint,
+                                     const BVH_VecNt& theMaxPoint)
+{
+  // Translate the line and box so that the box has center at the origin.
+  BVH_VecNt boxCenter, boxExtent;
+  boxCenter = (theMaxPoint + theMinPoint) * (T) 0.5;
+  boxExtent = (theMaxPoint - theMinPoint) * (T) 0.5;
+  BVH_VecNt point = theOrigin - boxCenter;
+  BVH_VecNt direction = theDirect;
+
+  Result result;
+  DoQuery (point, direction, boxExtent, result);
+
+  // Compute the closest point on the line.
+  result.closestPoint[0] =
+    theOrigin + theDirect * result.lineParameter;
+
+  // Compute the closest point on the box.
+  result.closestPoint[1] = boxCenter + point;
+  return result;
+}
+
+// =======================================================================
+// function : DoQuery
+// purpose  :
+// =======================================================================
+template<class T, int N>
+void LineToBoxDistance<T, N>::DoQuery (
+  BVH_VecNt& point, BVH_VecNt& direction,
+  BVH_VecNt const& boxExtent, Result& result)
+{
+  result.sqrDistance = (T)0;
+  result.lineParameter = (T)0;
+
+  // Apply reflections so that direction vector has nonnegative components.
+  bool reflect[3];
+
+  for (int i = 0; i < 3; ++i)
+  {
+    if (direction[i] < (T)0)
+    {
+      point[i] = -point[i];
+      direction[i] = -direction[i];
+      reflect[i] = true;
+    }
+    else
+    {
+      reflect[i] = false;
+    }
+  }
+
+  if (direction[0] > (T)0)
+  {
+    if (direction[1] > (T)0)
+    {
+      if (direction[2] > (T)0)  // (+,+,+)
+      {
+        CaseNoZeros (point, direction, boxExtent, result);
+      }
+      else  // (+,+,0)
+      {
+        Case0 (0, 1, 2, point, direction, boxExtent, result);
+      }
+    }
+    else
+    {
+      if (direction[2] > (T)0)  // (+,0,+)
+      {
+        Case0 (0, 2, 1, point, direction, boxExtent, result);
+      }
+      else  // (+,0,0)
+      {
+        Case00 (0, 1, 2, point, direction, boxExtent, result);
+      }
+    }
+  }
+  else
+  {
+    if (direction[1] > (T)0)
+    {
+      if (direction[2] > (T)0)  // (0,+,+)
+      {
+        Case0 (1, 2, 0, point, direction, boxExtent, result);
+      }
+      else  // (0,+,0)
+      {
+        Case00 (1, 0, 2, point, direction, boxExtent, result);
+      }
+    }
+    else
+    {
+      if (direction[2] > (T)0)  // (0,0,+)
+      {
+        Case00 (2, 0, 1, point, direction, boxExtent, result);
+      }
+      else  // (0,0,0)
+      {
+        Case000 (point, boxExtent, result);
+      }
+    }
+  }
+
+  // Undo the reflections applied previously.
+  for (int i = 0; i < 3; ++i)
+  {
+    if (reflect[i])
+    {
+      point[i] = -point[i];
+    }
+  }
+
+  result.distance = sqrt (result.sqrDistance);
+}
+
+// =======================================================================
+// function : Face
+// purpose  :
+// =======================================================================
+template<class T, int N>
+void LineToBoxDistance<T, N>::Face (int i0, int i1,
+                                    int i2, BVH_VecNt& pnt, BVH_VecNt const& dir,
+                                    BVH_VecNt const& PmE, BVH_VecNt const& boxExtent, Result& result)
+{
+  BVH_VecNt PpE;
+  T lenSqr, inv, tmp, param, t, delta;
+
+  PpE[i1] = pnt[i1] + boxExtent[i1];
+  PpE[i2] = pnt[i2] + boxExtent[i2];
+
+  if (dir[i0] * PpE[i1] >= dir[i1] * PmE[i0])
+  {
+    if (dir[i0] * PpE[i2] >= dir[i2] * PmE[i0])
+    {
+      // v[i1] >= -e[i1], v[i2] >= -e[i2] (distance = 0)
+      pnt[i0] = boxExtent[i0];
+      inv = ((T)1) / dir[i0];
+      pnt[i1] -= dir[i1] * PmE[i0] * inv;
+      pnt[i2] -= dir[i2] * PmE[i0] * inv;
+      result.lineParameter = -PmE[i0] * inv;
+    }
+    else
+    {
+      // v[i1] >= -e[i1], v[i2] < -e[i2]
+      lenSqr = dir[i0] * dir[i0] + dir[i2] * dir[i2];
+      tmp = lenSqr * PpE[i1] - dir[i1] * (dir[i0] * PmE[i0] +
+                                          dir[i2] * PpE[i2]);
+
+      if (tmp <= ((T)2)*lenSqr * boxExtent[i1])
+      {
+        t = tmp / lenSqr;
+        lenSqr += dir[i1] * dir[i1];
+        tmp = PpE[i1] - t;
+        delta = dir[i0] * PmE[i0] + dir[i1] * tmp + dir[i2] * PpE[i2];
+        param = -delta / lenSqr;
+        result.sqrDistance += PmE[i0] * PmE[i0] + tmp * tmp +
+                              PpE[i2] * PpE[i2] + delta * param;
+
+        result.lineParameter = param;
+        pnt[i0] = boxExtent[i0];
+        pnt[i1] = t - boxExtent[i1];
+        pnt[i2] = -boxExtent[i2];
+      }
+      else
+      {
+        lenSqr += dir[i1] * dir[i1];
+        delta = dir[i0] * PmE[i0] + dir[i1] * PmE[i1] + dir[i2] * PpE[i2];
+        param = -delta / lenSqr;
+        result.sqrDistance += PmE[i0] * PmE[i0] + PmE[i1] * PmE[i1]
+                              + PpE[i2] * PpE[i2] + delta * param;
+
+        result.lineParameter = param;
+        pnt[i0] = boxExtent[i0];
+        pnt[i1] = boxExtent[i1];
+        pnt[i2] = -boxExtent[i2];
+      }
+    }
+  }
+  else
+  {
+    if (dir[i0] * PpE[i2] >= dir[i2] * PmE[i0])
+    {
+      // v[i1] < -e[i1], v[i2] >= -e[i2]
+      lenSqr = dir[i0] * dir[i0] + dir[i1] * dir[i1];
+      tmp = lenSqr * PpE[i2] - dir[i2] * (dir[i0] * PmE[i0] +
+                                          dir[i1] * PpE[i1]);
+
+      if (tmp <= ((T)2)*lenSqr * boxExtent[i2])
+      {
+        t = tmp / lenSqr;
+        lenSqr += dir[i2] * dir[i2];
+        tmp = PpE[i2] - t;
+        delta = dir[i0] * PmE[i0] + dir[i1] * PpE[i1] + dir[i2] * tmp;
+        param = -delta / lenSqr;
+        result.sqrDistance += PmE[i0] * PmE[i0] + PpE[i1] * PpE[i1] +
+                              tmp * tmp + delta * param;
+
+        result.lineParameter = param;
+        pnt[i0] = boxExtent[i0];
+        pnt[i1] = -boxExtent[i1];
+        pnt[i2] = t - boxExtent[i2];
+      }
+      else
+      {
+        lenSqr += dir[i2] * dir[i2];
+        delta = dir[i0] * PmE[i0] + dir[i1] * PpE[i1] + dir[i2] * PmE[i2];
+        param = -delta / lenSqr;
+        result.sqrDistance += PmE[i0] * PmE[i0] + PpE[i1] * PpE[i1] +
+                              PmE[i2] * PmE[i2] + delta * param;
+
+        result.lineParameter = param;
+        pnt[i0] = boxExtent[i0];
+        pnt[i1] = -boxExtent[i1];
+        pnt[i2] = boxExtent[i2];
+      }
+    }
+    else
+    {
+      // v[i1] < -e[i1], v[i2] < -e[i2]
+      lenSqr = dir[i0] * dir[i0] + dir[i2] * dir[i2];
+      tmp = lenSqr * PpE[i1] - dir[i1] * (dir[i0] * PmE[i0] +
+                                          dir[i2] * PpE[i2]);
+
+      if (tmp >= (T)0)
+      {
+        // v[i1]-edge is closest
+        if (tmp <= ((T)2)*lenSqr * boxExtent[i1])
+        {
+          t = tmp / lenSqr;
+          lenSqr += dir[i1] * dir[i1];
+          tmp = PpE[i1] - t;
+          delta = dir[i0] * PmE[i0] + dir[i1] * tmp + dir[i2] * PpE[i2];
+          param = -delta / lenSqr;
+          result.sqrDistance += PmE[i0] * PmE[i0] + tmp * tmp +
+                                PpE[i2] * PpE[i2] + delta * param;
+
+          result.lineParameter = param;
+          pnt[i0] = boxExtent[i0];
+          pnt[i1] = t - boxExtent[i1];
+          pnt[i2] = -boxExtent[i2];
+        }
+        else
+        {
+          lenSqr += dir[i1] * dir[i1];
+          delta = dir[i0] * PmE[i0] + dir[i1] * PmE[i1]
+                  + dir[i2] * PpE[i2];
+          param = -delta / lenSqr;
+          result.sqrDistance += PmE[i0] * PmE[i0] + PmE[i1] * PmE[i1]
+                                + PpE[i2] * PpE[i2] + delta * param;
+
+          result.lineParameter = param;
+          pnt[i0] = boxExtent[i0];
+          pnt[i1] = boxExtent[i1];
+          pnt[i2] = -boxExtent[i2];
+        }
+
+        return;
+      }
+
+      lenSqr = dir[i0] * dir[i0] + dir[i1] * dir[i1];
+      tmp = lenSqr * PpE[i2] - dir[i2] * (dir[i0] * PmE[i0] +
+                                          dir[i1] * PpE[i1]);
+
+      if (tmp >= (T)0)
+      {
+        // v[i2]-edge is closest
+        if (tmp <= ((T)2)*lenSqr * boxExtent[i2])
+        {
+          t = tmp / lenSqr;
+          lenSqr += dir[i2] * dir[i2];
+          tmp = PpE[i2] - t;
+          delta = dir[i0] * PmE[i0] + dir[i1] * PpE[i1] + dir[i2] * tmp;
+          param = -delta / lenSqr;
+          result.sqrDistance += PmE[i0] * PmE[i0] + PpE[i1] * PpE[i1] +
+                                tmp * tmp + delta * param;
+
+          result.lineParameter = param;
+          pnt[i0] = boxExtent[i0];
+          pnt[i1] = -boxExtent[i1];
+          pnt[i2] = t - boxExtent[i2];
+        }
+        else
+        {
+          lenSqr += dir[i2] * dir[i2];
+          delta = dir[i0] * PmE[i0] + dir[i1] * PpE[i1]
+                  + dir[i2] * PmE[i2];
+          param = -delta / lenSqr;
+          result.sqrDistance += PmE[i0] * PmE[i0] + PpE[i1] * PpE[i1]
+                                + PmE[i2] * PmE[i2] + delta * param;
+
+          result.lineParameter = param;
+          pnt[i0] = boxExtent[i0];
+          pnt[i1] = -boxExtent[i1];
+          pnt[i2] = boxExtent[i2];
+        }
+
+        return;
+      }
+
+      // (v[i1],v[i2])-corner is closest
+      lenSqr += dir[i2] * dir[i2];
+      delta = dir[i0] * PmE[i0] + dir[i1] * PpE[i1] + dir[i2] * PpE[i2];
+      param = -delta / lenSqr;
+      result.sqrDistance += PmE[i0] * PmE[i0] + PpE[i1] * PpE[i1]
+                            + PpE[i2] * PpE[i2] + delta * param;
+
+      result.lineParameter = param;
+      pnt[i0] = boxExtent[i0];
+      pnt[i1] = -boxExtent[i1];
+      pnt[i2] = -boxExtent[i2];
+    }
+  }
+}
+
+// =======================================================================
+// function : CaseNoZeros
+// purpose  :
+// =======================================================================
+template<class T, int N>
+void LineToBoxDistance<T, N>::CaseNoZeros (
+  BVH_VecNt& pnt, BVH_VecNt const& dir,
+  BVH_VecNt const& boxExtent, Result& result)
+{
+  BVH_VecNt PmE = pnt - boxExtent;
+  T prodDxPy = dir[0] * PmE[1];
+  T prodDyPx = dir[1] * PmE[0];
+  T prodDzPx, prodDxPz, prodDzPy, prodDyPz;
+
+  if (prodDyPx >= prodDxPy)
+  {
+    prodDzPx = dir[2] * PmE[0];
+    prodDxPz = dir[0] * PmE[2];
+
+    if (prodDzPx >= prodDxPz)
+    {
+      // line intersects x = e0
+      Face (0, 1, 2, pnt, dir, PmE, boxExtent, result);
+    }
+    else
+    {
+      // line intersects z = e2
+      Face (2, 0, 1, pnt, dir, PmE, boxExtent, result);
+    }
+  }
+  else
+  {
+    prodDzPy = dir[2] * PmE[1];
+    prodDyPz = dir[1] * PmE[2];
+
+    if (prodDzPy >= prodDyPz)
+    {
+      // line intersects y = e1
+      Face (1, 2, 0, pnt, dir, PmE, boxExtent, result);
+    }
+    else
+    {
+      // line intersects z = e2
+      Face (2, 0, 1, pnt, dir, PmE, boxExtent, result);
+    }
+  }
+}
+
+// =======================================================================
+// function : Case0
+// purpose  :
+// =======================================================================
+template<class T, int N>
+void LineToBoxDistance<T, N>::Case0 (int i0, int i1,
+                                     int i2, BVH_VecNt& pnt, BVH_VecNt const& dir,
+                                     BVH_VecNt const& boxExtent, Result& result)
+{
+  T PmE0 = pnt[i0] - boxExtent[i0];
+  T PmE1 = pnt[i1] - boxExtent[i1];
+  T prod0 = dir[i1] * PmE0;
+  T prod1 = dir[i0] * PmE1;
+  T delta, invLSqr, inv;
+
+  if (prod0 >= prod1)
+  {
+    // line intersects P[i0] = e[i0]
+    pnt[i0] = boxExtent[i0];
+
+    T PpE1 = pnt[i1] + boxExtent[i1];
+    delta = prod0 - dir[i0] * PpE1;
+
+    if (delta >= (T)0)
+    {
+      invLSqr = ((T)1) / (dir[i0] * dir[i0] + dir[i1] * dir[i1]);
+      result.sqrDistance += delta * delta * invLSqr;
+      pnt[i1] = -boxExtent[i1];
+      result.lineParameter = - (dir[i0] * PmE0 + dir[i1] * PpE1) * invLSqr;
+    }
+    else
+    {
+      inv = ((T)1) / dir[i0];
+      pnt[i1] -= prod0 * inv;
+      result.lineParameter = -PmE0 * inv;
+    }
+  }
+  else
+  {
+    // line intersects P[i1] = e[i1]
+    pnt[i1] = boxExtent[i1];
+
+    T PpE0 = pnt[i0] + boxExtent[i0];
+    delta = prod1 - dir[i1] * PpE0;
+
+    if (delta >= (T)0)
+    {
+      invLSqr = ((T)1) / (dir[i0] * dir[i0] + dir[i1] * dir[i1]);
+      result.sqrDistance += delta * delta * invLSqr;
+      pnt[i0] = -boxExtent[i0];
+      result.lineParameter = - (dir[i0] * PpE0 + dir[i1] * PmE1) * invLSqr;
+    }
+    else
+    {
+      inv = ((T)1) / dir[i1];
+      pnt[i0] -= prod1 * inv;
+      result.lineParameter = -PmE1 * inv;
+    }
+  }
+
+  if (pnt[i2] < -boxExtent[i2])
+  {
+    delta = pnt[i2] + boxExtent[i2];
+    result.sqrDistance += delta * delta;
+    pnt[i2] = -boxExtent[i2];
+  }
+  else if (pnt[i2] > boxExtent[i2])
+  {
+    delta = pnt[i2] - boxExtent[i2];
+    result.sqrDistance += delta * delta;
+    pnt[i2] = boxExtent[i2];
+  }
+}
+
+// =======================================================================
+// function : Case00
+// purpose  :
+// =======================================================================
+template<class T, int N>
+void LineToBoxDistance<T, N>::Case00 (int i0, int i1,
+                                      int i2, BVH_VecNt& pnt, BVH_VecNt const& dir,
+                                      BVH_VecNt const& boxExtent, Result& result)
+{
+  T delta;
+
+  result.lineParameter = (boxExtent[i0] - pnt[i0]) / dir[i0];
+
+  pnt[i0] = boxExtent[i0];
+
+  if (pnt[i1] < -boxExtent[i1])
+  {
+    delta = pnt[i1] + boxExtent[i1];
+    result.sqrDistance += delta * delta;
+    pnt[i1] = -boxExtent[i1];
+  }
+  else if (pnt[i1] > boxExtent[i1])
+  {
+    delta = pnt[i1] - boxExtent[i1];
+    result.sqrDistance += delta * delta;
+    pnt[i1] = boxExtent[i1];
+  }
+
+  if (pnt[i2] < -boxExtent[i2])
+  {
+    delta = pnt[i2] + boxExtent[i2];
+    result.sqrDistance += delta * delta;
+    pnt[i2] = -boxExtent[i2];
+  }
+  else if (pnt[i2] > boxExtent[i2])
+  {
+    delta = pnt[i2] - boxExtent[i2];
+    result.sqrDistance += delta * delta;
+    pnt[i2] = boxExtent[i2];
+  }
+}
+
+// =======================================================================
+// function : Case000
+// purpose  :
+// =======================================================================
+template<class T, int N>
+void LineToBoxDistance<T, N>::Case000 (
+  BVH_VecNt& pnt, BVH_VecNt const& boxExtent, Result& result)
+{
+  T delta;
+
+  if (pnt[0] < -boxExtent[0])
+  {
+    delta = pnt[0] + boxExtent[0];
+    result.sqrDistance += delta * delta;
+    pnt[0] = -boxExtent[0];
+  }
+  else if (pnt[0] > boxExtent[0])
+  {
+    delta = pnt[0] - boxExtent[0];
+    result.sqrDistance += delta * delta;
+    pnt[0] = boxExtent[0];
+  }
+
+  if (pnt[1] < -boxExtent[1])
+  {
+    delta = pnt[1] + boxExtent[1];
+    result.sqrDistance += delta * delta;
+    pnt[1] = -boxExtent[1];
+  }
+  else if (pnt[1] > boxExtent[1])
+  {
+    delta = pnt[1] - boxExtent[1];
+    result.sqrDistance += delta * delta;
+    pnt[1] = boxExtent[1];
+  }
+
+  if (pnt[2] < -boxExtent[2])
+  {
+    delta = pnt[2] + boxExtent[2];
+    result.sqrDistance += delta * delta;
+    pnt[2] = -boxExtent[2];
+  }
+  else if (pnt[2] > boxExtent[2])
+  {
+    delta = pnt[2] - boxExtent[2];
+    result.sqrDistance += delta * delta;
+    pnt[2] = boxExtent[2];
+  }
+}
+
+#endif // _BRepMesh_LineBoxDistance_Header

src/BRepMesh/BRepMesh_MinStCut.cxx

@@ -0,0 +1,676 @@
+// Created on: 2015-11-27
+// Created by: Danila ULYANOV
+// Copyright (c) 2015 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 <stdio.h>
+#include <string.h>
+#include <stdexcept>
+
+#include <BRepMesh_MinStCut.hxx>
+
+GraphEval::GraphEval (const int theNodeNumMax)
+  : m_nodeNum (0), m_nodeNumMax (theNodeNumMax)
+{
+  m_nodes = (Node*) malloc (m_nodeNumMax * sizeof (Node));
+  m_edgeBlock  = new Storage<Edge> (EDGE_BLOCK_SIZE);
+  m_flow = 0;
+}
+
+GraphEval::~GraphEval()
+{
+  free (m_nodes);
+  delete m_edgeBlock;
+}
+
+GraphEval::NodeId GraphEval::AddNode (int theNum)
+{
+  NodeId i = m_nodeNum;
+  m_nodeNum += theNum;
+
+  if (m_nodeNum > m_nodeNumMax)
+  {
+    throw std::runtime_error ("Error: the number of nodes is exceeded!");
+  }
+
+  memset (m_nodes + i, 0, theNum * sizeof (NodeStruct));
+
+  return i;
+}
+
+void GraphEval::AddEdge (NodeId theNodeFrom, NodeId theNodeTo, CapacityType theCap, CapacityType theRevCap)
+{
+  Edge* a, *a_rev;
+  Node* from = m_nodes + theNodeFrom;
+  Node* to = m_nodes + theNodeTo;
+
+  a = m_edgeBlock->New (2);
+  a_rev = a + 1;
+
+  a->Sister = a_rev;
+  a_rev->Sister = a;
+  a->Next = from->First;
+  from->First = a;
+  a_rev->Next = ( (Node*) to)->First;
+  to->First = a_rev;
+  a->Head = to;
+  a_rev->Head = from;
+  a->ReverseCap = theCap;
+  a_rev->ReverseCap = theRevCap;
+}
+
+void GraphEval::AddTWeights (NodeId theNode, CapacityType theCapSource, CapacityType theCapSink)
+{
+  register CapacityType delta = m_nodes[theNode].TrCap;
+
+  if (delta > 0)
+  {
+    theCapSource += delta;
+  }
+
+  else
+  {
+    theCapSink   -= delta;
+  }
+
+  m_flow += (theCapSource < theCapSink) ? theCapSource : theCapSink;
+  m_nodes[theNode].TrCap = theCapSource - theCapSink;
+}
+
+#define TERM ( (Edge *) 1 )
+#define ORPH ( (Edge *) 2 )
+
+#define INFINITE 1000000000
+
+inline void GraphEval::setActive (Node* theNode)
+{
+  if (!theNode->Next)
+  {
+    if (m_queueLast[1])
+    {
+      m_queueLast[1]->Next = theNode;
+    }
+
+    else
+    {
+      m_queueFirst[1] = theNode;
+    }
+
+    m_queueLast[1] = theNode;
+    theNode->Next = theNode;
+  }
+}
+
+inline GraphEval::Node* GraphEval::nextActive()
+{
+  Node* aNode;
+
+  while (1)
+  {
+    if (! (aNode = m_queueFirst[0]))
+    {
+      m_queueFirst[0] = aNode = m_queueFirst[1];
+      m_queueLast[0]  = m_queueLast[1];
+      m_queueFirst[1] = NULL;
+      m_queueLast[1]  = NULL;
+
+      if (!aNode)
+      {
+        return NULL;
+      }
+    }
+
+    if (aNode->Next == aNode)
+    {
+      m_queueFirst[0] = m_queueLast[0] = NULL;
+    }
+
+    else
+    {
+      m_queueFirst[0] = aNode->Next;
+    }
+
+    aNode->Next = NULL;
+
+    if (aNode->Parent)
+    {
+      return aNode;
+    }
+  }
+}
+
+void GraphEval::maxflowInit()
+{
+  Node* aNode;
+
+  m_queueFirst[0] = m_queueLast[0] = NULL;
+  m_queueFirst[1] = m_queueLast[1] = NULL;
+  m_orphanFirst = NULL;
+
+  for (aNode = m_nodes; aNode < m_nodes + m_nodeNum; aNode++)
+  {
+    aNode->Next = NULL;
+    aNode->TS = 0;
+
+    if (aNode->TrCap > 0)
+    {
+      aNode->IsSink = 0;
+      aNode->Parent = TERM;
+      setActive (aNode);
+      aNode->TS = 0;
+      aNode->DIST = 1;
+    }
+
+    else if (aNode->TrCap < 0)
+    {
+      aNode->IsSink = 1;
+      aNode->Parent = TERM;
+      setActive (aNode);
+      aNode->TS = 0;
+      aNode->DIST = 1;
+    }
+
+    else
+    {
+      aNode->Parent = NULL;
+    }
+  }
+
+  m_time = 0;
+}
+
+void GraphEval::augment (Edge* theMiddleEdge)
+{
+  Node* aNode;
+  Edge* anEdge;
+  TCapacityType aBottleneck;
+  NodePtr* aNodePtr;
+
+  aBottleneck = theMiddleEdge->ReverseCap;
+
+  for (aNode = theMiddleEdge->Sister->Head; ; aNode = anEdge->Head)
+  {
+    anEdge = aNode->Parent;
+
+    if (anEdge == TERM)
+    {
+      break;
+    }
+
+    if (aBottleneck > anEdge->Sister->ReverseCap)
+    {
+      aBottleneck = anEdge->Sister->ReverseCap;
+    }
+  }
+
+  if (aBottleneck > aNode->TrCap)
+  {
+    aBottleneck = aNode->TrCap;
+  }
+
+  for (aNode = theMiddleEdge->Head; ; aNode = anEdge->Head)
+  {
+    anEdge = aNode->Parent;
+
+    if (anEdge == TERM)
+    {
+      break;
+    }
+
+    if (aBottleneck > anEdge->ReverseCap)
+    {
+      aBottleneck = anEdge->ReverseCap;
+    }
+  }
+
+  if (aBottleneck > - aNode->TrCap)
+  {
+    aBottleneck = - aNode->TrCap;
+  }
+
+
+  theMiddleEdge->Sister->ReverseCap += aBottleneck;
+  theMiddleEdge->ReverseCap -= aBottleneck;
+
+  for (aNode = theMiddleEdge->Sister->Head; ; aNode = anEdge->Head)
+  {
+    anEdge = aNode->Parent;
+
+    if (anEdge == TERM)
+    {
+      break;
+    }
+
+    anEdge->ReverseCap += aBottleneck;
+    anEdge->Sister->ReverseCap -= aBottleneck;
+
+    if (!anEdge->Sister->ReverseCap)
+    {
+      aNode->Parent = ORPH;
+      aNodePtr = m_nodePtrBlock->New();
+      aNodePtr->Ptr = aNode;
+      aNodePtr->Next = m_orphanFirst;
+      m_orphanFirst = aNodePtr;
+    }
+  }
+
+  aNode->TrCap -= aBottleneck;
+
+  if (!aNode->TrCap)
+  {
+    aNode->Parent = ORPH;
+    aNodePtr = m_nodePtrBlock->New();
+    aNodePtr->Ptr = aNode;
+    aNodePtr->Next = m_orphanFirst;
+    m_orphanFirst = aNodePtr;
+  }
+
+  for (aNode = theMiddleEdge->Head; ; aNode = anEdge->Head)
+  {
+    anEdge = aNode->Parent;
+
+    if (anEdge == TERM)
+    {
+      break;
+    }
+
+    anEdge->Sister->ReverseCap += aBottleneck;
+    anEdge->ReverseCap -= aBottleneck;
+
+    if (!anEdge->ReverseCap)
+    {
+      aNode->Parent = ORPH;
+      aNodePtr = m_nodePtrBlock->New();
+      aNodePtr->Ptr = aNode;
+      aNodePtr->Next = m_orphanFirst;
+      m_orphanFirst = aNodePtr;
+    }
+  }
+
+  aNode->TrCap += aBottleneck;
+
+  if (!aNode->TrCap)
+  {
+    aNode->Parent = ORPH;
+    aNodePtr = m_nodePtrBlock->New();
+    aNodePtr->Ptr = aNode;
+    aNodePtr->Next = m_orphanFirst;
+    m_orphanFirst = aNodePtr;
+  }
+
+
+  m_flow += aBottleneck;
+}
+
+void GraphEval::processSourceOrphan (Node* theNode)
+{
+  Node* aNode;
+  Edge* anEdge0, *anEdge0Min = NULL, *anEdge;
+  NodePtr* aNodePtr;
+  int d, d_min = INFINITE;
+
+  for (anEdge0 = theNode->First; anEdge0; anEdge0 = anEdge0->Next)
+    if (anEdge0->Sister->ReverseCap)
+    {
+      aNode = anEdge0->Head;
+
+      if (!aNode->IsSink && (anEdge = aNode->Parent))
+      {
+        d = 0;
+
+        while (1)
+        {
+          if (aNode->TS == m_time)
+          {
+            d += aNode->DIST;
+            break;
+          }
+
+          anEdge = aNode->Parent;
+          d ++;
+
+          if (anEdge == TERM)
+          {
+            aNode->TS = m_time;
+            aNode->DIST = 1;
+            break;
+          }
+
+          if (anEdge == ORPH)
+          {
+            d = INFINITE;
+            break;
+          }
+
+          aNode = anEdge->Head;
+        }
+
+        if (d < INFINITE)
+        {
+          if (d < d_min)
+          {
+            anEdge0Min = anEdge0;
+            d_min = d;
+          }
+
+          for (aNode = anEdge0->Head; aNode->TS != m_time; aNode = aNode->Parent->Head)
+          {
+            aNode->TS = m_time;
+            aNode->DIST = d --;
+          }
+        }
+      }
+    }
+
+  if (theNode->Parent = anEdge0Min)
+  {
+    theNode->TS = m_time;
+    theNode->DIST = d_min + 1;
+  }
+
+  else
+  {
+    theNode->TS = 0;
+
+    for (anEdge0 = theNode->First; anEdge0; anEdge0 = anEdge0->Next)
+    {
+      aNode = anEdge0->Head;
+
+      if (!aNode->IsSink && (anEdge = aNode->Parent))
+      {
+        if (anEdge0->Sister->ReverseCap)
+        {
+          setActive (aNode);
+        }
+
+        if (anEdge != TERM && anEdge != ORPH && anEdge->Head == theNode)
+        {
+          aNode->Parent = ORPH;
+          aNodePtr = m_nodePtrBlock->New();
+          aNodePtr->Ptr = aNode;
+
+          if (m_orphanLast)
+          {
+            m_orphanLast->Next = aNodePtr;
+          }
+
+          else
+          {
+            m_orphanFirst        = aNodePtr;
+          }
+
+          m_orphanLast = aNodePtr;
+          aNodePtr->Next = NULL;
+        }
+      }
+    }
+  }
+}
+
+void GraphEval::processSinkOrphan (Node* theNode)
+{
+  Node* aNode;
+  Edge* anEdge0, *anEdge0Min = NULL, *anEdge;
+  NodePtr* aNodePtr;
+  int d, d_min = INFINITE;
+
+  for (anEdge0 = theNode->First; anEdge0; anEdge0 = anEdge0->Next)
+    if (anEdge0->ReverseCap)
+    {
+      aNode = anEdge0->Head;
+
+      if (aNode->IsSink && (anEdge = aNode->Parent))
+      {
+        d = 0;
+
+        while (1)
+        {
+          if (aNode->TS == m_time)
+          {
+            d += aNode->DIST;
+            break;
+          }
+
+          anEdge = aNode->Parent;
+          d ++;
+
+          if (anEdge == TERM)
+          {
+            aNode->TS = m_time;
+            aNode->DIST = 1;
+            break;
+          }
+
+          if (anEdge == ORPH)
+          {
+            d = INFINITE;
+            break;
+          }
+
+          aNode = anEdge->Head;
+        }
+
+        if (d < INFINITE)
+        {
+          if (d < d_min)
+          {
+            anEdge0Min = anEdge0;
+            d_min = d;
+          }
+
+          for (aNode = anEdge0->Head; aNode->TS != m_time; aNode = aNode->Parent->Head)
+          {
+            aNode->TS = m_time;
+            aNode->DIST = d --;
+          }
+        }
+      }
+    }
+
+  if (theNode->Parent = anEdge0Min)
+  {
+    theNode->TS = m_time;
+    theNode->DIST = d_min + 1;
+  }
+
+  else
+  {
+    theNode->TS = 0;
+
+    for (anEdge0 = theNode->First; anEdge0; anEdge0 = anEdge0->Next)
+    {
+      aNode = anEdge0->Head;
+
+      if (aNode->IsSink && (anEdge = aNode->Parent))
+      {
+        if (anEdge0->ReverseCap)
+        {
+          setActive (aNode);
+        }
+
+        if (anEdge != TERM && anEdge != ORPH && anEdge->Head == theNode)
+        {
+          aNode->Parent = ORPH;
+          aNodePtr = m_nodePtrBlock->New();
+          aNodePtr->Ptr = aNode;
+
+          if (m_orphanLast)
+          {
+            m_orphanLast->Next = aNodePtr;
+          }
+
+          else
+          {
+            m_orphanFirst        = aNodePtr;
+          }
+
+          m_orphanLast = aNodePtr;
+          aNodePtr->Next = NULL;
+        }
+      }
+    }
+  }
+}
+
+GraphEval::FlowType GraphEval::MaximumFlow()
+{
+  Node* i, *j, *current_node = NULL;
+  Edge* a;
+  NodePtr* np, *np_next;
+
+  maxflowInit();
+  m_nodePtrBlock = new DataBlock<NodePtr> (NODEPTR_BLOCK_SIZE);
+
+  while (1)
+  {
+    if (i = current_node)
+    {
+      i->Next = NULL;
+
+      if (!i->Parent)
+      {
+        i = NULL;
+      }
+    }
+
+    if (!i)
+    {
+      if (! (i = nextActive()))
+      {
+        break;
+      }
+    }
+
+    if (!i->IsSink)
+    {
+      for (a = i->First; a; a = a->Next)
+        if (a->ReverseCap)
+        {
+          j = a->Head;
+
+          if (!j->Parent)
+          {
+            j->IsSink = 0;
+            j->Parent = a->Sister;
+            j->TS = i->TS;
+            j->DIST = i->DIST + 1;
+            setActive (j);
+          }
+
+          else if (j->IsSink)
+          {
+            break;
+          }
+
+          else if (j->TS <= i->TS &&
+                   j->DIST > i->DIST)
+          {
+            j->Parent = a->Sister;
+            j->TS = i->TS;
+            j->DIST = i->DIST + 1;
+          }
+        }
+    }
+
+    else
+    {
+      for (a = i->First; a; a = a->Next)
+        if (a->Sister->ReverseCap)
+        {
+          j = a->Head;
+
+          if (!j->Parent)
+          {
+            j->IsSink = 1;
+            j->Parent = a->Sister;
+            j->TS = i->TS;
+            j->DIST = i->DIST + 1;
+            setActive (j);
+          }
+
+          else if (!j->IsSink)
+          {
+            a = a->Sister;
+            break;
+          }
+
+          else if (j->TS <= i->TS &&
+                   j->DIST > i->DIST)
+          {
+            j->Parent = a->Sister;
+            j->TS = i->TS;
+            j->DIST = i->DIST + 1;
+          }
+        }
+    }
+
+    m_time ++;
+
+    if (a)
+    {
+      i->Next = i;
+      current_node = i;
+
+      augment (a);
+
+      while (np = m_orphanFirst)
+      {
+        np_next = np->Next;
+        np->Next = NULL;
+
+        while (np = m_orphanFirst)
+        {
+          m_orphanFirst = np->Next;
+          i = np->Ptr;
+          m_nodePtrBlock->Delete (np);
+
+          if (!m_orphanFirst)
+          {
+            m_orphanLast = NULL;
+          }
+
+          if (i->IsSink)
+          {
+            processSinkOrphan (i);
+          }
+
+          else
+          {
+            processSourceOrphan (i);
+          }
+        }
+
+        m_orphanFirst = np_next;
+      }
+    }
+
+    else
+    {
+      current_node = NULL;
+    }
+  }
+
+  delete m_nodePtrBlock;
+
+  return m_flow;
+}
+
+GraphEval::TerminalType GraphEval::Label (NodeId theNode)
+{
+  if (m_nodes[theNode].Parent && ! (m_nodes[theNode].IsSink))
+  {
+    return SOURCE;
+  }
+
+  return SINK;
+}

src/BRepMesh/BRepMesh_MinStCut.hxx

@@ -0,0 +1,182 @@
+// Created on: 2015-11-27
+// Created by: Danila ULYANOV
+// Copyright (c) 2015 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 _BRepMesh_MinStCut_HeaderFile
+#define _BRepMesh_MinStCut_HeaderFile
+
+#include <cassert>
+
+#include <BRepMesh_Block.hxx>
+
+#define EDGE_BLOCK_SIZE 1024
+#define NODEPTR_BLOCK_SIZE 128
+
+class GraphEval
+{
+public:
+  typedef enum
+  {
+    SOURCE = 0,
+    SINK = 1
+  } TerminalType;
+
+  typedef int NodeId;
+  typedef double CapacityType;
+  typedef double FlowType;
+  typedef double TCapacityType;
+
+  GraphEval (const int theNodeNumMax);
+
+  ~GraphEval();
+
+  NodeId AddNode (int theNum = 1);
+  void AddEdge (NodeId theNodeFrom, NodeId theNodeTo, CapacityType theCap, CapacityType theRevCap);
+  void AddTWeights (NodeId theNode, CapacityType theCapSource, CapacityType theCapSink);
+  TerminalType Label (NodeId theNode);
+  FlowType MaximumFlow();
+
+private:
+  struct EdgeStruct;
+
+  typedef struct NodeStruct
+  {
+    EdgeStruct* First;
+
+    EdgeStruct* Parent;
+    NodeStruct* Next;
+    int TS;
+    int DIST;
+    short IsSink;
+
+    CapacityType TrCap;
+  } Node;
+
+  typedef struct EdgeStruct
+  {
+    NodeStruct* Head;
+    EdgeStruct* Next;
+    EdgeStruct* Sister;
+
+    CapacityType ReverseCap;
+  } Edge;
+
+  typedef struct NodePtrStruct
+  {
+    NodeStruct* Ptr;
+    NodePtrStruct* Next;
+  } NodePtr;
+
+  NodeStruct* m_nodes;
+  int m_nodeNum, m_nodeNumMax;
+  Storage<Edge>* m_edgeBlock;
+  DataBlock<NodePtr>* m_nodePtrBlock;
+
+  FlowType m_flow;
+
+  Node* m_queueFirst[2], *m_queueLast[2];
+  NodePtr* m_orphanFirst, *m_orphanLast;
+  int m_time;
+
+  void setActive (Node* theNode);
+  Node* nextActive();
+
+  void maxflowInit();
+  void augment (Edge* theMiddleEdge);
+  void processSourceOrphan (Node* theNode);
+  void processSinkOrphan (Node* theNode);
+};
+
+class Energy : GraphEval
+{
+public:
+  typedef NodeId Var;
+
+  typedef CapacityType Value;
+  typedef FlowType TotalValue;
+
+  Energy (const int theMaxNodes);
+
+  ~Energy();
+
+  Var AddVariable();
+
+  void AddConstant (Value E);
+
+  void AddTerm (Var x,
+                Value E0, Value E1);
+
+  void AddPairwise (Var x, Var y,
+                    Value E00, Value E01,
+                    Value E10, Value E11);
+
+  TotalValue Minimize();
+
+  int Label (Var x);
+
+};
+
+inline Energy::Energy (const int theMaxNodes)
+  : GraphEval (theMaxNodes)
+{
+}
+
+inline Energy::~Energy() {}
+
+inline Energy::Var Energy::AddVariable()
+{
+  return AddNode();
+}
+
+inline void Energy::AddTerm (Var x,
+                             Value A, Value B)
+{
+  AddTWeights (x, B, A);
+}
+
+inline void Energy::AddPairwise (Var x, Var y,
+                                 Value A, Value B,
+                                 Value C, Value D)
+{
+  AddTWeights (x, D, A);
+  B -= A;
+  C -= D;
+
+  assert (B + C >= 0);
+
+  if (B < 0)
+  {
+    AddTWeights (x, 0, B);
+    AddTWeights (y, 0, -B);
+    AddEdge (x, y, 0, B + C);
+  }
+
+  else if (C < 0)
+  {
+    AddTWeights (x, 0, -C);
+    AddTWeights (y, 0, C);
+    AddEdge (x, y, B + C, 0);
+  }
+
+  else
+  {
+    AddEdge (x, y, B, C);
+  }
+}
+
+inline Energy::TotalValue Energy::Minimize() { return MaximumFlow(); }
+
+inline int Energy::Label (Var x) { return (int) GraphEval::Label (x); }
+
+#endif // _BRepMesh_MinStCut_HeaderFile

src/BRepMesh/BRepMesh_RestoreOrientationTool.cxx

@@ -0,0 +1,954 @@
+// Created on: 2015-08-27
+// Created by: Danila ULYANOV
+// Copyright (c) 2015 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 <BRepMesh_RestoreOrientationTool.hxx>
+
+#include <OSD_Timer.hxx>
+#include <TCollection_AsciiString.hxx>
+#include <TopExp_Explorer.hxx>
+#include <TopoDS_Compound.hxx>
+
+#include <BRep_Builder.hxx>
+#include <BRep_Tool.hxx>
+#include <BRepAdaptor_Surface.hxx>
+#include <Poly_Triangulation.hxx>
+
+#include <NCollection_Map.hxx>
+#include <NCollection_DataMap.hxx>
+#include <Message_ProgressSentry.hxx>
+#include <Standard_Integer.hxx>
+
+#include <Precision.hxx>
+
+#include <math_BullardGenerator.hxx>
+
+#include <set>
+#include <map>
+#include <memory>
+
+#include <BRepMesh_MinStCut.hxx>
+
+#include <Eigen/Core>
+#include <Eigen/Sparse>
+
+// #define MEASURE_PERFORMANCE
+
+IMPLEMENT_STANDARD_HANDLE (BRepMesh_RestoreOrientationTool, Standard_Transient)
+IMPLEMENT_STANDARD_RTTIEXT(BRepMesh_RestoreOrientationTool, Standard_Transient)
+
+namespace
+{
+
+// The maximum positive double value.
+const Standard_Real THE_MAX_VAL = std::numeric_limits<Standard_Real>::max();
+
+// The 'small' (safe for inversion) 4D vector of doubles.
+const BVH_Vec3d THE_MIN_VEC (1e-100, 1e-100, 1e-100);
+
+// Relative weight of coherency component for optimization.
+const Standard_Real COHERENCY_WEIGHT = 50.0;
+
+// Minimum number of traced rays per face.
+const Standard_Integer MIN_FACE_RAYS = 50;
+
+// Maximum number of traced rays per face.
+const Standard_Integer MAX_FACE_RAYS = 2000;
+
+// Indicates that ray does not hit any geometry.
+const Standard_Integer INVALID_HIT = -1;
+
+//! Samples value from array according to probabilities.
+Standard_Integer sampleValue (Standard_Real theKsi, const std::vector<Standard_Real>& theCDF)
+{
+  return static_cast<Standard_Integer> (std::lower_bound (theCDF.begin(), theCDF.end(), theKsi) - theCDF.begin());
+}
+
+//! Utility class to manage coherence table.
+class SparseCoherencyTable
+{
+public:
+
+  //! Allocates square coherency table with theMaxSize * theMaxSize maximum number of elements.
+  SparseCoherencyTable (Standard_Size theMaxSize)
+    : myMaxSize (theMaxSize),
+      mySize (0),
+      myTable ((Standard_Integer)theMaxSize,
+               (Standard_Integer)theMaxSize)
+  {
+    //
+  }
+
+  void Reserve (Standard_Size theNumElements)
+  {
+    myTable.reserve ((Standard_Integer)theNumElements);
+  }
+
+  //! Sets coherence value.
+  void setCoherence (Standard_Size theI, Standard_Size theJ, Standard_Real theValue)
+  {
+    if (theI < theJ) std::swap (theI, theJ);
+    myTable.coeffRef ((Standard_Integer)theI, (Standard_Integer)theJ) = static_cast<Standard_ShortReal> (theValue);
+  }
+
+  //! Returns coherence value.
+  Standard_Real getCoherence (Standard_Size theI, Standard_Size theJ)
+  {
+    if (theI < theJ) std::swap (theI, theJ);
+    return static_cast<Standard_Real> (myTable.coeff ((Standard_Integer)theI, (Standard_Integer)theJ));
+  }
+
+  //! Returns actual table size.
+  Standard_Size Size()
+  {
+    return mySize;
+  }
+
+  //! Sets actual table size.
+  void SetSize (Standard_Size theSize)
+  {
+    mySize = theSize;
+  }
+
+private:
+
+  Standard_Size myMaxSize;
+  Standard_Size mySize;
+
+  Eigen::SparseMatrix<Standard_ShortReal> myTable;
+};
+
+//! Utility class to manage coherence table.
+class CoherencyTable
+{
+public:
+
+  //! Allocates square coherency table with theMaxSize * theMaxSize maximum number of elements.
+  CoherencyTable (Standard_Size theMaxSize)
+    : myMaxSize (theMaxSize),
+      mySize (0),
+      myTable (theMaxSize * theMaxSize, 0.f)
+  {
+    //
+  }
+
+  //! Sets coherence value.
+  void setCoherence (Standard_Size theI, Standard_Size theJ, Standard_Real theValue)
+  {
+    if (theI < theJ) std::swap (theI, theJ);
+    myTable[theI * myMaxSize + theJ] = static_cast<Standard_ShortReal> (theValue);
+  }
+
+  //! Returns coherence value.
+  Standard_Real getCoherence (Standard_Size theI, Standard_Size theJ)
+  {
+    if (theI < theJ) std::swap (theI, theJ);
+    return static_cast<Standard_Real> (myTable[theI * myMaxSize + theJ]);
+  }
+
+  //! Returns actual table size.
+  Standard_Size Size()
+  {
+    return mySize;
+  }
+
+  //! Sets actual table size.
+  void SetSize (Standard_Size theSize)
+  {
+    mySize = theSize;
+  }
+
+private:
+
+  Standard_Size myMaxSize;
+  Standard_Size mySize;
+
+  std::vector<Standard_ShortReal> myTable;
+};
+
+//! Returns time of intersection of ray and triangle.
+Standard_Real rayTriangleIntersection (const BVH_Vec3d& theOrigin,
+                                       const BVH_Vec3d& theDirect,
+                                       const BVH_Vec3d& thePoint0,
+                                       const BVH_Vec3d& thePoint1,
+                                       const BVH_Vec3d& thePoint2)
+{
+  const BVH_Vec3d  aEdge0  = thePoint1 - thePoint0;
+  const BVH_Vec3d  aEdge1  = thePoint0 - thePoint2;
+
+  BVH_Vec3d        aNormal = BVH_Vec3d::Cross (aEdge1, aEdge0);
+  const Standard_Real aDot = aNormal.Dot (theDirect);
+  if (fabs (aDot) < 1.0e-16)
+  {
+    return THE_MAX_VAL; // ray in the plane of triangle
+  }
+
+  const BVH_Vec3d aEdge2 = (thePoint0 - theOrigin) * (1.0 / aDot);
+
+  const Standard_Real aTime = aNormal.Dot (aEdge2);
+  if (aTime < 0.0)
+  {
+    return THE_MAX_VAL;
+  }
+
+  const BVH_Vec3d theInc = BVH_Vec3d::Cross (theDirect, aEdge2);
+
+  Standard_Real U = theInc.Dot (aEdge1);
+  Standard_Real V = theInc.Dot (aEdge0);
+
+  if ((U >= 0.0) && (V >= 0.0) && (U + V <= 1.0))
+  {
+    return aTime;
+  }
+
+  return THE_MAX_VAL;
+}
+
+//! Checks if two boxes (AABBs) are overlapped.
+bool overlapBoxes (const BVH_Vec3d& theBoxMin1,
+                   const BVH_Vec3d& theBoxMax1,
+                   const BVH_Vec3d& theBoxMin2,
+                   const BVH_Vec3d& theBoxMax2) {
+
+  return !(theBoxMin1.x() > theBoxMax2.x()
+        || theBoxMax1.x() < theBoxMin2.x()
+        || theBoxMin1.y() > theBoxMax2.y()
+        || theBoxMax1.y() < theBoxMin2.y()
+        || theBoxMin1.z() > theBoxMax2.z()
+        || theBoxMax1.z() < theBoxMin2.z());
+}
+
+}
+
+// =======================================================================
+// function : BRepMesh_RestoreOrientationTool
+// purpose  :
+// =======================================================================
+BRepMesh_RestoreOrientationTool::BRepMesh_RestoreOrientationTool (const bool theVisibilityOnly /* false */)
+  : myIsDone (false),
+    myVisibilityOnly (theVisibilityOnly)
+{
+  //
+}
+
+// =======================================================================
+// function : BRepMesh_RestoreOrientationTool
+// purpose  :
+// =======================================================================
+BRepMesh_RestoreOrientationTool::BRepMesh_RestoreOrientationTool (const TopoDS_Shape& theShape,
+                                                                  const bool theVisibilityOnly /* false */)
+  : myIsDone (false),
+    myVisibilityOnly (theVisibilityOnly)
+{
+  Init (theShape);
+}
+
+// =======================================================================
+// function : Init
+// purpose  :
+// =======================================================================
+void BRepMesh_RestoreOrientationTool::Init (const TopoDS_Shape& theShape)
+{
+  if (theShape.IsNull())
+  {
+    Standard_NullObject::Raise();
+  }
+
+  myFaceList.clear();
+
+  for (TopExp_Explorer anIter (theShape, TopAbs_FACE); anIter.More(); anIter.Next())
+  {
+    myFaceList.push_back (static_cast<const TopoDS_Face&> (anIter.Current()));
+  }
+
+  myTriangulation.Elements.clear();
+  myTriangulation.Vertices.clear();
+
+  Standard_Integer aCurrentTriangle = 0;
+
+  for (BRepMesh_FaceList::iterator anIter = myFaceList.begin(); anIter != myFaceList.end(); ++anIter)
+  {
+    Handle (BRepMesh_TriangulatedPatch) aNewPatch (new BRepMesh_TriangulatedPatch(*anIter));
+
+    Standard_Integer aTriangleCount = extractTriangulation (*anIter, *aNewPatch);
+
+    if (aTriangleCount == 0)
+    {
+      continue;
+    }
+
+    aCurrentTriangle += aTriangleCount;
+    myFaceIdIntervals.push_back (aCurrentTriangle - 1);
+    aNewPatch->PreprocessData();
+
+    myPatches.push_back (aNewPatch);
+
+    extractTriangulation (*anIter, myTriangulation);
+  }
+}
+
+// =======================================================================
+// function : extractTriangulation
+// purpose  :
+// =======================================================================
+Standard_Integer BRepMesh_RestoreOrientationTool::extractTriangulation (const TopoDS_Face& theFace,
+                                                                        BVH_Triangulation<Standard_Real, 3>& theTriangulation)
+{
+  TopLoc_Location aLocation;
+
+  Handle (Poly_Triangulation) aPolyTriangulation = BRep_Tool::Triangulation (theFace, aLocation);
+
+  if (aPolyTriangulation.IsNull())
+  {
+    return 0;
+  }
+
+  bool isMirrored = aLocation.Transformation().VectorialPart().Determinant() < 0;
+  bool needToFlip = static_cast<bool> (isMirrored ^ (theFace.Orientation() == TopAbs_REVERSED));
+
+  BRepAdaptor_Surface aFaceAdaptor (theFace, Standard_False);
+
+  const Standard_Integer aVertOffset =
+    static_cast<Standard_Integer> (theTriangulation.Vertices.size()) - 1;
+
+  for (Standard_Integer aVertIdx = 1; aVertIdx <= aPolyTriangulation->NbNodes(); ++aVertIdx)
+  {
+    gp_Pnt aVertex = aPolyTriangulation->Nodes().Value (aVertIdx);
+
+    aVertex.Transform (aLocation.Transformation());
+
+    theTriangulation.Vertices.push_back (BVH_Vec3d (aVertex.X(),
+                                                    aVertex.Y(),
+                                                    aVertex.Z()));
+  }
+
+  for (Standard_Integer aTriIdx = 1; aTriIdx <= aPolyTriangulation->NbTriangles(); ++aTriIdx)
+  {
+    Standard_Integer aVertex1;
+    Standard_Integer aVertex2;
+    Standard_Integer aVertex3;
+
+    if (!needToFlip)
+    {
+      aPolyTriangulation->Triangles().Value (aTriIdx).Get (aVertex1,
+                                                           aVertex2,
+                                                           aVertex3);
+    }
+    else
+    {
+      aPolyTriangulation->Triangles().Value (aTriIdx).Get (aVertex3,
+                                                           aVertex2,
+                                                           aVertex1);
+    }
+
+    theTriangulation.Elements.push_back (BVH_Vec4i (aVertex1 + aVertOffset,
+                                                    aVertex2 + aVertOffset,
+                                                    aVertex3 + aVertOffset,
+                                                    1));
+  }
+
+  theTriangulation.MarkDirty(); // triangulation needs BVH rebuilding
+
+  return aPolyTriangulation->NbTriangles();
+}
+
+// =======================================================================
+// function : computeVisibility
+// purpose  :
+// =======================================================================
+void BRepMesh_RestoreOrientationTool::computeVisibility (BVH_Triangulation<Standard_Real, 3>& theTriangulation,
+                                                         Handle (BRepMesh_TriangulatedPatch)& thePatch,
+                                                         Standard_Integer theRaysNb)
+{
+  // deterministic sequence
+  math_BullardGenerator anRNG (0);
+
+  Standard_Real aFrontVisibility = 0.0;
+  Standard_Real aBackVisibility = 0.0;
+
+  Standard_Real aFrontDistance = 0.0;
+  Standard_Real aBackDistance = 0.0;
+
+  for (Standard_Integer anIdx = 0; anIdx < theRaysNb; ++anIdx)
+  {
+    // select triangle with probability proportional to triangle area
+    Standard_Integer aTrgIdx = sampleValue (anRNG.NextReal(), thePatch->TriangleCDF());
+
+    const BVH_Vec4i& aTriangle = thePatch->Elements[aTrgIdx];
+
+    // sample point in triangle uniformly
+    Standard_Real aKsi, aPsi;
+    do
+    {
+      aKsi = anRNG.NextReal();
+      aPsi = anRNG.NextReal();
+    } while (aKsi + aPsi > 1.0);
+
+    const BVH_Vec3d& aPoint0 = thePatch->Vertices[aTriangle.x()];
+    const BVH_Vec3d& aPoint1 = thePatch->Vertices[aTriangle.y()];
+    const BVH_Vec3d& aPoint2 = thePatch->Vertices[aTriangle.z()];
+
+    BVH_Vec3d aPoint (aPoint0 * (1 - aKsi - aPsi)
+                    + aPoint1 * aKsi
+                    + aPoint2 * aPsi);
+
+    const BVH_Vec3d aEdge0  = aPoint1 - aPoint0;
+    const BVH_Vec3d aEdge1  = aPoint0 - aPoint2;
+
+    const BVH_Vec3d aNormal = (BVH_Vec3d::Cross (aEdge1, aEdge0).Normalized());
+
+    // hemisphere direction
+    BVH_Vec3d aDirection;
+    Standard_Real aDot = 0.0;
+
+    do
+    {
+      aKsi = anRNG.NextReal();
+      aPsi = anRNG.NextReal();
+
+      Standard_Real anSqrtPsi = sqrt (aPsi);
+
+      aDirection = BVH_Vec3d (anSqrtPsi * cos (2.f * M_PI * aKsi),
+                              anSqrtPsi * sin (2.f * M_PI * aKsi),
+                              sqrt (1.f - aPsi));
+
+      aDot = aDirection.Dot (aNormal);
+
+    } while (std::abs (aDot) < 0.2);
+
+    if (aDot < 0.0)
+    {
+      aDirection = -aDirection;
+    }
+
+    const Standard_Real anEpsilon = std::max (1.0e-6, 1.0e-4 * theTriangulation.Box().Size().Modulus());
+
+    Standard_Real aDist = 0.0;
+    if (traceRay (aPoint + aNormal * anEpsilon, aDirection, theTriangulation, aDist) == INVALID_HIT)
+    {
+      aFrontVisibility += 1.0;
+    }
+    else
+    {
+      aFrontDistance += aDist;
+    }
+
+    if (traceRay (aPoint - aNormal * anEpsilon, -aDirection, theTriangulation, aDist) == INVALID_HIT)
+    {
+      aBackVisibility += 1.0;
+    }
+    else
+    {
+      aBackDistance += aDist;
+    }
+  }
+
+  Standard_Real aInvRayNb = 1.0 / theRaysNb;
+  thePatch->SetFrontVisibility (aFrontVisibility * aInvRayNb);
+  thePatch->SetBackVisibility (aBackVisibility * aInvRayNb);
+  thePatch->SetFrontDistance (aFrontDistance * aInvRayNb);
+  thePatch->SetBackDistance (aBackDistance * aInvRayNb);
+}
+
+// =======================================================================
+// function : traceRay
+// purpose  :
+// =======================================================================
+Standard_Integer BRepMesh_RestoreOrientationTool::traceRay (const BVH_Vec3d& theOrigin,
+                                                            const BVH_Vec3d& theDirect,
+                                                            BVH_Triangulation<Standard_Real, 3>& theTriangulation,
+                                                            Standard_Real& theDistance)
+{
+  Standard_Integer aStack[32];
+
+  const NCollection_Handle<BVH_Tree<Standard_Real, 3> >& aBVH = theTriangulation.BVH();
+
+  if (aBVH.IsNull())
+  {
+    return INVALID_HIT;
+  }
+
+  Standard_Integer aNode = 0;
+
+  Standard_Real aMinDist = std::numeric_limits<Standard_Real>::max();
+  Standard_Integer aHitTriangle = -1;
+
+  Standard_Integer aHead = -1; // stack head position
+
+  BVH_Vec3d anInvDirect = theDirect.cwiseAbs().cwiseMax (THE_MIN_VEC);
+  anInvDirect = BVH_Vec3d (1.0 / anInvDirect.x(), 1.0 / anInvDirect.y(), 1.0 / anInvDirect.z());
+
+  anInvDirect = BVH_Vec3d (_copysign (anInvDirect.x(), theDirect.x()),
+                           _copysign (anInvDirect.y(), theDirect.y()),
+                           _copysign (anInvDirect.z(), theDirect.z()));
+
+  for (;;)
+  {
+    if (aBVH->IsOuter (aNode))
+    {
+      for (Standard_Integer anElemNb = aBVH->BegPrimitive (aNode); anElemNb <= aBVH->EndPrimitive (aNode); ++anElemNb)
+      {
+        const BVH_Vec4i& aTriangle = theTriangulation.Elements[anElemNb];
+        const BVH_Vec3d& aVrt0 = theTriangulation.Vertices[aTriangle.x()];
+        const BVH_Vec3d& aVrt1 = theTriangulation.Vertices[aTriangle.y()];
+        const BVH_Vec3d& aVrt2 = theTriangulation.Vertices[aTriangle.z()];
+
+        Standard_Real aDist = rayTriangleIntersection (theOrigin, theDirect, aVrt0, aVrt1, aVrt2);
+
+
+        if (aDist < aMinDist)
+        {
+          aMinDist = aDist;
+          aHitTriangle = anElemNb;
+        }
+      }
+
+      if (aHead < 0)
+      {
+        break;
+      }
+
+      aNode = aStack[aHead--];
+    }
+    else
+    {
+      BVH_Vec3d aTime0 = (aBVH->MinPoint (aBVH->LeftChild (aNode)) - theOrigin) * anInvDirect;
+      BVH_Vec3d aTime1 = (aBVH->MaxPoint (aBVH->LeftChild (aNode)) - theOrigin) * anInvDirect;
+      
+      BVH_Vec3d aTimeMax = aTime0.cwiseMax (aTime1);
+      BVH_Vec3d aTimeMin = aTime0.cwiseMin (aTime1);
+
+      aTime0 = (aBVH->MinPoint (aBVH->RightChild (aNode)) - theOrigin) * anInvDirect;
+      aTime1 = (aBVH->MaxPoint (aBVH->RightChild (aNode)) - theOrigin) * anInvDirect;
+
+      Standard_Real aTimeFinal = aTimeMax.minComp();
+      Standard_Real aTimeStart = aTimeMin.maxComp();
+
+      aTimeMax = aTime0.cwiseMax (aTime1);
+      aTimeMin = aTime0.cwiseMin (aTime1);
+
+      Standard_Real aTimeMin1 = (aTimeStart <= aTimeFinal) && (aTimeFinal >= 0.0)
+                               ? aTimeStart : THE_MAX_VAL;
+
+      aTimeFinal = aTimeMax.minComp();
+      aTimeStart = aTimeMin.maxComp();
+
+      Standard_Real aTimeMin2 = (aTimeStart <= aTimeFinal) && (aTimeFinal >= 0.f)
+                               ? aTimeStart : THE_MAX_VAL;
+
+      const bool aHitLft = (aTimeMin1 != THE_MAX_VAL);
+      const bool aHitRgh = (aTimeMin2 != THE_MAX_VAL);
+
+      if (aHitLft && aHitRgh)
+      {
+        int aLeft = aBVH->LeftChild (aNode);
+        int aRight = aBVH->RightChild (aNode);
+
+        aNode = (aTimeMin1 < aTimeMin2) ? aLeft : aRight;
+        aStack[++aHead] = (aTimeMin2 <= aTimeMin1) ? aLeft : aRight;
+      }
+      else
+      {
+        if (aHitLft || aHitRgh)
+        {
+          aNode = aHitLft ? aBVH->LeftChild (aNode) : aBVH->RightChild (aNode);
+        }
+        else
+        {
+          if (aHead < 0)
+          {
+            break;
+          }
+
+          aNode = aStack[aHead--];
+        }
+      }
+    }
+  }
+
+  theDistance = aMinDist;
+
+  return aHitTriangle;
+}
+
+// =======================================================================
+// function : edgeCoherence
+// purpose  :
+// =======================================================================
+Standard_Real BRepMesh_RestoreOrientationTool::edgeCoherence (const BRepMesh_OrientedEdge& theEdge1,
+                                                              const BRepMesh_TriangulatedPatch&  thePatch1,
+                                                              const BRepMesh_OrientedEdge& theEdge2,
+                                                              const BRepMesh_TriangulatedPatch&  thePatch2)
+{
+  BVH_Vec3d anE1V1;
+  BVH_Vec3d anE1V2;
+  BVH_Vec3d anE2V1;
+  BVH_Vec3d anE2V2;
+
+  thePatch1.GetEdgeVertices (theEdge1, anE1V1, anE1V2);
+  thePatch2.GetEdgeVertices (theEdge2, anE2V1, anE2V2);
+
+  BVH_Vec3d anEdge1 = anE1V2 - anE1V1;
+  BVH_Vec3d anEdge2 = anE2V2 - anE2V1;
+
+  Standard_Real aScalarProduct = anEdge1.Dot (anEdge2);
+
+  Standard_Real aDist = BRepMesh_TriangulatedPatch::EdgeSet::EdgeToEdgeDistance (anE1V1, anE1V2, anE2V1, anE2V2);
+  Standard_Real aSign = -_copysign (1.0, aScalarProduct);
+
+  return aSign * sqrt (std::abs (aScalarProduct)) / (1.0 + aDist);
+}
+
+// =======================================================================
+// function : computeCoherence
+// purpose  :
+// =======================================================================
+Standard_Real BRepMesh_RestoreOrientationTool::computeCoherence (Handle (BRepMesh_TriangulatedPatch)& thePatch0,
+                                                                 Handle (BRepMesh_TriangulatedPatch)& thePatch1)
+{
+  Standard_Real aTotalCoherence = 0.0;
+
+  // Ensure that BVH already built before parallel section
+  thePatch0->BoundaryEdgeSet().BVH();
+  thePatch1->BoundaryEdgeSet().BVH();
+
+  #pragma omp parallel for
+  for (Standard_Integer anIdx = 0; anIdx < thePatch1->BoundaryEdgeSet().Size(); ++anIdx)
+  {
+    BRepMesh_OrientedEdge aCurrentEdge = thePatch1->BoundaryEdge (anIdx);
+    BVH_Vec3d aVertex1;
+    BVH_Vec3d aVertex2;
+    thePatch1->GetEdgeVertices (aCurrentEdge, aVertex1, aVertex2);
+    BRepMesh_OrientedEdge aClosestEdge = thePatch0->BoundaryEdge (
+      thePatch0->BoundaryEdgeSet().CoherentEdge (aVertex1, aVertex2));
+
+    #pragma omp critical
+    {
+      aTotalCoherence += edgeCoherence (aCurrentEdge, *thePatch1, aClosestEdge, *thePatch0);
+    }
+  }
+
+  return aTotalCoherence;
+}
+
+// =======================================================================
+// function : Perform
+// purpose  :
+// =======================================================================
+void BRepMesh_RestoreOrientationTool::Perform (const Handle(Message_ProgressIndicator) thePI)
+{
+  Standard_Real aMaxArea = std::numeric_limits<Standard_Real>::epsilon();
+
+#ifdef MEASURE_PERFORMANCE
+  OSD_Timer aTimer;
+  aTimer.Start();
+#endif
+
+  // Flipped flags for all patches.
+  std::vector<char> aFlipped (myPatches.size(), 0);
+  std::vector<char> aFlippedC (myPatches.size(), 0);
+
+  Standard_Real aMaxCoherence = -std::numeric_limits<Standard_Real>::max();
+
+  const Standard_Real anEpsilon = std::max (1.0e-4, 1.0e-2 * myTriangulation.Box().Size().Modulus());
+  BVH_Vec3d anEpsilonVec (anEpsilon, anEpsilon, anEpsilon);
+
+  std::unique_ptr<SparseCoherencyTable> aTable;
+
+  aTable->Reserve (myPatches.size());
+
+  if (!myVisibilityOnly)
+  {
+    aTable.reset (new SparseCoherencyTable (myPatches.size() * 3));
+    aTable->SetSize (myPatches.size());
+
+    if (!thePI.IsNull()) thePI->NewScope ( 80, "Coherence" );
+    Message_ProgressSentry aPS (thePI, "Compute initial coherence", 0.0, (Standard_Real) (myPatches.size() - 1), 1.0);
+
+    // Compute coherence
+    for (Standard_Size i = 0; i < myPatches.size() - 1; ++i, aPS.Next())
+    {
+      Handle (BRepMesh_TriangulatedPatch)& aTriPatch1 = myPatches[i];
+
+      for (Standard_Size j = i + 1; j < myPatches.size(); ++j)
+      {
+        Handle (BRepMesh_TriangulatedPatch)& aTriPatch2 = myPatches[j];
+
+        if (overlapBoxes (aTriPatch1->Box().CornerMin() - anEpsilonVec,
+                          aTriPatch1->Box().CornerMax() + anEpsilonVec,
+                          aTriPatch2->Box().CornerMin() - anEpsilonVec,
+                          aTriPatch2->Box().CornerMax() + anEpsilonVec))
+        {
+          Standard_Real aCoherence = computeCoherence (aTriPatch1, aTriPatch2) + computeCoherence (aTriPatch2, aTriPatch1);
+
+          aMaxCoherence = std::max (aMaxCoherence, aCoherence);
+
+          aTable->setCoherence (i, j, aCoherence);
+        }
+      }
+    }
+
+    if (!thePI.IsNull()) thePI->EndScope();
+
+    std::vector<BRepMesh_SuperPatch> aMetaPatches;
+    aMetaPatches.reserve (myPatches.size() * 3);
+
+    // Prepare meta patches
+    for (Standard_Size i = 0; i < myPatches.size(); ++i)
+    {
+      aMetaPatches.push_back (BRepMesh_SuperPatch (static_cast<Standard_Integer> (i)));
+    }
+
+    myPatchQueue.clear();
+    myPatchQueue.reserve (myPatches.size());
+
+    for (Standard_Size i = 0; i < myPatches.size() - 1; ++i)
+    {
+      for (Standard_Size j = i + 1; j < myPatches.size(); ++j)
+      {
+        Standard_Real aCoherence = aTable->getCoherence (i, j);
+        if (aCoherence != 0.0)
+        {
+          myPatchQueue.push_back (BRepMesh_SuperPatchPair (&aMetaPatches[i], &aMetaPatches[j], aTable->getCoherence (i, j)));
+        }
+      }
+    }
+
+    std::set<Standard_Integer> aRemovedSuperPatches;
+
+    // Greedy optimization of coherence
+    while (!myPatchQueue.empty())
+    {
+      std::sort (myPatchQueue.begin(), myPatchQueue.end());
+
+      const BRepMesh_SuperPatchPair& aBestPair = myPatchQueue.back();
+
+      if (aBestPair.Coherence < 0)
+      {
+        // Do flip of second patch
+        Standard_Integer anIndex = aBestPair.SecondPatch->Index;
+
+        for (Standard_Size j = 0; j < aBestPair.SecondPatch->PatchIds.size(); ++j)
+        {
+          Standard_Integer aPatchId = aBestPair.SecondPatch->PatchIds[j];
+
+          if (aPatchId == anIndex)
+          {
+            aFlippedC[aPatchId] = aFlippedC[aPatchId] == 0 ? 1 : 0;
+          }
+
+          aFlipped[aPatchId] = aFlipped[aPatchId] == 0 ? 1 : 0;
+        }
+
+        for (Standard_Size i = 0; i < aTable->Size(); ++i)
+        {
+          Standard_Real aCoherence = aTable->getCoherence (i, anIndex);
+
+          aTable->setCoherence (i, anIndex, -aCoherence);
+        }
+      }
+
+      aRemovedSuperPatches.insert (aBestPair.FirstPatch->Index);
+      aRemovedSuperPatches.insert (aBestPair.SecondPatch->Index);
+
+      aMetaPatches.push_back (BRepMesh_SuperPatch (*aBestPair.FirstPatch, *aBestPair.SecondPatch, (Standard_Integer)aTable->Size()));
+
+      const BRepMesh_SuperPatch& aNewPatch = aMetaPatches.back();
+
+      // Update coherence
+      for (Standard_Size i = 0; i < aTable->Size(); ++i)
+      {
+        if (aRemovedSuperPatches.find ((Standard_Integer)i) == aRemovedSuperPatches.end())
+        {
+          Standard_Real aNewCoherence = 0.0;
+          aNewCoherence += aTable->getCoherence (i, aBestPair.FirstPatch->Index);
+          aNewCoherence += aTable->getCoherence (i, aBestPair.SecondPatch->Index);
+
+          aTable->setCoherence (i, aNewPatch.Index, aNewCoherence);
+        }
+      }
+
+      std::vector<BRepMesh_SuperPatchPair> aNewPatchQueue;
+      aNewPatchQueue.reserve (myPatchQueue.size());
+
+      // Update pairs
+      for (Standard_Size i = 0; i < myPatchQueue.size(); ++i)
+      {
+        if (!myPatchQueue[i].HavePatch (aBestPair.FirstPatch->Index)
+         && !myPatchQueue[i].HavePatch (aBestPair.SecondPatch->Index))
+        {
+          // Accept a pair
+          aNewPatchQueue.push_back (myPatchQueue[i]);
+        }
+      }
+
+      myPatchQueue.swap (aNewPatchQueue);
+
+      // Add new pairs
+      for (Standard_Size i = 0; i < aTable->Size(); ++i)
+      {
+        if (aRemovedSuperPatches.find ((Standard_Integer)i) == aRemovedSuperPatches.end())
+        {
+          Standard_Real aCoherence = aTable->getCoherence (i, aNewPatch.Index);
+
+          if (aCoherence != 0.0)
+          {
+            myPatchQueue.push_back (BRepMesh_SuperPatchPair (&aNewPatch, &aMetaPatches[i], aCoherence));
+          }
+        }
+      }
+
+      aTable->SetSize (aTable->Size() + 1);
+    }
+
+#ifdef MEASURE_PERFORMANCE
+    aTimer.Stop();
+    std::cout << "Coherency time:  " << aTimer.ElapsedTime() << std::endl;
+    aTimer.Reset();
+    aTimer.Start();
+#endif
+  }
+
+  // Compute max area
+  for (Standard_Size i = 0; i < myPatches.size(); ++i)
+  {
+    aMaxArea = std::max (aMaxArea, myPatches[i]->TotalArea());
+  }
+
+  const Standard_Real anInvMaxArea = 1.0 / aMaxArea;
+
+  const Standard_Integer aPatchesNb = static_cast<Standard_Integer> (myPatches.size());
+
+  // Ensure that BVH already built before parallel section
+  myTriangulation.BVH();
+
+  if (!thePI.IsNull()) thePI->NewScope (myVisibilityOnly ? 100 : 20, "Visibility");
+  Message_ProgressSentry aPS (thePI, "Compute visibility", 0, aPatchesNb, 1);
+
+  // Compute visibility
+  #pragma omp parallel for
+  for (Standard_Integer i = 0; i < aPatchesNb; ++i)
+  {
+    Handle (BRepMesh_TriangulatedPatch)& aTriPatch = myPatches[i];
+    Standard_Integer aRaysNb = Max (MIN_FACE_RAYS, static_cast<Standard_Integer> (aTriPatch->TotalArea() * anInvMaxArea * MAX_FACE_RAYS));
+    computeVisibility (myTriangulation, aTriPatch, aRaysNb);
+
+    if (!myVisibilityOnly && aFlipped[i])
+    {
+      aTriPatch->FlipVisibility();
+    }
+
+    if (!thePI.IsNull())
+    {
+      #pragma omp critical
+      {
+        aPS.Next();
+      }
+    }
+  }
+
+  if (!thePI.IsNull()) thePI->EndScope();
+
+#ifdef MEASURE_PERFORMANCE
+  aTimer.Stop();
+  std::cout << "Visibility time:  " << aTimer.ElapsedTime() << std::endl;
+  aTimer.Reset();
+  aTimer.Start();
+#endif
+
+  // Optimization
+  Energy* aGraph = new Energy ((Standard_Integer)myPatches.size());
+
+  std::vector<Energy::Var> aVariables (myPatches.size());
+
+  for (Standard_Size i = 0; i < myPatches.size(); ++i)
+  {
+    Handle (BRepMesh_TriangulatedPatch)& aTriPatch = myPatches[i];
+    aVariables[i] = aGraph->AddVariable();
+
+    aGraph->AddTerm (aVariables[i], aTriPatch->BackVisibility() - aTriPatch->FrontVisibility(),
+                                    aTriPatch->FrontVisibility() - aTriPatch->BackVisibility());
+  }
+
+  if (!myVisibilityOnly)
+  {
+    const Standard_Real aInvMaxCoherence = 1.0 / aMaxCoherence;
+
+    for (Standard_Size i = 0; i < myPatches.size() - 1; ++i)
+    {
+      Handle (BRepMesh_TriangulatedPatch)& aTriPatch1 = myPatches[i];
+      for (Standard_Size j = i + 1; j < myPatches.size(); ++j)
+      {
+        Handle (BRepMesh_TriangulatedPatch)& aTriPatch2 = myPatches[j];
+
+        if (!overlapBoxes (aTriPatch1->Box().CornerMin() - anEpsilonVec,
+                           aTriPatch1->Box().CornerMax() + anEpsilonVec,
+                           aTriPatch2->Box().CornerMin() - anEpsilonVec,
+                           aTriPatch2->Box().CornerMax() + anEpsilonVec))
+        {
+          continue;
+        }
+
+        bool isCoherenceFlipped = (aFlipped[i] ^ aFlipped[j] ^ aFlippedC[i] ^ aFlippedC[j]) != 0;
+        Standard_Real aCoherence = aTable->getCoherence (i, j) * aInvMaxCoherence  * (isCoherenceFlipped ? -1.0 : 1.0);
+
+        if (aCoherence > 0.0)
+        {
+          aCoherence *= COHERENCY_WEIGHT;
+
+          // Prevent setting different labels for already coherent patches.
+          aGraph->AddPairwise (aVariables[i], aVariables[j], 0, aCoherence, aCoherence, 0);
+        }
+      }
+    }
+  }
+
+  try
+  {
+    aGraph->Minimize();
+  }
+  catch (std::exception& theError)
+  {
+    std::cout << theError.what() << std::endl;
+  }
+
+  for (Standard_Size i = 0; i < myPatches.size(); ++i)
+  {
+    if (aGraph->Label (aVariables[i]) == 1)
+    {
+      aFlipped[i] = aFlipped[i] == 0 ? 1 : 0;
+    }
+  }
+
+#ifdef MEASURE_PERFORMANCE
+  aTimer.Stop();
+  std::cout << "Optimization time:  " << aTimer.ElapsedTime() << std::endl;
+#endif
+
+  BRep_Builder aBuilder;
+  TopoDS_Compound aCompound;
+  aBuilder.MakeCompound (aCompound);
+  for (Standard_Size i = 0; i < myPatches.size(); ++i)
+  {
+    TopoDS_Face aFace = myPatches[i]->Face();
+    if (aFlipped[i])
+    {
+      aFace.Reverse();
+    }
+    aBuilder.Add (aCompound, aFace);
+  }
+
+  myShape = aCompound;
+  myIsDone = true;
+}

src/BRepMesh/BRepMesh_RestoreOrientationTool.hxx

@@ -0,0 +1,190 @@
+// Created on: 2015-08-27
+// Created by: Danila ULYANOV
+// Copyright (c) 2015 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 _BRepMesh_RestoreOrientationTool_HeaderFile
+#define _BRepMesh_RestoreOrientationTool_HeaderFile
+
+#include <Standard.hxx>
+#include <TopoDS_Shape.hxx>
+#include <TopoDS_Face.hxx>
+#include <Message_ProgressIndicator.hxx>
+
+#include <BRepMesh_TriangulatedPatch.hxx>
+
+//! Class representing a collection of simple patches.
+struct BRepMesh_SuperPatch
+{
+  //! Create super patch from single patch.
+  BRepMesh_SuperPatch (Standard_Integer thePatchId)
+    : Index (thePatchId)
+  {
+    PatchIds.push_back (thePatchId);
+  }
+
+  //! Create super patch from a pair of super patches.
+  BRepMesh_SuperPatch (const BRepMesh_SuperPatch& theFirst,
+                       const BRepMesh_SuperPatch& theSecond,
+                       Standard_Integer theIndex)
+    : Index (theIndex)
+  {
+    PatchIds.resize (theFirst.PatchIds.size() + theSecond.PatchIds.size());
+
+    std::vector<Standard_Integer>::iterator aNextPos =
+      std::copy (theFirst.PatchIds.begin(), theFirst.PatchIds.end(), PatchIds.begin());
+
+    std::copy (theSecond.PatchIds.begin(), theSecond.PatchIds.end(), aNextPos);
+  }
+
+  std::vector<Standard_Integer> PatchIds;
+  Standard_Integer Index;
+};
+
+//! Class representing a pair of super patches.
+struct BRepMesh_SuperPatchPair
+{
+  BRepMesh_SuperPatchPair (const BRepMesh_SuperPatch* theFirstPatch = NULL,
+                           const BRepMesh_SuperPatch* theSecondPatch  = NULL,
+                           Standard_Real theCoherency = 0.0)
+    : FirstPatch (theFirstPatch),
+      SecondPatch (theSecondPatch),
+      Coherence (theCoherency)
+  {
+    //
+  }
+
+  bool operator < (const BRepMesh_SuperPatchPair& theOther) const
+  {
+    return (std::abs (Coherence) < std::abs (theOther.Coherence));
+  }
+
+  bool HavePatch (Standard_Integer thePatch)
+  {
+    return FirstPatch->Index == thePatch || SecondPatch->Index == thePatch;
+  }
+
+  const BRepMesh_SuperPatch* FirstPatch;
+  const BRepMesh_SuperPatch* SecondPatch;
+  Standard_Real Coherence;
+};
+
+//! List of faces.
+typedef std::vector<TopoDS_Face> BRepMesh_FaceList;
+
+//! This tool intended to restore consistent orientation of surfaces.
+//! It takes as input OCCT shape and outputs a set of re-oriented faces.
+class BRepMesh_RestoreOrientationTool : public Standard_Transient
+{
+public:
+
+  DEFINE_STANDARD_ALLOC
+
+  //! Creates uninitialized orientation tool.
+  Standard_EXPORT BRepMesh_RestoreOrientationTool (const bool theVisibilityOnly = false);
+
+  //! Creates orientation tool from the given shape.
+  Standard_EXPORT BRepMesh_RestoreOrientationTool (const TopoDS_Shape& theShape,
+                                                   const bool theVisibilityOnly = false);
+
+public:
+
+  //! Loads the given shape into orientation tool.
+  Standard_EXPORT void Init (const TopoDS_Shape& theShape);
+
+  //! Performs restoring of consistent orientation.
+  Standard_EXPORT void Perform (const Handle(Message_ProgressIndicator) thePI = NULL);
+
+  //! Returns "Visibility only" mode.
+  bool VisibilityOnly() const
+  {
+    return myVisibilityOnly;
+  }
+
+  //! Set "Visibility only" mode.
+  void SetVisibilityOnly (bool isEnabled)
+  {
+    myVisibilityOnly = isEnabled;
+  }
+
+  //! True if the algorithm successfully performed.
+  bool IsDone() const
+  {
+    return myIsDone;
+  }
+
+  //! Returns shape with restored normal orientation.
+  TopoDS_Shape Shape() const
+  {
+    return myShape;
+  }
+
+protected:
+
+  //! Calculates the coherence value between two patches.
+  static Standard_Real computeCoherence (Handle (BRepMesh_TriangulatedPatch)& thePatch0,
+                                         Handle (BRepMesh_TriangulatedPatch)& thePatch1);
+
+  //! Calculates the visibility value for the patch.
+  void computeVisibility (BVH_Triangulation<Standard_Real, 3>& theTriangulation,
+                          Handle (BRepMesh_TriangulatedPatch)& thePatch,
+                          Standard_Integer theRaysNb);
+
+  //! Returns the closest triangle hit by the ray.
+  Standard_Integer traceRay (const BVH_Vec3d& theOrigin,
+                             const BVH_Vec3d& theDirect,
+                             BVH_Triangulation<Standard_Real, 3>& theTriangulation,
+                             Standard_Real& theDistance);
+
+  //! Extracts a triangulation from the face.
+  static Standard_Integer extractTriangulation (const TopoDS_Face& theShape,
+                                                BVH_Triangulation<Standard_Real, 3>& theTriangulation);
+
+  //! Calculates the coherence value between two edges.
+  static Standard_Real edgeCoherence (const BRepMesh_OrientedEdge& theEdge1,
+                                      const BRepMesh_TriangulatedPatch&  thePatch1,
+                                      const BRepMesh_OrientedEdge& theEdge2,
+                                      const BRepMesh_TriangulatedPatch&  thePatch2);
+
+protected:
+
+  //! Is output ready?
+  bool myIsDone;
+
+  //! Resulted shape.
+  TopoDS_Shape myShape;
+
+  //! List of triangulated faces of the shape.
+  BRepMesh_FaceList myFaceList;
+
+  //! List of patches.
+  std::vector<Handle (BRepMesh_TriangulatedPatch)> myPatches;
+
+  //! Vertex and normal data of triangulated shape.
+  BVH_Triangulation<Standard_Real, 3> myTriangulation;
+
+  std::vector<BRepMesh_SuperPatchPair> myPatchQueue;
+
+  //! Maps Id intervals to topological faces.
+  std::vector<Standard_Integer> myFaceIdIntervals;
+
+  //! Use only visibility metric?
+  bool myVisibilityOnly;
+
+  DEFINE_STANDARD_RTTI (BRepMesh_RestoreOrientationTool)
+};
+
+DEFINE_STANDARD_HANDLE(BRepMesh_RestoreOrientationTool, Standard_Transient)
+
+
+#endif

src/BRepMesh/BRepMesh_TriangulatedPatch.cxx

@@ -0,0 +1,233 @@
+// Created on: 2015-08-28
+// Created by: Danila ULYANOV
+// Copyright (c) 2015 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 <BRepMesh_TriangulatedPatch.hxx>
+#include <BVH_LinearBuilder.hxx>
+
+#include <set>
+#include <map>
+
+IMPLEMENT_STANDARD_HANDLE (BRepMesh_TriangulatedPatch, Standard_Transient)
+IMPLEMENT_STANDARD_RTTIEXT(BRepMesh_TriangulatedPatch, Standard_Transient)
+
+namespace {
+
+  //! Compare functor for symmetrical pairs of integers.
+  struct SymmetricalPairComparator
+  {
+    bool operator() (const std::pair<int, int>& theEdge1,
+                      const std::pair<int, int>& theEdge2) const
+    {
+      const int aMin1 = std::min (theEdge1.first, theEdge1.second);
+      const int aMin2 = std::min (theEdge2.first, theEdge2.second);
+
+      return (aMin1 < aMin2) ||
+        (aMin1 == aMin2 && std::max (theEdge1.first, theEdge1.second) < std::max (theEdge2.first, theEdge2.second));
+    }
+
+    bool equal (const std::pair<int, int>& theEdge1,
+                const std::pair<int, int>& theEdge2) const
+    {
+      return std::min (theEdge1.first, theEdge1.second) == std::min (theEdge2.first, theEdge2.second)
+          && std::max (theEdge1.first, theEdge1.second) == std::max (theEdge2.first, theEdge2.second);
+    }
+  };
+
+  // Returns the area of triangle
+  Standard_Real triangleArea (const BVH_Vec3d& theV1,
+                              const BVH_Vec3d& theV2,
+                              const BVH_Vec3d& theV3)
+  {
+    BVH_Vec3d aE1 = theV2 - theV1;
+    BVH_Vec3d aE2 = theV3 - theV1;
+
+    return BVH_Vec3d::Cross (aE1, aE2).Modulus() * 0.5;
+  }
+
+  //! Set of integers.
+  typedef std::set<int> SetOfInteger;
+
+  //! Type of edge-triangles connectivity map.
+  typedef std::map<std::pair<int, int>, SetOfInteger, SymmetricalPairComparator> EdgeFaceMap;
+}
+
+// =======================================================================
+// function : BRepMesh_TriangulatedPatch
+// purpose  :
+// =======================================================================
+BRepMesh_TriangulatedPatch::BRepMesh_TriangulatedPatch (TopoDS_Face theFace)
+  : BVH_Triangulation<Standard_Real, 3>(),
+    myBoundaryEdgeSet (Vertices),
+    myTotalArea (0.0),
+    myFrontVisibility (0.0),
+    myBackVisibility (0.0),
+    myFrontDistance (0.0),
+    myBackDistance (0.0),
+    myFace (theFace)
+{
+  myBuilder = new BVH_LinearBuilder<Standard_Real, 3> (5 /* leaf size */, 32 /* max height */);
+}
+
+// =======================================================================
+// function : Append
+// purpose  :
+// =======================================================================
+void BRepMesh_TriangulatedPatch::Append (const Handle (BRepMesh_TriangulatedPatch)& thePatch)
+{
+  // Vertices
+  Standard_Size aVertOffset = Vertices.size();
+  Vertices.resize (Vertices.size() + thePatch->Vertices.size());
+
+  std::copy (thePatch->Vertices.begin(), thePatch->Vertices.end(), Vertices.begin() + aVertOffset);
+
+  // Triangles
+  struct
+  {
+    Standard_Integer Offset;
+    BVH_Vec4i operator() (const BVH_Vec4i& theTriangle) { return theTriangle + BVH_Vec4i (Offset, Offset, Offset, 0); }
+  } AddOffsetHelper;
+
+  AddOffsetHelper.Offset = static_cast<Standard_Integer> (aVertOffset);
+
+  Standard_Size aFaceOffset = Elements.size();
+  Elements.resize (Elements.size() + thePatch->Elements.size());
+
+  std::transform (thePatch->Elements.begin(), thePatch->Elements.end(), Elements.begin() + aFaceOffset, AddOffsetHelper);
+
+  // Area
+  myTotalArea += thePatch->TotalArea();
+}
+
+// =======================================================================
+// function : computeBoundareEdges
+// purpose  :
+// =======================================================================
+void BRepMesh_TriangulatedPatch::computeBoundaryEdges()
+{
+  myBoundaryEdgeSet.MarkDirty();
+
+  if (!myBoundaryEdgeSet.Elements.empty())
+  {
+    return;
+  }
+
+  EdgeFaceMap anEdgeFaceMap;
+
+  // use triangle indices for connectivity extraction
+  for (Standard_Size aTrgIdx = 0; aTrgIdx < Elements.size(); ++aTrgIdx)
+  {
+    const BVH_Vec4i& aTriangle = Elements[aTrgIdx];
+
+    for (Standard_Integer anEdgeStart = 0, edgeEnd = 1; anEdgeStart < 3; ++anEdgeStart, edgeEnd = (anEdgeStart + 1) % 3)
+    {
+      std::pair<Standard_Integer, Standard_Integer> anEdge (aTriangle[anEdgeStart], aTriangle[edgeEnd]);
+
+      anEdgeFaceMap[anEdge].insert (static_cast<Standard_Integer> (aTrgIdx));
+    }
+  }
+
+  // extract boundary and non-manifold edges
+  myBoundaryEdgeSet.Elements.clear();
+
+  for (EdgeFaceMap::iterator anIter = anEdgeFaceMap.begin(); anIter != anEdgeFaceMap.end(); ++anIter)
+  {
+    const SetOfInteger& anAdjacentFaces = anIter->second;
+
+    const std::pair<Standard_Integer, Standard_Integer>& anEdge = anIter->first;
+
+    if (anAdjacentFaces.size() == 1) // boundary
+    {
+      myBoundaryEdgeSet.Elements.push_back (BVH_Vec4i (anEdge.first, anEdge.second, 0, 1));
+    }
+  }
+
+  myBoundaryEdgeSet.BVH();
+}
+
+// =======================================================================
+// function : PreprocessData
+// purpose  :
+// =======================================================================
+void BRepMesh_TriangulatedPatch::PreprocessData()
+{
+
+  computeBoundaryEdges();
+  updateArea();
+}
+
+// =======================================================================
+// function : Flip
+// purpose  :
+// =======================================================================
+void BRepMesh_TriangulatedPatch::Flip()
+{
+  // mark that triangles was flipped
+  struct
+  {
+    void operator() (BVH_Vec4i& theTriangle) { theTriangle.w() = -theTriangle.w(); }
+  } TriangleFlipHelper;
+
+  std::for_each (Elements.begin(), Elements.end(), TriangleFlipHelper);
+
+  // mark that edges was flipped
+  struct
+  {
+    void operator() (BVH_Vec4i& theEdge) { theEdge.w() = -theEdge.w(); }
+  } EdgeFlipHelper;
+
+  std::for_each (myBoundaryEdgeSet.Elements.begin(), myBoundaryEdgeSet.Elements.end(), EdgeFlipHelper);
+
+  FlipVisibility();
+}
+
+// =======================================================================
+// function : FlipVisibility
+// purpose  :
+// =======================================================================
+void BRepMesh_TriangulatedPatch::FlipVisibility()
+{
+  std::swap (myFrontVisibility, myBackVisibility);
+  std::swap (myFrontDistance, myBackDistance);
+}
+
+// =======================================================================
+// function : updateArea
+// purpose  :
+// =======================================================================
+void BRepMesh_TriangulatedPatch::updateArea()
+{
+  myTotalArea = 0.0;
+
+  myTriangleCDF.resize (Elements.size());
+
+  for (Standard_Size aTrgIdx = 0; aTrgIdx < Elements.size(); ++aTrgIdx)
+  {
+    const BVH_Vec4i& aTriangle = Elements[aTrgIdx];
+
+    const BVH_Vec3d& aV1 = Vertices[aTriangle.x()];
+    const BVH_Vec3d& aV2 = Vertices[aTriangle.y()];
+    const BVH_Vec3d& aV3 = Vertices[aTriangle.z()];
+
+    myTriangleCDF[aTrgIdx] = triangleArea (aV1, aV2, aV3) + (aTrgIdx == 0 ? 0.0 : myTriangleCDF[aTrgIdx - 1]);
+  }
+
+  myTotalArea = myTriangleCDF.back();
+  Standard_Real aInvArea = 1.0 / myTotalArea;
+
+  for (Standard_Size aTrgIdx = 0; aTrgIdx < myTriangleCDF.size(); ++aTrgIdx)
+  {
+    myTriangleCDF[aTrgIdx] *= aInvArea;
+  }
+}

src/BRepMesh/BRepMesh_TriangulatedPatch.hxx

@@ -0,0 +1,201 @@
+// Created on: 2015-08-28
+// Created by: Danila ULYANOV
+// Copyright (c) 2015 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 _BRepMesh_TriangulatedPatch_HeaderFile
+#define _BRepMesh_TriangulatedPatch_HeaderFile
+
+#include <BVH_Triangulation.hxx>
+
+#include <BRepMesh_OrientedEdge.hxx>
+#include <BRepMesh_EdgeSet.hxx>
+
+#include <TopoDS_Face.hxx>
+
+class Handle (BRepMesh_TriangulatedPatch);
+
+//! Class representing a patch composed from triangles.
+class BRepMesh_TriangulatedPatch : public BVH_Triangulation<Standard_Real, 3>, public Standard_Transient
+{
+
+public:
+
+  DEFINE_STANDARD_ALLOC
+
+  typedef BRepMesh_EdgeSet<Standard_Real, 3> EdgeSet;
+
+public:
+
+  //! Creates new triangulated patch.
+  Standard_EXPORT BRepMesh_TriangulatedPatch(TopoDS_Face theFace);
+
+  //! Updates edge and triangle data.
+  Standard_EXPORT void PreprocessData();
+
+  //! Flips the normals orientation of the patch.
+  Standard_EXPORT void Flip();
+
+  //! Flips only visibility coefficients of the patch.
+  Standard_EXPORT void FlipVisibility();
+
+  //! Returns triangle Id stored in 4th component of triangle record.
+  Standard_Integer TriangleId (Standard_Integer theIndex) const
+  {
+    return Elements[theIndex].w();
+  }
+
+  //! Returns set of patch boundary edges.
+  BRepMesh_EdgeSet<Standard_Real, 3>& BoundaryEdgeSet()
+  {
+    return myBoundaryEdgeSet;
+  }
+
+  //! Returns set of patch boundary edges.
+  const BRepMesh_EdgeSet<Standard_Real, 3>& BoundaryEdgeSet() const
+  {
+    return myBoundaryEdgeSet;
+  }
+
+  //! Returns boundary edge by index.
+  BRepMesh_OrientedEdge BoundaryEdge (Standard_Integer theIndex) const
+  {
+    BVH_Vec4i anEdge = myBoundaryEdgeSet.Elements[theIndex];
+    return BRepMesh_OrientedEdge (anEdge.w() > 0.0 ? anEdge.x() : anEdge.y(),
+                                  anEdge.w() > 0.0 ? anEdge.y() : anEdge.x());
+  }
+
+  //! Returns edge vertices.
+  void GetEdgeVertices (const BRepMesh_OrientedEdge& theEdge,
+                        BVH_Vec3d& theVertex1,
+                        BVH_Vec3d& theVertex2) const
+  {
+    theVertex1 = Vertices[theEdge.FirstNode()];
+    theVertex2 = Vertices[theEdge.LastNode()];
+  }
+
+  //! Appends another patch geometry.
+  Standard_EXPORT void Append (const Handle (BRepMesh_TriangulatedPatch)& thePatch);
+
+  //! Returns precomputed CDF for triangles.
+  const std::vector<Standard_Real>& TriangleCDF() const
+  {
+    return myTriangleCDF;
+  }
+
+  //! Returns total area of the patch.
+  Standard_Real TotalArea() const
+  {
+    return myTotalArea;
+  }
+
+  //! Sets total area of the patch.
+  void SetTotalArea (Standard_Real theVal)
+  {
+    myTotalArea = theVal;
+  }
+
+  //! Returns front visibility of the patch.
+  Standard_Real FrontVisibility() const
+  {
+    return myFrontVisibility;
+  }
+
+  //! Sets front visibility of the patch.
+  void SetFrontVisibility (Standard_Real theVal)
+  {
+    myFrontVisibility = theVal;
+  }
+
+  //! Returns back visibility of the patch.
+  Standard_Real BackVisibility() const
+  {
+    return myBackVisibility;
+  }
+
+  //! Sets back visibility of the patch.
+  void SetBackVisibility (Standard_Real theVal)
+  {
+    myBackVisibility = theVal;
+  }
+
+  //! Returns front distance for the patch.
+  Standard_Real FrontDistance() const
+  {
+    return myFrontDistance;
+  }
+
+  //! Sets front distance for the patch.
+  void SetFrontDistance (Standard_Real theVal)
+  {
+    myFrontDistance = theVal;
+  }
+
+  //! Returns back distance for the patch.
+  Standard_Real BackDistance() const
+  {
+    return myBackDistance;
+  }
+
+  //! Sets back distance for the patch.
+  void SetBackDistance (Standard_Real theVal)
+  {
+    myBackDistance = theVal;
+  }
+
+  //! Returns original topological face.
+  TopoDS_Face Face()
+  {
+    return myFace;
+  }
+
+protected:
+
+  //! Extracts edges from triangulation. Fills edge -> faces map. Extracts boundary edges. Fills edge set.
+  void computeBoundaryEdges();
+
+  //! Computes areas of triangles. Computes total patch area. Computes sampling probabilities.
+  void updateArea();
+
+private:
+
+  //! Set of boundary edges.
+  EdgeSet myBoundaryEdgeSet;
+
+  //! Total triangle area.
+  Standard_Real myTotalArea;
+
+  //! Total front visibility.
+  Standard_Real myFrontVisibility;
+
+  //! Total back visibility.
+  Standard_Real myBackVisibility;
+
+  //! Accumulated front distance of the patch.
+  Standard_Real myFrontDistance;
+
+  //! Accumulated back distance of the patch.
+  Standard_Real myBackDistance;
+
+  //! Cumulative distribution function.
+  std::vector<Standard_Real> myTriangleCDF;
+
+  //! Reference to original topological face.
+  TopoDS_Face myFace;
+
+  DEFINE_STANDARD_RTTI (BRepMesh_TriangulatedPatch)
+};
+
+DEFINE_STANDARD_HANDLE(BRepMesh_TriangulatedPatch, Standard_Transient)
+
+#endif

src/BRepMesh/FILES

@@ -54,4 +54,14 @@ BRepMesh_IncrementalMesh.cxx
 BRepMesh_GeomTool.hxx
 BRepMesh_GeomTool.cxx
 BRepMesh_PairOfPolygon.hxx
+BRepMesh_RestoreOrientationTool.cxx
+BRepMesh_RestoreOrientationTool.hxx
+BRepMesh_TriangulatedPatch.cxx
+BRepMesh_TriangulatedPatch.hxx
+BRepMesh_EdgeSet.hxx
+BRepMesh_EdgeSet.lxx
+BRepMesh_Block.hxx
+BRepMesh_MinStCut.cxx
+BRepMesh_MinStCut.hxx
+BRepMesh_LineBoxDistance.hxx
 EXTERNLIB

src/MeshTest/MeshTest.cxx

@@ -14,87 +14,77 @@
 // Alternatively, this file may be used under the terms of Open CASCADE
 // commercial license or contractual agreement.
 
-#include <Standard_Stream.hxx>
-
-#include <stdio.h>
-
-#include <MeshTest.ixx>
 
+#include <AppCont_ContMatrices.hxx>
+#include <Bnd_Box.hxx>
+#include <BRep_Builder.hxx>
+#include <BRep_Tool.hxx>
+#include <BRepAdaptor_Surface.hxx>
+#include <BRepBndLib.hxx>
+#include <BRepBuilderAPI_MakeFace.hxx>
+#include <BRepBuilderAPI_MakePolygon.hxx>
+#include <BRepBuilderAPI_MakeVertex.hxx>
+#include <BRepLib.hxx>
+#include <BRepMesh_DataStructureOfDelaun.hxx>
+#include <BRepMesh_Delaun.hxx>
+#include <BRepMesh_Edge.hxx>
+#include <BRepMesh_FastDiscret.hxx>
+#include <BRepMesh_IncrementalMesh.hxx>
+#include <BRepMesh_Triangle.hxx>
+#include <BRepMesh_Vertex.hxx>
+#include <BRepTest.hxx>
+#include <BRepTools.hxx>
+#include <CSLib.hxx>
+#include <CSLib_DerivativeStatus.hxx>
+#include <DBRep.hxx>
+#include <Draw.hxx>
+#include <Draw_Appli.hxx>
+#include <Draw_Interpretor.hxx>
+#include <Draw_Marker3D.hxx>
+#include <Draw_MarkerShape.hxx>
+#include <Draw_Segment2D.hxx>
+#include <DrawTrSurf.hxx>
+#include <Extrema_LocateExtPC.hxx>
+#include <GCPnts_UniformAbscissa.hxx>
+#include <Geom_Curve.hxx>
+#include <Geom_Plane.hxx>
+#include <Geom_Surface.hxx>
+#include <GeomAdaptor_Curve.hxx>
+#include <GeometryTest.hxx>
+#include <gp_Pln.hxx>
+#include <gp_Trsf.hxx>
+#include <math.hxx>
+#include <math_Matrix.hxx>
+#include <math_Vector.hxx>
+#include <MeshTest.hxx>
 #include <MeshTest_DrawableMesh.hxx>
+#include <PLib.hxx>
+#include <Poly_Connect.hxx>
+#include <Poly_PolygonOnTriangulation.hxx>
+#include <Poly_Triangulation.hxx>
+#include <Precision.hxx>
+#include <Standard_Stream.hxx>
+#include <TColgp_Array1OfPnt2d.hxx>
+#include <TCollection_AsciiString.hxx>
+#include <TColStd_HArray1OfInteger.hxx>
+#include <TColStd_ListIteratorOfListOfInteger.hxx>
+#include <TColStd_MapIteratorOfMapOfInteger.hxx>
 #include <TopAbs_ShapeEnum.hxx>
+#include <TopExp_Explorer.hxx>
+#include <TopLoc_Location.hxx>
 #include <TopoDS.hxx>
+#include <TopoDS_Compound.hxx>
 #include <TopoDS_Edge.hxx>
 #include <TopoDS_Face.hxx>
 #include <TopoDS_Shape.hxx>
-#include <TopoDS_Compound.hxx>
-#include <TopExp_Explorer.hxx>
-#include <TopTools_ListIteratorOfListOfShape.hxx>
-#include <DBRep.hxx>
-#include <BRepTest.hxx>
-#include <GeometryTest.hxx>
-#include <BRep_Tool.hxx>
-#include <BRep_Builder.hxx>
-#include <Draw_MarkerShape.hxx>
-#include <Draw_Appli.hxx>
-#include <Draw.hxx>
-#include <DrawTrSurf.hxx>
-#include <BRepMesh_Triangle.hxx>
-#include <BRepMesh_DataStructureOfDelaun.hxx>
-#include <BRepMesh_Delaun.hxx>
-#include <BRepMesh_FastDiscret.hxx>
-#include <BRepMesh_Vertex.hxx>
-#include <BRepMesh_Edge.hxx>
-#include <BRepMesh_IncrementalMesh.hxx>
-#include <TColStd_ListIteratorOfListOfInteger.hxx>
-#include <TColStd_MapIteratorOfMapOfInteger.hxx>
-#include <Bnd_Box.hxx>
-#include <Precision.hxx>
-#include <Draw_Interpretor.hxx>
-#include <Geom_Plane.hxx>
-#include <Geom_Surface.hxx>
-#include <Draw_Marker3D.hxx>
-#include <Draw_Segment2D.hxx>
-#include <TCollection_AsciiString.hxx>
-
-#include <GCPnts_UniformAbscissa.hxx>
-#include <GeomAdaptor_Curve.hxx>
-#include <Geom_Curve.hxx>
-#include <Extrema_LocateExtPC.hxx>
-
-#include <TopLoc_Location.hxx>
-#include <gp_Trsf.hxx>
-#include <Poly_Triangulation.hxx>
-#include <Poly_Connect.hxx>
-#include <TColgp_Array1OfPnt2d.hxx>
-#include <TColStd_HArray1OfInteger.hxx>
-#include <TopExp_Explorer.hxx>
-#include <gp_Pln.hxx>
-
-#include <PLib.hxx>
-#include <AppCont_ContMatrices.hxx>
-#include <math_Vector.hxx>
-#include <math_Matrix.hxx>
-#include <math.hxx>
-
-#include <CSLib_DerivativeStatus.hxx>
-#include <CSLib.hxx>
-#include <BRepAdaptor_Surface.hxx>
-#include <Bnd_Box.hxx>
-#include <BRepBndLib.hxx>
-#include <BRepLib.hxx>
-
-
-//epa Memory leaks test
-#include <BRepBuilderAPI_MakePolygon.hxx>
 #include <TopoDS_Wire.hxx>
-#include <BRepBuilderAPI_MakeFace.hxx>
-#include <BRepTools.hxx>
-
-//OAN: for triepoints
-#include <BRepBuilderAPI_MakeVertex.hxx>
-#include <Poly_PolygonOnTriangulation.hxx>
+#include <TopTools_ListIteratorOfListOfShape.hxx>
 #include <TopTools_MapIteratorOfMapOfShape.hxx>
+#include <BRepMesh_RestoreOrientationTool.hxx>
 
+#include <stdio.h>
+//epa Memory leaks test
+//OAN: for triepoints
 #ifdef WNT
 Standard_IMPORT Draw_Viewer dout;
 #endif
@@ -1615,6 +1605,45 @@ Standard_Integer correctnormals (Draw_Interpretor& theDI,
   return 0;
 }
 
+//=======================================================================
+//function : RestoreOrientation
+//purpose  : Makes normals of a shape consistently oriented
+//=======================================================================
+Standard_Integer RestoreOrientation (Draw_Interpretor& /*theDI*/,
+                                     Standard_Integer theNArg,
+                                     const char** theArgVal)
+{
+  bool aVisibilityOnly = false;
+
+  if (theNArg < 2)
+  {
+    std::cout << "restoreorient shape [-visibility]\n";
+    std::cout << "Please specify the shape\n";
+  }
+
+  if (theNArg > 2)
+  {
+    if (!strcmp (theArgVal[2], "-visibility"))
+    {
+      aVisibilityOnly = true;
+    }
+  }
+
+  TopoDS_Shape aShape = DBRep::Get (theArgVal[1]);
+
+  Handle(BRepMesh_RestoreOrientationTool) aTool = new BRepMesh_RestoreOrientationTool (aShape, aVisibilityOnly);
+  aTool->Perform();
+
+  if (aTool->IsDone())
+  {
+    std::string aShapeName = std::string (theArgVal[1]) + "_restored";
+    DBRep::Set (aShapeName.c_str(), aTool->Shape());
+    std::cout << "Resulting shape:  " << aShapeName << std::endl;
+  }
+
+  return 0;
+}
+
 //=======================================================================
 void  MeshTest::Commands(Draw_Interpretor& theCommands)
 //=======================================================================
@@ -1651,6 +1680,8 @@ void  MeshTest::Commands(Draw_Interpretor& theCommands)
 
   theCommands.Add("correctnormals", "correctnormals shape",__FILE__, correctnormals, g);
 
+  theCommands.Add("restoreorient", "restoreorient shape [-visibility]",__FILE__, RestoreOrientation, g);
+
 #if 0
   theCommands.Add("extrema","extrema ",__FILE__, extrema, g);
   theCommands.Add("vb","vb ",__FILE__, vb, g);