wassname/scikit-image
1
0
Fork 0
mirror of https://github.com/wassname/scikit-image.git synced 2026-07-19 11:27:45 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #780 from JDWarner/spacing_marching_cubes

Browse Source 
DOC: Change sampling kwarg name to spacing in marching_cubes
This commit is contained in:
Johannes Schönberger
2013-10-13 12:00:03 -07:00
parent 0a8493805e b9f5dd3ad5
commit d1ba9fa6e5
6 changed files with 59 additions and 61 deletions
Download Patch File Download Diff File Expand all files Collapse all files
doc/examples/plot_marching_cubes.py
+1 -1
View File
@@ -17,7 +17,7 @@ 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.
`spacing` kwarg.
"""
import numpy as np
skimage/draw/draw3d.py
+12 -14
View File
@@ -3,10 +3,10 @@ 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):
def ellipsoid(a, b, c, spacing=(1., 1., 1.), levelset=False):
"""
Generates ellipsoid with semimajor axes aligned with grid dimensions
on grid with specified `sampling`.
on grid with specified `spacing`.
Parameters
----------
@@ -16,8 +16,8 @@ def ellipsoid(a, b, c, sampling=(1., 1., 1.), levelset=False):
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.
spacing : tuple of floats, length 3
Spacing 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.
@@ -26,7 +26,7 @@ def ellipsoid(a, b, c, sampling=(1., 1., 1.), levelset=False):
Returns
-------
ellip : (N, M, P) array
Ellipsoid centered in a correctly sized array for given `sampling`.
Ellipsoid centered in a correctly sized array for given `spacing`.
Boolean dtype unless `levelset=True`, in which case a float array is
returned with the level set above 0.0 representing the ellipsoid.
@@ -34,7 +34,7 @@ def ellipsoid(a, b, c, sampling=(1., 1., 1.), levelset=False):
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]
offset = np.r_[1, 1, 1] * np.r_[spacing]
# Calculate limits, and ensure output volume is odd & symmetric
low = np.ceil((- np.r_[a, b, c] - offset))
@@ -43,14 +43,14 @@ def ellipsoid(a, b, c, sampling=(1., 1., 1.), levelset=False):
for dim in range(3):
if (high[dim] - low[dim]) % 2 == 0:
low[dim] -= 1
num = np.arange(low[dim], high[dim], sampling[dim])
num = np.arange(low[dim], high[dim], spacing[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]]
x, y, z = np.mgrid[low[0]:high[0]:spacing[0],
low[1]:high[1]:spacing[1],
low[2]:high[2]:spacing[2]]
if not levelset:
arr = ((x / float(a)) ** 2 +
@@ -64,10 +64,10 @@ def ellipsoid(a, b, c, sampling=(1., 1., 1.), levelset=False):
return arr
def ellipsoid_stats(a, b, c, sampling=(1., 1., 1.)):
def ellipsoid_stats(a, b, c):
"""
Calculates analytical surface area and volume for ellipsoid with
semimajor axes aligned with grid dimensions of specified `sampling`.
semimajor axes aligned with grid dimensions of specified `spacing`.
Parameters
----------
@@ -77,8 +77,6 @@ def ellipsoid_stats(a, b, c, sampling=(1., 1., 1.)):
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
-------
skimage/draw/tests/test_draw3d.py
+2 -2
View File
@@ -6,7 +6,7 @@ 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 = ellipsoid(2, 2, 4, spacing=(1., 1., 2.))
test_anisotropic = test_anisotropic[1:-1, 1:-1, 1:-1]
expected = np.array([[[0, 0, 0, 0, 0],
@@ -45,7 +45,7 @@ def test_ellipsoid_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.),
test_anisotropic = ellipsoid(2, 2, 4, spacing=(1., 1., 2.),
levelset=True)
test_anisotropic = test_anisotropic[1:-1, 1:-1, 1:-1]
skimage/measure/_marching_cubes.py
+3 -3
View File
@@ -2,7 +2,7 @@ import numpy as np
from . import _marching_cubes_cy
def marching_cubes(volume, level, sampling=(1., 1., 1.)):
def marching_cubes(volume, level, spacing=(1., 1., 1.)):
"""
Marching cubes algorithm to find iso-valued surfaces in 3d volumetric data
@@ -12,7 +12,7 @@ def marching_cubes(volume, level, sampling=(1., 1., 1.)):
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
spacing : length-3 tuple of floats
Voxel spacing in spatial dimensions corresponding to numpy array
indexing dimensions (M, N, P) as in `volume`.
@@ -107,7 +107,7 @@ def marching_cubes(volume, level, sampling=(1., 1., 1.)):
# 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)
spacing)
# Find and collect unique vertices, storing triangle verts as indices.
# Returns a true mesh with no degenerate faces.
skimage/measure/_marching_cubes_cy.pyx
+36 -36
View File
@@ -56,32 +56,32 @@ def unpack_unique_verts(list trilist):
def iterate_and_store_3d(double[:, :, ::1] arr, double level,
tuple sampling=(1., 1., 1.)):
tuple spacing=(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
If `spacing` 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
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(sampling) != 3:
raise ValueError("`sampling` must be (double, double, double)")
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_sampling, plus_z
cdef bint odd_spacing, plus_z
plus_z = False
if [float(i) for i in sampling] == [1.0, 1.0, 1.0]:
odd_sampling = False
if [float(i) for i in spacing] == [1.0, 1.0, 1.0]:
odd_spacing = False
else:
odd_sampling = True
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,11 +107,11 @@ def iterate_and_store_3d(double[:, :, ::1] arr, double level,
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]
# 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) *
@@ -138,15 +138,15 @@ 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_sampling:
if odd_spacing:
# 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]
# 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
@@ -193,11 +193,11 @@ def iterate_and_store_3d(double[:, :, ::1] arr, double level,
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
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
@@ -208,15 +208,15 @@ def iterate_and_store_3d(double[:, :, ::1] arr, double level,
# 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]
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
skimage/measure/tests/test_marching_cubes.py
+5 -5
View File
@@ -16,12 +16,12 @@ def test_marching_cubes_isotropic():
def test_marching_cubes_anisotropic():
sampling = (1., 10 / 6., 16 / 6.)
ellipsoid_anisotropic = ellipsoid(6, 10, 16, sampling=sampling,
spacing = (1., 10 / 6., 16 / 6.)
ellipsoid_anisotropic = ellipsoid(6, 10, 16, spacing=spacing,
levelset=True)
_, surf = ellipsoid_stats(6, 10, 16, sampling=sampling)
_, surf = ellipsoid_stats(6, 10, 16)
verts, faces = marching_cubes(ellipsoid_anisotropic, 0.,
sampling=sampling)
spacing=spacing)
surf_calc = mesh_surface_area(verts, faces)
# Test within 1.5% tolerance for anisotropic. Will always underestimate.
@@ -32,7 +32,7 @@ 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))
spacing=(1, 2))
assert_raises(ValueError, marching_cubes, np.zeros((20, 20)), 0)