0027792: TKMath, BVH -- Add tesselator API and implement base marching cube (MC) tessellator

The goal is to extend BVH package with complete set of algorithms allowing generation of distant fields from meshes (or another objects) and tessellation of distant fields. This allows to implement some geometrical algorithm in quite elegant way. E.g., the offset mesh for the given mesh can be calculated in the following way: the input mesh is converted into distant field that represents the mesh in form of implicit function D(x, y, z) = 0. After that, for generation of offset mesh we need to perform tessellation of implicit function D(x, y, z) - T = 0, where T is desired offset distance.
2025-08-09 13:22:24 +03:00 · 2016-08-18 16:15:26 +03:00
14 changed files with 1411 additions and 357 deletions
--- a/src/BVH/BVH.cxx
+++ b/src/BVH/BVH.cxx
@@ -22,6 +22,7 @@
 #include <BVH_BinnedBuilder.hxx>
 #include <BVH_SweepPlaneBuilder.hxx>
 #include <BVH_SpatialMedianBuilder.hxx>
+#include <BVH_MarchingCubes.hxx>

 #include <Standard_Real.hxx>

@@ -151,5 +152,11 @@ template class BVH_DistanceField<Standard_Real, 4>;
 template class BVH_DistanceField<Standard_ShortReal, 3>;
 template class BVH_DistanceField<Standard_ShortReal, 4>;

+template class BVH_Tessellator<Standard_ShortReal, 3>;
+template class BVH_Tessellator<Standard_ShortReal, 4>;
+
+template class BVH_MarchingCubes<Standard_ShortReal, 3>;
+template class BVH_MarchingCubes<Standard_ShortReal, 4>;
+
 template class BVH_Transform<Standard_Real, 4>;
 template class BVH_Transform<Standard_ShortReal, 4>;
--- a/src/BVH/BVH_DistanceField.hxx
+++ b/src/BVH/BVH_DistanceField.hxx
@@ -16,7 +16,11 @@
 #ifndef _BVH_DistanceField_Header
 #define _BVH_DistanceField_Header

+#include <map>
+#include <set>
+
 #include <BVH_Geometry.hxx>
+#include <BVH_ImplicitFunction.hxx>

 template<class T, int N> class BVH_ParallelDistanceFieldBuilder;

@@ -26,7 +30,7 @@ template<class T, int N> class BVH_ParallelDistanceFieldBuilder;
 //! the structure. Distance fields are used as alternative surface representations
 //! (like polygons or NURBS).
 template<class T, int N>
-class BVH_DistanceField
+class BVH_DistanceField : public BVH_ImplicitFunction<T, N>
 {
  friend class BVH_ParallelDistanceFieldBuilder<T, N>;

@@ -43,8 +47,11 @@ public:
  //! Releases resources of 3D distance field.
  virtual ~BVH_DistanceField();

+  //! Calculates distance at the given point.
+  virtual T Calculate (const BVH_VecNt& thePoint);
+
  //! Builds 3D distance field from BVH geometry.
-  Standard_Boolean Build (BVH_Geometry<T, N>& theGeometry);
+  Standard_Boolean Build (BVH_Triangulation<T, N>* theMesh);

 public:

@@ -59,7 +66,11 @@ public:
            const Standard_Integer theY,
            const Standard_Integer theZ)
  {
-    return myVoxelData[theX + (theY + theZ * myDimensionY) * myDimensionX];
+    const Standard_Integer aX = Max (Min (myDimensionX - 1, theX), 0);
+    const Standard_Integer aY = Max (Min (myDimensionY - 1, theY), 0);
+    const Standard_Integer aZ = Max (Min (myDimensionZ - 1, theZ), 0);
+
+    return myVoxelData[aX + (aY + aZ * myDimensionY) * myDimensionX];
  }

  //! Returns distance value for the given voxel.
@@ -67,7 +78,11 @@ public:
           const Standard_Integer theY,
           const Standard_Integer theZ) const
  {
-    return myVoxelData[theX + (theY + theZ * myDimensionY) * myDimensionX];
+    const Standard_Integer aX = Max (Min (myDimensionX - 1, theX), 0);
+    const Standard_Integer aY = Max (Min (myDimensionY - 1, theY), 0);
+    const Standard_Integer aZ = Max (Min (myDimensionZ - 1, theZ), 0);
+
+    return myVoxelData[aX + (aY + aZ * myDimensionY) * myDimensionX];
  }

  //! Returns size of voxel grid in X dimension.
@@ -94,23 +109,59 @@ public:
    return myVoxelSize;
  }

-  //! Returns minimum corner of voxel grid.
-  const BVH_VecNt& CornerMin() const
-  {
-    return myCornerMin;
-  }
+public:

-  //! Returns maximum corner of voxel grid.
-  const BVH_VecNt& CornerMax() const
-  {
-    return myCornerMax;
-  }
+  //! Type of vertex mesh feature.
+  static const Standard_Integer VERT_ID = -2;
+
+  //! Type of face (triangle) mesh feature.
+  static const Standard_Integer FACE_ID = -1;
+
+  //! Feature (vertex, edge, face) of the mesh.
+  typedef std::pair<Standard_Integer,
+                    Standard_Integer> Feature;

 protected:

-  //! Performs building of distance field for the given Z slices.
-  void BuildSlices (BVH_Geometry<T, N>& theGeometry,
-    const Standard_Integer theStartZ, const Standard_Integer theFinalZ);
+  //! Compare functor for symmetrical pairs of integers.
+  struct SymmetricalPairComparator
+  {
+    bool operator() (const Feature& thePair1,
+                     const Feature& thePair2) const
+    {
+
+      const Standard_Integer aMin1 = std::min (thePair1.first, thePair1.second);
+      const Standard_Integer aMax1 = std::max (thePair1.first, thePair1.second);
+      const Standard_Integer aMin2 = std::min (thePair2.first, thePair2.second);
+      const Standard_Integer aMax2 = std::max (thePair2.first, thePair2.second);
+
+      return (aMin1 < aMin2) || (aMin1 == aMin2 && aMax1 < aMax2);
+    }
+
+    bool equal (const Feature& thePair1,
+                const Feature& thePair2) const
+    {
+      return std::min (thePair1.first, thePair1.second) == std::min (thePair2.first, thePair2.second)
+          && std::max (thePair1.first, thePair1.second) == std::max (thePair2.first, thePair2.second);
+    }
+  };
+
+  //! Associates normal with each mesh feature (vertex, edge, face).
+  std::map<Feature, BVH_VecNt, SymmetricalPairComparator> myFeatureNormalMap;
+
+  //! Computes the normal to the given triangle feature (vertex, edge, face).
+  BVH_VecNt normalToFeature (const Feature& theFeature);
+
+protected:
+
+  //! Builds 3D distance field for the given Z slices.
+  void buildSlices (Standard_Integer theStartZ, Standard_Integer theFinalZ);
+
+  //! Computes squared distance from the point to BVH triangulation.
+  T squareDistanceToMesh (const BVH_VecNt& thePnt, Standard_Boolean& theIsOutside);
+
+  //! Returns the closest feature and direction to it from the given point.
+  std::pair<BVH_VecNt, Feature> directToTrg (const BVH_VecNt& thePnt, const Standard_Integer theTrg);

 protected:

@@ -120,12 +171,6 @@ protected:
  //! Size of single voxel.
  BVH_VecNt myVoxelSize;

-  //! Minimum corner of voxel grid.
-  BVH_VecNt myCornerMin;
-
-  //! Maximum corner of voxel grid.
-  BVH_VecNt myCornerMax;
-
  //! Size of voxel grid in X dimension.
  Standard_Integer myDimensionX;

@@ -141,6 +186,9 @@ protected:
  //! Enables/disables signing of distance field.
  Standard_Boolean myComputeSign;

+  //! Triangulation to compute signed distance field.
+  BVH_Triangulation<T, N>* myMesh;
+
 };

 #include <BVH_DistanceField.lxx>
--- a/src/BVH/BVH_DistanceField.lxx
+++ b/src/BVH/BVH_DistanceField.lxx
@@ -52,14 +52,6 @@ BVH_DistanceField<T, N>::~BVH_DistanceField()
  delete [] myVoxelData;
 }

-#if defined (_WIN32) && defined (max)
-  #undef max
-#endif
-
-#include <limits>
-
-#define BVH_DOT3(A, B) (A.x() * B.x() + A.y() * B.y() + A.z() * B.z())
-
 namespace BVH
 {
  //=======================================================================
@@ -92,322 +84,206 @@ namespace BVH
           aNearestY * aNearestY +
           aNearestZ * aNearestZ;
  }
+}

