From 5eba5a42a7e3f0c4d8ff3e50c9dff98a90d7b950 Mon Sep 17 00:00:00 2001
From: "Josh Warner (Mac)" <warner.joshua@mayo.edu>
Date: Sat, 16 Mar 2013 21:12:13 -0500
Subject: [PATCH 01/24] FEAT: Full 3d anisotropic implementation of marching
 cubes

Also includes surface area calculation algorithm from
generated mesh.
Convenient output to visualize with `mayavi.mlab`.
Efficient Cython implementation.
---
 skimage/measure/_marching_cubes.pyx | 924 ++++++++++++++++++++++++++++
 skimage/measure/marching_cubes.py   | 200 ++++++
 2 files changed, 1124 insertions(+)
 create mode 100644 skimage/measure/_marching_cubes.pyx
 create mode 100644 skimage/measure/marching_cubes.py

diff --git a/skimage/measure/_marching_cubes.pyx b/skimage/measure/_marching_cubes.pyx
new file mode 100644
index 00000000..8858acf4
--- /dev/null
+++ b/skimage/measure/_marching_cubes.pyx
@@ -0,0 +1,924 @@
+#cython: cdivision=True
+#cython: boundscheck=False
+#cython: nonecheck=False
+#cython: wraparound=False
+import numpy as np
+cimport numpy as cnp
+
+
+cdef inline double _get_fraction(double from_value, double to_value,
+                                 double level):
+    if (to_value == from_value):
+        return 0
+    return ((level - from_value) / (to_value - from_value))
+
+
+def iterate_and_store_3d(cnp.ndarray[double, ndim=3] arr,
+                         double level, tuple sampling=(1., 1., 1.)):
+    """Iterate across the given array in a marching-cubes fashion,
+    looking for volumes with edges that cross 'level'. If such a volume is
+    found, appropriate triangulations are added to a growing list of
+    triangles to be returned by the function.
+
+    If `sampling` is not provided, vertices are returned in the indexing
+    coordinate system (assuming all 3 spatial dimensions sampled equally).
+    If `sampling` is provided, vertices will be returned in volume coordinates
+    relative to the origin, regularly spaced as specified in each dimension.
+
+    """
+    if arr.shape[0] < 2 or arr.shape[1] < 2 or arr.shape[2] < 2:
+        raise ValueError("Input array must be at least 2x2x2.")
+    if len(sampling) != 3:
+        raise ValueError("`sampling` must be of form (double, double, double)")
+
+    cdef list tri_list = []
+    cdef list norm_list = []
+    cdef Py_ssize_t n
+    cdef bint odd_sampling, plus_z
+    plus_z = False
+    if ((sampling == (1., 1., 1.)) or
+         (sampling == (1., 1., 1)) or
+         (sampling == (1., 1, 1.)) or
+         (sampling == (1, 1., 1.)) or
+         (sampling == (1, 1, 1.)) or
+         (sampling == (1., 1, 1)) or
+         (sampling == (1, 1., 1)) or
+         (sampling == (1, 1, 1))):
+        odd_sampling = False
+    else:
+        odd_sampling = True
+
+    # The plan is to iterate a 2x2x2 cube across the input array. This means
+    # the upper-left corner of the cube needs to iterate across a sub-array
+    # of size one-less-large in each direction (so we can get away with no
+    # bounds checking in Cython). The cube is represented by eight vertices:
+    # v1, v2, ..., v8, oriented thus (see Lorensen, Figure 4):
+    #
+    #           v8 ------ v7
+    #          / |       / |        y
+    #         /  |      /  |        ^  z
+    #       v4 ------ v3   |        | /
+    #        |  v5 ----|- v6        |/          (note: NOT right handed!)
+    #        | /       |  /          ----> x
+    #        |/        | /
+    #       v1 ------ v2
+    #
+    # We also maintain the current 2D coordinates for v1, and ensure the array
+    # is of type 'double' and is C-contiguous (last index varies fastest).
+
+    # Coords start at (0, 0, 0).
+    cdef Py_ssize_t[3] coords
+    coords[0] = 0
+    coords[1] = 0
+    coords[2] = 0
+
+    # Extract doubles from `sampling` for speed
+    cdef double[3] sampling2
+    sampling2[0] = sampling[0]
+    sampling2[1] = sampling[1]
+    sampling2[2] = sampling[2]
+
+    # Calculate the number of iterations we'll need
+    cdef Py_ssize_t num_cube_steps = ((arr.shape[0] - 1) *
+                                      (arr.shape[1] - 1) *
+                                      (arr.shape[2] - 1))
+
+    cdef unsigned char cube_case = 0
+    cdef tuple e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, e11, e12
+    cdef double v1, v2, v3, v4, v5, v6, v7, v8, r0, r1, c0, c1, d0, d1
+    cdef Py_ssize_t x0, y0, z0, x1, y1, z1
+    e5, e6, e7, e8 = (0, 0, 0), (0, 0, 0), (0, 0, 0), (0, 0, 0)
+
+    for n in range(num_cube_steps):
+        # There are 255 unique values for `cube_case`. This algorithm follows
+        # the Lorensen paper in vertex and edge labeling, however, it should
+        # be noted that Lorensen used a left-handed coordinate system while
+        # NumPy uses a proper right handed system. Transforming between these
+        # coordinate systems was handled in the definitions of the cube
+        # vertices v1, v2, ..., v8.
+        #
+        # Refer to the paper, figure 4, for cube edge designations e1, ... e12
+
+        # Standard Py_ssize_t coordinates for indexing
+        x0, y0, z0 = coords[0], coords[1], coords[2]
+        x1, y1, z1 = x0 + 1, y0 + 1, z0 + 1
+
+        if odd_sampling:
+            # These doubles are the modified world coordinates; they are only
+            # calculated if non-default `sampling` provided.
+            r0 = coords[0] * sampling2[0]
+            c0 = coords[1] * sampling2[1]
+            d0 = coords[2] * sampling2[2]
+            r1 = r0 + sampling2[0]
+            c1 = c0 + sampling2[1]
+            d1 = d0 + sampling2[2]
+        else:
+            r0, c0, d0, r1, c1, d1 = x0, y0, z0, x1, y1, z1
+
+        # We use a right-handed coordinate system, UNlike the paper, but want
+        # to index in agreement - the coordinate adjustment takes place here.
+        v1 = arr[x0, y0, z0]
+        v2 = arr[x1, y0, z0]
+        v3 = arr[x1, y1, z0]
+        v4 = arr[x0, y1, z0]
+        v5 = arr[x0, y0, z1]
+        v6 = arr[x1, y0, z1]
+        v7 = arr[x1, y1, z1]
+        v8 = arr[x0, y1, z1]
+
+        # Unique triangulation cases
+        cube_case = 0
+        if (v1 > level): cube_case += 1
+        if (v2 > level): cube_case += 2
+        if (v3 > level): cube_case += 4
+        if (v4 > level): cube_case += 8
+        if (v5 > level): cube_case += 16
+        if (v6 > level): cube_case += 32
+        if (v7 > level): cube_case += 64
+        if (v8 > level): cube_case += 128
+
+        if (cube_case != 0 and cube_case != 255):
+            # Only do anything if there's a plane intersecting the cube.
+            # Cases 0 and 255 are entirely below/above the contour.
+
+            if cube_case > 127:
+                cube_case = 255 - cube_case
+
+            # Calculate cube edges, to become triangulation vertices.
+            # If we moved in a convenient direction, save 1/3 of the effort by
+            # re-assigning prior results.
+            if plus_z:
+                # Reassign prior calculated edges
+                e1 = e5
+                e2 = e6
+                e3 = e7
+                e4 = e8
+            else:
+                # Calculate edges normally
+                if odd_sampling:
+                    e1 = r0 + _get_fraction(v1, v2, level) * sampling2[0], c0, d0
+                    e2 = r1, c0 + _get_fraction(v2, v3, level) * sampling2[1], d0
+                    e3 = r0 + _get_fraction(v4, v3, level) * sampling2[0], c1, d0
+                    e4 = r0, c0 + _get_fraction(v1, v4, level) * sampling2[1], d0
+                else:
+                    e1 = r0 + _get_fraction(v1, v2, level), c0, d0
+                    e2 = r1, c0 + _get_fraction(v2, v3, level), d0
+                    e3 = r0 + _get_fraction(v4, v3, level), c1, d0
+                    e4 = r0, c0 + _get_fraction(v1, v4, level), d0
+
+            # These must be calculated at each point uunless we implemented a
+            # large, growing lookup table for all adjacent values; could save
+            # ~30% in terms of runtime at the expense of memory usage and
+            # much greater complexity.
+            if odd_sampling:
+                e5 = r0 + _get_fraction(v5, v6, level) * sampling2[0], c0, d1
+                e6 = r1, c0 + _get_fraction(v6, v7, level) * sampling2[1], d1
+                e7 = r0 + _get_fraction(v8, v7, level) * sampling2[0], c1, d1
+                e8 = r0, c0 + _get_fraction(v5, v8, level) * sampling2[1], d1
+                e9 = r0, c0, d0 + _get_fraction(v1, v5, level) * sampling2[2]
+                e10 = r1, c0, d0 + _get_fraction(v2, v6, level) * sampling2[2]
+                e11 = r0, c1, d0 + _get_fraction(v4, v8, level) * sampling2[2]
+                e12 = r1, c1, d0 + _get_fraction(v3, v7, level) * sampling2[2]
+            else:
+                e5 = r0 + _get_fraction(v5, v6, level), c0, d1
+                e6 = r1, c0 + _get_fraction(v6, v7, level), d1
+                e7 = r0 + _get_fraction(v8, v7, level), c1, d1
+                e8 = r0, c0 + _get_fraction(v5, v8, level), d1
+                e9 = r0, c0, d0 + _get_fraction(v1, v5, level)
+                e10 = r1, c0, d0 + _get_fraction(v2, v6, level)
+                e11 = r0, c1, d0 + _get_fraction(v4, v8, level)
+                e12 = r1, c1, d0 + _get_fraction(v3, v7, level)
+
+
+            # Append appropriate triangles to the growing output `tri_list`
+            _append_tris(tri_list, cube_case, e1, e2, e3, e4, e5,
+                           e6, e7, e8, e9, e10, e11, e12)
+
+        # Advance the coords indices
+        if coords[2] < arr.shape[2] - 2:
+            coords[2] += 1
+            plus_z = True
+        elif coords[1] < arr.shape[1] - 2:
+            coords[1] += 1
+            coords[2] = 0
+            plus_z = False
+        else:
+            coords[0] += 1
+            coords[1] = 0
+            coords[2] = 0
+            plus_z = False
+
+    return tri_list
+
+
+def _append_tris(list tri_list, unsigned char case, tuple e1, tuple e2,
+                 tuple e3, tuple e4, tuple e5, tuple e6, tuple e7, tuple e8,
+                 tuple e9, tuple e10, tuple e11, tuple e12):
+    # Permits recursive use for duplicated planes to conserve code - it's
+    # quite long enough as-is.
+
+    if (case == 1):
+        # front lower left corner
+        tri_list.append([e1, e4, e9])
+    elif (case == 2):
+        # front lower right corner
+        tri_list.append([e10, e2, e1])
+    elif (case == 3):
+        # front lower plane
+        tri_list.append([e2, e4, e9])
+        tri_list.append([e2, e9, e10])
+    elif (case == 4):
+        # front upper right corner
+        tri_list.append([e12, e3, e2])
+    elif (case == 5):
+        # lower left, upper right corners
+        _append_tris(tri_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 4, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 6):
+        # front right plane
+        tri_list.append([e12, e3, e1])
+        tri_list.append([e12, e1, e10])
+    elif (case == 7):
+        # Shelf including v1, v2, v3
+        tri_list.append([e3, e4, e12])
+        tri_list.append([e4, e9, e12])
+        tri_list.append([e12, e9, e10])
+    elif (case == 8):
+        # front upper left corner
+        tri_list.append([e3, e11, e4])
+    elif (case == 9):
+        # front left plane
+        tri_list.append([e3, e11, e9])
+        tri_list.append([e3, e9, e1])
+    elif (case == 10):
+        # upper left, lower right corners
+        _append_tris(tri_list, 2, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 11):
+        # Shelf including v4, v1, v2
+        tri_list.append([e3, e11, e2])
+        tri_list.append([e11, e10, e2])
+        tri_list.append([e11, e9, e10])
+    elif (case == 12):
+        # front upper plane
+        tri_list.append([e11, e4, e12])
+        tri_list.append([e2, e4, e12])
+    elif (case == 13):
+        # Shelf including v1, v4, v3
+        tri_list.append([e11, e9, e12])
+        tri_list.append([e12, e9, e1])
+        tri_list.append([e12, e1, e2])
+    elif (case == 14):
+        # Shelf including v2, v3, v4
+        tri_list.append([e11, e10, e12])
+        tri_list.append([e11, e4, e10])
+        tri_list.append([e4, e1, e10])
+    elif (case == 15):
+        # Plane parallel to x-axis through middle
+        tri_list.append([e11, e9, e12])
+        tri_list.append([e12, e9, e10])
+    elif (case == 16):
+        # back lower left corner
+        tri_list.append([e8, e9, e5])
+    elif (case == 17):
+        # lower left plane
+        tri_list.append([e4, e1, e8])
+        tri_list.append([e8, e1, e5])
+    elif (case == 18):
+        # lower left back, lower right front corners
+        _append_tris(tri_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 2, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 19):
+        # Shelf including v1, v2, v5
+        tri_list.append([e8, e4, e2])
+        tri_list.append([e8, e2, e10])
+        tri_list.append([e8, e10, e5])
+    elif (case == 20):
+        # lower left back, upper right front corners
+        _append_tris(tri_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 4, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 21):
+        # lower left plane + upper right front corner, v1, v3, v5
+        _append_tris(tri_list, 17, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 4, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 22):
+        # front right plane + lower left back corner, v2, v3, v5
+        _append_tris(tri_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 6, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 23):
+        # Rotated case 14 in the paper
+        tri_list.append([e3, e10, e8])
+        tri_list.append([e3, e10, e12])
+        tri_list.append([e8, e10, e5])
+        tri_list.append([e3, e4, e8])
+    elif (case == 24):
+        # upper front left, lower back left corners
+        _append_tris(tri_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 25):
+        # Shelf including v1, v4, v5
+        tri_list.append([e1, e5, e3])
+        tri_list.append([e3, e8, e11])
+        tri_list.append([e3, e5, e8])
+    elif (case == 26):
+        # Three isolated corners
+        _append_tris(tri_list, 2, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 27):
+        # Full corner v1, case 9 in paper: (v1, v2, v4, v5)
+        tri_list.append([e11, e3, e2])
+        tri_list.append([e11, e2, e10])
+        tri_list.append([e10, e11, e8])
+        tri_list.append([e8, e5, e10])
+    elif (case == 28):
+        # upper front plane + corner v5
+        _append_tris(tri_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 12, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 29):
+        # special case of 11 in the paper: (v1, v3, v4, v5)
+        tri_list.append([e11, e5, e2])
+        tri_list.append([e11, e12, e2])
+        tri_list.append([e11, e5, e8])
+        tri_list.append([e2, e1, e5])
+    elif (case == 30):
+        # Shelf (v2, v3, v4) and lower left back corner
+        _append_tris(tri_list, 14, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 31):
+        # Shelf: (v6, v7, v8) by inversion
+        tri_list.append([e11, e12, e10])
+        tri_list.append([e11, e8, e10])
+        tri_list.append([e8, e10, e5])
+    elif (case == 32):
+        # lower right back corner
+        tri_list.append([e6, e5, e10])
+    elif (case == 33):
+        # lower right back, lower left front corners
+        _append_tris(tri_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 34):
+        # lower right plane
+        tri_list.append([e1, e2, e5])
+        tri_list.append([e2, e6, e5])
+    elif (case == 35):
+        # Shelf: v1, v2, v6
+        tri_list.append([e4, e2, e6])
+        tri_list.append([e4, e9, e6])
+        tri_list.append([e6, e9, e5])
+    elif (case == 36):
+        # upper right front, lower right back corners
+        _append_tris(tri_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 4, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 37):
+        # lower left front, upper right front, lower right back corners
+        _append_tris(tri_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 4, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 38):
+        # Shelf: v2, v3, v6
+        tri_list.append([e3, e1, e5])
+        tri_list.append([e3, e5, e12])
+        tri_list.append([e12, e5, e6])
+    elif (case == 39):
+        # Full corner v2: (v1, v2, v3, v6)
+        tri_list.append([e3, e4, e5])
+        tri_list.append([e4, e9, e5])
+        tri_list.append([e3, e5, e6])
+        tri_list.append([e3, e12, e6])
+    elif (case == 40):
+        # upper left front, lower right back corners
+        _append_tris(tri_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 41):
+        # front left plane, lower right back corner
+        _append_tris(tri_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 9, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 42):
+        # lower right plane, upper front left corner
+        _append_tris(tri_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 34, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 43):
+        # Rotated case 11 in paper
+        tri_list.append([e11, e3, e9])
+        tri_list.append([e3, e9, e6])
+        tri_list.append([e3, e2, e6])
+        tri_list.append([e9, e5, e6])
+    elif (case == 44):
+        # upper front plane, lower right back corner
+        _append_tris(tri_list, 12, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 45):
+        # Shelf: (v1, v3, v4) + lower right back corner
+        _append_tris(tri_list, 13, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 46):
+        # Rotated case 14 in paper
+        tri_list.append([e4, e11, e12])
+        tri_list.append([e4, e12, e5])
+        tri_list.append([e12, e5, e6])
+        tri_list.append([e4, e5, e1])
+    elif (case == 47):
+        # Shelf: (v5, v8, v7) by inversion
+        tri_list.append([e11, e9, e12])
+        tri_list.append([e12, e9, e5])
+        tri_list.append([e12, e5, e6])
+    elif (case == 48):
+        # Back lower plane
+        tri_list.append([e9, e10, e6])
+        tri_list.append([e9, e6, e8])
+    elif (case == 49):
+        # Shelf: (v1, v5, v6)
+        tri_list.append([e4, e8, e6])
+        tri_list.append([e4, e6, e1])
+        tri_list.append([e6, e1, e10])
+    elif (case == 50):
+        # Shelf: (v2, v5, v6)
+        tri_list.append([e8, e6, e2])
+        tri_list.append([e8, e2, e1])
+        tri_list.append([e8, e9, e1])
+    elif (case == 51):
+        # Plane through middle of cube, parallel to x-z axis
+        tri_list.append([e4, e8, e2])
+        tri_list.append([e8, e2, e6])
+    elif (case == 52):
+        # Back lower plane, and front upper right corner
+        _append_tris(tri_list, 48, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 4, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 53):
+        # Shelf (v1, v5, v6) and front upper right corner
+        _append_tris(tri_list, 49, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 4, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 54):
+        # Rotated case 11 from paper (v2, v3, v5, v6)
+        tri_list.append([e1, e9, e3])
+        tri_list.append([e9, e3, e6])
+        tri_list.append([e9, e8, e6])
+        tri_list.append([e12, e3, e6])
+    elif (case == 55):
+        # Shelf: (v4, v8, v7) by inversion
+        tri_list.append([e4, e8, e6])
+        tri_list.append([e4, e6, e3])
+        tri_list.append([e6, e3, e12])
+    elif (case == 56):
+        # Back lower plane + upper left front corner
+        _append_tris(tri_list, 48, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 57):
+        # Rotated case 14 from paper (v4, v1, v5, v6)
+        tri_list.append([e3, e11, e8])
+        tri_list.append([e3, e8, e10])
+        tri_list.append([e10, e6, e8])
+        tri_list.append([e3, e1, e10])
+    elif (case == 58):
+        # Shelf: (v2, v6, v5) + upper left front corner
+        _append_tris(tri_list, 50, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 59):
+        # Shelf: (v3, v7, v8) by inversion
+        tri_list.append([e2, e6, e8])
+        tri_list.append([e8, e2, e3])
+        tri_list.append([e8, e3, e11])
+    elif (case == 60):
+        # AMBIGUOUS CASE: parallel planes (front upper, back lower)
+        _append_tris(tri_list, 48, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 12, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 61):
+        # Upper back plane + lower right front corner by inversion
+        _append_tris(tri_list, 63, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 2, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 62):
+        # Upper back plane + lower left front corner by inversion
+        _append_tris(tri_list, 63, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 63):
+        # Upper back plane
+        tri_list.append([e11, e12, e6])
+        tri_list.append([e11, e8, e6])
+    elif (case == 64):
+        # Upper right back corner
+        tri_list.append([e12, e7, e6])
+    elif (case == 65):
+        # upper right back, lower left front corners
+        _append_tris(tri_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 66):
+        # upper right back, lower right front corners
+        _append_tris(tri_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 2, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 67):
+        # lower front plane + upper right back corner
+        _append_tris(tri_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 3, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 68):
+        # upper right plane
+        tri_list.append([e3, e2, e6])
+        tri_list.append([e3, e7, e6])
+    elif (case == 69):
+        # Upper right plane, lower left front corner
+        _append_tris(tri_list, 68, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 70):
+        # Shelf: (v2, v3, v7)
+        tri_list.append([e1, e3, e7])
+        tri_list.append([e1, e10, e7])
+        tri_list.append([e7, e10, e6])
+    elif (case == 71):
+        # Rotated version of case 11 in paper (v1, v2, v3, v7)
+        tri_list.append([e10, e7, e4])
+        tri_list.append([e4, e3, e7])
+        tri_list.append([e10, e4, e9])
+        tri_list.append([e7, e10, e6])
+    elif (case == 72):
+        # upper left front, upper right back corners
+        _append_tris(tri_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 73):
+        # front left plane, upper right back corner
+        _append_tris(tri_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 9, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 74):
+        # Three isolated corners, exactly case 7 in paper
+        _append_tris(tri_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 2, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 75):
+        # Shelf: (v1, v2, v4) + upper right back corner
+        _append_tris(tri_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 11, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 76):
+        # Shelf: (v4, v3, v7)
+        tri_list.append([e4, e2, e6])
+        tri_list.append([e4, e11, e7])
+        tri_list.append([e4, e7, e6])
+    elif (case == 77):
+        # Rotated case 14 in paper (v1, v4, v3, v7)
+        tri_list.append([e11, e9, e1])
+        tri_list.append([e11, e1, e6])
+        tri_list.append([e1, e6, e2])
+        tri_list.append([e11, e6, e7])
+    elif (case == 78):
+        # Full corner v3: (v2, v3, v4, v7)
+        tri_list.append([e1, e4, e7])
+        tri_list.append([e1, e7, e6])
+        tri_list.append([e4, e11, e7])
+        tri_list.append([e1, e10, e6])
+    elif (case == 79):
+        # Shelf: (v6, v5, v8) by inversion
+        tri_list.append([e9, e11, e10])
+        tri_list.append([e11, e7, e10])
+        tri_list.append([e7, e10, e6])
+    elif (case == 80):
+        # lower left back, upper right back corners (v5, v7)
+        _append_tris(tri_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 81):
+        # lower left plane, upper right back corner
+        _append_tris(tri_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 17, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 82):
+        # isolated corners (v2, v5, v7)
+        _append_tris(tri_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 2, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 83):
+        # Shelf: (v1, v2, v5) + upper right back corner
+        _append_tris(tri_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 19, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 84):
+        # upper right plane, lower left back corner
+        _append_tris(tri_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 68, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 85):
+        # AMBIGUOUS CASE: upper right and lower left parallel planes
+        _append_tris(tri_list, 17, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 68, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 86):
+        # Shelf: (v2, v3, v7) + lower left back corner
+        _append_tris(tri_list, 70, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 87):
+        # Upper left plane + lower right back corner, by inversion
+        _append_tris(tri_list, 119, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 88):
+        # Isolated corners v4, v5, v7
+        _append_tris(tri_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 89):
+        # Shelf: (v1, v4, v5) + isolated corner v7
+        _append_tris(tri_list, 25, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 90):
+        # Four isolated corners v2, v4, v5, v7
+        _append_tris(tri_list, 2, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 91):
+        # Three isolated corners, v3, v6, v8 by inversion
+        _append_tris(tri_list, 4, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 127, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 92):
+        # Shelf (v4, v3, v7) + isolated corner v5
+        _append_tris(tri_list, 76, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 93):
+        # Lower right plane + isolated corner v8 by inversion
+        _append_tris(tri_list, 127, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 34, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 94):
+        # Isolated corners v1, v6, v8 by inversion
+        _append_tris(tri_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 127, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 95):
+        # Isolated corners v6, v8 by inversion
+        _append_tris(tri_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 127, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 96):
+        # back right plane
+        tri_list.append([e7, e12, e5])
+        tri_list.append([e5, e10, e12])
+    elif (case == 97):
+        # back right plane + isolated corner v1
+        _append_tris(tri_list, 96, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 98):
+        # Shelf: (v2, v6, v7)
+        tri_list.append([e1, e7, e5])
+        tri_list.append([e7, e1, e12])
+        tri_list.append([e1, e12, e2])
+    elif (case == 99):
+        # Rotated case 14 in paper: (v1, v2, v6, v7)
+        tri_list.append([e9, e2, e7])
+        tri_list.append([e9, e2, e4])
+        tri_list.append([e2, e7, e12])
+        tri_list.append([e7, e9, e5])
+    elif (case == 100):
+        # Shelf: (v3, v6, v7)
+        tri_list.append([e3, e7, e5])
+        tri_list.append([e3, e5, e2])
+        tri_list.append([e2, e5, e10])
+    elif (case == 101):
+        # Shelf: (v3, v6, v7) + isolated corner v1
+        _append_tris(tri_list, 100, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 102):
+        # Plane bisecting left-right halves of cube
+        tri_list.append([e1, e3, e7])
+        tri_list.append([e1, e7, e5])
+    elif (case == 103):
+        # Shelf: (v4, v5, v8) by inversion
+        tri_list.append([e3, e7, e5])
+        tri_list.append([e3, e5, e4])
+        tri_list.append([e4, e5, e9])
+    elif (case == 104):
+        # Back right plane + isolated corner v4
+        _append_tris(tri_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 96, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 105):
+        # AMBIGUOUS CASE: back right and front left planes
+        _append_tris(tri_list, 96, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 9, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 106):
+        # Shelf: (v2, v6, v7) + isolated corner v4
+        _append_tris(tri_list, 98, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 107):
+        # Back left plane + isolated corner v3 by inversion
+        _append_tris(tri_list, 4, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 111, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 108):
+        # Rotated case 11 from paper: (v4, v3, v7, v6)
+        tri_list.append([e4, e10, e7])
+        tri_list.append([e4, e10, e2])
+        tri_list.append([e4, e11, e7])
+        tri_list.append([e7, e10, e5])
+    elif (case == 109):
+        # Back left plane + isolated corner v2 by inversion
+        _append_tris(tri_list, 111, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 2, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 110):
+        # Shelf: (v1, v5, v8) by inversion
+        tri_list.append([e1, e5, e7])
+        tri_list.append([e1, e7, e11])
+        tri_list.append([e1, e11, e4])
+    elif (case == 111):
+        # Back left plane
+        tri_list.append([e11, e9, e7])
+        tri_list.append([e9, e7, e5])
+    elif (case == 112):
+        # Shelf: (v5, v6, v7)
+        tri_list.append([e9, e10, e12])
+        tri_list.append([e9, e12, e7])
+        tri_list.append([e9, e7, e8])
+    elif (case == 113):
+        # Exactly case 11 from paper: (v1, v5, v6, v7)
+        tri_list.append([e1, e8, e12])
+        tri_list.append([e1, e8, e4])
+        tri_list.append([e8, e7, e12])
+        tri_list.append([e12, e1, e10])
+    elif (case == 114):
+        # Full corner v6: (v2, v6, v7, v5)
+        tri_list.append([e1, e9, e7])
+        tri_list.append([e1, e7, e12])
+        tri_list.append([e1, e12, e2])
+        tri_list.append([e9, e8, e7])
+    elif (case == 115):
+        # Shelf: (v3, v4, v8)
+        tri_list.append([e2, e4, e8])
+        tri_list.append([e2, e12, e7])
+        tri_list.append([e2, e8, e7])
+    elif (case == 116):
+        # Rotated case 14 in paper: (v5, v6, v7, v3)
+        tri_list.append([e9, e2, e7])
+        tri_list.append([e9, e2, e10])
+        tri_list.append([e9, e8, e7])
+        tri_list.append([e2, e3, e7])
+    elif (case == 117):
+        # upper left plane + isolated corner v2 by inversion
+        _append_tris(tri_list, 2, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 119, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 118):
+        # Shelf: (v1, v4, v8)
+        tri_list.append([e1, e3, e7])
+        tri_list.append([e7, e1, e8])
+        tri_list.append([e1, e8, e9])
+    elif (case == 119):
+        # Upper left plane
+        tri_list.append([e4, e3, e7])
+        tri_list.append([e4, e8, e7])
+    elif (case == 120):
+        # Shelf: (v1, v2, v3) + isolated corner v8
+        _append_tris(tri_list, 7, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 127, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 121):
+        # Front right plane + isolated corner v8
+        _append_tris(tri_list, 6, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 127, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 122):
+        # Isolated corners v1, v3, v8
+        _append_tris(tri_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 4, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 127, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 123):
+        # Isolated corners v3, v8
+        _append_tris(tri_list, 4, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 127, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 124):
+        # Front lower plane + isolated corner v8
+        _append_tris(tri_list, 3, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 127, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 125):
+        # Isolated corners v2, v8
+        _append_tris(tri_list, 2, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 127)
+    elif (case == 126):
+        # Isolated corners v1, v8
+        _append_tris(tri_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 127, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 127):
+        # Isolated corner v8
+        tri_list.append([e11, e7, e8])
+
+    return
diff --git a/skimage/measure/marching_cubes.py b/skimage/measure/marching_cubes.py
new file mode 100644
index 00000000..10980580
--- /dev/null
+++ b/skimage/measure/marching_cubes.py
@@ -0,0 +1,200 @@
+import numpy as np
+from . import _marching_cubes
+
+
+def marching_cubes(volume, level, sampling=(1., 1., 1.)):
+    """
+    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` if
+        not provided in this format.
+    level : float
+        Contour value to search for isosurfaces in `volume`.
+    sampling : length-3 tuple of floats
+        Voxel spacing in spatial dimensions corresponding to numpy array
+        indexing dimensions (M, N, P) as in `volume`.
+
+    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`.
+
+    Notes
+    -----
+    The marching cubes algorithm is implemented as described in [1]_.
+    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 [1]_, 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 [1]_, 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.
+
+    To quantify the area of an isosurface generated by this algorithm, pass
+    the output 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, tris = 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],
+                             tris)
+    >>> mlab.show()
+
+    References
+    ----------
+    .. [1] 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).
+
+    See Also
+    --------
+    skimage.measure.mesh_surface_area
+
+    """
+    # Check inputs
+    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)
+
+    # Extract raw triangles using marching cubes in Cython
+    #   Returns a list of length-3 lists, each sub-list containing three
+    #   tuples. The tuples hold (x, y, z) coordinates for triangle vertices.
+    # Note: this algorithm is fast, but returns degenerate "triangles" which
+    #   have repeated vertices - and equivalent vertices are redundantly
+    #   placed in every triangle they connect with.
+    raw_tris = _marching_cubes.iterate_and_store_3d(volume, float(level),
+                                                    sampling)
+
+    # Find and collect unique vertices, storing triangle verts as indices.
+    # Removes much redundancy and eliminates degenerate "triangles".
+    vert_list, tri_list = _unpack_unique_verts(raw_tris)
+
+    return vert_list, tri_list
+
+
+def _unpack_unique_verts(trilist):
+    """
+    Converts a list of lists of tuples corresponding to triangle vertices into
+    a unique vertex list, and a list of triangle faces w/indices corresponding
+    to entries of the vertex list.
+
+    """
+    idx = 0
+    vert_index = {}
+    vert_list = []
+    tri_list = []
+
+    # Iterate over triangles
+    for i in range(len(trilist)):
+        templist = []
+
+        # Only parse vertices from non-degenerate triangles
+        if not ((trilist[i][0] == trilist[i][1]) or
+                (trilist[i][0] == trilist[i][2]) or
+                (trilist[i][1] == trilist[i][2])):
+
+            # Iterate over vertices within each triangle
+            for j in range(3):
+                vert = trilist[i][j]
+
+                # Check if a new unique vertex found
+                if vert not in vert_index:
+                    vert_index[vert] = idx
+                    templist.append(idx)
+                    vert_list.append(vert)
+                    idx += 1
+                else:
+                    templist.append(vert_index[vert])
+
+            tri_list.append(templist)
+
+    return vert_list, tri_list
+
+
+def mesh_surface_area(verts, tris):
+    """
+    Compute surface area, given vertices & triangular faces
+
+    Parameters
+    ----------
+    verts : list
+        List of length-3 NumPy arrays containing vertex coordinates.
+        Units in each dimension should be consistent.
+
+    tris : list
+        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.
+
+    Notes
+    -----
+    The arguments expected by this function are the exact outputs from
+    `skimage.measure.marching_cubes`. For unit correct output, ensure correct
+    `spacing` was passed to `skimage.measure.marching_cubes`.
+
+    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])
+    a = actual_verts[:, 0, :] - actual_verts[:, 1, :]
+    b = actual_verts[:, 0, :] - actual_verts[:, 2, :]
+    del actual_verts
+
+    # Area of triangle = 1/2 * Euclidean norm of cross product
+    return ((np.cross(a, b) ** 2).sum(axis=1) ** 0.5).sum() / 2.

From 84618af7c7fa945cca946522e2b5483a195cb762 Mon Sep 17 00:00:00 2001
From: "Josh Warner (Mac)" <warner.joshua@mayo.edu>
Date: Sat, 16 Mar 2013 21:14:13 -0500
Subject: [PATCH 02/24] import `marching_cubes` and `mesh_surface_area`

---
 skimage/measure/__init__.py | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/skimage/measure/__init__.py b/skimage/measure/__init__.py
index 616071d3..108cd7d9 100755
--- a/skimage/measure/__init__.py
+++ b/skimage/measure/__init__.py
@@ -1,4 +1,5 @@
 from .find_contours import find_contours
+from ._marching_cubes import marching_cubes, mesh_surface_area
 from ._regionprops import regionprops, perimeter
 from ._structural_similarity import structural_similarity
 from ._polygon import approximate_polygon, subdivide_polygon
@@ -21,4 +22,7 @@ __all__ = ['find_contours',
            'moments',
            'moments_central',
            'moments_normalized',
-           'moments_hu']
+           'moments_hu',
+           'sum_blocks',
+           'marching_cubes',
+           'mesh_surface_area']

From 855031ad477012a8ec4bda3c1cf504923f7fc4cc Mon Sep 17 00:00:00 2001
From: "Josh Warner (Mac)" <warner.joshua@mayo.edu>
Date: Sat, 16 Mar 2013 21:51:33 -0500
Subject: [PATCH 03/24] add _marching_cubes.pyx to list of Cython modules to
 build

---
 skimage/measure/setup.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/skimage/measure/setup.py b/skimage/measure/setup.py
index 21d9964e..6bc8212d 100644
--- a/skimage/measure/setup.py
+++ b/skimage/measure/setup.py
@@ -14,6 +14,7 @@ def configuration(parent_package='', top_path=None):
 
     cython(['_find_contours.pyx'], working_path=base_path)
     cython(['_moments.pyx'], working_path=base_path)
+    cython(['_marching_cubes.pyx'], working_path=base_path)
 
     config.add_extension('_find_contours', sources=['_find_contours.c'],
                          include_dirs=[get_numpy_include_dirs()])

From 5b990a0856b7568aeab6dd486b37cce0ea44b618 Mon Sep 17 00:00:00 2001
From: "Josh Warner (Mac)" <warner.joshua@mayo.edu>
Date: Sat, 16 Mar 2013 22:22:26 -0500
Subject: [PATCH 04/24] fix import of _marching_cubes

---
 skimage/measure/setup.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/skimage/measure/setup.py b/skimage/measure/setup.py
index 6bc8212d..4b18a3a0 100644
--- a/skimage/measure/setup.py
+++ b/skimage/measure/setup.py
@@ -20,6 +20,8 @@ def configuration(parent_package='', top_path=None):
                          include_dirs=[get_numpy_include_dirs()])
     config.add_extension('_moments', sources=['_moments.c'],
                          include_dirs=[get_numpy_include_dirs()])
+    config.add_extension('_marching_cubes', sources=['_marching_cubes.c'],
+                         include_dirs=[get_numpy_include_dirs()])
 
     return config
 

From 9779cd5e4a0017082c9ede6fba2b0d64ca456b99 Mon Sep 17 00:00:00 2001
From: "Josh Warner (Mac)" <warner.joshua@mayo.edu>
Date: Sun, 17 Mar 2013 00:46:39 -0500
Subject: [PATCH 05/24] fix: bug in second `_append_tris` call in case 125

---
 skimage/measure/_marching_cubes.pyx | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/skimage/measure/_marching_cubes.pyx b/skimage/measure/_marching_cubes.pyx
index 8858acf4..be7c601f 100644
--- a/skimage/measure/_marching_cubes.pyx
+++ b/skimage/measure/_marching_cubes.pyx
@@ -910,7 +910,8 @@ def _append_tris(list tri_list, unsigned char case, tuple e1, tuple e2,
         # Isolated corners v2, v8
         _append_tris(tri_list, 2, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 127)
+        _append_tris(tri_list, 127, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
     elif (case == 126):
         # Isolated corners v1, v8
         _append_tris(tri_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,

From 9a86aba2263e4729bf0b0c396b9fc514f0ae5b3a Mon Sep 17 00:00:00 2001
From: "Josh Warner (Mac)" <warner.joshua@mayo.edu>
Date: Sun, 17 Mar 2013 08:48:32 -0500
Subject: [PATCH 06/24] rename files to conform with naming conventions

---
 skimage/measure/{marching_cubes.py => _marching_cubes.py}  | 7 +++----
 .../{_marching_cubes.pyx => _marching_cubes_cy.pyx}        | 2 +-
 skimage/measure/setup.py                                   | 4 ++--
 3 files changed, 6 insertions(+), 7 deletions(-)
 rename skimage/measure/{marching_cubes.py => _marching_cubes.py} (97%)
 rename skimage/measure/{_marching_cubes.pyx => _marching_cubes_cy.pyx} (99%)

diff --git a/skimage/measure/marching_cubes.py b/skimage/measure/_marching_cubes.py
similarity index 97%
rename from skimage/measure/marching_cubes.py
rename to skimage/measure/_marching_cubes.py
index 10980580..d375c9fe 100644
--- a/skimage/measure/marching_cubes.py
+++ b/skimage/measure/_marching_cubes.py
@@ -1,5 +1,5 @@
 import numpy as np
-from . import _marching_cubes
+from . import _marching_cubes_cy
 
 
 def marching_cubes(volume, level, sampling=(1., 1., 1.)):
@@ -111,8 +111,8 @@ def marching_cubes(volume, level, sampling=(1., 1., 1.)):
     # Note: this algorithm is fast, but returns degenerate "triangles" which
     #   have repeated vertices - and equivalent vertices are redundantly
     #   placed in every triangle they connect with.
-    raw_tris = _marching_cubes.iterate_and_store_3d(volume, float(level),
-                                                    sampling)
+    raw_tris = _marching_cubes_cy.iterate_and_store_3d(volume, float(level),
+                                                       sampling)
 
     # Find and collect unique vertices, storing triangle verts as indices.
     # Removes much redundancy and eliminates degenerate "triangles".
@@ -169,7 +169,6 @@ def mesh_surface_area(verts, tris):
     verts : list
         List of length-3 NumPy arrays containing vertex coordinates.
         Units in each dimension should be consistent.
-
     tris : list
         List of length-3 lists of integers, referencing vertex coordinates as
         provided in `verts`
diff --git a/skimage/measure/_marching_cubes.pyx b/skimage/measure/_marching_cubes_cy.pyx
similarity index 99%
rename from skimage/measure/_marching_cubes.pyx
rename to skimage/measure/_marching_cubes_cy.pyx
index be7c601f..68faa8de 100644
--- a/skimage/measure/_marching_cubes.pyx
+++ b/skimage/measure/_marching_cubes_cy.pyx
@@ -166,7 +166,7 @@ def iterate_and_store_3d(cnp.ndarray[double, ndim=3] arr,
                     e3 = r0 + _get_fraction(v4, v3, level), c1, d0
                     e4 = r0, c0 + _get_fraction(v1, v4, level), d0
 
-            # These must be calculated at each point uunless we implemented a
+            # These must be calculated at each point unless we implemented a
             # large, growing lookup table for all adjacent values; could save
             # ~30% in terms of runtime at the expense of memory usage and
             # much greater complexity.
diff --git a/skimage/measure/setup.py b/skimage/measure/setup.py
index 4b18a3a0..19d69429 100644
--- a/skimage/measure/setup.py
+++ b/skimage/measure/setup.py
@@ -14,13 +14,13 @@ def configuration(parent_package='', top_path=None):
 
     cython(['_find_contours.pyx'], working_path=base_path)
     cython(['_moments.pyx'], working_path=base_path)
-    cython(['_marching_cubes.pyx'], working_path=base_path)
+    cython(['_marching_cubes_cy.pyx'], working_path=base_path)
 
     config.add_extension('_find_contours', sources=['_find_contours.c'],
                          include_dirs=[get_numpy_include_dirs()])
     config.add_extension('_moments', sources=['_moments.c'],
                          include_dirs=[get_numpy_include_dirs()])
-    config.add_extension('_marching_cubes', sources=['_marching_cubes.c'],
+    config.add_extension('_marching_cubes_cy', sources=['_marching_cubes.c'],
                          include_dirs=[get_numpy_include_dirs()])
 
     return config

From d2e51aa03533613c26c78d36c785af40bcceca4d Mon Sep 17 00:00:00 2001
From: "Josh Warner (Mac)" <warner.joshua@mayo.edu>
Date: Sun, 17 Mar 2013 09:07:56 -0500
Subject: [PATCH 07/24] fix import error from `_marching_cubes_cy`

---
 skimage/measure/setup.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/skimage/measure/setup.py b/skimage/measure/setup.py
index 19d69429..be57ca7b 100644
--- a/skimage/measure/setup.py
+++ b/skimage/measure/setup.py
@@ -20,7 +20,8 @@ def configuration(parent_package='', top_path=None):
                          include_dirs=[get_numpy_include_dirs()])
     config.add_extension('_moments', sources=['_moments.c'],
                          include_dirs=[get_numpy_include_dirs()])
-    config.add_extension('_marching_cubes_cy', sources=['_marching_cubes.c'],
+    config.add_extension('_marching_cubes_cy',
+                         sources=['_marching_cubes_cy.c'],
                          include_dirs=[get_numpy_include_dirs()])
 
     return config

From 44c04996c7c25a1d531744aec431a889a22aff12 Mon Sep 17 00:00:00 2001
From: JDWarner <warner.joshua@mayo.edu>
Date: Wed, 20 Mar 2013 15:42:57 -0500
Subject: [PATCH 08/24] move `unpack_unique_verts` to Cython for speed

---
 skimage/measure/_marching_cubes.py     | 41 +------------------------
 skimage/measure/_marching_cubes_cy.pyx | 42 ++++++++++++++++++++++++++
 2 files changed, 43 insertions(+), 40 deletions(-)

diff --git a/skimage/measure/_marching_cubes.py b/skimage/measure/_marching_cubes.py
index d375c9fe..e5c4bcff 100644
--- a/skimage/measure/_marching_cubes.py
+++ b/skimage/measure/_marching_cubes.py
@@ -116,46 +116,7 @@ def marching_cubes(volume, level, sampling=(1., 1., 1.)):
 
     # Find and collect unique vertices, storing triangle verts as indices.
     # Removes much redundancy and eliminates degenerate "triangles".
-    vert_list, tri_list = _unpack_unique_verts(raw_tris)
-
-    return vert_list, tri_list
-
-
-def _unpack_unique_verts(trilist):
-    """
-    Converts a list of lists of tuples corresponding to triangle vertices into
-    a unique vertex list, and a list of triangle faces w/indices corresponding
-    to entries of the vertex list.
-
-    """
-    idx = 0
-    vert_index = {}
-    vert_list = []
-    tri_list = []
-
-    # Iterate over triangles
-    for i in range(len(trilist)):
-        templist = []
-
-        # Only parse vertices from non-degenerate triangles
-        if not ((trilist[i][0] == trilist[i][1]) or
-                (trilist[i][0] == trilist[i][2]) or
-                (trilist[i][1] == trilist[i][2])):
-
-            # Iterate over vertices within each triangle
-            for j in range(3):
-                vert = trilist[i][j]
-
-                # Check if a new unique vertex found
-                if vert not in vert_index:
-                    vert_index[vert] = idx
-                    templist.append(idx)
-                    vert_list.append(vert)
-                    idx += 1
-                else:
-                    templist.append(vert_index[vert])
-
-            tri_list.append(templist)
+    vert_list, tri_list = _marching_cubes_cy.unpack_unique_verts(raw_tris)
 
     return vert_list, tri_list
 
diff --git a/skimage/measure/_marching_cubes_cy.pyx b/skimage/measure/_marching_cubes_cy.pyx
index 68faa8de..186ea8e8 100644
--- a/skimage/measure/_marching_cubes_cy.pyx
+++ b/skimage/measure/_marching_cubes_cy.pyx
@@ -13,6 +13,48 @@ cdef inline double _get_fraction(double from_value, double to_value,
     return ((level - from_value) / (to_value - from_value))
 
 
+def unpack_unique_verts(list trilist):
+    """
+    Converts a list of lists of tuples corresponding to triangle vertices into
+    a unique vertex list, and a list of triangle faces w/indices corresponding
+    to entries of the vertex list.
+
+    """
+    cdef Py_ssize_t idx = 0
+    cdef Py_ssize_t n_tris = len(trilist)
+    cdef Py_ssize_t i, j
+    cdef dict vert_index = {}
+    cdef list vert_list = []
+    cdef list tri_list = []
+    cdef list templist
+
+    # Iterate over triangles
+    for i in range(n_tris):
+        templist = []
+
+        # Only parse vertices from non-degenerate triangles
+        if not ((trilist[i][0] == trilist[i][1]) or
+                (trilist[i][0] == trilist[i][2]) or
+                (trilist[i][1] == trilist[i][2])):
+
+            # Iterate over vertices within each triangle
+            for j in range(3):
+                vert = trilist[i][j]
+
+                # Check if a new unique vertex found
+                if vert not in vert_index:
+                    vert_index[vert] = idx
+                    templist.append(idx)
+                    vert_list.append(vert)
+                    idx += 1
+                else:
+                    templist.append(vert_index[vert])
+
+            tri_list.append(templist)
+
+    return vert_list, tri_list
+
+
 def iterate_and_store_3d(cnp.ndarray[double, ndim=3] arr,
                          double level, tuple sampling=(1., 1., 1.)):
     """Iterate across the given array in a marching-cubes fashion,

