wassname/scikit-image
1
0
Fork 0
mirror of https://github.com/wassname/scikit-image.git synced 2026-07-18 12:40:14 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #1115 from JDWarner/enh_marching_cubes

Browse Source 
ENH: Simplify `marching_cubes` for maintainability
This commit is contained in:
Stefan van der Walt
2014-09-01 17:32:26 +01:00
parent 725c199018 af59d06461
commit 85439b173a
2 changed files with 22 additions and 65 deletions
Download Patch File Download Diff File Expand all files Collapse all files
skimage/measure/_marching_cubes.py
+6 -3
View File
@@ -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):
skimage/measure/_marching_cubes_cy.pyx
+16 -62
View File
@@ -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`