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

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.

diff --git a/src/BVH/BVH.cxx b/src/BVH/BVH.cxx
index 2301f3f68eb06d50dda16d27c040f5b914b5a8f0..209df7c9bad5da5721f2304c13fadced71014a6f 100644 (file)
--- 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>;

diff --git a/src/BVH/BVH_DistanceField.hxx b/src/BVH/BVH_DistanceField.hxx
index e17f0110815281a3e3e6db2467f36ea750843f66..f6991fd2b1ce853b38bcbe0428ad7213e5dc9cce 100644 (file)
--- 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:
+
+  //! 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:
 
-  //! Returns maximum corner of voxel grid.
-  const BVH_VecNt& CornerMax() const
+  //! Compare functor for symmetrical pairs of integers.
+  struct SymmetricalPairComparator
   {
-    return myCornerMax;
-  }
+    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:
 
-  //! 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);
+  //! 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>

diff --git a/src/BVH/BVH_DistanceField.lxx b/src/BVH/BVH_DistanceField.lxx
index c3678bcecf0c6fae1351eef635992ef19bb8abb0..ca9bcc35d7bdf2a5a2fb5e508fce2339869cb3d3 100644 (file)
--- 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,360 +84,371 @@ 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)
-  {
-    Standard_STATIC_ASSERT (N == 3 || N == 4);
+#define BVH_DOT3(A, B) (A.x() * B.x() + A.y() * B.y() + A.z() * B.z())
 
-    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;
+//=======================================================================
+//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 T aABdotAP = BVH_DOT3 (aAB, aAP);
-    const T aACdotAP = BVH_DOT3 (aAC, aAP);
+  const BVH_Vec4i aTrg = myMesh->Elements[theTrgIdx];
 
-    if (aABdotAP <= static_cast<T> (0) && aACdotAP <= static_cast<T> (0))
-    {
-      return aAP;
-    }
+  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 typename VectorType<T, N>::Type aBC = theVertC - theVertB;
-    const typename VectorType<T, N>::Type aBP = thePoint - theVertB;
+  const BVH_VecNt aAB = aVertB - aVertA;
+  const BVH_VecNt aAC = aVertC - aVertA;
+  const BVH_VecNt aAP = thePnt - aVertA;
 
-    const T aBAdotBP = -BVH_DOT3 (aAB, aBP);
-    const T aBCdotBP =  BVH_DOT3 (aBC, aBP);
+  const T aABdotAP = BVH_DOT3 (aAB, aAP);
+  const T aACdotAP = BVH_DOT3 (aAC, aAP);
 
-    if (aBAdotBP <= static_cast<T> (0) && aBCdotBP <= static_cast<T> (0))
-    {
-      return aBP;
-    }
+  if (aABdotAP <= static_cast<T> (0) && aACdotAP <= static_cast<T> (0))
+  {
+    return std::make_pair (aAP, Feature (VERT_ID, aTrg.x()));
+  }
 
-    const typename VectorType<T, N>::Type aCP = thePoint - theVertC;
+  const BVH_VecNt aBC = aVertC - aVertB;
+  const BVH_VecNt aBP = thePnt - aVertB;
 
-    const T aCBdotCP = -BVH_DOT3 (aBC, aCP);
-    const T aCAdotCP = -BVH_DOT3 (aAC, aCP);
+  const T aBAdotBP = -BVH_DOT3 (aAB, aBP);
+  const T aBCdotBP =  BVH_DOT3 (aBC, aBP);
 
-    if (aCAdotCP <= static_cast<T> (0) && aCBdotCP <= static_cast<T> (0))
-    {
-      return aCP;
-    }
+  if (aBAdotBP <= static_cast<T> (0) && aBCdotBP <= static_cast<T> (0))
+  {
+    return std::make_pair (aBP, Feature (VERT_ID, aTrg.y()));
+  }
 
-    const T aACdotBP = BVH_DOT3 (aAC, aBP);
+  const BVH_VecNt aCP = thePnt - aVertC;
 
-    const T aVC = aABdotAP * aACdotBP + aBAdotBP * aACdotAP;
+  const T aCBdotCP = -BVH_DOT3 (aBC, aCP);
+  const T aCAdotCP = -BVH_DOT3 (aAC, aCP);
 
-    if (aVC <= static_cast<T> (0) && aABdotAP >= static_cast<T> (0) && aBAdotBP >= static_cast<T> (0))
-    {
-      return aAP - aAB * (aABdotAP / (aABdotAP + aBAdotBP));
-    }
+  if (aCAdotCP <= static_cast<T> (0) && aCBdotCP <= static_cast<T> (0))
+  {
+    return std::make_pair (aCP, Feature (VERT_ID, aTrg.z()));
+  }
 
-    const T aABdotCP = BVH_DOT3 (aAB, aCP);
+  const T aACdotBP = BVH_DOT3 (aAC, aBP);
 
-    const T aVA = aBAdotBP * aCAdotCP - aABdotCP * aACdotBP;
+  const T aVC = aABdotAP * aACdotBP + aBAdotBP * aACdotAP;
 
-    if (aVA <= static_cast<T> (0) && aBCdotBP >= static_cast<T> (0) && aCBdotCP >= static_cast<T> (0))
-    {
-      return aBP - aBC * (aBCdotBP / (aBCdotBP + aCBdotCP));
-    }
+  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 aVB = aABdotCP * aACdotAP + aABdotAP * aCAdotCP;
+  const T aABdotCP = BVH_DOT3 (aAB, aCP);
 
-    if (aVB <= static_cast<T> (0) && aACdotAP >= static_cast<T> (0) && aCAdotCP >= static_cast<T> (0))
-    {
-      return aAP - aAC * (aACdotAP / (aACdotAP + aCAdotCP));
-    }
+  const T aVA = aBAdotBP * aCAdotCP - aABdotCP * aACdotBP;
 
-    const T aNorm = static_cast<T> (1.0) / (aVA + aVB + aVC);
+  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 aU = aVA * aNorm;
-    const T aV = aVB * aNorm;
+  const T aVB = aABdotCP * aACdotAP + aABdotAP * aCAdotCP;
 
-    return thePoint - (theVertA * aU + theVertB * aV + theVertC * (static_cast<T> (1.0) - aU - aV));
+  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()));
   }
 
