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. 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 diff --git a/doc/examples/plot_marching_cubes.py b/doc/examples/plot_marching_cubes.py new file mode 100644 index 00000000..a57a40a2 --- /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 +from mpl_toolkits.mplot3d.art3d import Poly3DCollection + +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 collection of triangles +mesh = 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() 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..0b6fcb2d --- /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 + Calculated volume of ellipsoid. + surf : float + 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/draw/tests/test_draw3d.py b/skimage/draw/tests/test_draw3d.py new file mode 100644 index 00000000..59a7f6b3 --- /dev/null +++ b/skimage/draw/tests/test_draw3d.py @@ -0,0 +1,104 @@ +import numpy as np +from numpy.testing import assert_array_equal, assert_allclose + +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.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, 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.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. , 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_allclose(test, expected) + assert_allclose(test_anisotropic, expected) + + +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() 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'] diff --git a/skimage/measure/_marching_cubes.py b/skimage/measure/_marching_cubes.py new file mode 100644 index 00000000..772c52f1 --- /dev/null +++ b/skimage/measure/_marching_cubes.py @@ -0,0 +1,157 @@ +import numpy as np +from . import _marching_cubes_cy + + +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`. + 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 + ------- + verts : (V, 3) array + Spatial coordinates for V unique mesh vertices. Coordinate order + matches input `volume` (M, N, P). + faces : (F, 3) array + Define triangular faces via referencing vertex indices from ``verts``. + This algorithm specifically outputs triangles, so each face has + exactly three indices. + + 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 outputs directly into `skimage.measure.mesh_surface_area`. + + Regarding visualization of algorithm output, the ``mayavi`` package + is recommended. To contour a volume named `myvolume` about the level 0.0: + + >>> from mayavi import mlab + >>> verts, 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 and ensure `volume` is C-contiguous for memoryviews + if volume.ndim != 3: + 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 + # 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_cy.iterate_and_store_3d(volume, float(level), + sampling) + + # Find and collect unique vertices, storing triangle verts as indices. + # Returns a true mesh with no degenerate faces. + verts, faces = _marching_cubes_cy.unpack_unique_verts(raw_tris) + + return np.asarray(verts), np.asarray(faces) + + +def mesh_surface_area(verts, tris): + """ + Compute surface area, given vertices & triangular faces + + Parameters + ---------- + 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 now [coordinate units] ** 2. + + 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`. + + This algorithm works properly only if the ``faces`` provided are all + triangles. + + See Also + -------- + skimage.measure.marching_cubes + + """ + # 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 + + # 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 new file mode 100644 index 00000000..a0f63d1c --- /dev/null +++ b/skimage/measure/_marching_cubes_cy.pyx @@ -0,0 +1,987 @@ +#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 unpack_unique_verts(list trilist): + """ + 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. + + """ + 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 face_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]) + + face_list.append(templist) + + return vert_list, face_list + + +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 + 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). + 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 (double, double, double)") + + cdef list face_list = [] + cdef list norm_list = [] + cdef Py_ssize_t n + cdef bint odd_sampling, plus_z + plus_z = False + if [float(i) for i in sampling] == [1.0, 1.0, 1.0]: + 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: + 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 + # 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 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. + 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 `face_list` + _append_tris(face_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 face_list + + +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 + # quite long enough as-is. + + if (case == 1): + # front lower left corner + face_list.append([e1, e4, e9]) + elif (case == 2): + # front lower right corner + face_list.append([e10, e2, e1]) + elif (case == 3): + # front lower plane + face_list.append([e2, e4, e9]) + face_list.append([e2, e9, e10]) + elif (case == 4): + # front upper right corner + face_list.append([e12, e3, e2]) + elif (case == 5): + # lower left, upper right corners + _append_tris(face_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _append_tris(face_list, 4, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + elif (case == 6): + # front right plane + face_list.append([e12, e3, e1]) + face_list.append([e12, e1, e10]) + elif (case == 7): + # Shelf including v1, v2, v3 + 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 + face_list.append([e3, e11, e4]) + elif (case == 9): + # front left plane + face_list.append([e3, e11, e9]) + face_list.append([e3, e9, e1]) + elif (case == 10): + # upper left, lower right corners + _append_tris(face_list, 2, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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 + face_list.append([e3, e11, e2]) + face_list.append([e11, e10, e2]) + face_list.append([e11, e9, e10]) + elif (case == 12): + # front upper plane + face_list.append([e11, e4, e12]) + face_list.append([e2, e4, e12]) + elif (case == 13): + # Shelf including v1, v4, v3 + 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 + 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 + face_list.append([e11, e9, e12]) + face_list.append([e12, e9, e10]) + elif (case == 16): + # back lower left corner + face_list.append([e8, e9, e5]) + elif (case == 17): + # lower left plane + face_list.append([e4, e1, e8]) + face_list.append([e8, e1, e5]) + elif (case == 18): + # lower left back, lower right front corners + _append_tris(face_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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 + 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(face_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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(face_list, 17, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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(face_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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 + 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(face_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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 + 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(face_list, 2, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _append_tris(face_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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) + 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(face_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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) + 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(face_list, 14, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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 + 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 + face_list.append([e6, e5, e10]) + elif (case == 33): + # lower right back, lower left front corners + _append_tris(face_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _append_tris(face_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + elif (case == 34): + # lower right plane + face_list.append([e1, e2, e5]) + face_list.append([e2, e6, e5]) + elif (case == 35): + # Shelf: v1, v2, v6 + 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(face_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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(face_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _append_tris(face_list, 4, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _append_tris(face_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + elif (case == 38): + # Shelf: v2, v3, v6 + 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) + 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(face_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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(face_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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(face_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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 + 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(face_list, 12, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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(face_list, 13, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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 + 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 + face_list.append([e11, e9, e12]) + face_list.append([e12, e9, e5]) + face_list.append([e12, e5, e6]) + elif (case == 48): + # Back lower plane + face_list.append([e9, e10, e6]) + face_list.append([e9, e6, e8]) + elif (case == 49): + # Shelf: (v1, v5, v6) + 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) + 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 + 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(face_list, 48, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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(face_list, 49, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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) + 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 + 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(face_list, 48, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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) + 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(face_list, 50, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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 + 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(face_list, 48, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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(face_list, 63, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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(face_list, 63, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _append_tris(face_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + elif (case == 63): + # Upper back plane + face_list.append([e11, e12, e6]) + face_list.append([e11, e8, e6]) + elif (case == 64): + # Upper right back corner + face_list.append([e12, e7, e6]) + elif (case == 65): + # upper right back, lower left front corners + _append_tris(face_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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(face_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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(face_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _append_tris(face_list, 3, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + elif (case == 68): + # upper right plane + face_list.append([e3, e2, e6]) + face_list.append([e3, e7, e6]) + elif (case == 69): + # Upper right plane, lower left front corner + _append_tris(face_list, 68, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _append_tris(face_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + elif (case == 70): + # Shelf: (v2, v3, v7) + 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) + 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(face_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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(face_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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(face_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _append_tris(face_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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(face_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _append_tris(face_list, 11, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + elif (case == 76): + # Shelf: (v4, v3, v7) + 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) + 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) + 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 + 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(face_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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(face_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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(face_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _append_tris(face_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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(face_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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(face_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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(face_list, 17, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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(face_list, 70, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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(face_list, 119, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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(face_list, 64, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _append_tris(face_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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(face_list, 25, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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(face_list, 2, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _append_tris(face_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _append_tris(face_list, 16, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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(face_list, 4, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _append_tris(face_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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(face_list, 76, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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(face_list, 127, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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(face_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _append_tris(face_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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(face_list, 32, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _append_tris(face_list, 127, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + elif (case == 96): + # back right plane + face_list.append([e7, e12, e5]) + face_list.append([e5, e10, e12]) + elif (case == 97): + # back right plane + isolated corner v1 + _append_tris(face_list, 96, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _append_tris(face_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + elif (case == 98): + # Shelf: (v2, v6, v7) + 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) + 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) + 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(face_list, 100, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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 + face_list.append([e1, e3, e7]) + face_list.append([e1, e7, e5]) + elif (case == 103): + # Shelf: (v4, v5, v8) by inversion + 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(face_list, 8, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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(face_list, 96, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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(face_list, 98, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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(face_list, 4, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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) + 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(face_list, 111, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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 + face_list.append([e1, e5, e7]) + face_list.append([e1, e7, e11]) + face_list.append([e1, e11, e4]) + elif (case == 111): + # Back left plane + face_list.append([e11, e9, e7]) + face_list.append([e9, e7, e5]) + elif (case == 112): + # Shelf: (v5, v6, v7) + 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) + 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) + 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) + 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) + 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(face_list, 2, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _append_tris(face_list, 119, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + elif (case == 118): + # Shelf: (v1, v4, v8) + face_list.append([e1, e3, e7]) + face_list.append([e7, e1, e8]) + face_list.append([e1, e8, e9]) + elif (case == 119): + # Upper left plane + face_list.append([e4, e3, e7]) + face_list.append([e4, e8, e7]) + elif (case == 120): + # Shelf: (v5, v6, v7) + isolated corner v4 + _append_tris(face_list, 112, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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(face_list, 6, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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(face_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _append_tris(face_list, 4, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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(face_list, 4, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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(face_list, 3, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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(face_list, 2, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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(face_list, 1, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _append_tris(face_list, 127, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + elif (case == 127): + # Isolated corner v8 + 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(face_list, 6, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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(face_list, 119, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _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(face_list, 63, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + _append_tris(face_list, 3, e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, + e11, e12) + + return diff --git a/skimage/measure/setup.py b/skimage/measure/setup.py index 21d9964e..be57ca7b 100644 --- a/skimage/measure/setup.py +++ b/skimage/measure/setup.py @@ -14,11 +14,15 @@ 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_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_cy', + sources=['_marching_cubes_cy.c'], + include_dirs=[get_numpy_include_dirs()]) return config diff --git a/skimage/measure/tests/test_marching_cubes.py b/skimage/measure/tests/test_marching_cubes.py new file mode 100644 index 00000000..7a1fd40a --- /dev/null +++ b/skimage/measure/tests/test_marching_cubes.py @@ -0,0 +1,40 @@ +import numpy as np +from numpy.testing import assert_raises + +from skimage.draw import ellipsoid, ellipsoid_stats +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) + 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_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) + + # 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()