From 296d1020474540f2397a111f815d189927c925c2 Mon Sep 17 00:00:00 2001
From: JDWarner <warner.joshua@mayo.edu>
Date: Wed, 20 Mar 2013 16:35:55 -0500
Subject: [PATCH 09/24] fix: symmetric ambiguous dual-plane cases require
 special treatment

---
 skimage/measure/_marching_cubes_cy.pyx | 29 +++++++++++++++++++++++++-
 1 file changed, 28 insertions(+), 1 deletion(-)

diff --git a/skimage/measure/_marching_cubes_cy.pyx b/skimage/measure/_marching_cubes_cy.pyx
index 186ea8e8..2fe1a1bd 100644
--- a/skimage/measure/_marching_cubes_cy.pyx
+++ b/skimage/measure/_marching_cubes_cy.pyx
@@ -184,7 +184,10 @@ def iterate_and_store_3d(cnp.ndarray[double, ndim=3] arr,
             # Cases 0 and 255 are entirely below/above the contour.
 
             if cube_case > 127:
-                cube_case = 255 - cube_case
+                if ((cube_case != 150) and
+                     (cube_case != 170) and
+                     (cube_case != 195)):
+                    cube_case = 255 - cube_case
 
             # Calculate cube edges, to become triangulation vertices.
             # If we moved in a convenient direction, save 1/3 of the effort by
@@ -963,5 +966,29 @@ def _append_tris(list tri_list, unsigned char case, tuple e1, tuple e2,
     elif (case == 127):
         # Isolated corner v8
         tri_list.append([e11, e7, e8])
+    elif (case == 150):
+        # AMBIGUOUS CASE: back right and front left planes
+        # In these cube_case > 127 cases, the vertices are identical BUT
+        # they are connected in the opposite fashion.
+        _append_tris(tri_list, 6, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 111, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 170):
+        # AMBIGUOUS CASE: upper left and lower right planes
+        # In these cube_case > 127 cases, the vertices are identical BUT
+        # they are connected in the opposite fashion.
+        _append_tris(tri_list, 119, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 34, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+    elif (case == 195):
+        # AMBIGUOUS CASE: back upper and front lower planes
+        # In these cube_case > 127 cases, the vertices are identical BUT
+        # they are connected in the opposite fashion.
+        _append_tris(tri_list, 63, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
+        _append_tris(tri_list, 3, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+                     e11, e12)
 
     return

From 165de870d0cbee07ee7547c4df3ea4f35c55aaba Mon Sep 17 00:00:00 2001
From: "Josh Warner (Mac)" <warner.joshua@mayo.edu>
Date: Thu, 28 Mar 2013 14:56:57 -0500
Subject: [PATCH 10/24] fix error in triangulation of case 120

---
 skimage/measure/_marching_cubes_cy.pyx | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/skimage/measure/_marching_cubes_cy.pyx b/skimage/measure/_marching_cubes_cy.pyx
index 2fe1a1bd..cf3ad93d 100644
--- a/skimage/measure/_marching_cubes_cy.pyx
+++ b/skimage/measure/_marching_cubes_cy.pyx
@@ -920,10 +920,10 @@ def _append_tris(list tri_list, unsigned char case, tuple e1, tuple e2,
         tri_list.append([e4, e3, e7])
         tri_list.append([e4, e8, e7])
     elif (case == 120):
-        # Shelf: (v1, v2, v3) + isolated corner v8
-        _append_tris(tri_list, 7, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        # Shelf: (v5, v6, v7) + isolated corner v4
+        _append_tris(tri_list, 112, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 127, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(tri_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 121):
         # Front right plane + isolated corner v8

From f8cc4799da581057849416ced72fc33198cae45b Mon Sep 17 00:00:00 2001
From: "Josh Warner (Mac)" <warner.joshua@mayo.edu>
Date: Fri, 31 May 2013 09:46:11 -0500
Subject: [PATCH 11/24] FIX: dimension ordering now correct (x, y, z)

---
 skimage/measure/_marching_cubes_cy.pyx | 27 ++++++++++----------------
 1 file changed, 10 insertions(+), 17 deletions(-)

diff --git a/skimage/measure/_marching_cubes_cy.pyx b/skimage/measure/_marching_cubes_cy.pyx
index cf3ad93d..18e14ddb 100644
--- a/skimage/measure/_marching_cubes_cy.pyx
+++ b/skimage/measure/_marching_cubes_cy.pyx
@@ -78,14 +78,7 @@ def iterate_and_store_3d(cnp.ndarray[double, ndim=3] arr,
     cdef Py_ssize_t n
     cdef bint odd_sampling, plus_z
     plus_z = False
-    if ((sampling == (1., 1., 1.)) or
-         (sampling == (1., 1., 1)) or
-         (sampling == (1., 1, 1.)) or
-         (sampling == (1, 1., 1.)) or
-         (sampling == (1, 1, 1.)) or
-         (sampling == (1., 1, 1)) or
-         (sampling == (1, 1., 1)) or
-         (sampling == (1, 1, 1))):
+    if [float(i) for i in sampling] == [1.0, 1.0, 1.0]:
         odd_sampling = False
     else:
         odd_sampling = True
@@ -159,14 +152,14 @@ def iterate_and_store_3d(cnp.ndarray[double, ndim=3] arr,
 
         # We use a right-handed coordinate system, UNlike the paper, but want
         # to index in agreement - the coordinate adjustment takes place here.
-        v1 = arr[x0, y0, z0]
-        v2 = arr[x1, y0, z0]
-        v3 = arr[x1, y1, z0]
-        v4 = arr[x0, y1, z0]
-        v5 = arr[x0, y0, z1]
-        v6 = arr[x1, y0, z1]
-        v7 = arr[x1, y1, z1]
-        v8 = arr[x0, y1, z1]
+        v1 = arr[z0, y0, x0]
+        v2 = arr[z0, y0, x1]
+        v3 = arr[z0, y1, x1]
+        v4 = arr[z0, y1, x0]
+        v5 = arr[z1, y0, x0]
+        v6 = arr[z1, y0, x1]
+        v7 = arr[z1, y1, x1]
+        v8 = arr[z1, y1, x0]
 
         # Unique triangulation cases
         cube_case = 0
@@ -240,7 +233,7 @@ def iterate_and_store_3d(cnp.ndarray[double, ndim=3] arr,
                            e6, e7, e8, e9, e10, e11, e12)
 
         # Advance the coords indices
-        if coords[2] < arr.shape[2] - 2:
+        if coords[2] < arr.shape[0] - 2:
             coords[2] += 1
             plus_z = True
         elif coords[1] < arr.shape[1] - 2:

From 164b416bd9abd4a13fee7bc3dc70a8f9fd17ea20 Mon Sep 17 00:00:00 2001
From: "Josh Warner (Mac)" <warner.joshua@mayo.edu>
Date: Fri, 31 May 2013 10:27:54 -0500
Subject: [PATCH 12/24] FIX: volume now directly cast, docfixes, C-contiguous
 ordering

---
 skimage/measure/_marching_cubes.py | 51 ++++++++++++++----------------
 1 file changed, 24 insertions(+), 27 deletions(-)

diff --git a/skimage/measure/_marching_cubes.py b/skimage/measure/_marching_cubes.py
index e5c4bcff..418f945e 100644
--- a/skimage/measure/_marching_cubes.py
+++ b/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

From fcd20e5660742aa3d5d085ff79249caefe8a8bb3 Mon Sep 17 00:00:00 2001
From: "Josh Warner (Mac)" <warner.joshua@mayo.edu>
Date: Fri, 31 May 2013 10:29:09 -0500
Subject: [PATCH 13/24] FIX: use memoryviews, docfix, and change tri_list to
 face_list

---
 skimage/measure/_marching_cubes_cy.pyx | 666 ++++++++++++-------------
 1 file changed, 333 insertions(+), 333 deletions(-)

diff --git a/skimage/measure/_marching_cubes_cy.pyx b/skimage/measure/_marching_cubes_cy.pyx
index 18e14ddb..7a280366 100644
--- a/skimage/measure/_marching_cubes_cy.pyx
+++ b/skimage/measure/_marching_cubes_cy.pyx
@@ -13,11 +13,11 @@ cdef inline double _get_fraction(double from_value, double to_value,
     return ((level - from_value) / (to_value - from_value))
 
 
-def unpack_unique_verts(list trilist):
+def unpack_unique_vert_list(list trilist):
     """
-    Converts a list of lists of tuples corresponding to triangle vertices into
-    a unique vertex list, and a list of triangle faces w/indices corresponding
-    to entries of the vertex list.
+    Convert_list a list of lists of tuples corresponding to triangle vertices
+    into a unique vertex list, and a list of triangle faces w/indices
+    corresponding to entries of the vertex list.
 
     """
     cdef Py_ssize_t idx = 0
@@ -25,7 +25,7 @@ def unpack_unique_verts(list trilist):
     cdef Py_ssize_t i, j
     cdef dict vert_index = {}
     cdef list vert_list = []
-    cdef list tri_list = []
+    cdef list face_list = []
     cdef list templist
 
     # Iterate over triangles
@@ -50,13 +50,13 @@ def unpack_unique_verts(list trilist):
                 else:
                     templist.append(vert_index[vert])
 
-            tri_list.append(templist)
+            face_list.append(templist)
 
-    return vert_list, tri_list
+    return vert_list, face_list
 
 
-def iterate_and_store_3d(cnp.ndarray[double, ndim=3] arr,
-                         double level, tuple sampling=(1., 1., 1.)):
+def iterate_and_store_3d(double[:, :, ::1] arr, double level,
+                         tuple sampling=(1., 1., 1.)):
     """Iterate across the given array in a marching-cubes fashion,
     looking for volumes with edges that cross 'level'. If such a volume is
     found, appropriate triangulations are added to a growing list of
@@ -71,9 +71,9 @@ def iterate_and_store_3d(cnp.ndarray[double, ndim=3] arr,
     if arr.shape[0] < 2 or arr.shape[1] < 2 or arr.shape[2] < 2:
         raise ValueError("Input array must be at least 2x2x2.")
     if len(sampling) != 3:
-        raise ValueError("`sampling` must be of form (double, double, double)")
+        raise ValueError("`sampling` must be (double, double, double)")
 
-    cdef list tri_list = []
+    cdef list face_list = []
     cdef list norm_list = []
     cdef Py_ssize_t n
     cdef bint odd_sampling, plus_z
@@ -228,8 +228,8 @@ def iterate_and_store_3d(cnp.ndarray[double, ndim=3] arr,
                 e12 = r1, c1, d0 + _get_fraction(v3, v7, level)
 
 
-            # Append appropriate triangles to the growing output `tri_list`
-            _append_tris(tri_list, cube_case, e1, e2, e3, e4, e5,
+            # Append appropriate triangles to the growing output `face_list`
+            _append_tris(face_list, cube_case, e1, e2, e3, e4, e5,
                            e6, e7, e8, e9, e10, e11, e12)
 
         # Advance the coords indices
@@ -246,10 +246,10 @@ def iterate_and_store_3d(cnp.ndarray[double, ndim=3] arr,
             coords[2] = 0
             plus_z = False
 
-    return tri_list
+    return face_list
 
 
-def _append_tris(list tri_list, unsigned char case, tuple e1, tuple e2,
+def _append_tris(list face_list, unsigned char case, tuple e1, tuple e2,
                  tuple e3, tuple e4, tuple e5, tuple e6, tuple e7, tuple e8,
                  tuple e9, tuple e10, tuple e11, tuple e12):
     # Permits recursive use for duplicated planes to conserve code - it's
@@ -257,731 +257,731 @@ def _append_tris(list tri_list, unsigned char case, tuple e1, tuple e2,
 
     if (case == 1):
         # front lower left corner
-        tri_list.append([e1, e4, e9])
+        face_list.append([e1, e4, e9])
     elif (case == 2):
         # front lower right corner
-        tri_list.append([e10, e2, e1])
+        face_list.append([e10, e2, e1])
     elif (case == 3):
         # front lower plane
-        tri_list.append([e2, e4, e9])
-        tri_list.append([e2, e9, e10])
+        face_list.append([e2, e4, e9])
+        face_list.append([e2, e9, e10])
     elif (case == 4):
         # front upper right corner
-        tri_list.append([e12, e3, e2])
+        face_list.append([e12, e3, e2])
     elif (case == 5):
         # lower left, upper right corners
-        _append_tris(tri_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 4, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 4, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 6):
         # front right plane
-        tri_list.append([e12, e3, e1])
-        tri_list.append([e12, e1, e10])
+        face_list.append([e12, e3, e1])
+        face_list.append([e12, e1, e10])
     elif (case == 7):
         # Shelf including v1, v2, v3
-        tri_list.append([e3, e4, e12])
-        tri_list.append([e4, e9, e12])
-        tri_list.append([e12, e9, e10])
+        face_list.append([e3, e4, e12])
+        face_list.append([e4, e9, e12])
+        face_list.append([e12, e9, e10])
     elif (case == 8):
         # front upper left corner
-        tri_list.append([e3, e11, e4])
+        face_list.append([e3, e11, e4])
     elif (case == 9):
         # front left plane
-        tri_list.append([e3, e11, e9])
-        tri_list.append([e3, e9, e1])
+        face_list.append([e3, e11, e9])
+        face_list.append([e3, e9, e1])
     elif (case == 10):
         # upper left, lower right corners
-        _append_tris(tri_list, 2, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 2, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 11):
         # Shelf including v4, v1, v2
-        tri_list.append([e3, e11, e2])
-        tri_list.append([e11, e10, e2])
-        tri_list.append([e11, e9, e10])
+        face_list.append([e3, e11, e2])
+        face_list.append([e11, e10, e2])
+        face_list.append([e11, e9, e10])
     elif (case == 12):
         # front upper plane
-        tri_list.append([e11, e4, e12])
-        tri_list.append([e2, e4, e12])
+        face_list.append([e11, e4, e12])
+        face_list.append([e2, e4, e12])
     elif (case == 13):
         # Shelf including v1, v4, v3
-        tri_list.append([e11, e9, e12])
-        tri_list.append([e12, e9, e1])
-        tri_list.append([e12, e1, e2])
+        face_list.append([e11, e9, e12])
+        face_list.append([e12, e9, e1])
+        face_list.append([e12, e1, e2])
     elif (case == 14):
         # Shelf including v2, v3, v4
-        tri_list.append([e11, e10, e12])
-        tri_list.append([e11, e4, e10])
-        tri_list.append([e4, e1, e10])
+        face_list.append([e11, e10, e12])
+        face_list.append([e11, e4, e10])
+        face_list.append([e4, e1, e10])
     elif (case == 15):
         # Plane parallel to x-axis through middle
-        tri_list.append([e11, e9, e12])
-        tri_list.append([e12, e9, e10])
+        face_list.append([e11, e9, e12])
+        face_list.append([e12, e9, e10])
     elif (case == 16):
         # back lower left corner
-        tri_list.append([e8, e9, e5])
+        face_list.append([e8, e9, e5])
     elif (case == 17):
         # lower left plane
-        tri_list.append([e4, e1, e8])
-        tri_list.append([e8, e1, e5])
+        face_list.append([e4, e1, e8])
+        face_list.append([e8, e1, e5])
     elif (case == 18):
         # lower left back, lower right front corners
-        _append_tris(tri_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 2, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 2, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 19):
         # Shelf including v1, v2, v5
-        tri_list.append([e8, e4, e2])
-        tri_list.append([e8, e2, e10])
-        tri_list.append([e8, e10, e5])
+        face_list.append([e8, e4, e2])
+        face_list.append([e8, e2, e10])
+        face_list.append([e8, e10, e5])
     elif (case == 20):
         # lower left back, upper right front corners
-        _append_tris(tri_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 4, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 4, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 21):
         # lower left plane + upper right front corner, v1, v3, v5
-        _append_tris(tri_list, 17, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 17, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 4, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 4, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 22):
         # front right plane + lower left back corner, v2, v3, v5
-        _append_tris(tri_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 6, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 6, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 23):
         # Rotated case 14 in the paper
-        tri_list.append([e3, e10, e8])
-        tri_list.append([e3, e10, e12])
-        tri_list.append([e8, e10, e5])
-        tri_list.append([e3, e4, e8])
+        face_list.append([e3, e10, e8])
+        face_list.append([e3, e10, e12])
+        face_list.append([e8, e10, e5])
+        face_list.append([e3, e4, e8])
     elif (case == 24):
         # upper front left, lower back left corners
-        _append_tris(tri_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 25):
         # Shelf including v1, v4, v5
-        tri_list.append([e1, e5, e3])
-        tri_list.append([e3, e8, e11])
-        tri_list.append([e3, e5, e8])
+        face_list.append([e1, e5, e3])
+        face_list.append([e3, e8, e11])
+        face_list.append([e3, e5, e8])
     elif (case == 26):
         # Three isolated corners
-        _append_tris(tri_list, 2, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 2, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 27):
         # Full corner v1, case 9 in paper: (v1, v2, v4, v5)
-        tri_list.append([e11, e3, e2])
-        tri_list.append([e11, e2, e10])
-        tri_list.append([e10, e11, e8])
-        tri_list.append([e8, e5, e10])
+        face_list.append([e11, e3, e2])
+        face_list.append([e11, e2, e10])
+        face_list.append([e10, e11, e8])
+        face_list.append([e8, e5, e10])
     elif (case == 28):
         # upper front plane + corner v5
-        _append_tris(tri_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 12, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 12, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 29):
         # special case of 11 in the paper: (v1, v3, v4, v5)
-        tri_list.append([e11, e5, e2])
-        tri_list.append([e11, e12, e2])
-        tri_list.append([e11, e5, e8])
-        tri_list.append([e2, e1, e5])
+        face_list.append([e11, e5, e2])
+        face_list.append([e11, e12, e2])
+        face_list.append([e11, e5, e8])
+        face_list.append([e2, e1, e5])
     elif (case == 30):
         # Shelf (v2, v3, v4) and lower left back corner
-        _append_tris(tri_list, 14, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 14, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 31):
         # Shelf: (v6, v7, v8) by inversion
-        tri_list.append([e11, e12, e10])
-        tri_list.append([e11, e8, e10])
-        tri_list.append([e8, e10, e5])
+        face_list.append([e11, e12, e10])
+        face_list.append([e11, e8, e10])
+        face_list.append([e8, e10, e5])
     elif (case == 32):
         # lower right back corner
-        tri_list.append([e6, e5, e10])
+        face_list.append([e6, e5, e10])
     elif (case == 33):
         # lower right back, lower left front corners
-        _append_tris(tri_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 34):
         # lower right plane
-        tri_list.append([e1, e2, e5])
-        tri_list.append([e2, e6, e5])
+        face_list.append([e1, e2, e5])
+        face_list.append([e2, e6, e5])
     elif (case == 35):
         # Shelf: v1, v2, v6
-        tri_list.append([e4, e2, e6])
-        tri_list.append([e4, e9, e6])
-        tri_list.append([e6, e9, e5])
+        face_list.append([e4, e2, e6])
+        face_list.append([e4, e9, e6])
+        face_list.append([e6, e9, e5])
     elif (case == 36):
         # upper right front, lower right back corners
-        _append_tris(tri_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 4, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 4, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 37):
         # lower left front, upper right front, lower right back corners
-        _append_tris(tri_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 4, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 4, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 38):
         # Shelf: v2, v3, v6
-        tri_list.append([e3, e1, e5])
-        tri_list.append([e3, e5, e12])
-        tri_list.append([e12, e5, e6])
+        face_list.append([e3, e1, e5])
+        face_list.append([e3, e5, e12])
+        face_list.append([e12, e5, e6])
     elif (case == 39):
         # Full corner v2: (v1, v2, v3, v6)
-        tri_list.append([e3, e4, e5])
-        tri_list.append([e4, e9, e5])
-        tri_list.append([e3, e5, e6])
-        tri_list.append([e3, e12, e6])
+        face_list.append([e3, e4, e5])
+        face_list.append([e4, e9, e5])
+        face_list.append([e3, e5, e6])
+        face_list.append([e3, e12, e6])
     elif (case == 40):
         # upper left front, lower right back corners
-        _append_tris(tri_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 41):
         # front left plane, lower right back corner
-        _append_tris(tri_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 9, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 9, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 42):
         # lower right plane, upper front left corner
-        _append_tris(tri_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 34, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 34, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 43):
         # Rotated case 11 in paper
-        tri_list.append([e11, e3, e9])
-        tri_list.append([e3, e9, e6])
-        tri_list.append([e3, e2, e6])
-        tri_list.append([e9, e5, e6])
+        face_list.append([e11, e3, e9])
+        face_list.append([e3, e9, e6])
+        face_list.append([e3, e2, e6])
+        face_list.append([e9, e5, e6])
     elif (case == 44):
         # upper front plane, lower right back corner
-        _append_tris(tri_list, 12, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 12, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 45):
         # Shelf: (v1, v3, v4) + lower right back corner
-        _append_tris(tri_list, 13, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 13, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 46):
         # Rotated case 14 in paper
-        tri_list.append([e4, e11, e12])
-        tri_list.append([e4, e12, e5])
-        tri_list.append([e12, e5, e6])
-        tri_list.append([e4, e5, e1])
+        face_list.append([e4, e11, e12])
+        face_list.append([e4, e12, e5])
+        face_list.append([e12, e5, e6])
+        face_list.append([e4, e5, e1])
     elif (case == 47):
         # Shelf: (v5, v8, v7) by inversion
-        tri_list.append([e11, e9, e12])
-        tri_list.append([e12, e9, e5])
-        tri_list.append([e12, e5, e6])
+        face_list.append([e11, e9, e12])
+        face_list.append([e12, e9, e5])
+        face_list.append([e12, e5, e6])
     elif (case == 48):
         # Back lower plane
-        tri_list.append([e9, e10, e6])
-        tri_list.append([e9, e6, e8])
+        face_list.append([e9, e10, e6])
+        face_list.append([e9, e6, e8])
     elif (case == 49):
         # Shelf: (v1, v5, v6)
-        tri_list.append([e4, e8, e6])
-        tri_list.append([e4, e6, e1])
-        tri_list.append([e6, e1, e10])
+        face_list.append([e4, e8, e6])
+        face_list.append([e4, e6, e1])
+        face_list.append([e6, e1, e10])
     elif (case == 50):
         # Shelf: (v2, v5, v6)
-        tri_list.append([e8, e6, e2])
-        tri_list.append([e8, e2, e1])
-        tri_list.append([e8, e9, e1])
+        face_list.append([e8, e6, e2])
+        face_list.append([e8, e2, e1])
+        face_list.append([e8, e9, e1])
     elif (case == 51):
         # Plane through middle of cube, parallel to x-z axis
-        tri_list.append([e4, e8, e2])
-        tri_list.append([e8, e2, e6])
+        face_list.append([e4, e8, e2])
+        face_list.append([e8, e2, e6])
     elif (case == 52):
         # Back lower plane, and front upper right corner
-        _append_tris(tri_list, 48, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 48, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 4, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 4, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 53):
         # Shelf (v1, v5, v6) and front upper right corner
-        _append_tris(tri_list, 49, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 49, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 4, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 4, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 54):
         # Rotated case 11 from paper (v2, v3, v5, v6)
-        tri_list.append([e1, e9, e3])
-        tri_list.append([e9, e3, e6])
-        tri_list.append([e9, e8, e6])
-        tri_list.append([e12, e3, e6])
+        face_list.append([e1, e9, e3])
+        face_list.append([e9, e3, e6])
+        face_list.append([e9, e8, e6])
+        face_list.append([e12, e3, e6])
     elif (case == 55):
         # Shelf: (v4, v8, v7) by inversion
-        tri_list.append([e4, e8, e6])
-        tri_list.append([e4, e6, e3])
-        tri_list.append([e6, e3, e12])
+        face_list.append([e4, e8, e6])
+        face_list.append([e4, e6, e3])
+        face_list.append([e6, e3, e12])
     elif (case == 56):
         # Back lower plane + upper left front corner
-        _append_tris(tri_list, 48, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 48, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 57):
         # Rotated case 14 from paper (v4, v1, v5, v6)
-        tri_list.append([e3, e11, e8])
-        tri_list.append([e3, e8, e10])
-        tri_list.append([e10, e6, e8])
-        tri_list.append([e3, e1, e10])
+        face_list.append([e3, e11, e8])
+        face_list.append([e3, e8, e10])
+        face_list.append([e10, e6, e8])
+        face_list.append([e3, e1, e10])
     elif (case == 58):
         # Shelf: (v2, v6, v5) + upper left front corner
-        _append_tris(tri_list, 50, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 50, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 59):
         # Shelf: (v3, v7, v8) by inversion
-        tri_list.append([e2, e6, e8])
-        tri_list.append([e8, e2, e3])
-        tri_list.append([e8, e3, e11])
+        face_list.append([e2, e6, e8])
+        face_list.append([e8, e2, e3])
+        face_list.append([e8, e3, e11])
     elif (case == 60):
         # AMBIGUOUS CASE: parallel planes (front upper, back lower)
-        _append_tris(tri_list, 48, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 48, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 12, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 12, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 61):
         # Upper back plane + lower right front corner by inversion
-        _append_tris(tri_list, 63, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 63, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 2, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 2, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 62):
         # Upper back plane + lower left front corner by inversion
-        _append_tris(tri_list, 63, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 63, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 63):
         # Upper back plane
-        tri_list.append([e11, e12, e6])
-        tri_list.append([e11, e8, e6])
+        face_list.append([e11, e12, e6])
+        face_list.append([e11, e8, e6])
     elif (case == 64):
         # Upper right back corner
-        tri_list.append([e12, e7, e6])
+        face_list.append([e12, e7, e6])
     elif (case == 65):
         # upper right back, lower left front corners
-        _append_tris(tri_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 66):
         # upper right back, lower right front corners
-        _append_tris(tri_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 2, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 2, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 67):
         # lower front plane + upper right back corner
-        _append_tris(tri_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 3, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 3, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 68):
         # upper right plane
-        tri_list.append([e3, e2, e6])
-        tri_list.append([e3, e7, e6])
+        face_list.append([e3, e2, e6])
+        face_list.append([e3, e7, e6])
     elif (case == 69):
         # Upper right plane, lower left front corner
-        _append_tris(tri_list, 68, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 68, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 70):
         # Shelf: (v2, v3, v7)
-        tri_list.append([e1, e3, e7])
-        tri_list.append([e1, e10, e7])
-        tri_list.append([e7, e10, e6])
+        face_list.append([e1, e3, e7])
+        face_list.append([e1, e10, e7])
+        face_list.append([e7, e10, e6])
     elif (case == 71):
         # Rotated version of case 11 in paper (v1, v2, v3, v7)
-        tri_list.append([e10, e7, e4])
-        tri_list.append([e4, e3, e7])
-        tri_list.append([e10, e4, e9])
-        tri_list.append([e7, e10, e6])
+        face_list.append([e10, e7, e4])
+        face_list.append([e4, e3, e7])
+        face_list.append([e10, e4, e9])
+        face_list.append([e7, e10, e6])
     elif (case == 72):
         # upper left front, upper right back corners
-        _append_tris(tri_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 73):
         # front left plane, upper right back corner
-        _append_tris(tri_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 9, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 9, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 74):
         # Three isolated corners, exactly case 7 in paper
-        _append_tris(tri_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 2, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 2, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 75):
         # Shelf: (v1, v2, v4) + upper right back corner
-        _append_tris(tri_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 11, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 11, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 76):
         # Shelf: (v4, v3, v7)
-        tri_list.append([e4, e2, e6])
-        tri_list.append([e4, e11, e7])
-        tri_list.append([e4, e7, e6])
+        face_list.append([e4, e2, e6])
+        face_list.append([e4, e11, e7])
+        face_list.append([e4, e7, e6])
     elif (case == 77):
         # Rotated case 14 in paper (v1, v4, v3, v7)
-        tri_list.append([e11, e9, e1])
-        tri_list.append([e11, e1, e6])
-        tri_list.append([e1, e6, e2])
-        tri_list.append([e11, e6, e7])
+        face_list.append([e11, e9, e1])
+        face_list.append([e11, e1, e6])
+        face_list.append([e1, e6, e2])
+        face_list.append([e11, e6, e7])
     elif (case == 78):
         # Full corner v3: (v2, v3, v4, v7)
-        tri_list.append([e1, e4, e7])
-        tri_list.append([e1, e7, e6])
-        tri_list.append([e4, e11, e7])
-        tri_list.append([e1, e10, e6])
+        face_list.append([e1, e4, e7])
+        face_list.append([e1, e7, e6])
+        face_list.append([e4, e11, e7])
+        face_list.append([e1, e10, e6])
     elif (case == 79):
         # Shelf: (v6, v5, v8) by inversion
-        tri_list.append([e9, e11, e10])
-        tri_list.append([e11, e7, e10])
-        tri_list.append([e7, e10, e6])
+        face_list.append([e9, e11, e10])
+        face_list.append([e11, e7, e10])
+        face_list.append([e7, e10, e6])
     elif (case == 80):
         # lower left back, upper right back corners (v5, v7)
-        _append_tris(tri_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 81):
         # lower left plane, upper right back corner
-        _append_tris(tri_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 17, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 17, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 82):
         # isolated corners (v2, v5, v7)
-        _append_tris(tri_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 2, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 2, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 83):
         # Shelf: (v1, v2, v5) + upper right back corner
-        _append_tris(tri_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 19, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 19, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 84):
         # upper right plane, lower left back corner
-        _append_tris(tri_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 68, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 68, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 85):
         # AMBIGUOUS CASE: upper right and lower left parallel planes
-        _append_tris(tri_list, 17, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 17, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 68, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 68, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 86):
         # Shelf: (v2, v3, v7) + lower left back corner
-        _append_tris(tri_list, 70, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 70, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 87):
         # Upper left plane + lower right back corner, by inversion
-        _append_tris(tri_list, 119, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 119, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 88):
         # Isolated corners v4, v5, v7
-        _append_tris(tri_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 89):
         # Shelf: (v1, v4, v5) + isolated corner v7
-        _append_tris(tri_list, 25, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 25, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 90):
         # Four isolated corners v2, v4, v5, v7
-        _append_tris(tri_list, 2, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 2, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 91):
         # Three isolated corners, v3, v6, v8 by inversion
-        _append_tris(tri_list, 4, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 4, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 127, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 127, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 92):
         # Shelf (v4, v3, v7) + isolated corner v5
-        _append_tris(tri_list, 76, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 76, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 93):
         # Lower right plane + isolated corner v8 by inversion
-        _append_tris(tri_list, 127, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 127, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 34, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 34, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 94):
         # Isolated corners v1, v6, v8 by inversion
-        _append_tris(tri_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 127, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 127, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 95):
         # Isolated corners v6, v8 by inversion
-        _append_tris(tri_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 127, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 127, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 96):
         # back right plane
-        tri_list.append([e7, e12, e5])
-        tri_list.append([e5, e10, e12])
+        face_list.append([e7, e12, e5])
+        face_list.append([e5, e10, e12])
     elif (case == 97):
         # back right plane + isolated corner v1
-        _append_tris(tri_list, 96, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 96, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 98):
         # Shelf: (v2, v6, v7)
-        tri_list.append([e1, e7, e5])
-        tri_list.append([e7, e1, e12])
-        tri_list.append([e1, e12, e2])
+        face_list.append([e1, e7, e5])
+        face_list.append([e7, e1, e12])
+        face_list.append([e1, e12, e2])
     elif (case == 99):
         # Rotated case 14 in paper: (v1, v2, v6, v7)
-        tri_list.append([e9, e2, e7])
-        tri_list.append([e9, e2, e4])
-        tri_list.append([e2, e7, e12])
-        tri_list.append([e7, e9, e5])
+        face_list.append([e9, e2, e7])
+        face_list.append([e9, e2, e4])
+        face_list.append([e2, e7, e12])
+        face_list.append([e7, e9, e5])
     elif (case == 100):
         # Shelf: (v3, v6, v7)
-        tri_list.append([e3, e7, e5])
-        tri_list.append([e3, e5, e2])
-        tri_list.append([e2, e5, e10])
+        face_list.append([e3, e7, e5])
+        face_list.append([e3, e5, e2])
+        face_list.append([e2, e5, e10])
     elif (case == 101):
         # Shelf: (v3, v6, v7) + isolated corner v1
-        _append_tris(tri_list, 100, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 100, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 102):
         # Plane bisecting left-right halves of cube
-        tri_list.append([e1, e3, e7])
-        tri_list.append([e1, e7, e5])
+        face_list.append([e1, e3, e7])
+        face_list.append([e1, e7, e5])
     elif (case == 103):
         # Shelf: (v4, v5, v8) by inversion
-        tri_list.append([e3, e7, e5])
-        tri_list.append([e3, e5, e4])
-        tri_list.append([e4, e5, e9])
+        face_list.append([e3, e7, e5])
+        face_list.append([e3, e5, e4])
+        face_list.append([e4, e5, e9])
     elif (case == 104):
         # Back right plane + isolated corner v4
-        _append_tris(tri_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 96, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 96, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 105):
         # AMBIGUOUS CASE: back right and front left planes
-        _append_tris(tri_list, 96, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 96, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 9, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 9, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 106):
         # Shelf: (v2, v6, v7) + isolated corner v4
-        _append_tris(tri_list, 98, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 98, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 107):
         # Back left plane + isolated corner v3 by inversion
-        _append_tris(tri_list, 4, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 4, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 111, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 111, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 108):
         # Rotated case 11 from paper: (v4, v3, v7, v6)
-        tri_list.append([e4, e10, e7])
-        tri_list.append([e4, e10, e2])
-        tri_list.append([e4, e11, e7])
-        tri_list.append([e7, e10, e5])
+        face_list.append([e4, e10, e7])
+        face_list.append([e4, e10, e2])
+        face_list.append([e4, e11, e7])
+        face_list.append([e7, e10, e5])
     elif (case == 109):
         # Back left plane + isolated corner v2 by inversion
-        _append_tris(tri_list, 111, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 111, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 2, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 2, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 110):
         # Shelf: (v1, v5, v8) by inversion
-        tri_list.append([e1, e5, e7])
-        tri_list.append([e1, e7, e11])
-        tri_list.append([e1, e11, e4])
+        face_list.append([e1, e5, e7])
+        face_list.append([e1, e7, e11])
+        face_list.append([e1, e11, e4])
     elif (case == 111):
         # Back left plane
-        tri_list.append([e11, e9, e7])
-        tri_list.append([e9, e7, e5])
+        face_list.append([e11, e9, e7])
+        face_list.append([e9, e7, e5])
     elif (case == 112):
         # Shelf: (v5, v6, v7)
-        tri_list.append([e9, e10, e12])
-        tri_list.append([e9, e12, e7])
-        tri_list.append([e9, e7, e8])
+        face_list.append([e9, e10, e12])
+        face_list.append([e9, e12, e7])
+        face_list.append([e9, e7, e8])
     elif (case == 113):
         # Exactly case 11 from paper: (v1, v5, v6, v7)
-        tri_list.append([e1, e8, e12])
-        tri_list.append([e1, e8, e4])
-        tri_list.append([e8, e7, e12])
-        tri_list.append([e12, e1, e10])
+        face_list.append([e1, e8, e12])
+        face_list.append([e1, e8, e4])
+        face_list.append([e8, e7, e12])
+        face_list.append([e12, e1, e10])
     elif (case == 114):
         # Full corner v6: (v2, v6, v7, v5)
-        tri_list.append([e1, e9, e7])
-        tri_list.append([e1, e7, e12])
-        tri_list.append([e1, e12, e2])
-        tri_list.append([e9, e8, e7])
+        face_list.append([e1, e9, e7])
+        face_list.append([e1, e7, e12])
+        face_list.append([e1, e12, e2])
+        face_list.append([e9, e8, e7])
     elif (case == 115):
         # Shelf: (v3, v4, v8)
-        tri_list.append([e2, e4, e8])
-        tri_list.append([e2, e12, e7])
-        tri_list.append([e2, e8, e7])
+        face_list.append([e2, e4, e8])
+        face_list.append([e2, e12, e7])
+        face_list.append([e2, e8, e7])
     elif (case == 116):
         # Rotated case 14 in paper: (v5, v6, v7, v3)
-        tri_list.append([e9, e2, e7])
-        tri_list.append([e9, e2, e10])
-        tri_list.append([e9, e8, e7])
-        tri_list.append([e2, e3, e7])
+        face_list.append([e9, e2, e7])
+        face_list.append([e9, e2, e10])
+        face_list.append([e9, e8, e7])
+        face_list.append([e2, e3, e7])
     elif (case == 117):
         # upper left plane + isolated corner v2 by inversion
-        _append_tris(tri_list, 2, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 2, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 119, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 119, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 118):
         # Shelf: (v1, v4, v8)
-        tri_list.append([e1, e3, e7])
-        tri_list.append([e7, e1, e8])
-        tri_list.append([e1, e8, e9])
+        face_list.append([e1, e3, e7])
+        face_list.append([e7, e1, e8])
+        face_list.append([e1, e8, e9])
     elif (case == 119):
         # Upper left plane
-        tri_list.append([e4, e3, e7])
-        tri_list.append([e4, e8, e7])
+        face_list.append([e4, e3, e7])
+        face_list.append([e4, e8, e7])
     elif (case == 120):
         # Shelf: (v5, v6, v7) + isolated corner v4
-        _append_tris(tri_list, 112, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 112, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 121):
         # Front right plane + isolated corner v8
-        _append_tris(tri_list, 6, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 6, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 127, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 127, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 122):
         # Isolated corners v1, v3, v8
-        _append_tris(tri_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 4, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 4, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 127, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 127, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 123):
         # Isolated corners v3, v8
-        _append_tris(tri_list, 4, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 4, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 127, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 127, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 124):
         # Front lower plane + isolated corner v8
-        _append_tris(tri_list, 3, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 3, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 127, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 127, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 125):
         # Isolated corners v2, v8
-        _append_tris(tri_list, 2, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 2, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 127, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 127, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 126):
         # Isolated corners v1, v8
-        _append_tris(tri_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 127, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 127, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 127):
         # Isolated corner v8
-        tri_list.append([e11, e7, e8])
+        face_list.append([e11, e7, e8])
     elif (case == 150):
         # AMBIGUOUS CASE: back right and front left planes
         # In these cube_case > 127 cases, the vertices are identical BUT
         # they are connected in the opposite fashion.
-        _append_tris(tri_list, 6, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 6, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 111, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 111, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 170):
         # AMBIGUOUS CASE: upper left and lower right planes
         # In these cube_case > 127 cases, the vertices are identical BUT
         # they are connected in the opposite fashion.
-        _append_tris(tri_list, 119, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 119, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 34, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 34, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
     elif (case == 195):
         # AMBIGUOUS CASE: back upper and front lower planes
         # In these cube_case > 127 cases, the vertices are identical BUT
         # they are connected in the opposite fashion.
-        _append_tris(tri_list, 63, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 63, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
-        _append_tris(tri_list, 3, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
+        _append_tris(face_list, 3, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10,
                      e11, e12)
 
     return

From 19b8831a0d2cf3325cbfa2868bf8ad9432db0ccd Mon Sep 17 00:00:00 2001
From: "Josh Warner (Mac)" <warner.joshua@mayo.edu>
Date: Sat, 1 Jun 2013 21:26:16 -0500
Subject: [PATCH 14/24] FIX: incorrect name for `unpack_unique_verts`, revert
 axis ordering

---
 skimage/measure/_marching_cubes_cy.pyx | 22 +++++++++++-----------
 1 file changed, 11 insertions(+), 11 deletions(-)

diff --git a/skimage/measure/_marching_cubes_cy.pyx b/skimage/measure/_marching_cubes_cy.pyx
index 7a280366..ec817b15 100644
--- a/skimage/measure/_marching_cubes_cy.pyx
+++ b/skimage/measure/_marching_cubes_cy.pyx
@@ -13,9 +13,9 @@ cdef inline double _get_fraction(double from_value, double to_value,
     return ((level - from_value) / (to_value - from_value))
 
 
-def unpack_unique_vert_list(list trilist):
+def unpack_unique_verts(list trilist):
     """
-    Convert_list a list of lists of tuples corresponding to triangle vertices
+    Convert a list of lists of tuples corresponding to triangle vertices
     into a unique vertex list, and a list of triangle faces w/indices
     corresponding to entries of the vertex list.
 
@@ -152,14 +152,14 @@ def iterate_and_store_3d(double[:, :, ::1] arr, double level,
 
         # We use a right-handed coordinate system, UNlike the paper, but want
         # to index in agreement - the coordinate adjustment takes place here.
-        v1 = arr[z0, y0, x0]
-        v2 = arr[z0, y0, x1]
-        v3 = arr[z0, y1, x1]
-        v4 = arr[z0, y1, x0]
-        v5 = arr[z1, y0, x0]
-        v6 = arr[z1, y0, x1]
-        v7 = arr[z1, y1, x1]
-        v8 = arr[z1, y1, x0]
+        v1 = arr[x0, y0, z0]
+        v2 = arr[x1, y0, z0]
+        v3 = arr[x1, y1, z0]
+        v4 = arr[x0, y1, z0]
+        v5 = arr[x0, y0, z1]
+        v6 = arr[x1, y0, z1]
+        v7 = arr[x1, y1, z1]
+        v8 = arr[x0, y1, z1]
 
         # Unique triangulation cases
         cube_case = 0
@@ -233,7 +233,7 @@ def iterate_and_store_3d(double[:, :, ::1] arr, double level,
                            e6, e7, e8, e9, e10, e11, e12)
 
         # Advance the coords indices
-        if coords[2] < arr.shape[0] - 2:
+        if coords[2] < arr.shape[2] - 2:
             coords[2] += 1
             plus_z = True
         elif coords[1] < arr.shape[1] - 2:

From b6f25906d652f60631ac769f551785e0582d78e4 Mon Sep 17 00:00:00 2001
From: "Josh Warner (Mac)" <warner.joshua@mayo.edu>
Date: Thu, 6 Jun 2013 15:56:29 -0500
Subject: [PATCH 15/24] DOCFIX: clarify docstrings and output vertex dimension
 ordering

---
 skimage/measure/_marching_cubes.py     | 6 +++---
 skimage/measure/_marching_cubes_cy.pyx | 2 +-
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/skimage/measure/_marching_cubes.py b/skimage/measure/_marching_cubes.py
index 418f945e..6af60672 100644
--- a/skimage/measure/_marching_cubes.py
+++ b/skimage/measure/_marching_cubes.py
@@ -9,8 +9,7 @@ def marching_cubes(volume, level, sampling=(1., 1., 1.)):
     Parameters
     ----------
     volume : (M, N, P) array of doubles
-        Input data volume to find isosurfaces. Will be cast to `np.float64` if
-        not provided in this format.
+        Input data volume to find isosurfaces. Will be cast to `np.float64`.
     level : float
         Contour value to search for isosurfaces in `volume`.
     sampling : length-3 tuple of floats
@@ -20,7 +19,8 @@ def marching_cubes(volume, level, sampling=(1., 1., 1.)):
     Returns
     -------
     verts : (V, 3) array
-        Spatial (x, y, z) coordinates for V unique vertices in mesh.
+        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
diff --git a/skimage/measure/_marching_cubes_cy.pyx b/skimage/measure/_marching_cubes_cy.pyx
index ec817b15..4343c1e5 100644
--- a/skimage/measure/_marching_cubes_cy.pyx
+++ b/skimage/measure/_marching_cubes_cy.pyx
@@ -60,7 +60,7 @@ def iterate_and_store_3d(double[:, :, ::1] arr, double level,
     """Iterate across the given array in a marching-cubes fashion,
     looking for volumes with edges that cross 'level'. If such a volume is
     found, appropriate triangulations are added to a growing list of
-    triangles to be returned by the function.
+    faces to be returned by this function.
 
     If `sampling` is not provided, vertices are returned in the indexing
     coordinate system (assuming all 3 spatial dimensions sampled equally).

From 997339beae952b0a3d63644546d02b257164a1a2 Mon Sep 17 00:00:00 2001
From: "Josh Warner (Mac)" <warner.joshua@mayo.edu>
Date: Tue, 18 Jun 2013 23:09:43 -0500
Subject: [PATCH 16/24] FEAT: add tests for marching cubes and mesh surface
 area

---
 skimage/measure/tests/test_marching_cubes.py | 141 +++++++++++++++++++
 1 file changed, 141 insertions(+)
 create mode 100644 skimage/measure/tests/test_marching_cubes.py

diff --git a/skimage/measure/tests/test_marching_cubes.py b/skimage/measure/tests/test_marching_cubes.py
new file mode 100644
index 00000000..e8135788
--- /dev/null
+++ b/skimage/measure/tests/test_marching_cubes.py
@@ -0,0 +1,141 @@
+import numpy as np
+from numpy.testing import assert_raises
+from scipy.special import (ellipkinc as ellip_F, ellipeinc as ellip_E)
+
+from skimage.measure import marching_cubes, mesh_surface_area
+
+
+def _ellipsoid(a, b, c, sampling=(1., 1., 1.), info=False, tight=False,
+               levelset=False):
+    """
+    Generates ellipsoid with semimajor axes aligned with grid dimensions,
+    on grid with specified `sampling`.
+
+    Parameters
+    ----------
+    a : float
+        Length of semimajor axis aligned with x-axis
+    b : float
+        Length of semimajor axis aligned with y-axis
+    c : float
+        Length of semimajor axis aligned with z-axis
+    sampling : tuple of floats, length 3
+        Sampling in each spatial dimension
+    info : bool
+        If False, only `bool_arr` returned.
+        If True, (`bool_arr`, `vol`, `surf`) returned; the additional
+        values are analytical volume and surface area calculated for
+        this ellipsoid.
+    tight : bool
+        Controls if the ellipsoid will precisely be contained within
+        the returned volume (tight=True) or if each dimension will be
+        2 longer than necessary (tight=False). For algorithms which
+        need both sides of a contour, use False.
+    levelset : bool
+        If True, returns the level set for this ellipsoid (signed level
+        set about zero, with positive denoting interior) as np.float64.
+        False returns a binarized version of said level set.
+
+    Returns
+    -------
+    bool_arr : (N, M, P) array
+        Sphere in an appropriately sized boolean array.
+    vol : float
+        Analytically calculated volume of ellipsoid. Only returned if
+        `info` is True.
+    surf : float
+        Analytically calculated surface area of ellipsoid. Only returned
+        if `info` is True.
+
+    """
+    if not tight:
+        offset = np.r_[1, 1, 1] * np.r_[sampling]
+    else:
+        offset = np.r_[0, 0, 0]
+
+    # Calculate limits, and ensure output volume is odd & symmetric
+    low = np.ceil((-np.r_[a, b, c] - offset))
+    high = np.floor((np.r_[a, b, c] + offset + 1))
+    for dim in range(3):
+        if (high[dim] - low[dim]) % 2 == 0:
+            low[dim] -= 1
+        num = np.arange(low[dim], high[dim], sampling[dim])
+        if 0 not in num:
+            low[dim] -= np.max(num[num < 0])
+
+    # Generate (anisotropic) spatial grid
+    x, y, z = np.mgrid[low[0]:high[0]:sampling[0],
+                       low[1]:high[1]:sampling[1],
+                       low[2]:high[2]:sampling[2]]
+
+    if not levelset:
+        arr = ((x / float(a)) ** 2 +
+               (y / float(b)) ** 2 +
+               (z / float(c)) ** 2) <= 1
+    else:
+        arr = ((x / float(a)) ** 2 +
+               (y / float(b)) ** 2 +
+               (z / float(c)) ** 2) - 1
+
+    if not info:
+        return arr
+    else:
+        # Surface calculation requires a >= b >= c and a != c.
+        abc = [a, b, c]
+        abc.sort(reverse=True)
+        a = abc[0]
+        b = abc[1]
+        c = abc[2]
+
+        # Volume
+        vol = 4 / 3. * np.pi * a * b * c
+
+        # Analytical ellipsoid surface area
+        phi = np.arcsin((1. - (c ** 2 / (a ** 2.))) ** 0.5)
+        d = float((a ** 2 - c ** 2) ** 0.5)
+        m = (a ** 2 * (b ** 2 - c ** 2) /
+             float(b ** 2 * (a ** 2 - c ** 2)))
+        F = ellip_F(phi, m)
+        E = ellip_E(phi, m)
+
+        surf = 2 * np.pi * (c ** 2 +
+                            b * c ** 2 / d * F +
+                            b * d * E)
+
+        return arr, vol, surf
+
+
+def test_marching_cubes_isotropic():
+    ellipsoid_isotropic, _, surf = _ellipsoid(6, 10, 16,
+                                              levelset=True,
+                                              info=True)
+    verts, faces = marching_cubes(ellipsoid_isotropic, 0.)
+    surf_calc = mesh_surface_area(verts, faces)
+
+    # Test within 1% tolerance for isotropic. Will always underestimate.
+    assert surf > surf_calc and surf_calc > surf * 0.99
+
+
+def test_marching_cubes_anisotropic():
+    sampling = (1., 10 / 6., 16 / 6.)
+    ellipsoid_isotropic, _, surf = _ellipsoid(6, 10, 16,
+                                              sampling=sampling,
+                                              levelset=True,
+                                              info=True)
+    verts, faces = marching_cubes(ellipsoid_isotropic, 0.,
+                                  sampling=sampling)
+    surf_calc = mesh_surface_area(verts, faces)
+    # Test within 1.5% tolerance for anisotropic. Will always underestimate.
+    assert surf > surf_calc and surf_calc > surf * 0.985
+
+
+def test_invalid_input():
+    assert_raises(ValueError, marching_cubes, np.zeros((2, 2, 1)), 0)
+    assert_raises(ValueError, marching_cubes, np.zeros((2, 2, 1)), 1)
+    assert_raises(ValueError, marching_cubes, np.ones((3, 3, 3)), 1,
+                  sampling=(1, 2))
+    assert_raises(ValueError, marching_cubes, np.zeros((20, 20)), 0)
+
+
+if __name__ == '__main__':
+    np.testing.run_module_suite()

From d0d9fee36e997d1578a424e49edacb8db202f863 Mon Sep 17 00:00:00 2001
From: "Josh Warner (Mac)" <warner.joshua@mayo.edu>
Date: Sun, 1 Sep 2013 20:49:58 -0500
Subject: [PATCH 17/24] DOC: Minor documentation formatting fixes in marching
 cubes

---
 skimage/measure/_marching_cubes.py     | 6 +++---
 skimage/measure/_marching_cubes_cy.pyx | 4 ++--
 2 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/skimage/measure/_marching_cubes.py b/skimage/measure/_marching_cubes.py
index 6af60672..772c52f1 100644
--- a/skimage/measure/_marching_cubes.py
+++ b/skimage/measure/_marching_cubes.py
@@ -41,8 +41,8 @@ def marching_cubes(volume, level, sampling=(1., 1., 1.)):
          /  |      /  |        ^  z
        v4 ------ v3   |        | /
         |  v5 ----|- v6        |/          (note: NOT right handed!)
-        | /       |  /          ----> x
-        |/        | /
+        |  /      |  /          ----> x
+        | /       | /
        v1 ------ v2
 
     Most notably, if v4, v8, v2, and v6 are all >= `level` (or any
@@ -153,5 +153,5 @@ def mesh_surface_area(verts, tris):
     b = actual_verts[:, 0, :] - actual_verts[:, 2, :]
     del actual_verts
 
-    # Area of triangle = 1/2 * Euclidean norm of cross product
+    # Area of triangle in 3D = 1/2 * Euclidean norm of cross product
     return ((np.cross(a, b) ** 2).sum(axis=1) ** 0.5).sum() / 2.
diff --git a/skimage/measure/_marching_cubes_cy.pyx b/skimage/measure/_marching_cubes_cy.pyx
index 4343c1e5..a0f63d1c 100644
--- a/skimage/measure/_marching_cubes_cy.pyx
+++ b/skimage/measure/_marching_cubes_cy.pyx
@@ -94,8 +94,8 @@ def iterate_and_store_3d(double[:, :, ::1] arr, double level,
     #         /  |      /  |        ^  z
     #       v4 ------ v3   |        | /
     #        |  v5 ----|- v6        |/          (note: NOT right handed!)
-    #        | /       |  /          ----> x
-    #        |/        | /
+    #        |  /      |  /          ----> x
+    #        | /       | /
     #       v1 ------ v2
     #
     # We also maintain the current 2D coordinates for v1, and ensure the array

From 14a0685838b53ce70d7674625b127e33c7a05625 Mon Sep 17 00:00:00 2001
From: "Josh Warner (Mac)" <warner.joshua@mayo.edu>
Date: Sun, 1 Sep 2013 20:51:20 -0500
Subject: [PATCH 18/24] MAINT: Refactor ellipsoid generator into skimage.draw

---
 skimage/draw/__init__.py                     |   3 +
 skimage/draw/draw3d.py                       | 117 +++++++++++++++++++
 skimage/measure/tests/test_marching_cubes.py | 117 ++-----------------
 3 files changed, 128 insertions(+), 109 deletions(-)
 create mode 100644 skimage/draw/draw3d.py

diff --git a/skimage/draw/__init__.py b/skimage/draw/__init__.py
index 38c114f7..4fb222d1 100644
--- a/skimage/draw/__init__.py
+++ b/skimage/draw/__init__.py
@@ -1,11 +1,14 @@
 from .draw import circle, ellipse, set_color
 from ._draw import line, polygon, ellipse_perimeter, circle_perimeter, \
                    bezier_segment
+from .draw3d import ellipsoid, ellipsoid_stats
 
 __all__ = ['line',
            'polygon',
            'ellipse',
            'ellipse_perimeter',
+           'ellipsoid',
+           'ellipsoid_stats',
            'circle',
            'circle_perimeter',
            'set_color']
diff --git a/skimage/draw/draw3d.py b/skimage/draw/draw3d.py
new file mode 100644
index 00000000..e247e360
--- /dev/null
+++ b/skimage/draw/draw3d.py
@@ -0,0 +1,117 @@
+# coding: utf-8
+import numpy as np
+from scipy.special import (ellipkinc as ellip_F, ellipeinc as ellip_E)
+
+
+def ellipsoid(a, b, c, sampling=(1., 1., 1.), levelset=False):
+    """
+    Generates ellipsoid with semimajor axes aligned with grid dimensions
+    on grid with specified `sampling`.
+
+    Parameters
+    ----------
+    a : float
+        Length of semimajor axis aligned with x-axis.
+    b : float
+        Length of semimajor axis aligned with y-axis.
+    c : float
+        Length of semimajor axis aligned with z-axis.
+    sampling : tuple of floats, length 3
+        Sampling in (x, y, z) spatial dimensions.
+    levelset : bool
+        If True, returns the level set for this ellipsoid (signed level
+        set about zero, with positive denoting interior) as np.float64.
+        False returns a binarized version of said level set.
+
+    Returns
+    -------
+    ellip : (N, M, P) array
+        Ellipsoid centered in a correctly sized array for given `sampling`.
+        Boolean dtype unless `levelset=True`, in which case a float array is
+        returned with the level set above 0.0 representing the ellipsoid.
+
+    """
+    if (a <= 0) or (b <= 0) or (c <= 0):
+        raise ValueError('Parameters a, b, and c must all be > 0')
+
+    offset = np.r_[1, 1, 1] * np.r_[sampling]
+
+    # Calculate limits, and ensure output volume is odd & symmetric
+    low = np.ceil((- np.r_[a, b, c] - offset))
+    high = np.floor((np.r_[a, b, c] + offset + 1))
+
+    for dim in range(3):
+        if (high[dim] - low[dim]) % 2 == 0:
+            low[dim] -= 1
+        num = np.arange(low[dim], high[dim], sampling[dim])
+        if 0 not in num:
+            low[dim] -= np.max(num[num < 0])
+
+    # Generate (anisotropic) spatial grid
+    x, y, z = np.mgrid[low[0]:high[0]:sampling[0],
+                       low[1]:high[1]:sampling[1],
+                       low[2]:high[2]:sampling[2]]
+
+    if not levelset:
+        arr = ((x / float(a)) ** 2 +
+               (y / float(b)) ** 2 +
+               (z / float(c)) ** 2) <= 1
+    else:
+        arr = ((x / float(a)) ** 2 +
+               (y / float(b)) ** 2 +
+               (z / float(c)) ** 2) - 1
+
+    return arr
+
+
+def ellipsoid_stats(a, b, c, sampling=(1., 1., 1.)):
+    """
+    Calculates analytical surface area and volume for ellipsoid with
+    semimajor axes aligned with grid dimensions of specified `sampling`.
+
+    Parameters
+    ----------
+    a : float
+        Length of semimajor axis aligned with x-axis.
+    b : float
+        Length of semimajor axis aligned with y-axis.
+    c : float
+        Length of semimajor axis aligned with z-axis.
+    sampling : tuple of floats, length 3
+        Sampling in (x, y, z) spatial dimensions.
+
+    Returns
+    -------
+    vol : float
+        Analytically calculated volume of ellipsoid.
+    surf : float
+        Analytically calculated surface area of ellipsoid.
+
+    """
+    if (a <= 0) or (b <= 0) or (c <= 0):
+        raise ValueError('Parameters a, b, and c must all be > 0')
+
+    # Calculate volume & surface area
+    # Surface calculation requires a >= b >= c and a != c.
+    abc = [a, b, c]
+    abc.sort(reverse=True)
+    a = abc[0]
+    b = abc[1]
+    c = abc[2]
+
+    # Volume
+    vol = 4 / 3. * np.pi * a * b * c
+
+    # Analytical ellipsoid surface area
+    phi = np.arcsin((1. - (c ** 2 / (a ** 2.))) ** 0.5)
+    d = float((a ** 2 - c ** 2) ** 0.5)
+    m = (a ** 2 * (b ** 2 - c ** 2) /
+         float(b ** 2 * (a ** 2 - c ** 2)))
+    F = ellip_F(phi, m)
+    E = ellip_E(phi, m)
+
+    surf = 2 * np.pi * (c ** 2 +
+                        b * c ** 2 / d * F +
+                        b * d * E)
+
+    return vol, surf
diff --git a/skimage/measure/tests/test_marching_cubes.py b/skimage/measure/tests/test_marching_cubes.py
index e8135788..991da060 100644
--- a/skimage/measure/tests/test_marching_cubes.py
+++ b/skimage/measure/tests/test_marching_cubes.py
@@ -1,114 +1,13 @@
 import numpy as np
 from numpy.testing import assert_raises
-from scipy.special import (ellipkinc as ellip_F, ellipeinc as ellip_E)
 
+from skimage.draw import ellipsoid, ellipsoid_stats
 from skimage.measure import marching_cubes, mesh_surface_area
 
 
-def _ellipsoid(a, b, c, sampling=(1., 1., 1.), info=False, tight=False,
-               levelset=False):
-    """
-    Generates ellipsoid with semimajor axes aligned with grid dimensions,
-    on grid with specified `sampling`.
-
-    Parameters
-    ----------
-    a : float
-        Length of semimajor axis aligned with x-axis
-    b : float
-        Length of semimajor axis aligned with y-axis
-    c : float
-        Length of semimajor axis aligned with z-axis
-    sampling : tuple of floats, length 3
-        Sampling in each spatial dimension
-    info : bool
-        If False, only `bool_arr` returned.
-        If True, (`bool_arr`, `vol`, `surf`) returned; the additional
-        values are analytical volume and surface area calculated for
-        this ellipsoid.
-    tight : bool
-        Controls if the ellipsoid will precisely be contained within
-        the returned volume (tight=True) or if each dimension will be
-        2 longer than necessary (tight=False). For algorithms which
-        need both sides of a contour, use False.
-    levelset : bool
-        If True, returns the level set for this ellipsoid (signed level
-        set about zero, with positive denoting interior) as np.float64.
-        False returns a binarized version of said level set.
-
-    Returns
-    -------
-    bool_arr : (N, M, P) array
-        Sphere in an appropriately sized boolean array.
-    vol : float
-        Analytically calculated volume of ellipsoid. Only returned if
-        `info` is True.
-    surf : float
-        Analytically calculated surface area of ellipsoid. Only returned
-        if `info` is True.
-
-    """
-    if not tight:
-        offset = np.r_[1, 1, 1] * np.r_[sampling]
-    else:
-        offset = np.r_[0, 0, 0]
-
-    # Calculate limits, and ensure output volume is odd & symmetric
-    low = np.ceil((-np.r_[a, b, c] - offset))
-    high = np.floor((np.r_[a, b, c] + offset + 1))
-    for dim in range(3):
-        if (high[dim] - low[dim]) % 2 == 0:
-            low[dim] -= 1
-        num = np.arange(low[dim], high[dim], sampling[dim])
-        if 0 not in num:
-            low[dim] -= np.max(num[num < 0])
-
-    # Generate (anisotropic) spatial grid
-    x, y, z = np.mgrid[low[0]:high[0]:sampling[0],
-                       low[1]:high[1]:sampling[1],
-                       low[2]:high[2]:sampling[2]]
-
-    if not levelset:
-        arr = ((x / float(a)) ** 2 +
-               (y / float(b)) ** 2 +
-               (z / float(c)) ** 2) <= 1
-    else:
-        arr = ((x / float(a)) ** 2 +
-               (y / float(b)) ** 2 +
-               (z / float(c)) ** 2) - 1
-
-    if not info:
-        return arr
-    else:
-        # Surface calculation requires a >= b >= c and a != c.
-        abc = [a, b, c]
-        abc.sort(reverse=True)
-        a = abc[0]
-        b = abc[1]
-        c = abc[2]
-
-        # Volume
-        vol = 4 / 3. * np.pi * a * b * c
-
-        # Analytical ellipsoid surface area
-        phi = np.arcsin((1. - (c ** 2 / (a ** 2.))) ** 0.5)
-        d = float((a ** 2 - c ** 2) ** 0.5)
-        m = (a ** 2 * (b ** 2 - c ** 2) /
-             float(b ** 2 * (a ** 2 - c ** 2)))
-        F = ellip_F(phi, m)
-        E = ellip_E(phi, m)
-
-        surf = 2 * np.pi * (c ** 2 +
-                            b * c ** 2 / d * F +
-                            b * d * E)
-
-        return arr, vol, surf
-
-
 def test_marching_cubes_isotropic():
-    ellipsoid_isotropic, _, surf = _ellipsoid(6, 10, 16,
-                                              levelset=True,
-                                              info=True)
+    ellipsoid_isotropic = ellipsoid(6, 10, 16, levelset=True)
+    _, surf = ellipsoid_stats(6, 10, 16, levelset=True)
     verts, faces = marching_cubes(ellipsoid_isotropic, 0.)
     surf_calc = mesh_surface_area(verts, faces)
 
@@ -118,13 +17,13 @@ def test_marching_cubes_isotropic():
 
 def test_marching_cubes_anisotropic():
     sampling = (1., 10 / 6., 16 / 6.)
-    ellipsoid_isotropic, _, surf = _ellipsoid(6, 10, 16,
-                                              sampling=sampling,
-                                              levelset=True,
-                                              info=True)
-    verts, faces = marching_cubes(ellipsoid_isotropic, 0.,
+    ellipsoid_anisotropic, _, surf = ellipsoid(6, 10, 16, sampling=sampling,
+                                               levelset=True)
+    _, surf = ellipsoid_stats(6, 10, 16, sampling=sampling, levelset=True)
+    verts, faces = marching_cubes(ellipsoid_anisotropic, 0.,
                                   sampling=sampling)
     surf_calc = mesh_surface_area(verts, faces)
+
     # Test within 1.5% tolerance for anisotropic. Will always underestimate.
     assert surf > surf_calc and surf_calc > surf * 0.985
 

From c678fa55f44872c470ddc7830775b0f54007325b Mon Sep 17 00:00:00 2001
From: "Josh Warner (Mac)" <warner.joshua@mayo.edu>
Date: Sun, 1 Sep 2013 20:52:06 -0500
Subject: [PATCH 19/24] DOC: Add concise marching cubes example to gallery
 using Matplotlib

---
 doc/examples/plot_marching_cubes.py | 56 +++++++++++++++++++++++++++++
 1 file changed, 56 insertions(+)
 create mode 100644 doc/examples/plot_marching_cubes.py

diff --git a/doc/examples/plot_marching_cubes.py b/doc/examples/plot_marching_cubes.py
new file mode 100644
index 00000000..0b321bfb
--- /dev/null
+++ b/doc/examples/plot_marching_cubes.py
@@ -0,0 +1,56 @@
+"""
+==============
+Marching Cubes
+==============
+
+Marching cubes is an algorithm to extract a 2D surface mesh from a 3D volume.
+This can be conceptualized as a 3D generalization of isolines on topographical
+or weather maps. It works by iterating across the volume, looking for regions
+which cross the level of interest. If such regions are found, triangulations
+are generated and added to an output mesh. The final result is a set of
+vertices and a set of triangular faces.
+
+The algorithm requires a data volume and an isosurface value. For example, in
+CT imaging Hounsfield units of +700 to +3000 represent bone. So, one potential
+input would be a reconstructed CT set of data and the value +700, to extract
+a mesh for regions of bone or bone-like density.
+
+This implementation also works correctly on anisotropic datasets, where the
+voxel spacing is not equal for every spatial dimension, through use of the
+`sampling` kwarg.
+
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+from mpl_toolkits.mplot3d import Axes3D
+import mpl_toolkits.mplot3d as a3
+
+from skimage import measure
+from skimage.draw import ellipsoid
+
+# Generate a level set about zero of two identical ellipsoids in 3D
+ellip_base = ellipsoid(6, 10, 16, levelset=True)
+ellip_double = np.concatenate((ellip_base[:-1, ...],
+                               ellip_base[2:, ...]), axis=0)
+
+# Use marching cubes to obtain the surface mesh of these ellipsoids
+verts, faces = measure.marching_cubes(ellip_double, 0)
+
+# Display resulting triangular mesh using Matplotlib. This can also be done
+# with mayavi (see skimage.measure.marching_cubes docstring).
+fig = plt.figure(figsize=(10, 12))
+ax = fig.add_subplot(111, projection='3d')
+
+# Fancy indexing: `verts[faces]` to generate a Poly3DCollection
+mesh = a3.art3d.Poly3DCollection(verts[faces])
+ax.add_collection3d(mesh)
+
+ax.set_xlabel("x-axis: a = 6 per ellipsoid")
+ax.set_ylabel("y-axis: b = 10")
+ax.set_zlabel("z-axis: c = 16")
+
+ax.set_xlim(0, 24)  # a = 6 (times two for 2nd ellipsoid)
+ax.set_ylim(0, 20)  # b = 10
+ax.set_zlim(0, 32)  # c = 16
+
+plt.show()

From 50d73fdd62503143f56b9c800b34a72d4e72b5ad Mon Sep 17 00:00:00 2001
From: "Josh Warner (Mac)" <warner.joshua@mayo.edu>
Date: Sun, 1 Sep 2013 21:31:53 -0500
Subject: [PATCH 20/24] TEST: Add test suite for draw3d
 (ellipsoid/ellipsoid_stats)

---
 skimage/draw/draw3d.py            |   4 +-
 skimage/draw/tests/test_draw3d.py | 104 ++++++++++++++++++++++++++++++
 2 files changed, 106 insertions(+), 2 deletions(-)
 create mode 100644 skimage/draw/tests/test_draw3d.py

diff --git a/skimage/draw/draw3d.py b/skimage/draw/draw3d.py
index e247e360..0b6fcb2d 100644
--- a/skimage/draw/draw3d.py
+++ b/skimage/draw/draw3d.py
@@ -83,9 +83,9 @@ def ellipsoid_stats(a, b, c, sampling=(1., 1., 1.)):
     Returns
     -------
     vol : float
-        Analytically calculated volume of ellipsoid.
+        Calculated volume of ellipsoid.
     surf : float
-        Analytically calculated surface area of ellipsoid.
+        Calculated surface area of ellipsoid.
 
     """
     if (a <= 0) or (b <= 0) or (c <= 0):
diff --git a/skimage/draw/tests/test_draw3d.py b/skimage/draw/tests/test_draw3d.py
new file mode 100644
index 00000000..26c63cde
--- /dev/null
+++ b/skimage/draw/tests/test_draw3d.py
@@ -0,0 +1,104 @@
+from numpy.testing import assert_array_equal
+import numpy as np
+
+from skimage.draw import ellipsoid, ellipsoid_stats
+
+
+def test_ellipsoid_bool():
+    test = ellipsoid(2, 2, 2)[1:-1, 1:-1, 1:-1]
+    test_anisotropic = ellipsoid(2, 2, 4, sampling=(1., 1., 2.))
+    test_anisotropic = test_anisotropic[1:-1, 1:-1, 1:-1]
+
+    expected = np.r_[[[[0, 0, 0, 0, 0],
+                       [0, 0, 0, 0, 0],
+                       [0, 0, 1, 0, 0],
+                       [0, 0, 0, 0, 0],
+                       [0, 0, 0, 0, 0]],
+
+                      [[0, 0, 0, 0, 0],
+                       [0, 1, 1, 1, 0],
+                       [0, 1, 1, 1, 0],
+                       [0, 1, 1, 1, 0],
+                       [0, 0, 0, 0, 0]],
+
+                      [[0, 0, 1, 0, 0],
+                       [0, 1, 1, 1, 0],
+                       [1, 1, 1, 1, 1],
+                       [0, 1, 1, 1, 0],
+                       [0, 0, 1, 0, 0]],
+
+                      [[0, 0, 0, 0, 0],
+                       [0, 1, 1, 1, 0],
+                       [0, 1, 1, 1, 0],
+                       [0, 1, 1, 1, 0],
+                       [0, 0, 0, 0, 0]],
+
+                      [[0, 0, 0, 0, 0],
+                       [0, 0, 0, 0, 0],
+                       [0, 0, 1, 0, 0],
+                       [0, 0, 0, 0, 0],
+                       [0, 0, 0, 0, 0]]]]
+
+    assert_array_equal(test, expected.astype(bool))
+    assert_array_equal(test_anisotropic, expected.astype(bool))
+
+
+def test_ellipsoid_levelset():
+    test = ellipsoid(2, 2, 2, levelset=True)[1:-1, 1:-1, 1:-1]
+    test_anisotropic = ellipsoid(2, 2, 4, sampling=(1., 1., 2.),
+                                 levelset=True)
+    test_anisotropic = test_anisotropic[1:-1, 1:-1, 1:-1]
+
+    expected = np.r_[[[[ 2.  ,  1.25,  1.  ,  1.25,  2.  ],
+                       [ 1.25,  0.5 ,  0.25,  0.5 ,  1.25],
+                       [ 1.  ,  0.25,  0.  ,  0.25,  1.  ],
+                       [ 1.25,  0.5 ,  0.25,  0.5 ,  1.25],
+                       [ 2.  ,  1.25,  1.  ,  1.25,  2.  ]],
+
+                      [[ 1.25,  0.5 ,  0.25,  0.5 ,  1.25],
+                       [ 0.5 , -0.25, -0.5 , -0.25,  0.5 ],
+                       [ 0.25, -0.5 , -0.75, -0.5 ,  0.25],
+                       [ 0.5 , -0.25, -0.5 , -0.25,  0.5 ],
+                       [ 1.25,  0.5 ,  0.25,  0.5 ,  1.25]],
+
+                      [[ 1.  ,  0.25,  0.  ,  0.25,  1.  ],
+                       [ 0.25, -0.5 , -0.75, -0.5 ,  0.25],
+                       [ 0.  , -0.75, -1.  , -0.75,  0.  ],
+                       [ 0.25, -0.5 , -0.75, -0.5 ,  0.25],
+                       [ 1.  ,  0.25,  0.  ,  0.25,  1.  ]],
+
+                      [[ 1.25,  0.5 ,  0.25,  0.5 ,  1.25],
+                       [ 0.5 , -0.25, -0.5 , -0.25,  0.5 ],
+                       [ 0.25, -0.5 , -0.75, -0.5 ,  0.25],
+                       [ 0.5 , -0.25, -0.5 , -0.25,  0.5 ],
+                       [ 1.25,  0.5 ,  0.25,  0.5 ,  1.25]],
+
+                      [[ 2.  ,  1.25,  1.  ,  1.25,  2.  ],
+                       [ 1.25,  0.5 ,  0.25,  0.5 ,  1.25],
+                       [ 1.  ,  0.25,  0.  ,  0.25,  1.  ],
+                       [ 1.25,  0.5 ,  0.25,  0.5 ,  1.25],
+                       [ 2.  ,  1.25,  1.  ,  1.25,  2.  ]]]]
+
+    assert_array_equal(test, expected.astype(bool))
+    assert_array_equal(test_anisotropic, expected.astype(bool))
+
+
+def test_ellipsoid_stats():
+    # Test comparison values generated by Wolfram Alpha
+    vol, surf = ellipsoid_stats(6, 10, 16)
+    assert(round(1280 * np.pi, 4) == round(vol, 4))
+    assert(1383.28 == round(surf, 2))
+
+    # Test when a <= b <= c does not hold
+    vol, surf = ellipsoid_stats(16, 6, 10)
+    assert(round(1280 * np.pi, 4) == round(vol, 4))
+    assert(1383.28 == round(surf, 2))
+
+    # Larger test to ensure reliability over broad range
+    vol, surf = ellipsoid_stats(17, 27, 169)
+    assert(round(103428 * np.pi, 4) == round(vol, 4))
+    assert(37426.3 == round(surf, 1))
+
+
+if __name__ == "__main__":
+    np.testing.run_module_suite()

From cd306bad5704ac66c3ec8594ef30e93f11e4af3e Mon Sep 17 00:00:00 2001
From: "Josh Warner (Mac)" <warner.joshua@mayo.edu>
Date: Sun, 1 Sep 2013 21:41:50 -0500
Subject: [PATCH 21/24] FIX: Add _marching_cubes_cy.pyx to bento.info

---
 bento.info | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/bento.info b/bento.info
index acdb155d..b1169244 100644
--- a/bento.info
+++ b/bento.info
@@ -51,6 +51,9 @@ Library:
     Extension: skimage.measure._moments
         Sources:
             skimage/measure/_moments.pyx
+    Extension: skimage.measure._marching_cubes_cy
+        Sources:
+            skimage/measure/_marching_cubes_cy.pyx
     Extension: skimage.graph._mcp
         Sources:
             skimage/graph/_mcp.pyx

From 1174d6fbc5ec7b9799e243ca20540d8effa39831 Mon Sep 17 00:00:00 2001
From: "Josh Warner (Mac)" <warner.joshua@mayo.edu>
Date: Sun, 1 Sep 2013 22:05:39 -0500
Subject: [PATCH 22/24] FIX: Errors in unit tests

---
 skimage/draw/tests/test_draw3d.py            | 4 ++--
 skimage/measure/tests/test_marching_cubes.py | 8 ++++----
 2 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/skimage/draw/tests/test_draw3d.py b/skimage/draw/tests/test_draw3d.py
index 26c63cde..14bf140f 100644
--- a/skimage/draw/tests/test_draw3d.py
+++ b/skimage/draw/tests/test_draw3d.py
@@ -79,8 +79,8 @@ def test_ellipsoid_levelset():
                        [ 1.25,  0.5 ,  0.25,  0.5 ,  1.25],
                        [ 2.  ,  1.25,  1.  ,  1.25,  2.  ]]]]
 
-    assert_array_equal(test, expected.astype(bool))
-    assert_array_equal(test_anisotropic, expected.astype(bool))
+    assert_array_equal(test, expected)
+    assert_array_equal(test_anisotropic, expected)
 
 
 def test_ellipsoid_stats():
diff --git a/skimage/measure/tests/test_marching_cubes.py b/skimage/measure/tests/test_marching_cubes.py
index 991da060..7a1fd40a 100644
--- a/skimage/measure/tests/test_marching_cubes.py
+++ b/skimage/measure/tests/test_marching_cubes.py
@@ -7,7 +7,7 @@ from skimage.measure import marching_cubes, mesh_surface_area
 
 def test_marching_cubes_isotropic():
     ellipsoid_isotropic = ellipsoid(6, 10, 16, levelset=True)
-    _, surf = ellipsoid_stats(6, 10, 16, levelset=True)
+    _, surf = ellipsoid_stats(6, 10, 16)
     verts, faces = marching_cubes(ellipsoid_isotropic, 0.)
     surf_calc = mesh_surface_area(verts, faces)
 
@@ -17,9 +17,9 @@ def test_marching_cubes_isotropic():
 
 def test_marching_cubes_anisotropic():
     sampling = (1., 10 / 6., 16 / 6.)
-    ellipsoid_anisotropic, _, surf = ellipsoid(6, 10, 16, sampling=sampling,
-                                               levelset=True)
-    _, surf = ellipsoid_stats(6, 10, 16, sampling=sampling, levelset=True)
+    ellipsoid_anisotropic = ellipsoid(6, 10, 16, sampling=sampling,
+                                      levelset=True)
+    _, surf = ellipsoid_stats(6, 10, 16, sampling=sampling)
     verts, faces = marching_cubes(ellipsoid_anisotropic, 0.,
                                   sampling=sampling)
     surf_calc = mesh_surface_area(verts, faces)

From abc7a16a8099ce0110b02b2dd934ef4873048943 Mon Sep 17 00:00:00 2001
From: "Josh Warner (Mac)" <warner.joshua@mayo.edu>
Date: Mon, 2 Sep 2013 15:31:51 -0500
Subject: [PATCH 23/24] STYLE: Use `np.array` instead of `np.r_` and import
 Poly3DCollection

---
 doc/examples/plot_marching_cubes.py |   6 +-
 skimage/draw/tests/test_draw3d.py   | 106 ++++++++++++++--------------
 2 files changed, 56 insertions(+), 56 deletions(-)

diff --git a/doc/examples/plot_marching_cubes.py b/doc/examples/plot_marching_cubes.py
index 0b321bfb..a57a40a2 100644
--- a/doc/examples/plot_marching_cubes.py
+++ b/doc/examples/plot_marching_cubes.py
@@ -23,7 +23,7 @@ voxel spacing is not equal for every spatial dimension, through use of the
 import numpy as np
 import matplotlib.pyplot as plt
 from mpl_toolkits.mplot3d import Axes3D
-import mpl_toolkits.mplot3d as a3
+from mpl_toolkits.mplot3d.art3d import Poly3DCollection
 
 from skimage import measure
 from skimage.draw import ellipsoid
@@ -41,8 +41,8 @@ verts, faces = measure.marching_cubes(ellip_double, 0)
 fig = plt.figure(figsize=(10, 12))
 ax = fig.add_subplot(111, projection='3d')
 
-# Fancy indexing: `verts[faces]` to generate a Poly3DCollection
-mesh = a3.art3d.Poly3DCollection(verts[faces])
+# Fancy indexing: `verts[faces]` to generate a collection of triangles
+mesh = Poly3DCollection(verts[faces])
 ax.add_collection3d(mesh)
 
 ax.set_xlabel("x-axis: a = 6 per ellipsoid")
diff --git a/skimage/draw/tests/test_draw3d.py b/skimage/draw/tests/test_draw3d.py
index 14bf140f..59a7f6b3 100644
--- a/skimage/draw/tests/test_draw3d.py
+++ b/skimage/draw/tests/test_draw3d.py
@@ -1,5 +1,5 @@
-from numpy.testing import assert_array_equal
 import numpy as np
+from numpy.testing import assert_array_equal, assert_allclose
 
 from skimage.draw import ellipsoid, ellipsoid_stats
 
@@ -9,35 +9,35 @@ def test_ellipsoid_bool():
     test_anisotropic = ellipsoid(2, 2, 4, sampling=(1., 1., 2.))
     test_anisotropic = test_anisotropic[1:-1, 1:-1, 1:-1]
 
-    expected = np.r_[[[[0, 0, 0, 0, 0],
-                       [0, 0, 0, 0, 0],
-                       [0, 0, 1, 0, 0],
-                       [0, 0, 0, 0, 0],
-                       [0, 0, 0, 0, 0]],
+    expected = np.array([[[0, 0, 0, 0, 0],
+                          [0, 0, 0, 0, 0],
+                          [0, 0, 1, 0, 0],
+                          [0, 0, 0, 0, 0],
+                          [0, 0, 0, 0, 0]],
 
-                      [[0, 0, 0, 0, 0],
-                       [0, 1, 1, 1, 0],
-                       [0, 1, 1, 1, 0],
-                       [0, 1, 1, 1, 0],
-                       [0, 0, 0, 0, 0]],
+                         [[0, 0, 0, 0, 0],
+                          [0, 1, 1, 1, 0],
+                          [0, 1, 1, 1, 0],
+                          [0, 1, 1, 1, 0],
+                          [0, 0, 0, 0, 0]],
 
-                      [[0, 0, 1, 0, 0],
-                       [0, 1, 1, 1, 0],
-                       [1, 1, 1, 1, 1],
-                       [0, 1, 1, 1, 0],
-                       [0, 0, 1, 0, 0]],
+                         [[0, 0, 1, 0, 0],
+                          [0, 1, 1, 1, 0],
+                          [1, 1, 1, 1, 1],
+                          [0, 1, 1, 1, 0],
+                          [0, 0, 1, 0, 0]],
 
-                      [[0, 0, 0, 0, 0],
-                       [0, 1, 1, 1, 0],
-                       [0, 1, 1, 1, 0],
-                       [0, 1, 1, 1, 0],
-                       [0, 0, 0, 0, 0]],
+                         [[0, 0, 0, 0, 0],
+                          [0, 1, 1, 1, 0],
+                          [0, 1, 1, 1, 0],
+                          [0, 1, 1, 1, 0],
+                          [0, 0, 0, 0, 0]],
 
-                      [[0, 0, 0, 0, 0],
-                       [0, 0, 0, 0, 0],
-                       [0, 0, 1, 0, 0],
-                       [0, 0, 0, 0, 0],
-                       [0, 0, 0, 0, 0]]]]
+                         [[0, 0, 0, 0, 0],
+                          [0, 0, 0, 0, 0],
+                          [0, 0, 1, 0, 0],
+                          [0, 0, 0, 0, 0],
+                          [0, 0, 0, 0, 0]]])
 
     assert_array_equal(test, expected.astype(bool))
     assert_array_equal(test_anisotropic, expected.astype(bool))
@@ -49,38 +49,38 @@ def test_ellipsoid_levelset():
                                  levelset=True)
     test_anisotropic = test_anisotropic[1:-1, 1:-1, 1:-1]
 
-    expected = np.r_[[[[ 2.  ,  1.25,  1.  ,  1.25,  2.  ],
-                       [ 1.25,  0.5 ,  0.25,  0.5 ,  1.25],
-                       [ 1.  ,  0.25,  0.  ,  0.25,  1.  ],
-                       [ 1.25,  0.5 ,  0.25,  0.5 ,  1.25],
-                       [ 2.  ,  1.25,  1.  ,  1.25,  2.  ]],
+    expected = np.array([[[ 2.  ,  1.25,  1.  ,  1.25,  2.  ],
+                          [ 1.25,  0.5 ,  0.25,  0.5 ,  1.25],
+                          [ 1.  ,  0.25,  0.  ,  0.25,  1.  ],
+                          [ 1.25,  0.5 ,  0.25,  0.5 ,  1.25],
+                          [ 2.  ,  1.25,  1.  ,  1.25,  2.  ]],
 
-                      [[ 1.25,  0.5 ,  0.25,  0.5 ,  1.25],
-                       [ 0.5 , -0.25, -0.5 , -0.25,  0.5 ],
-                       [ 0.25, -0.5 , -0.75, -0.5 ,  0.25],
-                       [ 0.5 , -0.25, -0.5 , -0.25,  0.5 ],
-                       [ 1.25,  0.5 ,  0.25,  0.5 ,  1.25]],
+                         [[ 1.25,  0.5 ,  0.25,  0.5 ,  1.25],
+                          [ 0.5 , -0.25, -0.5 , -0.25,  0.5 ],
+                          [ 0.25, -0.5 , -0.75, -0.5 ,  0.25],
+                          [ 0.5 , -0.25, -0.5 , -0.25,  0.5 ],
+                          [ 1.25,  0.5 ,  0.25,  0.5 ,  1.25]],
 
-                      [[ 1.  ,  0.25,  0.  ,  0.25,  1.  ],
-                       [ 0.25, -0.5 , -0.75, -0.5 ,  0.25],
-                       [ 0.  , -0.75, -1.  , -0.75,  0.  ],
-                       [ 0.25, -0.5 , -0.75, -0.5 ,  0.25],
-                       [ 1.  ,  0.25,  0.  ,  0.25,  1.  ]],
+                         [[ 1.  ,  0.25,  0.  ,  0.25,  1.  ],
+                          [ 0.25, -0.5 , -0.75, -0.5 ,  0.25],
+                          [ 0.  , -0.75, -1.  , -0.75,  0.  ],
+                          [ 0.25, -0.5 , -0.75, -0.5 ,  0.25],
+                          [ 1.  ,  0.25,  0.  ,  0.25,  1.  ]],
 
-                      [[ 1.25,  0.5 ,  0.25,  0.5 ,  1.25],
-                       [ 0.5 , -0.25, -0.5 , -0.25,  0.5 ],
-                       [ 0.25, -0.5 , -0.75, -0.5 ,  0.25],
-                       [ 0.5 , -0.25, -0.5 , -0.25,  0.5 ],
-                       [ 1.25,  0.5 ,  0.25,  0.5 ,  1.25]],
+                         [[ 1.25,  0.5 ,  0.25,  0.5 ,  1.25],
+                          [ 0.5 , -0.25, -0.5 , -0.25,  0.5 ],
+                          [ 0.25, -0.5 , -0.75, -0.5 ,  0.25],
+                          [ 0.5 , -0.25, -0.5 , -0.25,  0.5 ],
+                          [ 1.25,  0.5 ,  0.25,  0.5 ,  1.25]],
 
-                      [[ 2.  ,  1.25,  1.  ,  1.25,  2.  ],
-                       [ 1.25,  0.5 ,  0.25,  0.5 ,  1.25],
-                       [ 1.  ,  0.25,  0.  ,  0.25,  1.  ],
-                       [ 1.25,  0.5 ,  0.25,  0.5 ,  1.25],
-                       [ 2.  ,  1.25,  1.  ,  1.25,  2.  ]]]]
+                         [[ 2.  ,  1.25,  1.  ,  1.25,  2.  ],
+                          [ 1.25,  0.5 ,  0.25,  0.5 ,  1.25],
+                          [ 1.  ,  0.25,  0.  ,  0.25,  1.  ],
+                          [ 1.25,  0.5 ,  0.25,  0.5 ,  1.25],
+                          [ 2.  ,  1.25,  1.  ,  1.25,  2.  ]]])
 
-    assert_array_equal(test, expected)
-    assert_array_equal(test_anisotropic, expected)
+    assert_allclose(test, expected)
+    assert_allclose(test_anisotropic, expected)
 
 
 def test_ellipsoid_stats():

From d2bc02c10bf4446f5b8a8d678641ea8330800dd4 Mon Sep 17 00:00:00 2001
From: "Josh Warner (Mac)" <warner.joshua@mayo.edu>
Date: Mon, 2 Sep 2013 20:23:00 -0500
Subject: [PATCH 24/24] Credit myself for Marching Cubes

---
 CONTRIBUTORS.txt | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt
index b304646b..28203676 100644
--- a/CONTRIBUTORS.txt
+++ b/CONTRIBUTORS.txt
@@ -114,7 +114,7 @@
 
 - Joshua Warner
   Multichannel random walker segmentation, unified peak finder backend,
-  n-dimensional array padding, bug and doc fixes.
+  n-dimensional array padding, marching cubes, bug and doc fixes.
 
 - Petter Strandmark
   Perimeter calculation in regionprops.