-  //=======================================================================
-  //function : DirectionToNearestPoint
-  //purpose  : Computes squared distance from point to triangle
-  // ======================================================================
-  template<class T, int N>
-  typename VectorType<T, N>::Type DirectionToNearestPoint (
-    const typename VectorType<T, N>::Type& thePoint,
-    const typename VectorType<T, N>::Type& theVertA,
-    const typename VectorType<T, N>::Type& theVertB,
-    const typename VectorType<T, N>::Type& theVertC)
+#define BVH_DOT3(A, B) (A.x() * B.x() + A.y() * B.y() + A.z() * B.z())
+
+//=======================================================================
+//function : directToTrg
+//purpose  : Returns closest feature and direction to it
+// ======================================================================
+template<class T, int N>
+std::pair<typename BVH_DistanceField<T, N>::BVH_VecNt, typename BVH_DistanceField<T, N>::Feature>
+  BVH_DistanceField<T, N>::directToTrg (const BVH_VecNt& thePnt, const Standard_Integer theTrgIdx)
+{
+  Standard_STATIC_ASSERT (N == 3 || N == 4);
+
+  const BVH_Vec4i aTrg = myMesh->Elements[theTrgIdx];
+
+  const BVH_VecNt aVertA = myMesh->Vertices[aTrg.x()];
+  const BVH_VecNt aVertB = myMesh->Vertices[aTrg.y()];
+  const BVH_VecNt aVertC = myMesh->Vertices[aTrg.z()];
+
+  const BVH_VecNt aAB = aVertB - aVertA;
+  const BVH_VecNt aAC = aVertC - aVertA;
+  const BVH_VecNt aAP = thePnt - aVertA;
+
+  const T aABdotAP = BVH_DOT3 (aAB, aAP);
+  const T aACdotAP = BVH_DOT3 (aAC, aAP);
+
+  if (aABdotAP <= static_cast<T> (0) && aACdotAP <= static_cast<T> (0))
  {
-    Standard_STATIC_ASSERT (N == 3 || N == 4);
-
-    const typename VectorType<T, N>::Type aAB = theVertB - theVertA;
-    const typename VectorType<T, N>::Type aAC = theVertC - theVertA;
-    const typename VectorType<T, N>::Type aAP = thePoint - theVertA;
-
-    const T aABdotAP = BVH_DOT3 (aAB, aAP);
-    const T aACdotAP = BVH_DOT3 (aAC, aAP);
-
-    if (aABdotAP <= static_cast<T> (0) && aACdotAP <= static_cast<T> (0))
-    {
-      return aAP;
-    }
-
-    const typename VectorType<T, N>::Type aBC = theVertC - theVertB;
-    const typename VectorType<T, N>::Type aBP = thePoint - theVertB;
-
-    const T aBAdotBP = -BVH_DOT3 (aAB, aBP);
-    const T aBCdotBP =  BVH_DOT3 (aBC, aBP);
-
-    if (aBAdotBP <= static_cast<T> (0) && aBCdotBP <= static_cast<T> (0))
-    {
-      return aBP;
-    }
-
-    const typename VectorType<T, N>::Type aCP = thePoint - theVertC;
-
-    const T aCBdotCP = -BVH_DOT3 (aBC, aCP);
-    const T aCAdotCP = -BVH_DOT3 (aAC, aCP);
-
-    if (aCAdotCP <= static_cast<T> (0) && aCBdotCP <= static_cast<T> (0))
-    {
-      return aCP;
-    }
-
-    const T aACdotBP = BVH_DOT3 (aAC, aBP);
-
-    const T aVC = aABdotAP * aACdotBP + aBAdotBP * aACdotAP;
-
-    if (aVC <= static_cast<T> (0) && aABdotAP >= static_cast<T> (0) && aBAdotBP >= static_cast<T> (0))
-    {
-      return aAP - aAB * (aABdotAP / (aABdotAP + aBAdotBP));
-    }
-
-    const T aABdotCP = BVH_DOT3 (aAB, aCP);
-
-    const T aVA = aBAdotBP * aCAdotCP - aABdotCP * aACdotBP;
-
-    if (aVA <= static_cast<T> (0) && aBCdotBP >= static_cast<T> (0) && aCBdotCP >= static_cast<T> (0))
-    {
-      return aBP - aBC * (aBCdotBP / (aBCdotBP + aCBdotCP));
-    }
-
-    const T aVB = aABdotCP * aACdotAP + aABdotAP * aCAdotCP;
-
-    if (aVB <= static_cast<T> (0) && aACdotAP >= static_cast<T> (0) && aCAdotCP >= static_cast<T> (0))
-    {
-      return aAP - aAC * (aACdotAP / (aACdotAP + aCAdotCP));
-    }
-
-    const T aNorm = static_cast<T> (1.0) / (aVA + aVB + aVC);
-
-    const T aU = aVA * aNorm;
-    const T aV = aVB * aNorm;
-
-    return thePoint - (theVertA * aU + theVertB * aV + theVertC * (static_cast<T> (1.0) - aU - aV));
+    return std::make_pair (aAP, Feature (VERT_ID, aTrg.x()));
  }

-  //=======================================================================
-  //function : SquareDistanceToObject
-  //purpose  : Computes squared distance from point to BVH triangulation
-  //=======================================================================
-  template<class T, int N>
-  T SquareDistanceToObject (BVH_Object<T, N>* theObject,
-    const typename VectorType<T, N>::Type& thePnt, Standard_Boolean& theIsOutside)
+  const BVH_VecNt aBC = aVertC - aVertB;
+  const BVH_VecNt aBP = thePnt - aVertB;
+
+  const T aBAdotBP = -BVH_DOT3 (aAB, aBP);
+  const T aBCdotBP =  BVH_DOT3 (aBC, aBP);
+
+  if (aBAdotBP <= static_cast<T> (0) && aBCdotBP <= static_cast<T> (0))
  {
-    Standard_STATIC_ASSERT (N == 3 || N == 4);
+    return std::make_pair (aBP, Feature (VERT_ID, aTrg.y()));
+  }

-    T aMinDistance = std::numeric_limits<T>::max();
+  const BVH_VecNt aCP = thePnt - aVertC;

-    BVH_Triangulation<T, N>* aTriangulation =
-      dynamic_cast<BVH_Triangulation<T, N>*> (theObject);
+  const T aCBdotCP = -BVH_DOT3 (aBC, aCP);
+  const T aCAdotCP = -BVH_DOT3 (aAC, aCP);

-    Standard_ASSERT_RETURN (aTriangulation != NULL,
-      "Error: Unsupported BVH object (non triangulation)", aMinDistance);
+  if (aCAdotCP <= static_cast<T> (0) && aCBdotCP <= static_cast<T> (0))
+  {
+    return std::make_pair (aCP, Feature (VERT_ID, aTrg.z()));
+  }

-    const NCollection_Handle<BVH_Tree<T, N> >& aBVH = aTriangulation->BVH();
+  const T aACdotBP = BVH_DOT3 (aAC, aBP);

-    if (aBVH.IsNull())
+  const T aVC = aABdotAP * aACdotBP + aBAdotBP * aACdotAP;
+
+  if (aVC <= static_cast<T> (0) && aABdotAP >= static_cast<T> (0) && aBAdotBP >= static_cast<T> (0))
+  {
+    return std::make_pair (aAP - aAB * (aABdotAP / (aABdotAP + aBAdotBP)), Feature (aTrg.x(), aTrg.y()));
+  }
+
+  const T aABdotCP = BVH_DOT3 (aAB, aCP);
+
+  const T aVA = aBAdotBP * aCAdotCP - aABdotCP * aACdotBP;
+
+  if (aVA <= static_cast<T> (0) && aBCdotBP >= static_cast<T> (0) && aCBdotCP >= static_cast<T> (0))
+  {
+    return std::make_pair (aBP - aBC * (aBCdotBP / (aBCdotBP + aCBdotCP)), Feature (aTrg.y(), aTrg.z()));
+  }
+
+  const T aVB = aABdotCP * aACdotAP + aABdotAP * aCAdotCP;
+
+  if (aVB <= static_cast<T> (0) && aACdotAP >= static_cast<T> (0) && aCAdotCP >= static_cast<T> (0))
+  {
+    return std::make_pair (aAP - aAC * (aACdotAP / (aACdotAP + aCAdotCP)), Feature (aTrg.x(), aTrg.z()));
+  }
+
+  const T aNorm = static_cast<T> (1.0) / (aVA + aVB + aVC);
+
+  const T aU = aVA * aNorm;
+  const T aV = aVB * aNorm;
+
+  const BVH_VecNt aDirect = thePnt - (aVertA * aU + aVertB * aV + aVertC * (static_cast<T> (1.0) - aU - aV));
+
+  return std::make_pair (aDirect, Feature (FACE_ID, theTrgIdx));
+}
+
+//=======================================================================
+//function : squareDistanceToMesh
+//purpose  : Computes squared distance from point to BVH triangulation
+//=======================================================================
+template<class T, int N>
+T BVH_DistanceField<T, N>::squareDistanceToMesh (const BVH_VecNt& thePnt, Standard_Boolean& theIsOut)
+{
+  Standard_STATIC_ASSERT (N == 3 || N == 4);
+
+  T aMinDist = std::numeric_limits<T>::max();
+
+  const NCollection_Handle<BVH_Tree<T, N> >& aBVH = myMesh->BVH();
+
+  if (aBVH.IsNull())
+  {
+    return Standard_False;
+  }
+
+  std::pair<Standard_Integer, T> aStack[32];
+
+  Standard_Integer aHead = -1;
+  Standard_Integer aNode =  0; // root node
+
+  for (;;)
+  {
+    BVH_Vec4i aData = aBVH->NodeInfoBuffer()[aNode];
+
+    if (aData.x() == 0) // if inner node
    {
-      return Standard_False;
-    }
+      const T aDistToLft = BVH::DistanceToBox<T, N> (thePnt,
+                                                     aBVH->MinPoint (aData.y()),
+                                                     aBVH->MaxPoint (aData.y()));

-    std::pair<Standard_Integer, T> aStack[32];
+      const T aDistToRgh = BVH::DistanceToBox<T, N> (thePnt,
+                                                     aBVH->MinPoint (aData.z()),
+                                                     aBVH->MaxPoint (aData.z()));

-    Standard_Integer aHead = -1;
-    Standard_Integer aNode =  0; // root node
+      const Standard_Boolean aHitLft = aDistToLft <= aMinDist;
+      const Standard_Boolean aHitRgh = aDistToRgh <= aMinDist;

-    for (;;)
-    {
-      BVH_Vec4i aData = aBVH->NodeInfoBuffer()[aNode];
-
-      if (aData.x() == 0) // if inner node
+      if (aHitLft & aHitRgh)
      {
-        const T aDistToLft = DistanceToBox<T, N> (thePnt,
-                                                  aBVH->MinPoint (aData.y()),
-                                                  aBVH->MaxPoint (aData.y()));
+        aNode = (aDistToLft < aDistToRgh) ? aData.y() : aData.z();

-        const T aDistToRgh = DistanceToBox<T, N> (thePnt,
-                                                  aBVH->MinPoint (aData.z()),
-                                                  aBVH->MaxPoint (aData.z()));
-
-        const Standard_Boolean aHitLft = aDistToLft <= aMinDistance;
-        const Standard_Boolean aHitRgh = aDistToRgh <= aMinDistance;
-
-        if (aHitLft & aHitRgh)
+        aStack[++aHead] = std::pair<Standard_Integer, T> (
+          aDistToLft < aDistToRgh ? aData.z() : aData.y(), Max (aDistToLft, aDistToRgh));
+      }
+      else
+      {
+        if (aHitLft | aHitRgh)
        {
-          aNode = (aDistToLft < aDistToRgh) ? aData.y() : aData.z();
-
-          aStack[++aHead] = std::pair<Standard_Integer, T> (
-            aDistToLft < aDistToRgh ? aData.z() : aData.y(), Max (aDistToLft, aDistToRgh));
+          aNode = aHitLft ? aData.y() : aData.z();
        }
        else
        {
-          if (aHitLft | aHitRgh)
-          {
-            aNode = aHitLft ? aData.y() : aData.z();
-          }
-          else
+          if (aHead < 0)
+            return aMinDist;
+
+          std::pair<Standard_Integer, T>& anInfo = aStack[aHead--];
+
+          while (anInfo.second > aMinDist)
          {
            if (aHead < 0)
-              return aMinDistance;
+              return aMinDist;

-            std::pair<Standard_Integer, T>& anInfo = aStack[aHead--];
-
-            while (anInfo.second > aMinDistance)
-            {
-              if (aHead < 0)
-                return aMinDistance;
-
-              anInfo = aStack[aHead--];
-            }
-
-            aNode = anInfo.first;
+            anInfo = aStack[aHead--];
          }
+
+          aNode = anInfo.first;
        }
      }
-      else  // if leaf node
-      {
-        for (Standard_Integer aTrgIdx = aData.y(); aTrgIdx <= aData.z(); ++aTrgIdx)
-        {
-          const BVH_Vec4i aTriangle = aTriangulation->Elements[aTrgIdx];
-
-          const typename VectorType<T, N>::Type aVertex0 = aTriangulation->Vertices[aTriangle.x()];
-          const typename VectorType<T, N>::Type aVertex1 = aTriangulation->Vertices[aTriangle.y()];
-          const typename VectorType<T, N>::Type aVertex2 = aTriangulation->Vertices[aTriangle.z()];
-
-          const typename VectorType<T, N>::Type aDirection =
-            DirectionToNearestPoint<T, N> (thePnt, aVertex0, aVertex1, aVertex2);
-
-          const T aDistance = BVH_DOT3 (aDirection, aDirection);
-
-          if (aDistance < aMinDistance)
-          {
-            aMinDistance = aDistance;
-
-            typename VectorType<T, N>::Type aTrgEdges[] = { aVertex1 - aVertex0,
-                                                            aVertex2 - aVertex0 };
-
-            typename VectorType<T, N>::Type aTrgNormal;
-            
-            aTrgNormal.x() = aTrgEdges[0].y() * aTrgEdges[1].z() - aTrgEdges[0].z() * aTrgEdges[1].y();
-            aTrgNormal.y() = aTrgEdges[0].z() * aTrgEdges[1].x() - aTrgEdges[0].x() * aTrgEdges[1].z();
-            aTrgNormal.z() = aTrgEdges[0].x() * aTrgEdges[1].y() - aTrgEdges[0].y() * aTrgEdges[1].x();
-
-            theIsOutside = BVH_DOT3 (aTrgNormal, aDirection) > 0;
-          }
-        }
-
-        if (aHead < 0)
-          return aMinDistance;
-
-        std::pair<Standard_Integer, T>& anInfo = aStack[aHead--];
-
-        while (anInfo.second > aMinDistance)
-        {
-          if (aHead < 0)
-            return aMinDistance;
-
-          anInfo = aStack[aHead--];
-        }
-
-        aNode = anInfo.first;
-      }
    }
-  }
-
-  //=======================================================================
-  //function : SquareDistanceToGeomerty
-  //purpose  : Computes squared distance from point to BVH geometry
-  //=======================================================================
-  template<class T, int N>
-  T SquareDistanceToGeomerty (BVH_Geometry<T, N>& theGeometry,
-    const typename VectorType<T, N>::Type& thePnt, Standard_Boolean& theIsOutside)
-  {
-    Standard_STATIC_ASSERT (N == 3 || N == 4);
-
-    const BVH_Tree<T, N, BVH_BinaryTree>* aBVH = theGeometry.BVH().get();
-
-    if (aBVH == NULL)
+    else  // if leaf node
    {
-      return Standard_False;
-    }
-
-    std::pair<Standard_Integer, T> aStack[32];
-
-    Standard_Integer aHead = -1;
-    Standard_Integer aNode =  0; // root node
-
-    T aMinDistance = std::numeric_limits<T>::max();
-
-    for (;;)
-    {
-      BVH_Vec4i aData = aBVH->NodeInfoBuffer()[aNode];
-
-      if (aData.x() == 0) // if inner node
+      for (Standard_Integer aTrgIdx = aData.y(); aTrgIdx <= aData.z(); ++aTrgIdx)
      {
-        const T aDistToLft = DistanceToBox<T, N> (thePnt,
-                                                  aBVH->MinPoint (aData.y()),
-                                                  aBVH->MaxPoint (aData.y()));
+        std::pair<BVH_VecNt, Feature> aDirToFeature = directToTrg (thePnt, aTrgIdx);

-        const T aDistToRgh = DistanceToBox<T, N> (thePnt,
-                                                  aBVH->MinPoint (aData.z()),
-                                                  aBVH->MaxPoint (aData.z()));
+        const T aDistance = BVH_DOT3 (aDirToFeature.first,
+                                      aDirToFeature.first);

-        const Standard_Boolean aHitLft = aDistToLft <= aMinDistance;
-        const Standard_Boolean aHitRgh = aDistToRgh <= aMinDistance;
-
-        if (aHitLft & aHitRgh)
+        if (aDistance < aMinDist)
        {
-          aNode = (aDistToLft < aDistToRgh) ? aData.y() : aData.z();
+          aMinDist = aDistance;

-          aStack[++aHead] = std::pair<Standard_Integer, T> (
-            aDistToLft < aDistToRgh ? aData.z() : aData.y(), Max (aDistToLft, aDistToRgh));
-        }
-        else
-        {
-          if (aHitLft | aHitRgh)
+          if (myComputeSign)
          {
-            aNode = aHitLft ? aData.y() : aData.z();
-          }
-          else
-          {
-            if (aHead < 0)
-              return aMinDistance;
+            const BVH_VecNt aNormal = normalToFeature (aDirToFeature.second);

-            std::pair<Standard_Integer, T>& anInfo = aStack[aHead--];
-
-            while (anInfo.second > aMinDistance)
-            {
-              if (aHead < 0)
-                return aMinDistance;
-
-              anInfo = aStack[aHead--];
-            }
-
-            aNode = anInfo.first;
+            theIsOut = BVH_DOT3 (aNormal, aDirToFeature.first) > T(0.0);
          }
        }
      }
-      else  // if leaf node
+
+      if (aHead < 0)
+        return aMinDist;
+
+      std::pair<Standard_Integer, T>& anInfo = aStack[aHead--];
+
+      while (anInfo.second > aMinDist)
      {
-        Standard_Boolean isOutside = Standard_True;
-
-        const T aDistance = SquareDistanceToObject (
-          theGeometry.Objects()(aNode).operator->(), thePnt, isOutside);
-
-        if (aDistance < aMinDistance)
-        {
-          aMinDistance = aDistance;
-          theIsOutside = isOutside;
-        }
-
        if (aHead < 0)
-          return aMinDistance;
+          return aMinDist;

-        std::pair<Standard_Integer, T>& anInfo = aStack[aHead--];
-
-        while (anInfo.second > aMinDistance)
-        {
-          if (aHead < 0)
-            return aMinDistance;
-
-          anInfo = aStack[aHead--];
-        }
-
-        aNode = anInfo.first;
+        anInfo = aStack[aHead--];
      }
+
+      aNode = anInfo.first;
    }
  }
 }

