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..e247e360 --- /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 + Analytically calculated volume of ellipsoid. + surf : float + Analytically 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/measure/tests/test_marching_cubes.py b/skimage/measure/tests/test_marching_cubes.py index e8135788..991da060 100644 --- a/skimage/measure/tests/test_marching_cubes.py +++ b/skimage/measure/tests/test_marching_cubes.py @@ -1,114 +1,13 @@ import numpy as np from numpy.testing import assert_raises -from scipy.special import (ellipkinc as ellip_F, ellipeinc as ellip_E) +from skimage.draw import ellipsoid, ellipsoid_stats from skimage.measure import marching_cubes, mesh_surface_area -def _ellipsoid(a, b, c, sampling=(1., 1., 1.), info=False, tight=False, - 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 each spatial dimension - info : bool - If False, only `bool_arr` returned. - If True, (`bool_arr`, `vol`, `surf`) returned; the additional - values are analytical volume and surface area calculated for - this ellipsoid. - tight : bool - Controls if the ellipsoid will precisely be contained within - the returned volume (tight=True) or if each dimension will be - 2 longer than necessary (tight=False). For algorithms which - need both sides of a contour, use False. - 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 - ------- - bool_arr : (N, M, P) array - Sphere in an appropriately sized boolean array. - vol : float - Analytically calculated volume of ellipsoid. Only returned if - `info` is True. - surf : float - Analytically calculated surface area of ellipsoid. Only returned - if `info` is True. - - """ - if not tight: - offset = np.r_[1, 1, 1] * np.r_[sampling] - else: - offset = np.r_[0, 0, 0] - - # 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 - - if not info: - return arr - else: - # 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 arr, vol, surf - - def test_marching_cubes_isotropic(): - ellipsoid_isotropic, _, surf = _ellipsoid(6, 10, 16, - levelset=True, - info=True) + ellipsoid_isotropic = ellipsoid(6, 10, 16, levelset=True) + _, surf = ellipsoid_stats(6, 10, 16, levelset=True) verts, faces = marching_cubes(ellipsoid_isotropic, 0.) surf_calc = mesh_surface_area(verts, faces) @@ -118,13 +17,13 @@ def test_marching_cubes_isotropic(): def test_marching_cubes_anisotropic(): sampling = (1., 10 / 6., 16 / 6.) - ellipsoid_isotropic, _, surf = _ellipsoid(6, 10, 16, - sampling=sampling, - levelset=True, - info=True) - verts, faces = marching_cubes(ellipsoid_isotropic, 0., + ellipsoid_anisotropic, _, surf = ellipsoid(6, 10, 16, sampling=sampling, + levelset=True) + _, surf = ellipsoid_stats(6, 10, 16, sampling=sampling, levelset=True) + 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