diff --git a/doc/examples/plot_marching_cubes.py b/doc/examples/plot_marching_cubes.py index a57a40a2..36685434 100644 --- a/doc/examples/plot_marching_cubes.py +++ b/doc/examples/plot_marching_cubes.py @@ -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 diff --git a/skimage/draw/draw3d.py b/skimage/draw/draw3d.py index 0b6fcb2d..db2f0d39 100644 --- a/skimage/draw/draw3d.py +++ b/skimage/draw/draw3d.py @@ -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 ------- diff --git a/skimage/draw/tests/test_draw3d.py b/skimage/draw/tests/test_draw3d.py index 59a7f6b3..2e1198eb 100644 --- a/skimage/draw/tests/test_draw3d.py +++ b/skimage/draw/tests/test_draw3d.py @@ -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] diff --git a/skimage/measure/_marching_cubes.py b/skimage/measure/_marching_cubes.py index 772c52f1..37dab1c0 100644 --- a/skimage/measure/_marching_cubes.py +++ b/skimage/measure/_marching_cubes.py @@ -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. diff --git a/skimage/measure/_marching_cubes_cy.pyx b/skimage/measure/_marching_cubes_cy.pyx index a0f63d1c..085108ab 100644 --- a/skimage/measure/_marching_cubes_cy.pyx +++ b/skimage/measure/_marching_cubes_cy.pyx @@ -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 diff --git a/skimage/measure/tests/test_marching_cubes.py b/skimage/measure/tests/test_marching_cubes.py index 7a1fd40a..b3c2ddc1 100644 --- a/skimage/measure/tests/test_marching_cubes.py +++ b/skimage/measure/tests/test_marching_cubes.py @@ -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)