-#undef BVH_DOT3
-
 #ifdef HAVE_TBB

 //! Tool object for parallel construction of distance field (uses Intel TBB).
@@ -416,36 +292,163 @@ class BVH_ParallelDistanceFieldBuilder
 {
 private:

-  //! Input BVH geometry.
-  BVH_Geometry<T, N>* myGeometry;
-
  //! Output distance field.
  BVH_DistanceField<T, N>* myOutField;

 public:

-  BVH_ParallelDistanceFieldBuilder (BVH_DistanceField<T, N>* theOutField, BVH_Geometry<T, N>* theGeometry)
-  : myGeometry (theGeometry),
-    myOutField (theOutField)
+  BVH_ParallelDistanceFieldBuilder (BVH_DistanceField<T, N>* theOutField)
+  : myOutField (theOutField)
  {
    //
  }

  void operator() (const tbb::blocked_range<Standard_Integer>& theRange) const
  {
-    myOutField->BuildSlices (*myGeometry, theRange.begin(), theRange.end());
+    myOutField->buildSlices (theRange.begin(), theRange.end());
  }
 };

 #endif

+//=======================================================================
+//function : normalToFeature
+//purpose  : Computes the normal to the given triangle feature
+//=======================================================================
+template<class T, int N>
+typename BVH_DistanceField<T, N>::BVH_VecNt BVH_DistanceField<T, N>::normalToFeature (const Feature& theFeature)
+{
+  std::map<Feature, BVH_VecNt, SymmetricalPairComparator>::iterator aFeatureIt = myFeatureNormalMap.find (theFeature);
+
+  Standard_ASSERT_RAISE (aFeatureIt != myFeatureNormalMap.end(), "Failed to retrieve feature normal");
+
+  return aFeatureIt->second;
+}
+
 // =======================================================================
-// function : BuildSlices
+// function : Build
+// purpose  : Builds 3D distance field from BVH geometry
+// =======================================================================
+template<class T, int N>
+Standard_Boolean BVH_DistanceField<T, N>::Build (BVH_Triangulation<T, N>* theMesh)
+{
+  myMesh = theMesh;
+
+  if (myMesh->Size() == 0)
+  {
+    return Standard_False;
+  }
+
+  //----------------------------------------------------------------------
+  // Compute pseudo-normals
+  //----------------------------------------------------------------------
+
+  myFeatureNormalMap.clear();
+
+  if (myComputeSign)
+  {
+    std::map<Feature, std::set<Standard_Integer>, SymmetricalPairComparator> aConnectivityMap;
+
+    for (size_t aTrgIdx = 0; aTrgIdx < myMesh->Elements.size(); ++aTrgIdx) // build connectivity map
+    {
+      const BVH_Vec4i aTrg = myMesh->Elements[aTrgIdx];
+
+      const BVH_VecNt aTrgVrt0 = myMesh->Vertices[aTrg.x()];
+      const BVH_VecNt aTrgVrt1 = myMesh->Vertices[aTrg.y()];
+      const BVH_VecNt aTrgVrt2 = myMesh->Vertices[aTrg.z()];
+
+      const BVH_VecNt aTrgEdges[] = { aTrgVrt1 - aTrgVrt0,
+                                      aTrgVrt2 - aTrgVrt0 };
+
+      BVH_VecNt aNormal;
+
+      aNormal.x() = aTrgEdges[0].y() * aTrgEdges[1].z() - aTrgEdges[0].z() * aTrgEdges[1].y();
+      aNormal.y() = aTrgEdges[0].z() * aTrgEdges[1].x() - aTrgEdges[0].x() * aTrgEdges[1].z();
+      aNormal.z() = aTrgEdges[0].x() * aTrgEdges[1].y() - aTrgEdges[0].y() * aTrgEdges[1].x();
+
+      myFeatureNormalMap[std::make_pair (FACE_ID, aTrgIdx)] = aNormal / std::sqrt (BVH_DOT3 (aNormal, aNormal));
+
+      for (Standard_Integer aVrtStart = 0, aVrtEnd = 1; aVrtStart < 3; ++aVrtStart, aVrtEnd = (aVrtEnd + 1) % 3)
+      {
+        aConnectivityMap[Feature (aTrg[aVrtStart], aTrg[aVrtEnd])].insert (aTrgIdx);
+      }
+
+      for (Standard_Integer aVrtIdx = 0; aVrtIdx < 3; ++aVrtIdx)
+      {
+        aConnectivityMap[Feature (VERT_ID, aTrg[aVrtIdx])].insert (aTrgIdx);
+      }
+    }
+
+    std::map<Feature, std::set<Standard_Integer>, SymmetricalPairComparator>::iterator aFeatureIter;
+
+    for (aFeatureIter = aConnectivityMap.begin(); aFeatureIter != aConnectivityMap.end(); ++aFeatureIter)
+    {
+      BVH_VecNt aNormal;
+
+      for (std::set<Standard_Integer>::iterator aTrgIt = aFeatureIter->second.begin(); aTrgIt != aFeatureIter->second.end(); ++aTrgIt)
+      {
+        aNormal += myFeatureNormalMap[std::make_pair (FACE_ID, *aTrgIt)];
+      }
+
+      myFeatureNormalMap[aFeatureIter->first] = aNormal;
+    }
+  }
+
+  //----------------------------------------------------------------------
+  // Adjust parameters of voxel set
+  //----------------------------------------------------------------------
+
+  const BVH_VecNt aBoxSize = myMesh->Box().Size();
+
+  const T aMaxBoxSide = Max (Max (aBoxSize.x(), aBoxSize.y()), aBoxSize.z());
+
+  myDimensionX = static_cast<Standard_Integer> (myMaximumSize * aBoxSize.x() / aMaxBoxSide);
+  myDimensionY = static_cast<Standard_Integer> (myMaximumSize * aBoxSize.y() / aMaxBoxSide);
+  myDimensionZ = static_cast<Standard_Integer> (myMaximumSize * aBoxSize.z() / aMaxBoxSide);
+
+  myDimensionX = Min (myMaximumSize, Max (myDimensionX, 16));
+  myDimensionY = Min (myMaximumSize, Max (myDimensionY, 16));
+  myDimensionZ = Min (myMaximumSize, Max (myDimensionZ, 16));
+
+  const BVH_VecNt aGlobalBoxMin = myMesh->Box().CornerMin();
+  const BVH_VecNt aGlobalBoxMax = myMesh->Box().CornerMax();
+
+  const Standard_Integer aVoxelOffset = 2;
+
+  myCornerMin.x() = aGlobalBoxMin.x() - aVoxelOffset * aBoxSize.x() / (myDimensionX - 2 * aVoxelOffset);
+  myCornerMin.y() = aGlobalBoxMin.y() - aVoxelOffset * aBoxSize.y() / (myDimensionY - 2 * aVoxelOffset);
+  myCornerMin.z() = aGlobalBoxMin.z() - aVoxelOffset * aBoxSize.z() / (myDimensionZ - 2 * aVoxelOffset);
+
+  myCornerMax.x() = aGlobalBoxMax.x() + aVoxelOffset * aBoxSize.x() / (myDimensionX - 2 * aVoxelOffset);
+  myCornerMax.y() = aGlobalBoxMax.y() + aVoxelOffset * aBoxSize.y() / (myDimensionY - 2 * aVoxelOffset);
+  myCornerMax.z() = aGlobalBoxMax.z() + aVoxelOffset * aBoxSize.z() / (myDimensionZ - 2 * aVoxelOffset);
+
+  myVoxelSize.x() = (myCornerMax.x() - myCornerMin.x()) / myDimensionX;
+  myVoxelSize.y() = (myCornerMax.y() - myCornerMin.y()) / myDimensionY;
+  myVoxelSize.z() = (myCornerMax.z() - myCornerMin.z()) / myDimensionZ;
+
+  //----------------------------------------------------------------------
+  // Perform calculations
+  //----------------------------------------------------------------------
+
+#ifdef HAVE_TBB
+  tbb::parallel_for (
+    tbb::blocked_range<Standard_Integer> (0, myDimensionZ), BVH_ParallelDistanceFieldBuilder<T, N> (this));
+#else
+  buildSlices (0, myDimensionZ);
+#endif
+
+  return Standard_True;
+}
+
+#undef BVH_DOT3
+
+// =======================================================================
+// function : buildSlices
 // purpose  : Performs building of distance field for the given Z slices
 // =======================================================================
 template<class T, int N>
-void BVH_DistanceField<T, N>::BuildSlices (BVH_Geometry<T, N>& theGeometry,
-  const Standard_Integer theStartSlice, const Standard_Integer theFinalSlice)
+void BVH_DistanceField<T, N>::buildSlices (Standard_Integer theStartSlice, Standard_Integer theFinalSlice)
 {
  for (Standard_Integer aZ = theStartSlice; aZ < theFinalSlice; ++aZ)
  {
@@ -461,66 +464,48 @@ void BVH_DistanceField<T, N>::BuildSlices (BVH_Geometry<T, N>& theGeometry,

        Standard_Boolean isOutside = Standard_True;

-        const T aDistance = sqrt (
-          BVH::SquareDistanceToGeomerty<T, N> (theGeometry, aCenter, isOutside));
+        const T aDistance = std::sqrt (squareDistanceToMesh (aCenter, isOutside));

        Voxel (aX, aY, aZ) = (!myComputeSign || isOutside) ? aDistance : -aDistance;
      }
    }
+
+    std::cout << "*";
  }
 }

 // =======================================================================
-// function : Build
-// purpose  : Builds 3D distance field from BVH geometry
+// function : Calculate
+// purpose  : Calculates distance at the given point
 // =======================================================================
 template<class T, int N>
