FIX: volume now directly cast, docfixes, C-contiguous ordering

skimage/measure/_marching_cubes.py

@@ -19,12 +19,12 @@ def marching_cubes(volume, level, sampling=(1., 1., 1.)):
 
     Returns
     -------
-    vert_list : list
-        Every entry in this list is a unique vertex on the isosurface.
-    tri_list : list
-        Every entry in this list is a length-3 list of integers. These
-        represent triangular faces; the integers in each sub-list correspond
-        to vertices held in `vert_list`.
+    verts : (V, 3) array
+        Spatial (x, y, z) coordinates for V unique vertices in mesh.
+    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.
 
     Notes
     -----
@@ -69,7 +69,7 @@ def marching_cubes(volume, level, sampling=(1., 1., 1.)):
     of how the input array is traversed, but can be relied upon.
 
     To quantify the area of an isosurface generated by this algorithm, pass
-    the output directly into `skimage.measure.mesh_surface_area`.
+    the 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:
@@ -93,17 +93,12 @@ def marching_cubes(volume, level, sampling=(1., 1., 1.)):
     skimage.measure.mesh_surface_area
 
     """
-    # Check inputs
+    # Check inputs and ensure `volume` is C-contiguous for memoryviews
     if volume.ndim != 3:
-        raise ValueError("Input volume must be 3d.")
-    if volume.dtype.kind == 'f':
-        volume = volume.astype(np.float)
-    else:
-        from skimage.util import img_as_float
-        # If incorrect type provided, convert BOTH contour value and input
-        # volume using same method
-        level = img_as_float(np.array(level, dtype=volume.dtype))[0]
-        volume = img_as_float(volume)
+        raise ValueError("Input volume must have 3 dimensions.")
+    if level < volume.min() or level > volume.max():
+        raise ValueError("Contour level must be within volume data range.")
+    volume = np.array(volume, dtype=np.float64, order="C")
 
     # Extract raw triangles using marching cubes in Cython
     #   Returns a list of length-3 lists, each sub-list containing three
@@ -115,10 +110,10 @@ def marching_cubes(volume, level, sampling=(1., 1., 1.)):
                                                        sampling)
 
     # Find and collect unique vertices, storing triangle verts as indices.
-    # Removes much redundancy and eliminates degenerate "triangles".
-    vert_list, tri_list = _marching_cubes_cy.unpack_unique_verts(raw_tris)
+    # Returns a true mesh with no degenerate faces.
+    verts, faces = _marching_cubes_cy.unpack_unique_verts(raw_tris)
 
-    return vert_list, tri_list
+    return np.asarray(verts), np.asarray(faces)
 
 
 def mesh_surface_area(verts, tris):
@@ -127,17 +122,16 @@ def mesh_surface_area(verts, tris):
 
     Parameters
     ----------
-    verts : list
-        List of length-3 NumPy arrays containing vertex coordinates.
-        Units in each dimension should be consistent.
-    tris : list
+    verts : (V, 3) array of floats
+        Array containing (x, y, z) coordinates for V unique mesh vertices.
+    faces : (F, 3) array of ints
         List of length-3 lists of integers, referencing vertex coordinates as
         provided in `verts`
 
     Returns
     -------
     area : float
-        Surface area of mesh. Units in coordinates maintained, but squared.
+        Surface area of mesh. Units now [coordinate units] ** 2.
 
     Notes
     -----
@@ -145,13 +139,16 @@ def mesh_surface_area(verts, tris):
     `skimage.measure.marching_cubes`. For unit correct output, ensure correct
     `spacing` was passed to `skimage.measure.marching_cubes`.
 
+    This algorithm works properly only if the ``faces`` provided are all
+    triangles.
+
     See Also
     --------
     skimage.measure.marching_cubes
 
     """
-    # Define two vector arrays `a` and `b` from triangle vertices
-    actual_verts = np.array([[verts[i] for i in tri] for tri in tris])
+    # Fancy indexing to define two vector arrays from triangle vertices
+    actual_verts = verts[tris]
     a = actual_verts[:, 0, :] - actual_verts[:, 1, :]
     b = actual_verts[:, 0, :] - actual_verts[:, 2, :]
     del actual_verts