-  //=======================================================================
-  //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)
-  {
-    Standard_STATIC_ASSERT (N == 3 || N == 4);
+  const T aNorm = static_cast<T> (1.0) / (aVA + aVB + aVC);
 
-    T aMinDistance = std::numeric_limits<T>::max();
+  const T aU = aVA * aNorm;
+  const T aV = aVB * aNorm;
 
-    BVH_Triangulation<T, N>* aTriangulation =
-      dynamic_cast<BVH_Triangulation<T, N>*> (theObject);
+  const BVH_VecNt aDirect = thePnt - (aVertA * aU + aVertB * aV + aVertC * (static_cast<T> (1.0) - aU - aV));
 
-    Standard_ASSERT_RETURN (aTriangulation != NULL,
-      "Error: Unsupported BVH object (non triangulation)", aMinDistance);
+  return std::make_pair (aDirect, Feature (FACE_ID, theTrgIdx));
+}
 
-    const NCollection_Handle<BVH_Tree<T, N> >& aBVH = aTriangulation->BVH();
+//=======================================================================
+//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);
 
-    if (aBVH.IsNull())
-    {
-      return Standard_False;
-    }
+  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];
 
-    std::pair<Standard_Integer, T> aStack[32];
+  Standard_Integer aHead = -1;
+  Standard_Integer aNode =  0; // root node
 