-Standard_Boolean BVH_DistanceField<T, N>::Build (BVH_Geometry<T, N>& theGeometry)
+T BVH_DistanceField<T, N>::Calculate (const BVH_VecNt& thePoint)
 {
-  if (theGeometry.Size() == 0)
-  {
-    return Standard_False;
-  }
+  const BVH_VecNt aLocalPnt = thePoint.cwiseMin (myCornerMax).cwiseMax (myCornerMin) - myCornerMin;

-  const BVH_VecNt aGlobalBoxSize = theGeometry.Box().Size();
+  const T aGridPntX = aLocalPnt.x() / myVoxelSize.x() - static_cast<T> (0.5);
+  const T aGridPntY = aLocalPnt.y() / myVoxelSize.y() - static_cast<T> (0.5);
+  const T aGridPntZ = aLocalPnt.z() / myVoxelSize.z() - static_cast<T> (0.5);

-  const T aMaxBoxSide = Max (Max (aGlobalBoxSize.x(), aGlobalBoxSize.y()), aGlobalBoxSize.z());
+  const Standard_Integer aCellX = static_cast<Standard_Integer> (std::floor (aGridPntX));
+  const Standard_Integer aCellY = static_cast<Standard_Integer> (std::floor (aGridPntY));
+  const Standard_Integer aCellZ = static_cast<Standard_Integer> (std::floor (aGridPntZ));

-  myDimensionX = static_cast<Standard_Integer> (myMaximumSize * aGlobalBoxSize.x() / aMaxBoxSide);
-  myDimensionY = static_cast<Standard_Integer> (myMaximumSize * aGlobalBoxSize.y() / aMaxBoxSide);
-  myDimensionZ = static_cast<Standard_Integer> (myMaximumSize * aGlobalBoxSize.z() / aMaxBoxSide);
+  const T aRatioX = aGridPntX - aCellX;
+  const T aRatioY = aGridPntY - aCellY;
+  const T aRatioZ = aGridPntZ - aCellZ;

-  myDimensionX = Min (myMaximumSize, Max (myDimensionX, 16));
-  myDimensionY = Min (myMaximumSize, Max (myDimensionY, 16));
-  myDimensionZ = Min (myMaximumSize, Max (myDimensionZ, 16));
+  const T aInterpValueMinYMinZ =
+    Voxel (aCellX + 0, aCellY + 0, aCellZ + 0) * (T(1.0) - aRatioX) + Voxel (aCellX + 1, aCellY + 0, aCellZ + 0) * aRatioX;
+  const T aInterpValueMaxYMinZ =
+    Voxel (aCellX + 0, aCellY + 1, aCellZ + 0) * (T(1.0) - aRatioX) + Voxel (aCellX + 1, aCellY + 1, aCellZ + 0) * aRatioX;
+  const T aInterpValueMinYMaxZ =
+    Voxel (aCellX + 0, aCellY + 0, aCellZ + 1) * (T(1.0) - aRatioX) + Voxel (aCellX + 1, aCellY + 0, aCellZ + 1) * aRatioX;
+  const T aInterpValueMaxYMaxZ =
+    Voxel (aCellX + 0, aCellY + 1, aCellZ + 1) * (T(1.0) - aRatioX) + Voxel (aCellX + 1, aCellY + 1, aCellZ + 1) * aRatioX;

-  const BVH_VecNt aGlobalBoxMin = theGeometry.Box().CornerMin();
-  const BVH_VecNt aGlobalBoxMax = theGeometry.Box().CornerMax();
+  const T aInterpValueMinZ = aInterpValueMinYMinZ * (T(1.0) - aRatioY) + aInterpValueMaxYMinZ * aRatioY;
+  const T aInterpValueMaxZ = aInterpValueMinYMaxZ * (T(1.0) - aRatioY) + aInterpValueMaxYMaxZ * aRatioY;

-  const Standard_Integer aVoxelOffset = 2;
-
-  myCornerMin.x() = aGlobalBoxMin.x() - aVoxelOffset * aGlobalBoxSize.x() / (myDimensionX - 2 * aVoxelOffset);
-  myCornerMin.y() = aGlobalBoxMin.y() - aVoxelOffset * aGlobalBoxSize.y() / (myDimensionY - 2 * aVoxelOffset);
-  myCornerMin.z() = aGlobalBoxMin.z() - aVoxelOffset * aGlobalBoxSize.z() / (myDimensionZ - 2 * aVoxelOffset);
-
-  myCornerMax.x() = aGlobalBoxMax.x() + aVoxelOffset * aGlobalBoxSize.x() / (myDimensionX - 2 * aVoxelOffset);
-  myCornerMax.y() = aGlobalBoxMax.y() + aVoxelOffset * aGlobalBoxSize.y() / (myDimensionY - 2 * aVoxelOffset);
-  myCornerMax.z() = aGlobalBoxMax.z() + aVoxelOffset * aGlobalBoxSize.z() / (myDimensionZ - 2 * aVoxelOffset);
-  
-  myVoxelSize.x() = (myCornerMax.x() - myCornerMin.x()) / myDimensionX;
-  myVoxelSize.y() = (myCornerMax.y() - myCornerMin.y()) / myDimensionY;
-  myVoxelSize.z() = (myCornerMax.z() - myCornerMin.z()) / myDimensionZ;
-
-#ifdef HAVE_TBB
-  
-  tbb::parallel_for (tbb::blocked_range<Standard_Integer> (0, myDimensionZ),
-    BVH_ParallelDistanceFieldBuilder<T, N> (this, &theGeometry));
-
-#else
-
-  BuildSlices (theGeometry, 0, myDimensionZ);
-
-#endif
-
-  return Standard_True;
+  return aInterpValueMinZ * (T(1.0) - aRatioZ) + aInterpValueMaxZ * aRatioZ;
 }
