Using Voxel Centers for Marching Cubes

README.rst

@@ -8,6 +8,8 @@ three-dimensional NumPy array or as a Python function ``f(x, y, z)``. The first
 option is much faster, but it requires more memory and becomes unfeasible for
 very large volumes.
 
+**In this branch, the marching cubes algorithm uses the voxel centers to extract isosurfaces.**
+
 PyMCubes also provides a function to export the results of the marching cubes as
 COLLADA ``(.dae)`` files. This requires the
 `PyCollada <https://github.com/pycollada/pycollada>`_ library.

examples/spheres.py

@@ -12,9 +12,14 @@
 vertices1, triangles1 = mcubes.marching_cubes(u, 0)
 
 # Export the result to sphere.dae
-mcubes.export_mesh(vertices1, triangles1, "sphere1.dae", "MySphere")
-
-print("Done. Result saved in 'sphere1.dae'.")
+try:
+    import collada
+    mcubes.export_mesh(vertices1, triangles1, "sphere1.dae", "MySphere")
+    print("Done. Result saved in 'sphere1.dae'.")
+except ImportError:
+    print("Could not import collada. Saving as OFF instead.")
+    mcubes.export_off(vertices1, triangles1, "sphere1.off")
+    print("Done. Result saved in 'sphere1.off'.")
 
 print("Example 2: Isosurface in Python function...")
 print("(this might take a while...)")
@@ -31,8 +36,14 @@ def f(x, y, z):
         16)                         # Isosurface value
 
 # Export the result to sphere2.dae
-mcubes.export_mesh(vertices2, triangles2, "sphere2.dae", "MySphere")
-print("Done. Result saved in 'sphere2.dae'.")
+try:
+    import collada
+    mcubes.export_mesh(vertices2, triangles2, "sphere2.dae", "MySphere")
+    print("Done. Result saved in 'sphere2.dae'.")
+except ImportError:
+    print("Could not import collada. Saving as OFF instead.")
+    mcubes.export_off(vertices2, triangles2, "sphere2.off")
+    print("Done. Result saved in 'sphere2.off'.")
 
 try:
     print("Plotting mesh...")

mcubes/__init__.py

@@ -1,3 +1,3 @@
 
 from ._mcubes import marching_cubes, marching_cubes_func
-from .exporter import export_mesh, export_obj
+from .exporter import export_mesh, export_obj, export_off

mcubes/exporter.py

@@ -15,6 +15,23 @@ def export_obj(vertices, triangles, filename):
         for f in triangles:
             fh.write("f {} {} {}\n".format(*(f + 1)))
 
+
+def export_off(vertices, triangles, filename):
+    """
+    Exports a mesh in the (.off) format.
+    """
+
+    with open(filename, 'w') as fh:
+        fh.write('OFF\n')
+        fh.write('{} {} 0\n'.format(len(vertices), len(triangles)))
+
+        for v in vertices:
+            fh.write("{} {} {}\n".format(*v))
+
+        for f in triangles:
+            fh.write("3 {} {} {}\n".format(*f))
+
+
 def export_mesh(vertices, triangles, filename, mesh_name="mcubes_mesh"):
     """
     Exports a mesh in the COLLADA (.dae) format.

mcubes/src/marchingcubes.h

@@ -28,47 +28,47 @@ void marching_cubes(const vector3& lower, const vector3& upper,
     using namespace private_;
 
     // typedef decltype(lower[0]) coord_type;
-    
+
     // numx, numy and numz are the numbers of evaluations in each direction
     --numx; --numy; --numz;
-    
+
     coord_type dx = (upper[0] - lower[0])/static_cast<coord_type>(numx);
     coord_type dy = (upper[1] - lower[1])/static_cast<coord_type>(numy);
     coord_type dz = (upper[2] - lower[2])/static_cast<coord_type>(numz);
-    
+
     size_t* shared_indices = new size_t[2*numy*numz*3];
     const int z3 = numz*3;
     const int yz3 = numy*z3;
-    
+
     for(int i=0; i<numx; ++i)
     {
-        coord_type x = lower[0] + dx*i;
-        coord_type x_dx = lower[0] + dx*(i+1);
+        coord_type x = lower[0] + dx*i + dx/2;
+        coord_type x_dx = lower[0] + dx*(i+1) + dx/2;
         const int i_mod_2 = i % 2;
-        const int i_mod_2_inv = (i_mod_2 ? 0 : 1); 
+        const int i_mod_2_inv = (i_mod_2 ? 0 : 1);
 
         for(int j=0; j<numy; ++j)
         {
-            coord_type y = lower[1] + dy*j;
-            coord_type y_dy = lower[1] + dy*(j+1);
+            coord_type y = lower[1] + dy*j + dy/2;
+            coord_type y_dy = lower[1] + dy*(j+1) + dy/2;
             for(int k=0; k<numz; ++k)
             {
-                coord_type z = lower[2] + dz*k;
-                coord_type z_dz = lower[2] + dz*(k+1);
-                
+                coord_type z = lower[2] + dz*k + dz/2;
+                coord_type z_dz = lower[2] + dz*(k+1) + dz/2;
+
                 double v[8];
                 v[0] = f(x,y,z); v[1] = f(x_dx,y,z);
                 v[2] = f(x_dx,y_dy,z); v[3] = f(x, y_dy, z);
                 v[4] = f(x,y,z_dz); v[5] = f(x_dx,y,z_dz);
                 v[6] = f(x_dx,y_dy,z_dz); v[7] = f(x, y_dy, z_dz);
-                
+
                 unsigned int cubeindex = 0;
                 for(int m=0; m<8; ++m)
                     if(v[m] <= isovalue)
                         cubeindex |= 1<<m;
-                
+
                 // Generate vertices AVOIDING DUPLICATES.
-                
+
                 int edges = edge_table[cubeindex];
                 std::vector<size_t> indices(12, -1);
                 if(edges & 0x040)
@@ -86,10 +86,10 @@ void marching_cubes(const vector3& lower, const vector3& upper,
                 if(edges & 0x400)
                 {
                     indices[10] = vertices.size() / 3;
-                    shared_indices[i_mod_2*yz3 + j*z3 + k*3 + 2] = indices[10];                    
+                    shared_indices[i_mod_2*yz3 + j*z3 + k*3 + 2] = indices[10];
                     mc_add_vertex(x_dx, y+dx, z, z_dz, 2, v[2], v[6], isovalue, &vertices);
                 }
-                
+
                 if(edges & 0x001)
                 {
                     if(j == 0 || k == 0)
@@ -180,15 +180,15 @@ void marching_cubes(const vector3& lower, const vector3& upper,
                     else
                         indices[11] = shared_indices[i_mod_2_inv*yz3 + j*z3 + k*3 + 2];
                 }
-                
+
                 int tri;
                 int* triangle_table_ptr = triangle_table[cubeindex];
                 for(int m=0; tri = triangle_table_ptr[m], tri != -1; ++m)
                     polygons.push_back(indices[tri]);
             }
         }
     }
-    
+
     delete [] shared_indices;
 }