-    Standard_Integer aHead = -1;
-    Standard_Integer aNode =  0; // root node
+  for (;;)
+  {
+    BVH_Vec4i aData = aBVH->NodeInfoBuffer()[aNode];
 
-    for (;;)
+    if (aData.x() == 0) // if inner node
     {
-      BVH_Vec4i aData = aBVH->NodeInfoBuffer()[aNode];
+      const T aDistToLft = BVH::DistanceToBox<T, N> (thePnt,
+                                                     aBVH->MinPoint (aData.y()),
+                                                     aBVH->MaxPoint (aData.y()));
 
-      if (aData.x() == 0) // if inner node
-      {
-        const T aDistToLft = DistanceToBox<T, N> (thePnt,
-                                                  aBVH->MinPoint (aData.y()),
-                                                  aBVH->MaxPoint (aData.y()));
+      const T aDistToRgh = BVH::DistanceToBox<T, N> (thePnt,
+                                                     aBVH->MinPoint (aData.z()),
+                                                     aBVH->MaxPoint (aData.z()));
 
-        const T aDistToRgh = DistanceToBox<T, N> (thePnt,
-                                                  aBVH->MinPoint (aData.z()),
-                                                  aBVH->MaxPoint (aData.z()));
+      const Standard_Boolean aHitLft = aDistToLft <= aMinDist;
+      const Standard_Boolean aHitRgh = aDistToRgh <= aMinDist;
 
-        const Standard_Boolean aHitLft = aDistToLft <= aMinDistance;
-        const Standard_Boolean aHitRgh = aDistToRgh <= aMinDistance;
+      if (aHitLft & aHitRgh)
+      {
+        aNode = (aDistToLft < aDistToRgh) ? aData.y() : aData.z();
 
-        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 aMinDistance;
-
-            std::pair<Standard_Integer, T>& anInfo = aStack[aHead--];
+          if (aHead < 0)
+            return aMinDist;
 
-            while (anInfo.second > aMinDistance)
-            {
-              if (aHead < 0)
-                return aMinDistance;
+          std::pair<Standard_Integer, T>& anInfo = aStack[aHead--];
 
-              anInfo = aStack[aHead--];
-            }
+          while (anInfo.second > aMinDist)
+          {
+            if (aHead < 0)
+              return aMinDist;
 
-            aNode = anInfo.first;
+            anInfo = aStack[aHead--];
           }
+
+          aNode = anInfo.first;
         }
       }
-      else  // if leaf node
+    }
+    else  // if leaf node
+    {
+      for (Standard_Integer aTrgIdx = aData.y(); aTrgIdx <= aData.z(); ++aTrgIdx)
       {
-        for (Standard_Integer aTrgIdx = aData.y(); aTrgIdx <= aData.z(); ++aTrgIdx)
-        {
-          const BVH_Vec4i aTriangle = aTriangulation->Elements[aTrgIdx];
+        std::pair<BVH_VecNt, Feature> aDirToFeature = directToTrg (thePnt, 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 T aDistance = BVH_DOT3 (aDirToFeature.first,
+                                      aDirToFeature.first);
 
-          const typename VectorType<T, N>::Type aDirection =
-            DirectionToNearestPoint<T, N> (thePnt, aVertex0, aVertex1, aVertex2);
-
-          const T aDistance = BVH_DOT3 (aDirection, aDirection);
+        if (aDistance < aMinDist)
+        {
+          aMinDist = aDistance;
 
-          if (aDistance < aMinDistance)
+          if (myComputeSign)
           {
-            aMinDistance = aDistance;
-
-            typename VectorType<T, N>::Type aTrgEdges[] = { aVertex1 - aVertex0,
-                                                            aVertex2 - aVertex0 };
+            const BVH_VecNt aNormal = normalToFeature (aDirToFeature.second);
 
-            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;
+            theIsOut = BVH_DOT3 (aNormal, aDirToFeature.first) > T(0.0);
           }
         }
+      }
 
-        if (aHead < 0)
-          return aMinDistance;
-
-        std::pair<Standard_Integer, T>& anInfo = aStack[aHead--];
+      if (aHead < 0)
+        return aMinDist;
 
-        while (anInfo.second > aMinDistance)
-        {
-          if (aHead < 0)
-            return aMinDistance;
+      std::pair<Standard_Integer, T>& anInfo = aStack[aHead--];
 
-          anInfo = aStack[aHead--];
-        }
+      while (anInfo.second > aMinDist)
+      {
+        if (aHead < 0)
+          return aMinDist;
 
-        aNode = anInfo.first;
+        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);
+#ifdef HAVE_TBB
 
-    const BVH_Tree<T, N, BVH_BinaryTree>* aBVH = theGeometry.BVH().get();
+//! Tool object for parallel construction of distance field (uses Intel TBB).
+template<class T, int N>
+class BVH_ParallelDistanceFieldBuilder
+{
+private:
 
-    if (aBVH == NULL)
-    {
-      return Standard_False;
-    }
+  //! Output distance field.
+  BVH_DistanceField<T, N>* myOutField;
 
-    std::pair<Standard_Integer, T> aStack[32];
+public:
 
-    Standard_Integer aHead = -1;
-    Standard_Integer aNode =  0; // root node
+  BVH_ParallelDistanceFieldBuilder (BVH_DistanceField<T, N>* theOutField)
+  : myOutField (theOutField)
+  {
+    //
+  }
 
-    T aMinDistance = std::numeric_limits<T>::max();
+  void operator() (const tbb::blocked_range<Standard_Integer>& theRange) const
+  {
+    myOutField->buildSlices (theRange.begin(), theRange.end());
+  }
+};
 
-    for (;;)
-    {
-      BVH_Vec4i aData = aBVH->NodeInfoBuffer()[aNode];
+#endif
 
-      if (aData.x() == 0) // if inner node
-      {
-        const T aDistToLft = DistanceToBox<T, N> (thePnt,
-                                                  aBVH->MinPoint (aData.y()),
-                                                  aBVH->MaxPoint (aData.y()));
+//=======================================================================
+//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);
 
-        const T aDistToRgh = DistanceToBox<T, N> (thePnt,
-                                                  aBVH->MinPoint (aData.z()),
-                                                  aBVH->MaxPoint (aData.z()));
+  Standard_ASSERT_RAISE (aFeatureIt != myFeatureNormalMap.end(), "Failed to retrieve feature normal");
 
-        const Standard_Boolean aHitLft = aDistToLft <= aMinDistance;
-        const Standard_Boolean aHitRgh = aDistToRgh <= aMinDistance;
+  return aFeatureIt->second;
+}
 
-        if (aHitLft & aHitRgh)
-        {
-          aNode = (aDistToLft < aDistToRgh) ? aData.y() : aData.z();
+// =======================================================================
+// 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;
 
-          aStack[++aHead] = std::pair<Standard_Integer, T> (
-            aDistToLft < aDistToRgh ? aData.z() : aData.y(), Max (aDistToLft, aDistToRgh));
-        }
-        else
-        {
-          if (aHitLft | aHitRgh)
-          {
-            aNode = aHitLft ? aData.y() : aData.z();
-          }
-          else
-          {
-            if (aHead < 0)
-              return aMinDistance;
+  if (myMesh->Size() == 0)
+  {
+    return Standard_False;
+  }
 
-            std::pair<Standard_Integer, T>& anInfo = aStack[aHead--];
+  //----------------------------------------------------------------------
+  // Compute pseudo-normals
+  //----------------------------------------------------------------------
 
-            while (anInfo.second > aMinDistance)
-            {
-              if (aHead < 0)
-                return aMinDistance;
+  myFeatureNormalMap.clear();
 
-              anInfo = aStack[aHead--];
-            }
+  if (myComputeSign)
+  {
+    std::map<Feature, std::set<Standard_Integer>, SymmetricalPairComparator> aConnectivityMap;
 
-            aNode = anInfo.first;
-          }
-        }
-      }
-      else  // if leaf node
-      {
-        Standard_Boolean isOutside = Standard_True;
+    for (size_t aTrgIdx = 0; aTrgIdx < myMesh->Elements.size(); ++aTrgIdx) // build connectivity map
+    {
+      const BVH_Vec4i aTrg = myMesh->Elements[aTrgIdx];
 
-        const T aDistance = SquareDistanceToObject (
-          theGeometry.Objects()(aNode).operator->(), thePnt, isOutside);
+      const BVH_VecNt aTrgVrt0 = myMesh->Vertices[aTrg.x()];
+      const BVH_VecNt aTrgVrt1 = myMesh->Vertices[aTrg.y()];
+      const BVH_VecNt aTrgVrt2 = myMesh->Vertices[aTrg.z()];
 
-        if (aDistance < aMinDistance)
-        {
-          aMinDistance = aDistance;
-          theIsOutside = isOutside;
-        }
+      const BVH_VecNt aTrgEdges[] = { aTrgVrt1 - aTrgVrt0,
+                                      aTrgVrt2 - aTrgVrt0 };
 
-        if (aHead < 0)
-          return aMinDistance;
+      BVH_VecNt aNormal;
 
-        std::pair<Standard_Integer, T>& anInfo = aStack[aHead--];
+      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();
 
-        while (anInfo.second > aMinDistance)
-        {
-          if (aHead < 0)
-            return aMinDistance;
+      myFeatureNormalMap[std::make_pair (FACE_ID, aTrgIdx)] = aNormal / std::sqrt (BVH_DOT3 (aNormal, aNormal));
 
-          anInfo = aStack[aHead--];
-        }
+      for (Standard_Integer aVrtStart = 0, aVrtEnd = 1; aVrtStart < 3; ++aVrtStart, aVrtEnd = (aVrtEnd + 1) % 3)
+      {
+        aConnectivityMap[Feature (aTrg[aVrtStart], aTrg[aVrtEnd])].insert (aTrgIdx);
+      }
 
-        aNode = anInfo.first;
+      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;
+    }
   }
-}
 
-#undef BVH_DOT3
+  //----------------------------------------------------------------------
+  // Adjust parameters of voxel set
+  //----------------------------------------------------------------------
 
-#ifdef HAVE_TBB
+  const BVH_VecNt aBoxSize = myMesh->Box().Size();
 
-//! Tool object for parallel construction of distance field (uses Intel TBB).
-template<class T, int N>
-class BVH_ParallelDistanceFieldBuilder
-{
-private:
+  const T aMaxBoxSide = Max (Max (aBoxSize.x(), aBoxSize.y()), aBoxSize.z());
 
-  //! Input BVH geometry.
-  BVH_Geometry<T, N>* myGeometry;
+  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);
 
-  //! Output distance field.
-  BVH_DistanceField<T, N>* myOutField;
+  myDimensionX = Min (myMaximumSize, Max (myDimensionX, 16));
+  myDimensionY = Min (myMaximumSize, Max (myDimensionY, 16));
+  myDimensionZ = Min (myMaximumSize, Max (myDimensionZ, 16));
 
-public:
+  const BVH_VecNt aGlobalBoxMin = myMesh->Box().CornerMin();
+  const BVH_VecNt aGlobalBoxMax = myMesh->Box().CornerMax();
 
-  BVH_ParallelDistanceFieldBuilder (BVH_DistanceField<T, N>* theOutField, BVH_Geometry<T, N>* theGeometry)
-  : myGeometry (theGeometry),
-    myOutField (theOutField)
-  {
-    //
-  }
+  const Standard_Integer aVoxelOffset = 2;
 
-  void operator() (const tbb::blocked_range<Standard_Integer>& theRange) const
-  {
-    myOutField->BuildSlices (*myGeometry, theRange.begin(), theRange.end());
-  }
-};
+  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
+// 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 aGlobalBoxSize = theGeometry.Box().Size();
-
-  const T aMaxBoxSide = Max (Max (aGlobalBoxSize.x(), aGlobalBoxSize.y()), aGlobalBoxSize.z());
-
-  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);
-
-  myDimensionX = Min (myMaximumSize, Max (myDimensionX, 16));
-  myDimensionY = Min (myMaximumSize, Max (myDimensionY, 16));
-  myDimensionZ = Min (myMaximumSize, Max (myDimensionZ, 16));
+  const BVH_VecNt aLocalPnt = thePoint.cwiseMin (myCornerMax).cwiseMax (myCornerMin) - myCornerMin;
 
-  const BVH_VecNt aGlobalBoxMin = theGeometry.Box().CornerMin();
-  const BVH_VecNt aGlobalBoxMax = theGeometry.Box().CornerMax();
+  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 Standard_Integer aVoxelOffset = 2;
+  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));
 
