From af59d06461b999bcf9bd881eabd0da5fd5686d98 Mon Sep 17 00:00:00 2001 From: "Josh Warner (Mac)" Date: Sun, 31 Aug 2014 10:57:51 +0100 Subject: [PATCH] ENH: Simplify marching_cubes for maintainability. This does not change execution time, but the entire algorithm is much simpler to read and understand without the branching Cython code paths for anisotropic inputs. Memory usage decreases a minor amount. --- skimage/measure/_marching_cubes.py | 9 ++- skimage/measure/_marching_cubes_cy.pyx | 78 ++++++-------------------- 2 files changed, 22 insertions(+), 65 deletions(-) diff --git a/skimage/measure/_marching_cubes.py b/skimage/measure/_marching_cubes.py index 8f9dfcbc..a334fbbd 100644 --- a/skimage/measure/_marching_cubes.py +++ b/skimage/measure/_marching_cubes.py @@ -105,6 +105,9 @@ def marching_cubes(volume, level, spacing=(1., 1., 1.)): 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.") + if len(spacing) != 3: + raise ValueError("`spacing` must consist of three floats.") + volume = np.array(volume, dtype=np.float64, order="C") # Extract raw triangles using marching cubes in Cython @@ -113,14 +116,14 @@ def marching_cubes(volume, level, spacing=(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_faces = _marching_cubes_cy.iterate_and_store_3d(volume, float(level), - spacing) + raw_faces = _marching_cubes_cy.iterate_and_store_3d(volume, float(level)) # 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_faces) - return np.asarray(verts), np.asarray(faces) + # Adjust for non-isotropic spacing in `verts` at time of return + return np.asarray(verts) * np.r_[spacing], np.asarray(faces) def mesh_surface_area(verts, faces): diff --git a/skimage/measure/_marching_cubes_cy.pyx b/skimage/measure/_marching_cubes_cy.pyx index 085108ab..096d502d 100644 --- a/skimage/measure/_marching_cubes_cy.pyx +++ b/skimage/measure/_marching_cubes_cy.pyx @@ -55,33 +55,21 @@ def unpack_unique_verts(list trilist): return vert_list, face_list -def iterate_and_store_3d(double[:, :, ::1] arr, double level, - tuple spacing=(1., 1., 1.)): +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 faces to be returned by this function. - If `spacing` is not provided, vertices are returned in the indexing - coordinate system (assuming all 3 spatial dimensions sampled equally). - If `spacing` 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(spacing) != 3: - raise ValueError("`spacing` must be (double, double, double)") cdef list face_list = [] cdef list norm_list = [] cdef Py_ssize_t n - cdef bint odd_spacing, plus_z + cdef bint plus_z plus_z = False - if [float(i) for i in spacing] == [1.0, 1.0, 1.0]: - odd_spacing = False - else: - odd_spacing = 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 @@ -107,12 +95,6 @@ def iterate_and_store_3d(double[:, :, ::1] arr, double level, coords[1] = 0 coords[2] = 0 - # Extract doubles from `spacing` for speed - cdef double[3] spacing2 - spacing2[0] = spacing[0] - spacing2[1] = spacing[1] - spacing2[2] = spacing[2] - # Calculate the number of iterations we'll need cdef Py_ssize_t num_cube_steps = ((arr.shape[0] - 1) * (arr.shape[1] - 1) * @@ -120,7 +102,7 @@ def iterate_and_store_3d(double[:, :, ::1] arr, double level, 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 double v1, v2, v3, v4, v5, v6, v7, v8 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) @@ -138,18 +120,6 @@ def iterate_and_store_3d(double[:, :, ::1] arr, double level, x0, y0, z0 = coords[0], coords[1], coords[2] x1, y1, z1 = x0 + 1, y0 + 1, z0 + 1 - if odd_spacing: - # These doubles are the modified world coordinates; they are only - # calculated if non-default `spacing` provided. - r0 = coords[0] * spacing2[0] - c0 = coords[1] * spacing2[1] - d0 = coords[2] * spacing2[2] - r1 = r0 + spacing2[0] - c1 = c0 + spacing2[1] - d1 = d0 + spacing2[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] @@ -192,40 +162,24 @@ def iterate_and_store_3d(double[:, :, ::1] arr, double level, e3 = e7 e4 = e8 else: - # Calculate edges normally - if odd_spacing: - e1 = r0 + _get_fraction(v1, v2, level) * spacing2[0], c0, d0 - e2 = r1, c0 + _get_fraction(v2, v3, level) * spacing2[1], d0 - e3 = r0 + _get_fraction(v4, v3, level) * spacing2[0], c1, d0 - e4 = r0, c0 + _get_fraction(v1, v4, level) * spacing2[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 + # Calculate these edges normally + e1 = x0 + _get_fraction(v1, v2, level), y0, z0 + e2 = x1, y0 + _get_fraction(v2, v3, level), z0 + e3 = x0 + _get_fraction(v4, v3, level), y1, z0 + e4 = x0, y0 + _get_fraction(v1, v4, level), z0 # 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_spacing: - e5 = r0 + _get_fraction(v5, v6, level) * spacing2[0], c0, d1 - e6 = r1, c0 + _get_fraction(v6, v7, level) * spacing2[1], d1 - e7 = r0 + _get_fraction(v8, v7, level) * spacing2[0], c1, d1 - e8 = r0, c0 + _get_fraction(v5, v8, level) * spacing2[1], d1 - e9 = r0, c0, d0 + _get_fraction(v1, v5, level) * spacing2[2] - e10 = r1, c0, d0 + _get_fraction(v2, v6, level) * spacing2[2] - e11 = r0, c1, d0 + _get_fraction(v4, v8, level) * spacing2[2] - e12 = r1, c1, d0 + _get_fraction(v3, v7, level) * spacing2[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) + e5 = x0 + _get_fraction(v5, v6, level), y0, z1 + e6 = x1, y0 + _get_fraction(v6, v7, level), z1 + e7 = x0 + _get_fraction(v8, v7, level), y1, z1 + e8 = x0, y0 + _get_fraction(v5, v8, level), z1 + e9 = x0, y0, z0 + _get_fraction(v1, v5, level) + e10 = x1, y0, z0 + _get_fraction(v2, v6, level) + e11 = x0, y1, z0 + _get_fraction(v4, v8, level) + e12 = x1, y1, z0 + _get_fraction(v3, v7, level) # Append appropriate triangles to the growing output `face_list`