--- a/src/BVH/BVH_ImplicitFunction.hxx
+++ b/src/BVH/BVH_ImplicitFunction.hxx
@@ -0,0 +1,83 @@
+// Created on: 2015-04-07
+// Created by: Denis BOGOLEPOV
+// Copyright (c) 2013-2014 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 _BVH_ImplicitFunction_Header
+#define _BVH_ImplicitFunction_Header
+
+#include <BVH_Types.hxx>
+
+//! Represents implicit function defined over 2D/3D space.
+template<class T, int N>
+class BVH_ImplicitFunction
+{
+public:
+
+  typedef typename BVH::VectorType<T, N>::Type BVH_VecNt;
+
+public:
+
+  //! Creates new implicit function.
+  BVH_ImplicitFunction() : myCornerMin (static_cast<T> (1.0)),
+                           myCornerMax (static_cast<T> (1.0))
+  {
+    for (int aDim = 0; aDim < std::min (N, 3); ++aDim)
+    {
+      myCornerMin[aDim] = -std::numeric_limits<T>::max();
+      myCornerMax[aDim] =  std::numeric_limits<T>::max();
+    }
+  }
+
+  //! Creates new implicit function with the given domain.
+  BVH_ImplicitFunction (const BVH_VecNt& theCornerMin,
+                        const BVH_VecNt& theCornerMax) : myCornerMin (theCornerMin),
+                                                         myCornerMax (theCornerMax)
+  {
+    Standard_ASSERT_INVOKE (N == 2 || N == 3 || N == 4);
+  }
+
+  //! Releases resources of implicit function.
+  virtual ~BVH_ImplicitFunction() { }
+
+  //! Returns minimum corner of function domain.
+  const BVH_VecNt& CornerMin() const
+  {
+    return myCornerMin;
+  }
+
+  //! Returns maximum corner of function domain.
+  const BVH_VecNt& CornerMax() const
+  {
+    return myCornerMax;
+  }
+
+  //! Calculates function value at the given point.
+  virtual T Calculate (const BVH_VecNt& thePoint) = 0;
+
+  //! Estimates surface normal at the given point.
+  BVH_VecNt EstimateNormal (const BVH_VecNt& thePoint);
+
+protected:
+
+  //! Minimum corner of function domain.
+  BVH_VecNt myCornerMin;
+
+  //! Maximum corner of function domain.
+  BVH_VecNt myCornerMax;
+
+};
+
+#include <BVH_ImplicitFunction.lxx>
+
+#endif // _BVH_ImplicitFunction_Header
--- a/src/BVH/BVH_ImplicitFunction.lxx
+++ b/src/BVH/BVH_ImplicitFunction.lxx
@@ -0,0 +1,37 @@
+// Created on: 2016-08-18
+// Created by: Denis BOGOLEPOV
+// Copyright (c) 2013-2016 OPEN CASCADE SAS
+//
+// This file is part of Open CASCADE Technology software library.
+//
+// This library is free software; you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License version 2.1 as published
+// by the Free Software Foundation, with special exception defined in the file
+// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
+// distribution for complete text of the license and disclaimer of any warranty.
+//
+// Alternatively, this file may be used under the terms of Open CASCADE
+// commercial license or contractual agreement.
+
+// =======================================================================
+// function : EstimateNormal
+// purpose  : Estimates surface normal at the given point
+// =======================================================================
+template<class T, int N>
+typename BVH_ImplicitFunction<T, N>::BVH_VecNt BVH_ImplicitFunction<T, N>::EstimateNormal (const BVH_VecNt& thePoint)
+{
+  BVH_VecNt aNormal;
+
+  for (int aDim = 0; aDim < std::min (N, 3); ++aDim)
+  {
+    BVH_VecNt aLftOffset (thePoint);
+    BVH_VecNt aRghOffset (thePoint);
+
+    aLftOffset[aDim] -= BVH::Precision<T>::Confusion();
+    aRghOffset[aDim] += BVH::Precision<T>::Confusion();
+
+    aNormal[aDim] = Calculate (aRghOffset) - Calculate (aLftOffset);
+  }
+
+  return aNormal;
+}
--- a/src/BVH/BVH_MarchingCubes.cxx
+++ b/src/BVH/BVH_MarchingCubes.cxx
@@ -0,0 +1,316 @@
+#include <BVH_MarchingCubes.hxx>
+
+namespace BVH
+{
+  int MarchingCubeTables::myEdgeTable[256] = {
+    0x0  , 0x109, 0x203, 0x30a, 0x406, 0x50f, 0x605, 0x70c,
+    0x80c, 0x905, 0xa0f, 0xb06, 0xc0a, 0xd03, 0xe09, 0xf00,
+    0x190, 0x99 , 0x393, 0x29a, 0x596, 0x49f, 0x795, 0x69c,
+    0x99c, 0x895, 0xb9f, 0xa96, 0xd9a, 0xc93, 0xf99, 0xe90,
+    0x230, 0x339, 0x33 , 0x13a, 0x636, 0x73f, 0x435, 0x53c,
+    0xa3c, 0xb35, 0x83f, 0x936, 0xe3a, 0xf33, 0xc39, 0xd30,
+    0x3a0, 0x2a9, 0x1a3, 0xaa , 0x7a6, 0x6af, 0x5a5, 0x4ac,
+    0xbac, 0xaa5, 0x9af, 0x8a6, 0xfaa, 0xea3, 0xda9, 0xca0,
+    0x460, 0x569, 0x663, 0x76a, 0x66 , 0x16f, 0x265, 0x36c,
+    0xc6c, 0xd65, 0xe6f, 0xf66, 0x86a, 0x963, 0xa69, 0xb60,
+    0x5f0, 0x4f9, 0x7f3, 0x6fa, 0x1f6, 0xff , 0x3f5, 0x2fc,
+    0xdfc, 0xcf5, 0xfff, 0xef6, 0x9fa, 0x8f3, 0xbf9, 0xaf0,
+    0x650, 0x759, 0x453, 0x55a, 0x256, 0x35f, 0x55 , 0x15c,
+    0xe5c, 0xf55, 0xc5f, 0xd56, 0xa5a, 0xb53, 0x859, 0x950,
+    0x7c0, 0x6c9, 0x5c3, 0x4ca, 0x3c6, 0x2cf, 0x1c5, 0xcc ,
+    0xfcc, 0xec5, 0xdcf, 0xcc6, 0xbca, 0xac3, 0x9c9, 0x8c0,
+    0x8c0, 0x9c9, 0xac3, 0xbca, 0xcc6, 0xdcf, 0xec5, 0xfcc,
+    0xcc , 0x1c5, 0x2cf, 0x3c6, 0x4ca, 0x5c3, 0x6c9, 0x7c0,
+    0x950, 0x859, 0xb53, 0xa5a, 0xd56, 0xc5f, 0xf55, 0xe5c,
+    0x15c, 0x55 , 0x35f, 0x256, 0x55a, 0x453, 0x759, 0x650,
+    0xaf0, 0xbf9, 0x8f3, 0x9fa, 0xef6, 0xfff, 0xcf5, 0xdfc,
+    0x2fc, 0x3f5, 0xff , 0x1f6, 0x6fa, 0x7f3, 0x4f9, 0x5f0,
+    0xb60, 0xa69, 0x963, 0x86a, 0xf66, 0xe6f, 0xd65, 0xc6c,
+    0x36c, 0x265, 0x16f, 0x66 , 0x76a, 0x663, 0x569, 0x460,
+    0xca0, 0xda9, 0xea3, 0xfaa, 0x8a6, 0x9af, 0xaa5, 0xbac,
+    0x4ac, 0x5a5, 0x6af, 0x7a6, 0xaa , 0x1a3, 0x2a9, 0x3a0,
+    0xd30, 0xc39, 0xf33, 0xe3a, 0x936, 0x83f, 0xb35, 0xa3c,
+    0x53c, 0x435, 0x73f, 0x636, 0x13a, 0x33 , 0x339, 0x230,
+    0xe90, 0xf99, 0xc93, 0xd9a, 0xa96, 0xb9f, 0x895, 0x99c,
+    0x69c, 0x795, 0x49f, 0x596, 0x29a, 0x393, 0x99 , 0x190,
+    0xf00, 0xe09, 0xd03, 0xc0a, 0xb06, 0xa0f, 0x905, 0x80c,
+    0x70c, 0x605, 0x50f, 0x406, 0x30a, 0x203, 0x109, 0x0
+  };
+
+  char MarchingCubeTables::myTriTable[256][16] = {
+    {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 0,  8,  3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 0,  1,  9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 1,  8,  3,  9,  8,  1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 1,  2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 0,  8,  3,  1,  2, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 9,  2, 10,  0,  2,  9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 2,  8,  3,  2, 10,  8, 10,  9,  8, -1, -1, -1, -1, -1, -1, -1},
+    { 3, 11,  2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 0, 11,  2,  8, 11,  0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 1,  9,  0,  2,  3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 1, 11,  2,  1,  9, 11,  9,  8, 11, -1, -1, -1, -1, -1, -1, -1},
+    { 3, 10,  1, 11, 10,  3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 0, 10,  1,  0,  8, 10,  8, 11, 10, -1, -1, -1, -1, -1, -1, -1},
+    { 3,  9,  0,  3, 11,  9, 11, 10,  9, -1, -1, -1, -1, -1, -1, -1},
+    { 9,  8, 10, 10,  8, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 4,  7,  8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 4,  3,  0,  7,  3,  4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 0,  1,  9,  8,  4,  7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 4,  1,  9,  4,  7,  1,  7,  3,  1, -1, -1, -1, -1, -1, -1, -1},
+    { 1,  2, 10,  8,  4,  7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 3,  4,  7,  3,  0,  4,  1,  2, 10, -1, -1, -1, -1, -1, -1, -1},
+    { 9,  2, 10,  9,  0,  2,  8,  4,  7, -1, -1, -1, -1, -1, -1, -1},
+    { 2, 10,  9,  2,  9,  7,  2,  7,  3,  7,  9,  4, -1, -1, -1, -1},
+    { 8,  4,  7,  3, 11,  2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {11,  4,  7, 11,  2,  4,  2,  0,  4, -1, -1, -1, -1, -1, -1, -1},
+    { 9,  0,  1,  8,  4,  7,  2,  3, 11, -1, -1, -1, -1, -1, -1, -1},
+    { 4,  7, 11,  9,  4, 11,  9, 11,  2,  9,  2,  1, -1, -1, -1, -1},
+    { 3, 10,  1,  3, 11, 10,  7,  8,  4, -1, -1, -1, -1, -1, -1, -1},
+    { 1, 11, 10,  1,  4, 11,  1,  0,  4,  7, 11,  4, -1, -1, -1, -1},
+    { 4,  7,  8,  9,  0, 11,  9, 11, 10, 11,  0,  3, -1, -1, -1, -1},
+    { 4,  7, 11,  4, 11,  9,  9, 11, 10, -1, -1, -1, -1, -1, -1, -1},
+    { 9,  5,  4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 9,  5,  4,  0,  8,  3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 0,  5,  4,  1,  5,  0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 8,  5,  4,  8,  3,  5,  3,  1,  5, -1, -1, -1, -1, -1, -1, -1},
+    { 1,  2, 10,  9,  5,  4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 3,  0,  8,  1,  2, 10,  4,  9,  5, -1, -1, -1, -1, -1, -1, -1},
+    { 5,  2, 10,  5,  4,  2,  4,  0,  2, -1, -1, -1, -1, -1, -1, -1},
+    { 2, 10,  5,  3,  2,  5,  3,  5,  4,  3,  4,  8, -1, -1, -1, -1},
+    { 9,  5,  4,  2,  3, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 0, 11,  2,  0,  8, 11,  4,  9,  5, -1, -1, -1, -1, -1, -1, -1},
+    { 0,  5,  4,  0,  1,  5,  2,  3, 11, -1, -1, -1, -1, -1, -1, -1},
+    { 2,  1,  5,  2,  5,  8,  2,  8, 11,  4,  8,  5, -1, -1, -1, -1},
+    {10,  3, 11, 10,  1,  3,  9,  5,  4, -1, -1, -1, -1, -1, -1, -1},
+    { 4,  9,  5,  0,  8,  1,  8, 10,  1,  8, 11, 10, -1, -1, -1, -1},
+    { 5,  4,  0,  5,  0, 11,  5, 11, 10, 11,  0,  3, -1, -1, -1, -1},
+    { 5,  4,  8,  5,  8, 10, 10,  8, 11, -1, -1, -1, -1, -1, -1, -1},
+    { 9,  7,  8,  5,  7,  9, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 9,  3,  0,  9,  5,  3,  5,  7,  3, -1, -1, -1, -1, -1, -1, -1},
+    { 0,  7,  8,  0,  1,  7,  1,  5,  7, -1, -1, -1, -1, -1, -1, -1},
+    { 1,  5,  3,  3,  5,  7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 9,  7,  8,  9,  5,  7, 10,  1,  2, -1, -1, -1, -1, -1, -1, -1},
+    {10,  1,  2,  9,  5,  0,  5,  3,  0,  5,  7,  3, -1, -1, -1, -1},
+    { 8,  0,  2,  8,  2,  5,  8,  5,  7, 10,  5,  2, -1, -1, -1, -1},
+    { 2, 10,  5,  2,  5,  3,  3,  5,  7, -1, -1, -1, -1, -1, -1, -1},
+    { 7,  9,  5,  7,  8,  9,  3, 11,  2, -1, -1, -1, -1, -1, -1, -1},
+    { 9,  5,  7,  9,  7,  2,  9,  2,  0,  2,  7, 11, -1, -1, -1, -1},
+    { 2,  3, 11,  0,  1,  8,  1,  7,  8,  1,  5,  7, -1, -1, -1, -1},
+    {11,  2,  1, 11,  1,  7,  7,  1,  5, -1, -1, -1, -1, -1, -1, -1},
+    { 9,  5,  8,  8,  5,  7, 10,  1,  3, 10,  3, 11, -1, -1, -1, -1},
+    { 5,  7,  0,  5,  0,  9,  7, 11,  0,  1,  0, 10, 11, 10,  0, -1},
+    {11, 10,  0, 11,  0,  3, 10,  5,  0,  8,  0,  7,  5,  7,  0, -1},
+    {11, 10,  5,  7, 11,  5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {10,  6,  5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 0,  8,  3,  5, 10,  6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 9,  0,  1,  5, 10,  6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 1,  8,  3,  1,  9,  8,  5, 10,  6, -1, -1, -1, -1, -1, -1, -1},
+    { 1,  6,  5,  2,  6,  1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 1,  6,  5,  1,  2,  6,  3,  0,  8, -1, -1, -1, -1, -1, -1, -1},
+    { 9,  6,  5,  9,  0,  6,  0,  2,  6, -1, -1, -1, -1, -1, -1, -1},
+    { 5,  9,  8,  5,  8,  2,  5,  2,  6,  3,  2,  8, -1, -1, -1, -1},
+    { 2,  3, 11, 10,  6,  5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {11,  0,  8, 11,  2,  0, 10,  6,  5, -1, -1, -1, -1, -1, -1, -1},
+    { 0,  1,  9,  2,  3, 11,  5, 10,  6, -1, -1, -1, -1, -1, -1, -1},
+    { 5, 10,  6,  1,  9,  2,  9, 11,  2,  9,  8, 11, -1, -1, -1, -1},
+    { 6,  3, 11,  6,  5,  3,  5,  1,  3, -1, -1, -1, -1, -1, -1, -1},
+    { 0,  8, 11,  0, 11,  5,  0,  5,  1,  5, 11,  6, -1, -1, -1, -1},
+    { 3, 11,  6,  0,  3,  6,  0,  6,  5,  0,  5,  9, -1, -1, -1, -1},
+    { 6,  5,  9,  6,  9, 11, 11,  9,  8, -1, -1, -1, -1, -1, -1, -1},
+    { 5, 10,  6,  4,  7,  8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 4,  3,  0,  4,  7,  3,  6,  5, 10, -1, -1, -1, -1, -1, -1, -1},
+    { 1,  9,  0,  5, 10,  6,  8,  4,  7, -1, -1, -1, -1, -1, -1, -1},
+    {10,  6,  5,  1,  9,  7,  1,  7,  3,  7,  9,  4, -1, -1, -1, -1},
+    { 6,  1,  2,  6,  5,  1,  4,  7,  8, -1, -1, -1, -1, -1, -1, -1},
+    { 1,  2,  5,  5,  2,  6,  3,  0,  4,  3,  4,  7, -1, -1, -1, -1},
+    { 8,  4,  7,  9,  0,  5,  0,  6,  5,  0,  2,  6, -1, -1, -1, -1},
+    { 7,  3,  9,  7,  9,  4,  3,  2,  9,  5,  9,  6,  2,  6,  9, -1},
+    { 3, 11,  2,  7,  8,  4, 10,  6,  5, -1, -1, -1, -1, -1, -1, -1},
+    { 5, 10,  6,  4,  7,  2,  4,  2,  0,  2,  7, 11, -1, -1, -1, -1},
+    { 0,  1,  9,  4,  7,  8,  2,  3, 11,  5, 10,  6, -1, -1, -1, -1},
+    { 9,  2,  1,  9, 11,  2,  9,  4, 11,  7, 11,  4,  5, 10,  6, -1},
+    { 8,  4,  7,  3, 11,  5,  3,  5,  1,  5, 11,  6, -1, -1, -1, -1},
+    { 5,  1, 11,  5, 11,  6,  1,  0, 11,  7, 11,  4,  0,  4, 11, -1},
+    { 0,  5,  9,  0,  6,  5,  0,  3,  6, 11,  6,  3,  8,  4,  7, -1},
+    { 6,  5,  9,  6,  9, 11,  4,  7,  9,  7, 11,  9, -1, -1, -1, -1},
+    {10,  4,  9,  6,  4, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 4, 10,  6,  4,  9, 10,  0,  8,  3, -1, -1, -1, -1, -1, -1, -1},
+    {10,  0,  1, 10,  6,  0,  6,  4,  0, -1, -1, -1, -1, -1, -1, -1},
+    { 8,  3,  1,  8,  1,  6,  8,  6,  4,  6,  1, 10, -1, -1, -1, -1},
+    { 1,  4,  9,  1,  2,  4,  2,  6,  4, -1, -1, -1, -1, -1, -1, -1},
+    { 3,  0,  8,  1,  2,  9,  2,  4,  9,  2,  6,  4, -1, -1, -1, -1},
+    { 0,  2,  4,  4,  2,  6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 8,  3,  2,  8,  2,  4,  4,  2,  6, -1, -1, -1, -1, -1, -1, -1},
+    {10,  4,  9, 10,  6,  4, 11,  2,  3, -1, -1, -1, -1, -1, -1, -1},
+    { 0,  8,  2,  2,  8, 11,  4,  9, 10,  4, 10,  6, -1, -1, -1, -1},
+    { 3, 11,  2,  0,  1,  6,  0,  6,  4,  6,  1, 10, -1, -1, -1, -1},
+    { 6,  4,  1,  6,  1, 10,  4,  8,  1,  2,  1, 11,  8, 11,  1, -1},
+    { 9,  6,  4,  9,  3,  6,  9,  1,  3, 11,  6,  3, -1, -1, -1, -1},
+    { 8, 11,  1,  8,  1,  0, 11,  6,  1,  9,  1,  4,  6,  4,  1, -1},
+    { 3, 11,  6,  3,  6,  0,  0,  6,  4, -1, -1, -1, -1, -1, -1, -1},
+    { 6,  4,  8, 11,  6,  8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 7, 10,  6,  7,  8, 10,  8,  9, 10, -1, -1, -1, -1, -1, -1, -1},
+    { 0,  7,  3,  0, 10,  7,  0,  9, 10,  6,  7, 10, -1, -1, -1, -1},
+    {10,  6,  7,  1, 10,  7,  1,  7,  8,  1,  8,  0, -1, -1, -1, -1},
+    {10,  6,  7, 10,  7,  1,  1,  7,  3, -1, -1, -1, -1, -1, -1, -1},
+    { 1,  2,  6,  1,  6,  8,  1,  8,  9,  8,  6,  7, -1, -1, -1, -1},
+    { 2,  6,  9,  2,  9,  1,  6,  7,  9,  0,  9,  3,  7,  3,  9, -1},
+    { 7,  8,  0,  7,  0,  6,  6,  0,  2, -1, -1, -1, -1, -1, -1, -1},
+    { 7,  3,  2,  6,  7,  2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 2,  3, 11, 10,  6,  8, 10,  8,  9,  8,  6,  7, -1, -1, -1, -1},
+    { 2,  0,  7,  2,  7, 11,  0,  9,  7,  6,  7, 10,  9, 10,  7, -1},
+    { 1,  8,  0,  1,  7,  8,  1, 10,  7,  6,  7, 10,  2,  3, 11, -1},
+    {11,  2,  1, 11,  1,  7, 10,  6,  1,  6,  7,  1, -1, -1, -1, -1},
+    { 8,  9,  6,  8,  6,  7,  9,  1,  6, 11,  6,  3,  1,  3,  6, -1},
+    { 0,  9,  1, 11,  6,  7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 7,  8,  0,  7,  0,  6,  3, 11,  0, 11,  6,  0, -1, -1, -1, -1},
+    { 7, 11,  6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 7,  6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 3,  0,  8, 11,  7,  6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 0,  1,  9, 11,  7,  6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 8,  1,  9,  8,  3,  1, 11,  7,  6, -1, -1, -1, -1, -1, -1, -1},
+    {10,  1,  2,  6, 11,  7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 1,  2, 10,  3,  0,  8,  6, 11,  7, -1, -1, -1, -1, -1, -1, -1},
+    { 2,  9,  0,  2, 10,  9,  6, 11,  7, -1, -1, -1, -1, -1, -1, -1},
+    { 6, 11,  7,  2, 10,  3, 10,  8,  3, 10,  9,  8, -1, -1, -1, -1},
+    { 7,  2,  3,  6,  2,  7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 7,  0,  8,  7,  6,  0,  6,  2,  0, -1, -1, -1, -1, -1, -1, -1},
+    { 2,  7,  6,  2,  3,  7,  0,  1,  9, -1, -1, -1, -1, -1, -1, -1},
+    { 1,  6,  2,  1,  8,  6,  1,  9,  8,  8,  7,  6, -1, -1, -1, -1},
+    {10,  7,  6, 10,  1,  7,  1,  3,  7, -1, -1, -1, -1, -1, -1, -1},
+    {10,  7,  6,  1,  7, 10,  1,  8,  7,  1,  0,  8, -1, -1, -1, -1},
+    { 0,  3,  7,  0,  7, 10,  0, 10,  9,  6, 10,  7, -1, -1, -1, -1},
+    { 7,  6, 10,  7, 10,  8,  8, 10,  9, -1, -1, -1, -1, -1, -1, -1},
+    { 6,  8,  4, 11,  8,  6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 3,  6, 11,  3,  0,  6,  0,  4,  6, -1, -1, -1, -1, -1, -1, -1},
+    { 8,  6, 11,  8,  4,  6,  9,  0,  1, -1, -1, -1, -1, -1, -1, -1},
+    { 9,  4,  6,  9,  6,  3,  9,  3,  1, 11,  3,  6, -1, -1, -1, -1},
+    { 6,  8,  4,  6, 11,  8,  2, 10,  1, -1, -1, -1, -1, -1, -1, -1},
+    { 1,  2, 10,  3,  0, 11,  0,  6, 11,  0,  4,  6, -1, -1, -1, -1},
+    { 4, 11,  8,  4,  6, 11,  0,  2,  9,  2, 10,  9, -1, -1, -1, -1},
+    {10,  9,  3, 10,  3,  2,  9,  4,  3, 11,  3,  6,  4,  6,  3, -1},
+    { 8,  2,  3,  8,  4,  2,  4,  6,  2, -1, -1, -1, -1, -1, -1, -1},
+    { 0,  4,  2,  4,  6,  2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 1,  9,  0,  2,  3,  4,  2,  4,  6,  4,  3,  8, -1, -1, -1, -1},
+    { 1,  9,  4,  1,  4,  2,  2,  4,  6, -1, -1, -1, -1, -1, -1, -1},
+    { 8,  1,  3,  8,  6,  1,  8,  4,  6,  6, 10,  1, -1, -1, -1, -1},
+    {10,  1,  0, 10,  0,  6,  6,  0,  4, -1, -1, -1, -1, -1, -1, -1},
+    { 4,  6,  3,  4,  3,  8,  6, 10,  3,  0,  3,  9, 10,  9,  3, -1},
+    {10,  9,  4,  6, 10,  4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 4,  9,  5,  7,  6, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 0,  8,  3,  4,  9,  5, 11,  7,  6, -1, -1, -1, -1, -1, -1, -1},
+    { 5,  0,  1,  5,  4,  0,  7,  6, 11, -1, -1, -1, -1, -1, -1, -1},
+    {11,  7,  6,  8,  3,  4,  3,  5,  4,  3,  1,  5, -1, -1, -1, -1},
+    { 9,  5,  4, 10,  1,  2,  7,  6, 11, -1, -1, -1, -1, -1, -1, -1},
+    { 6, 11,  7,  1,  2, 10,  0,  8,  3,  4,  9,  5, -1, -1, -1, -1},
+    { 7,  6, 11,  5,  4, 10,  4,  2, 10,  4,  0,  2, -1, -1, -1, -1},
+    { 3,  4,  8,  3,  5,  4,  3,  2,  5, 10,  5,  2, 11,  7,  6, -1},
+    { 7,  2,  3,  7,  6,  2,  5,  4,  9, -1, -1, -1, -1, -1, -1, -1},
+    { 9,  5,  4,  0,  8,  6,  0,  6,  2,  6,  8,  7, -1, -1, -1, -1},
+    { 3,  6,  2,  3,  7,  6,  1,  5,  0,  5,  4,  0, -1, -1, -1, -1},
+    { 6,  2,  8,  6,  8,  7,  2,  1,  8,  4,  8,  5,  1,  5,  8, -1},
+    { 9,  5,  4, 10,  1,  6,  1,  7,  6,  1,  3,  7, -1, -1, -1, -1},
+    { 1,  6, 10,  1,  7,  6,  1,  0,  7,  8,  7,  0,  9,  5,  4, -1},
+    { 4,  0, 10,  4, 10,  5,  0,  3, 10,  6, 10,  7,  3,  7, 10, -1},
+    { 7,  6, 10,  7, 10,  8,  5,  4, 10,  4,  8, 10, -1, -1, -1, -1},
+    { 6,  9,  5,  6, 11,  9, 11,  8,  9, -1, -1, -1, -1, -1, -1, -1},
+    { 3,  6, 11,  0,  6,  3,  0,  5,  6,  0,  9,  5, -1, -1, -1, -1},
+    { 0, 11,  8,  0,  5, 11,  0,  1,  5,  5,  6, 11, -1, -1, -1, -1},
+    { 6, 11,  3,  6,  3,  5,  5,  3,  1, -1, -1, -1, -1, -1, -1, -1},
+    { 1,  2, 10,  9,  5, 11,  9, 11,  8, 11,  5,  6, -1, -1, -1, -1},
+    { 0, 11,  3,  0,  6, 11,  0,  9,  6,  5,  6,  9,  1,  2, 10, -1},
+    {11,  8,  5, 11,  5,  6,  8,  0,  5, 10,  5,  2,  0,  2,  5, -1},
+    { 6, 11,  3,  6,  3,  5,  2, 10,  3, 10,  5,  3, -1, -1, -1, -1},
+    { 5,  8,  9,  5,  2,  8,  5,  6,  2,  3,  8,  2, -1, -1, -1, -1},
+    { 9,  5,  6,  9,  6,  0,  0,  6,  2, -1, -1, -1, -1, -1, -1, -1},
+    { 1,  5,  8,  1,  8,  0,  5,  6,  8,  3,  8,  2,  6,  2,  8, -1},
+    { 1,  5,  6,  2,  1,  6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 1,  3,  6,  1,  6, 10,  3,  8,  6,  5,  6,  9,  8,  9,  6, -1},
+    {10,  1,  0, 10,  0,  6,  9,  5,  0,  5,  6,  0, -1, -1, -1, -1},
+    { 0,  3,  8,  5,  6, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {10,  5,  6, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {11,  5, 10,  7,  5, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {11,  5, 10, 11,  7,  5,  8,  3,  0, -1, -1, -1, -1, -1, -1, -1},
+    { 5, 11,  7,  5, 10, 11,  1,  9,  0, -1, -1, -1, -1, -1, -1, -1},
+    {10,  7,  5, 10, 11,  7,  9,  8,  1,  8,  3,  1, -1, -1, -1, -1},
+    {11,  1,  2, 11,  7,  1,  7,  5,  1, -1, -1, -1, -1, -1, -1, -1},
+    { 0,  8,  3,  1,  2,  7,  1,  7,  5,  7,  2, 11, -1, -1, -1, -1},
+    { 9,  7,  5,  9,  2,  7,  9,  0,  2,  2, 11,  7, -1, -1, -1, -1},
+    { 7,  5,  2,  7,  2, 11,  5,  9,  2,  3,  2,  8,  9,  8,  2, -1},
+    { 2,  5, 10,  2,  3,  5,  3,  7,  5, -1, -1, -1, -1, -1, -1, -1},
+    { 8,  2,  0,  8,  5,  2,  8,  7,  5, 10,  2,  5, -1, -1, -1, -1},
+    { 9,  0,  1,  5, 10,  3,  5,  3,  7,  3, 10,  2, -1, -1, -1, -1},
+    { 9,  8,  2,  9,  2,  1,  8,  7,  2, 10,  2,  5,  7,  5,  2, -1},
+    { 1,  3,  5,  3,  7,  5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 0,  8,  7,  0,  7,  1,  1,  7,  5, -1, -1, -1, -1, -1, -1, -1},
+    { 9,  0,  3,  9,  3,  5,  5,  3,  7, -1, -1, -1, -1, -1, -1, -1},
+    { 9,  8,  7,  5,  9,  7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 5,  8,  4,  5, 10,  8, 10, 11,  8, -1, -1, -1, -1, -1, -1, -1},
+    { 5,  0,  4,  5, 11,  0,  5, 10, 11, 11,  3,  0, -1, -1, -1, -1},
+    { 0,  1,  9,  8,  4, 10,  8, 10, 11, 10,  4,  5, -1, -1, -1, -1},
+    {10, 11,  4, 10,  4,  5, 11,  3,  4,  9,  4,  1,  3,  1,  4, -1},
+    { 2,  5,  1,  2,  8,  5,  2, 11,  8,  4,  5,  8, -1, -1, -1, -1},
+    { 0,  4, 11,  0, 11,  3,  4,  5, 11,  2, 11,  1,  5,  1, 11, -1},
+    { 0,  2,  5,  0,  5,  9,  2, 11,  5,  4,  5,  8, 11,  8,  5, -1},
+    { 9,  4,  5,  2, 11,  3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 2,  5, 10,  3,  5,  2,  3,  4,  5,  3,  8,  4, -1, -1, -1, -1},
+    { 5, 10,  2,  5,  2,  4,  4,  2,  0, -1, -1, -1, -1, -1, -1, -1},
+    { 3, 10,  2,  3,  5, 10,  3,  8,  5,  4,  5,  8,  0,  1,  9, -1},
+    { 5, 10,  2,  5,  2,  4,  1,  9,  2,  9,  4,  2, -1, -1, -1, -1},
+    { 8,  4,  5,  8,  5,  3,  3,  5,  1, -1, -1, -1, -1, -1, -1, -1},
+    { 0,  4,  5,  1,  0,  5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 8,  4,  5,  8,  5,  3,  9,  0,  5,  0,  3,  5, -1, -1, -1, -1},
+    { 9,  4,  5, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 4, 11,  7,  4,  9, 11,  9, 10, 11, -1, -1, -1, -1, -1, -1, -1},
+    { 0,  8,  3,  4,  9,  7,  9, 11,  7,  9, 10, 11, -1, -1, -1, -1},
+    { 1, 10, 11,  1, 11,  4,  1,  4,  0,  7,  4, 11, -1, -1, -1, -1},
+    { 3,  1,  4,  3,  4,  8,  1, 10,  4,  7,  4, 11, 10, 11,  4, -1},
+    { 4, 11,  7,  9, 11,  4,  9,  2, 11,  9,  1,  2, -1, -1, -1, -1},
+    { 9,  7,  4,  9, 11,  7,  9,  1, 11,  2, 11,  1,  0,  8,  3, -1},
+    {11,  7,  4, 11,  4,  2,  2,  4,  0, -1, -1, -1, -1, -1, -1, -1},
+    {11,  7,  4, 11,  4,  2,  8,  3,  4,  3,  2,  4, -1, -1, -1, -1},
+    { 2,  9, 10,  2,  7,  9,  2,  3,  7,  7,  4,  9, -1, -1, -1, -1},
+    { 9, 10,  7,  9,  7,  4, 10,  2,  7,  8,  7,  0,  2,  0,  7, -1},
+    { 3,  7, 10,  3, 10,  2,  7,  4, 10,  1, 10,  0,  4,  0, 10, -1},
+    { 1, 10,  2,  8,  7,  4, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 4,  9,  1,  4,  1,  7,  7,  1,  3, -1, -1, -1, -1, -1, -1, -1},
+    { 4,  9,  1,  4,  1,  7,  0,  8,  1,  8,  7,  1, -1, -1, -1, -1},
+    { 4,  0,  3,  7,  4,  3, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 4,  8,  7, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 9, 10,  8, 10, 11,  8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 3,  0,  9,  3,  9, 11, 11,  9, 10, -1, -1, -1, -1, -1, -1, -1},
+    { 0,  1, 10,  0, 10,  8,  8, 10, 11, -1, -1, -1, -1, -1, -1, -1},
+    { 3,  1, 10, 11,  3, 10, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 1,  2, 11,  1, 11,  9,  9, 11,  8, -1, -1, -1, -1, -1, -1, -1},
+    { 3,  0,  9,  3,  9, 11,  1,  2,  9,  2, 11,  9, -1, -1, -1, -1},
+    { 0,  2, 11,  8,  0, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 3,  2, 11, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 2,  3,  8,  2,  8, 10, 10,  8,  9, -1, -1, -1, -1, -1, -1, -1},
+    { 9, 10,  2,  0,  9,  2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 2,  3,  8,  2,  8, 10,  0,  1,  8,  1, 10,  8, -1, -1, -1, -1},
+    { 1, 10,  2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 1,  3,  8,  9,  1,  8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 0,  9,  1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    { 0,  3,  8, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
+    {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}
+  };
+
+  //==============================================================================
+  //function : EdgeTable
+  //purpose  :
+  //==============================================================================
+  int (&MarchingCubeTables::EdgeTable())[256]
+  {
+    return myEdgeTable;
+  }
+
+  //==============================================================================
+  //function : TriTable
+  //purpose  :
+  //==============================================================================
+  char (&MarchingCubeTables::TriTable())[256][16]
+  {
+    return myTriTable;
+  }
+}
--- a/src/BVH/BVH_MarchingCubes.hxx
+++ b/src/BVH/BVH_MarchingCubes.hxx
@@ -0,0 +1,76 @@
+// Created on: 2016-04-07
+// Created by: Denis BOGOLEPOV
+// Copyright (c) 2013-2016 OPEN CASCADE SAS
+//
+// This file is part of Open CASCADE Technology software library.
+//
+// This library is free software; you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License version 2.1 as published
+// by the Free Software Foundation, with special exception defined in the file
+// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
+// distribution for complete text of the license and disclaimer of any warranty.
+//
+// Alternatively, this file may be used under the terms of Open CASCADE
+// commercial license or contractual agreement.
+
+#ifndef _BVH_MarchingCubes_Header
+#define _BVH_MarchingCubes_Header
+
+#include <BVH_Tessellator.hxx>
+
+//! Implements marching cubes (MC) tessellation algorithm.
+template<class T, int N>
+class BVH_MarchingCubes : public BVH_Tessellator<T, N>
+{
+public:
+
+  typedef typename BVH::VectorType<T, N>::Type BVH_VecNt;
+
+public:
+
+  //! Creates marching cubes (MC) tessellator.
+  BVH_MarchingCubes (const BVH_Box<T, N>& theBox, const Standard_Integer theMaxSlices = 128);
+
+  //! Generates tessellation for the given value of implicit function.
+  virtual Standard_Boolean Perform (BVH_ImplicitFunction<T, N>& theFunction, const T theValue);
+
+protected:
+
+  //! Returns minimum corner of the given voxel.
+  BVH_VecNt VoxelCorner (const int theX,
+                         const int theY,
+                         const int theZ) const;
+
+  //! Interpolates the intersection point based function values.
+  BVH_VecNt VertexInter (const BVH_VecNt& thePnt1,
+                         const BVH_VecNt& thePnt2,
+                         const T          theVal1,
+                         const T          theVal2) const;
+};
+
+namespace BVH
+{
+  //! Tool object providing pre-computed tables for MC tessellator.
+  class MarchingCubeTables
+  {
+  public:
+
+    //! Returns table mapping the vertices under the surface to the intersecting edges.
+    Standard_EXPORT static int (&EdgeTable())[256];
+
+    //! Returns table used for generation of triangles representing the surface within a cell.
+    Standard_EXPORT static char (&TriTable())[256][16];
+
+  private:
+
+    //! Table that maps the vertices under the surface to the intersecting edges.
+    static int myEdgeTable[256];
+
+    //! Table used for generation of triangles representing the surface within a cell.
+    static char myTriTable[256][16];
+  };
+}
+
+#include <BVH_MarchingCubes.lxx>
+
+#endif // _BVH_MarchingCubes_Header
--- a/src/BVH/BVH_MarchingCubes.lxx
+++ b/src/BVH/BVH_MarchingCubes.lxx
@@ -0,0 +1,254 @@
+// Created on: 2016-04-07
+// Created by: Denis BOGOLEPOV
+// Copyright (c) 2013-2016 OPEN CASCADE SAS
+//
+// This file is part of Open CASCADE Technology software library.
+//
+// This library is free software; you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License version 2.1 as published
+// by the Free Software Foundation, with special exception defined in the file
+// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
+// distribution for complete text of the license and disclaimer of any warranty.
+//
+// Alternatively, this file may be used under the terms of Open CASCADE
+// commercial license or contractual agreement.
+
+// =======================================================================
+// function : BVH_MarchingCubes
+// purpose  : Creates marching cubes (MC) tessellator
+// =======================================================================
+template<class T, int N>
+BVH_MarchingCubes<T, N>::BVH_MarchingCubes (const BVH_Box<T, N>& theBox, const Standard_Integer theMaxSlices)
+: BVH_Tessellator<T, N> (theBox, theMaxSlices)
+{
+  //
+}
+
+// =======================================================================
+// function : VoxelCorner
+// purpose  :
+// =======================================================================
+template<class T, int N>
+typename BVH_MarchingCubes<T, N>::BVH_VecNt BVH_MarchingCubes<T, N>::VoxelCorner (const int theX,
+                                                                                  const int theY,
+                                                                                  const int theZ) const
+{
+  BVH_VecNt aCorner = myCornerMin;
+
+  aCorner.x() += myVoxelSize.x() * theX;
+  aCorner.y() += myVoxelSize.y() * theY;
+  aCorner.z() += myVoxelSize.z() * theZ;
+
+  return aCorner;
+}
+
+// =======================================================================
+// function : VertexInter
+// purpose  :
+// =======================================================================
+template<class T, int N>
+typename BVH_MarchingCubes<T, N>::BVH_VecNt BVH_MarchingCubes<T, N>::VertexInter (const BVH_VecNt& thePnt1,
+                                                                                  const BVH_VecNt& thePnt2,
+                                                                                  const T          theVal1,
+                                                                                  const T          theVal2) const
+{
+  return thePnt1 - (thePnt2 - thePnt1) * (theVal1 / (theVal2 - theVal1));
+}
+
+// =======================================================================
+// function : Perform
+// purpose  : Generates tessellation for the given value
+// =======================================================================
+template<class T, int N>
+Standard_Boolean BVH_MarchingCubes<T, N>::Perform (BVH_ImplicitFunction<T, N>& theFunction, const T theValue)
+{
+  T* aGridValues = new T[(mySlicesX + 1) *
+                         (mySlicesY + 1) *
+                         (mySlicesZ + 1)];
+
+  const int aSX = mySlicesX + 1;
+  const int aSY = mySlicesY + 1;
+
+  for (Standard_Integer aZ = 0; aZ <= mySlicesZ; ++aZ)
+  {
+    for (Standard_Integer aY = 0; aY <= mySlicesY; ++aY)
+    {
+      for (Standard_Integer aX = 0; aX <= mySlicesX; ++aX)
+      {
+        aGridValues[aX + (aY + aZ * aSY) * aSX] = theFunction.Calculate (VoxelCorner (aX, aY, aZ)) - theValue;
+      }
+    }
+  }
+
+  myTessellation = new BVH_Triangulation<T, N>;
+
+  IndexMap aVertexMap; // used to merge vertices
+
+  for (Standard_Integer aZ = 0; aZ < mySlicesZ; ++aZ)
+  {
+    for (Standard_Integer aY = 0; aY < mySlicesY; ++aY)
+    {
+      for (Standard_Integer aX = 0; aX < mySlicesX; ++aX)
+      {
+        // determine the index into the edge table
+
+        int aCubeIndex = 0;
+
+        if (aGridValues[(aX + 0) + ((aY + 0) + (aZ + 0) * aSY) * aSX] < T(0.0))
+          aCubeIndex |= 1;
+
+        if (aGridValues[(aX + 1) + ((aY + 0) + (aZ + 0) * aSY) * aSX] < T(0.0))
+          aCubeIndex |= 2;
+
+        if (aGridValues[(aX + 1) + ((aY + 1) + (aZ + 0) * aSY) * aSX] < T(0.0))
+          aCubeIndex |= 4;
+
+        if (aGridValues[(aX + 0) + ((aY + 1) + (aZ + 0) * aSY) * aSX] < T(0.0))
+          aCubeIndex |= 8;
+
+        if (aGridValues[(aX + 0) + ((aY + 0) + (aZ + 1) * aSY) * aSX] < T(0.0))
+          aCubeIndex |= 16;
+
+        if (aGridValues[(aX + 1) + ((aY + 0) + (aZ + 1) * aSY) * aSX] < T(0.0))
+          aCubeIndex |= 32;
+
+        if (aGridValues[(aX + 1) + ((aY + 1) + (aZ + 1) * aSY) * aSX] < T(0.0))
+          aCubeIndex |= 64;
+
+        if (aGridValues[(aX + 0) + ((aY + 1) + (aZ + 1) * aSY) * aSX] < T(0.0))
+          aCubeIndex |= 128;
+
+        // cube is entirely in/out of the surface
+        if (BVH::MarchingCubeTables::EdgeTable()[aCubeIndex] != 0)
+        {
+          // find the vertices where the surface intersects the cube
+          BVH_VecNt aVertexList[12];
+          
+          if (BVH::MarchingCubeTables::EdgeTable()[aCubeIndex] & 1)
+          {
+            aVertexList[0] = VertexInter (VoxelCorner (aX + 0, aY + 0, aZ + 0),
+                                          VoxelCorner (aX + 1, aY + 0, aZ + 0),
+                                          aGridValues[(aX + 0) + ((aY + 0) + (aZ + 0) * aSY) * aSX],
+                                          aGridValues[(aX + 1) + ((aY + 0) + (aZ + 0) * aSY) * aSX]);
+          }
+
+          if (BVH::MarchingCubeTables::EdgeTable()[aCubeIndex] & 2)
+          {
+            aVertexList[1] = VertexInter (VoxelCorner (aX + 1, aY + 0, aZ + 0),
+                                          VoxelCorner (aX + 1, aY + 1, aZ + 0),
+                                          aGridValues[(aX + 1) + ((aY + 0) + (aZ + 0) * aSY) * aSX],
+                                          aGridValues[(aX + 1) + ((aY + 1) + (aZ + 0) * aSY) * aSX]);
+          }
+
+          if (BVH::MarchingCubeTables::EdgeTable()[aCubeIndex] & 4)
+          {
+            aVertexList[2] = VertexInter (VoxelCorner (aX + 1, aY + 1, aZ + 0),
+                                          VoxelCorner (aX + 0, aY + 1, aZ + 0),
+                                          aGridValues[(aX + 1) + ((aY + 1) + (aZ + 0) * aSY) * aSX],
+                                          aGridValues[(aX + 0) + ((aY + 1) + (aZ + 0) * aSY) * aSX]);
+          }
+
+          if (BVH::MarchingCubeTables::EdgeTable()[aCubeIndex] & 8)
+          {
+            aVertexList[3] = VertexInter (VoxelCorner (aX + 0, aY + 1, aZ + 0),
+                                          VoxelCorner (aX + 0, aY + 0, aZ + 0),
+                                          aGridValues[(aX + 0) + ((aY + 1) + (aZ + 0) * aSY) * aSX],
+                                          aGridValues[(aX + 0) + ((aY + 0) + (aZ + 0) * aSY) * aSX]);
+          }
+
+          if (BVH::MarchingCubeTables::EdgeTable()[aCubeIndex] & 16)
+          {
+            aVertexList[4] = VertexInter (VoxelCorner (aX + 0, aY + 0, aZ + 1),
+                                          VoxelCorner (aX + 1, aY + 0, aZ + 1),
+                                          aGridValues[(aX + 0) + ((aY + 0) + (aZ + 1) * aSY) * aSX],
+                                          aGridValues[(aX + 1) + ((aY + 0) + (aZ + 1) * aSY) * aSX]);
+          }
+
+          if (BVH::MarchingCubeTables::EdgeTable()[aCubeIndex] & 32)
+          {
+            aVertexList[5] = VertexInter (VoxelCorner (aX + 1, aY + 0, aZ + 1),
+                                          VoxelCorner (aX + 1, aY + 1, aZ + 1),
+                                          aGridValues[(aX + 1) + ((aY + 0) + (aZ + 1) * aSY) * aSX],
+                                          aGridValues[(aX + 1) + ((aY + 1) + (aZ + 1) * aSY) * aSX]);
+          }
+
+          if (BVH::MarchingCubeTables::EdgeTable()[aCubeIndex] & 64)
+          {
+            aVertexList[6] = VertexInter (VoxelCorner (aX + 1, aY + 1, aZ + 1),
+                                          VoxelCorner (aX + 0, aY + 1, aZ + 1),
+                                          aGridValues[(aX + 1) + ((aY + 1) + (aZ + 1) * aSY) * aSX],
+                                          aGridValues[(aX + 0) + ((aY + 1) + (aZ + 1) * aSY) * aSX]);
+          }
+
+          if (BVH::MarchingCubeTables::EdgeTable()[aCubeIndex] & 128)
+          {
+            aVertexList[7] = VertexInter (VoxelCorner (aX + 0, aY + 1, aZ + 1),
+                                          VoxelCorner (aX + 0, aY + 0, aZ + 1),
+                                          aGridValues[(aX + 0) + ((aY + 1) + (aZ + 1) * aSY) * aSX],
+                                          aGridValues[(aX + 0) + ((aY + 0) + (aZ + 1) * aSY) * aSX]);
+          }
+
+          if (BVH::MarchingCubeTables::EdgeTable()[aCubeIndex] & 256)
+          {
+            aVertexList[8] = VertexInter (VoxelCorner (aX + 0, aY + 0, aZ + 0),
+                                          VoxelCorner (aX + 0, aY + 0, aZ + 1),
+                                          aGridValues[(aX + 0) + ((aY + 0) + (aZ + 0) * aSY) * aSX],
+                                          aGridValues[(aX + 0) + ((aY + 0) + (aZ + 1) * aSY) * aSX]);
+          }
+
+          if (BVH::MarchingCubeTables::EdgeTable()[aCubeIndex] & 512)
+          {
+            aVertexList[9] = VertexInter (VoxelCorner (aX + 1, aY + 0, aZ + 0),
+                                          VoxelCorner (aX + 1, aY + 0, aZ + 1),
+                                          aGridValues[(aX + 1) + ((aY + 0) + (aZ + 0) * aSY) * aSX],
+                                          aGridValues[(aX + 1) + ((aY + 0) + (aZ + 1) * aSY) * aSX]);
+          }
+
+          if (BVH::MarchingCubeTables::EdgeTable()[aCubeIndex] & 1024)
+          {
+            aVertexList[10] = VertexInter (VoxelCorner (aX + 1, aY + 1, aZ + 0),
+                                           VoxelCorner (aX + 1, aY + 1, aZ + 1),
+                                           aGridValues[(aX + 1) + ((aY + 1) + (aZ + 0) * aSY) * aSX],
+                                           aGridValues[(aX + 1) + ((aY + 1) + (aZ + 1) * aSY) * aSX]);
+          }
+
+          if (BVH::MarchingCubeTables::EdgeTable()[aCubeIndex] & 2048)
+          {
+            aVertexList[11] = VertexInter (VoxelCorner (aX + 0, aY + 1, aZ + 0),
+                                           VoxelCorner (aX + 0, aY + 1, aZ + 1),
+                                           aGridValues[(aX + 0) + ((aY + 1) + (aZ + 0) * aSY) * aSX],
+                                           aGridValues[(aX + 0) + ((aY + 1) + (aZ + 1) * aSY) * aSX]);
+          }
+
+          const char* aTrgIdxs = BVH::MarchingCubeTables::TriTable()[aCubeIndex];
+
+          for (Standard_Integer i = 0; aTrgIdxs[i] != -1; i += 3)
+          {
+            BVH_Vec4i aTrg (0, 0, 0, 0);
+
+            for (Standard_Integer k = 0; k < 3; ++k)
+            {
+              IndexMap::iterator aVrtIter = aVertexMap.find (aVertexList [aTrgIdxs[i + 2 - k]]);
+
+              if (aVrtIter != aVertexMap.end())
+              {
+                aTrg[k] = aVrtIter->second;
+              }
+              else
+              {
+                aVertexMap[aVertexList [aTrgIdxs[i + 2 - k]]] = aTrg[k] =
+                  static_cast<Standard_Integer> (myTessellation->Vertices.size());
+
+                myTessellation->Vertices.push_back (aVertexList [aTrgIdxs[i + 2 - k]]);
+              }
+            }
+
+            myTessellation->Elements.push_back (aTrg);
+          }
+        }
+      }
+    }
+  }
+
+  return Standard_True;
+}
--- a/src/BVH/BVH_Tessellator.hxx
+++ b/src/BVH/BVH_Tessellator.hxx
@@ -0,0 +1,145 @@
+// Created on: 2016-04-07
+// Created by: Denis BOGOLEPOV
+// Copyright (c) 2013-2016 OPEN CASCADE SAS
+//
+// This file is part of Open CASCADE Technology software library.
+//
+// This library is free software; you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License version 2.1 as published
+// by the Free Software Foundation, with special exception defined in the file
+// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
+// distribution for complete text of the license and disclaimer of any warranty.
+//
+// Alternatively, this file may be used under the terms of Open CASCADE
+// commercial license or contractual agreement.
+
+#ifndef _BVH_Tessellator_Header
+#define _BVH_Tessellator_Header
+
+#include <map>
+#include <BVH_Triangulation.hxx>
+#include <BVH_ImplicitFunction.hxx>
+
+//! Tool class for generating tessellation of specific value
+//! of 3D/4D implicit function. Tessellation is generated in
+//! the specified bounding box, and its precision is defined
+//! by the resolution of uniform grid used to split this box
+//! into sub-volumes (voxels).
+template<class T, int N>
+class BVH_Tessellator
+{
+public:
+
+  typedef typename BVH::VectorType<T, N>::Type BVH_VecNt;
+
+public:
+
+  //! Creates tessellator with the given resolution.
+  BVH_Tessellator (const BVH_Box<T, N>& theBox, const Standard_Integer theMaxSlices = 128);
+
+  //! Releases resources of the tessellator.
+  virtual ~BVH_Tessellator() { }
+
+  //! Returns number of slices in X dimension.
+  Standard_Integer SlicesX() const
+  {
+    return mySlicesX;
+  }
+
+  //! Returns number of slices in Y dimension.
+  Standard_Integer SlicesY() const
+  {
+    return mySlicesY;
+  }
+
+  //! Returns number of slices in Z dimension.
+  Standard_Integer SlicesZ() const
+  {
+    return mySlicesZ;
+  }
+
+  //! Returns minimum corner of tessellation AABB.
+  const BVH_VecNt& CornerMin() const
+  {
+    return myCornerMin;
+  }
+
+  //! Returns maximum corner of tessellation AABB.
+  const BVH_VecNt& CornerMax() const
+  {
+    return myCornerMax;
+  }
+
+  //! Returns resulting tessellation (triangulation) of the function.
+  NCollection_Handle<BVH_Triangulation<T, N> > Tessellation() const
+  {
+    return myTessellation;
+  }
+
+  //! Generates tessellation for the given value of implicit function.
+  virtual Standard_Boolean Perform (BVH_ImplicitFunction<T, N>& theFunction, const T theValue) = 0;
+
+protected:
+
+  //! Minimum corner of AABB.
+  BVH_VecNt myCornerMin;
+
+  //! Maximum corner of AABB.
+  BVH_VecNt myCornerMax;
+
+  //! Size of single grid voxel.
+  BVH_VecNt myVoxelSize;
+
+  //! Number of slices in X dimension.
+  Standard_Integer mySlicesX;
+
+  //! Number of slices in Y dimension.
+  Standard_Integer mySlicesY;
+
+  //! Number of slices in Z dimension.
+  Standard_Integer mySlicesZ;
+
+  //! Resulting tessellation (triangulation) of the function.
+  NCollection_Handle<BVH_Triangulation<T, N> > myTessellation;
+
+protected:
+
+  //! Compare functor for 3D/4D vectors.
+  struct VertexComparator
+  {
+    bool operator() (const BVH_VecNt& thePoint1,
+                     const BVH_VecNt& thePoint2) const
+    {
+      if (thePoint1.x() < thePoint2.x() - BVH::Precision<T>::Confusion())
+      {
+        return true;
+      }
+      if (thePoint1.x() > thePoint2.x() + BVH::Precision<T>::Confusion())
+      {
+        return false;
+      }
+      if (thePoint1.y() < thePoint2.y() - BVH::Precision<T>::Confusion())
+      {
+        return true;
+      }
+      if (thePoint1.y() > thePoint2.y() + BVH::Precision<T>::Confusion())
+      {
+        return false;
+      }
+      if (thePoint1.z() < thePoint2.z() - BVH::Precision<T>::Confusion())
+      {
+        return true;
+      }
+
+      return false;
+    }
+  };
+
+  //! Map to store the set of indices for 3D mesh point.
+  typedef std::map<BVH_VecNt, int, VertexComparator> IndexMap;
+
+};
+
+#include <BVH_Tessellator.lxx>
+
+#endif // _BVH_Tessellator_Header
--- a/src/BVH/BVH_Tessellator.lxx
+++ b/src/BVH/BVH_Tessellator.lxx
@@ -0,0 +1,40 @@
+// Created on: 2016-04-07
+// Created by: Denis BOGOLEPOV
+// Copyright (c) 2013-2016 OPEN CASCADE SAS
+//
+// This file is part of Open CASCADE Technology software library.
+//
+// This library is free software; you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License version 2.1 as published
+// by the Free Software Foundation, with special exception defined in the file
+// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
+// distribution for complete text of the license and disclaimer of any warranty.
+//
+// Alternatively, this file may be used under the terms of Open CASCADE
+// commercial license or contractual agreement.
+
+// =======================================================================
+// function : BVH_Tessellator
+// purpose  : Creates base tessellator with the given resolution
+// =======================================================================
+template<class T, int N>
+BVH_Tessellator<T, N>::BVH_Tessellator (const BVH_Box<T, N>& theBox, const Standard_Integer theMaxSlices)
+: myCornerMin (theBox.CornerMin()),
+  myCornerMax (theBox.CornerMax())
+{
+  Standard_STATIC_ASSERT (N == 3 || N == 4);
+
+  const BVH_VecNt aBoxSize = myCornerMax - myCornerMin;
+
+  mySlicesX = static_cast<Standard_Integer> (theMaxSlices * aBoxSize.x() / aBoxSize.maxComp());
+  mySlicesY = static_cast<Standard_Integer> (theMaxSlices * aBoxSize.y() / aBoxSize.maxComp());
+  mySlicesZ = static_cast<Standard_Integer> (theMaxSlices * aBoxSize.z() / aBoxSize.maxComp());
+
+  mySlicesX = Min (theMaxSlices, Max (mySlicesX, 8));
+  mySlicesY = Min (theMaxSlices, Max (mySlicesY, 8));
+  mySlicesZ = Min (theMaxSlices, Max (mySlicesZ, 8));
+
+  myVoxelSize.x() = aBoxSize.x() / mySlicesX;
+  myVoxelSize.y() = aBoxSize.y() / mySlicesY;
+  myVoxelSize.z() = aBoxSize.z() / mySlicesZ;
+}
--- a/src/BVH/BVH_Triangulation.hxx
+++ b/src/BVH/BVH_Triangulation.hxx
@@ -63,6 +63,11 @@ public:
  virtual void Swap (const Standard_Integer theIndex1,
                     const Standard_Integer theIndex2);

+public:
+
+  //! Saves triangulation to OBJ file (for debug).
+  void SaveToOBJ (const char* theFileName);
+
 };

 #include <BVH_Triangulation.lxx>
--- a/src/BVH/BVH_Triangulation.lxx
+++ b/src/BVH/BVH_Triangulation.lxx
@@ -18,7 +18,7 @@
 // purpose  :
 // =======================================================================
 template<class T, int N>
-BVH_Triangulation<T, N>::BVH_Triangulation()
+BVH_Triangulation<T, N>::BVH_Triangulation() : BVH_PrimitiveSet<T, N>()
 {
  //
 }
@@ -99,3 +99,37 @@ void BVH_Triangulation<T, N>::Swap (const Standard_Integer theIndex1,

  std::swap (anIndices1, anIndices2);
 }
+
+#include <fstream>
+
+// =======================================================================
+// function : SaveToOBJ
+// purpose  :
+// =======================================================================
+template<class T, int N>
+void BVH_Triangulation<T, N>::SaveToOBJ (const char* theFileName)
+{
+  if (N > 2)
+  {
+    std::ofstream aFile (theFileName);
+
+    if (aFile.is_open())
+    {
+      for (size_t aVrtIdx = 0; aVrtIdx < Vertices.size(); ++aVrtIdx)
+      {
+        aFile << "v " << Vertices[aVrtIdx][0] << " "
+                      << Vertices[aVrtIdx][1] << " "
+                      << Vertices[aVrtIdx][2] << "\n";
+      }
+
+      for (size_t aFaceIdx = 0; aFaceIdx < Elements.size(); ++aFaceIdx)
+      {
+        aFile << "f " << (Elements[aFaceIdx][0] + 1) << " "
+                      << (Elements[aFaceIdx][1] + 1) << " "
+                      << (Elements[aFaceIdx][2] + 1) << "\n";
+      }
+
+      aFile.close();
+    }
+  }
+}
--- a/src/BVH/BVH_Types.hxx
+++ b/src/BVH/BVH_Types.hxx
@@ -266,6 +266,23 @@ namespace BVH

    return aRes - static_cast<Standard_Integer> (aRes > theValue);
  }
+
+  //! Defines precision for numeric type.
+  template<class T> struct Precision { };
+
+  //! Defines precision for 64-bit FP type.
+  template<> struct Precision<Standard_Real>
+  {
+    //! Returns precision constant.
+    static Standard_Real Confusion() { return 1.0e-7; }
+  };
+
+  //! Defines precision for 32-bit FP type.
+  template<> struct Precision<Standard_ShortReal>
+  {
+    //! Returns precision constant.
+    static Standard_ShortReal Confusion() { return 1.0e-5f; }
+  };
 }

 #endif // _BVH_Types_Header
--- a/src/BVH/FILES
+++ b/src/BVH/FILES
@@ -46,3 +46,10 @@ BVH_QuadTree.lxx
 BVH_Triangulation.hxx
 BVH_Triangulation.lxx
 BVH_Types.hxx
+BVH_ImplicitFunction.hxx
+BVH_ImplicitFunction.lxx
+BVH_Tessellator.hxx
+BVH_Tessellator.lxx
+BVH_MarchingCubes.hxx
+BVH_MarchingCubes.lxx
+BVH_MarchingCubes.cxx