-  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);
+  const T aRatioX = aGridPntX - aCellX;
+  const T aRatioY = aGridPntY - aCellY;
+  const T aRatioZ = aGridPntZ - aCellZ;
 
-  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;
+  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;
 
-#ifdef HAVE_TBB
-  
-  tbb::parallel_for (tbb::blocked_range<Standard_Integer> (0, myDimensionZ),
-    BVH_ParallelDistanceFieldBuilder<T, N> (this, &theGeometry));
+  const T aInterpValueMinZ = aInterpValueMinYMinZ * (T(1.0) - aRatioY) + aInterpValueMaxYMinZ * aRatioY;
+  const T aInterpValueMaxZ = aInterpValueMinYMaxZ * (T(1.0) - aRatioY) + aInterpValueMaxYMaxZ * aRatioY;
 
-#else
-
-  BuildSlices (theGeometry, 0, myDimensionZ);
-
-#endif
-
-  return Standard_True;
+  return aInterpValueMinZ * (T(1.0) - aRatioZ) + aInterpValueMaxZ * aRatioZ;
 }

diff --git a/src/BVH/BVH_ImplicitFunction.hxx b/src/BVH/BVH_ImplicitFunction.hxx
new file mode 100644 (file)
index 0000000..fd1ed9d
--- /dev/null
+++ 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

