marching-cubes

marching_cubes

skimage.measure.marching_cubes(volume, level, spacing=(1.0, 1.0, 1.0), gradient_direction='descent') [source]

Marching cubes algorithm to find iso-valued surfaces in 3d volumetric data

Parameters:

volume : (M, N, P) array of doubles Input data volume to find isosurfaces. Will be cast to np.float64. level : float Contour value to search for isosurfaces in volume. spacing : length-3 tuple of floats Voxel spacing in spatial dimensions corresponding to numpy array indexing dimensions (M, N, P) as in volume. gradient_direction : string Controls if the mesh was generated from an isosurface with gradient descent toward objects of interest (the default), or the opposite. The two options are: * descent : Object was greater than exterior * ascent : Exterior was greater than object

Returns:

verts : (V, 3) array Spatial coordinates for V unique mesh vertices. Coordinate order matches input volume (M, N, P). faces : (F, 3) array Define triangular faces via referencing vertex indices from verts. This algorithm specifically outputs triangles, so each face has exactly three indices.

See also

skimage.measure.correct_mesh_orientation, skimage.measure.mesh_surface_area

Notes

The marching cubes algorithm is implemented as described in [R264]. A simple explanation is available here:

http://www.essi.fr/~lingrand/MarchingCubes/algo.html

There are several known ambiguous cases in the marching cubes algorithm. Using point labeling as in [R264], Figure 4, as shown:

    v8 ------ v7
   / |       / |        y
  /  |      /  |        ^  z
v4 ------ v3   |        | /
 |  v5 ----|- v6        |/          (note: NOT right handed!)
 |  /      |  /          ----> x
 | /       | /
v1 ------ v2

Most notably, if v4, v8, v2, and v6 are all >= level (or any generalization of this case) two parallel planes are generated by this algorithm, separating v4 and v8 from v2 and v6. An equally valid interpretation would be a single connected thin surface enclosing all four points. This is the best known ambiguity, though there are others.

This algorithm does not attempt to resolve such ambiguities; it is a naive implementation of marching cubes as in [R264], but may be a good beginning for work with more recent techniques (Dual Marching Cubes, Extended Marching Cubes, Cubic Marching Squares, etc.).

Because of interactions between neighboring cubes, the isosurface(s) generated by this algorithm are NOT guaranteed to be closed, particularly for complicated contours. Furthermore, this algorithm does not guarantee a single contour will be returned. Indeed, ALL isosurfaces which cross level will be found, regardless of connectivity.

The output is a triangular mesh consisting of a set of unique vertices and connecting triangles. The order of these vertices and triangles in the output list is determined by the position of the smallest x,y,z (in lexicographical order) coordinate in the contour. This is a side-effect of how the input array is traversed, but can be relied upon.

The generated mesh guarantees coherent orientation as of version 0.12.

To quantify the area of an isosurface generated by this algorithm, pass outputs directly into skimage.measure.mesh_surface_area.

Regarding visualization of algorithm output, the mayavi package is recommended. To contour a volume named myvolume about the level 0.0:

>>> from mayavi import mlab 
>>> verts, faces = marching_cubes(myvolume, 0.0, (1., 1., 2.)) 
>>> mlab.triangular_mesh([vert[0] for vert in verts],
...                      [vert[1] for vert in verts],
...                      [vert[2] for vert in verts],
...                      faces) 
>>> mlab.show() 

References

[R264]

(1, 2, 3, 4) Lorensen, William and Harvey E. Cline. Marching Cubes: A High Resolution 3D Surface Construction Algorithm. Computer Graphics (SIGGRAPH 87 Proceedings) 21(4) July 1987, p. 163-170).