diff --git a/src/BVH/BVH_ImplicitFunction.lxx b/src/BVH/BVH_ImplicitFunction.lxx
new file mode 100644 (file)
index 0000000..4abe9da
--- /dev/null
+++ 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;
+}

diff --git a/src/BVH/BVH_MarchingCubes.cxx b/src/BVH/BVH_MarchingCubes.cxx
new file mode 100644 (file)
index 0000000..4a4478e
--- /dev/null
+++ 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;
+  }
+}
\ No newline at end of file

diff --git a/src/BVH/BVH_MarchingCubes.hxx b/src/BVH/BVH_MarchingCubes.hxx
new file mode 100644 (file)
index 0000000..466335a
--- /dev/null
+++ 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

diff --git a/src/BVH/BVH_MarchingCubes.lxx b/src/BVH/BVH_MarchingCubes.lxx
new file mode 100644 (file)
index 0000000..a0934fa
--- /dev/null
+++ 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;
+}
\ No newline at end of file

diff --git a/src/BVH/BVH_Tessellator.hxx b/src/BVH/BVH_Tessellator.hxx
new file mode 100644 (file)
index 0000000..9f1c05d
--- /dev/null
+++ 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

diff --git a/src/BVH/BVH_Tessellator.lxx b/src/BVH/BVH_Tessellator.lxx
new file mode 100644 (file)
index 0000000..e87ce4d
--- /dev/null
+++ 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;
+}
\ No newline at end of file

diff --git a/src/BVH/BVH_Triangulation.hxx b/src/BVH/BVH_Triangulation.hxx
index ea6bce2b01458195bdaab67ab2fbb6316d4ba3d6..b049d0cca399d638f25f3a47df28ae5a3a32ffac 100644 (file)
--- 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>

diff --git a/src/BVH/BVH_Triangulation.lxx b/src/BVH/BVH_Triangulation.lxx
index 19ff30f49c87c93ad9b3a6a205d9e14363ba4a63..7735070da6a678625b1b299ded4788de0e25656e 100644 (file)
--- 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();
+    }
+  }
+}

diff --git a/src/BVH/BVH_Types.hxx b/src/BVH/BVH_Types.hxx
index a556672efdc574ea31bdea04a746f873097684cc..64446a60e05ef0ccf16b010bf86707b440eb4400 100644 (file)
--- 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

diff --git a/src/BVH/FILES b/src/BVH/FILES
index 88ff44dca158a7ad729ad94ff80141cb4e9ceefa..7cd0b09c763f45cc1c4c4879758b873328027b30 100644 (file)
--- 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
\ No newline at end of file