diff --git a/skimage/measure/_marching_cubes.py b/skimage/measure/_marching_cubes.py index 96fb1102..cd0e1b60 100644 --- a/skimage/measure/_marching_cubes.py +++ b/skimage/measure/_marching_cubes.py @@ -1,8 +1,10 @@ import numpy as np +import scipy.ndimage as ndi from . import _marching_cubes_cy -def marching_cubes(volume, level, spacing=(1., 1., 1.)): +def marching_cubes(volume, level, spacing=(1., 1., 1.), + gradient_direction='descent'): """ Marching cubes algorithm to find iso-valued surfaces in 3d volumetric data @@ -15,6 +17,12 @@ def marching_cubes(volume, level, spacing=(1., 1., 1.)): spacing : length-3 tuple of floats Voxel spacing in spatial dimensions corresponding to numpy array indexing dimensions (M, N, P) as in `volume`. + gradient_direction : string + Controls if the mesh was generated from an isosurface with gradient + descent toward objects of interest (the default), or the opposite. + The two options are: + * descent : Object was greater than exterior + * ascent : Exterior was greater than object Returns ------- @@ -68,14 +76,10 @@ def marching_cubes(volume, level, spacing=(1., 1., 1.)): lexicographical order) coordinate in the contour. This is a side-effect of how the input array is traversed, but can be relied upon. - The generated mesh does not guarantee coherent orientation because of how - symmetry is used in the algorithm. If this is required, e.g. due to a - particular visualization package or for generating 3D printing STL files, - the utility ``skimage.measure.correct_mesh_orientation`` is available to - fix this in post-processing. + The generated mesh guarantees coherent orientation as of version 0.12. To quantify the area of an isosurface generated by this algorithm, pass - the outputs directly into `skimage.measure.mesh_surface_area`. + outputs directly into `skimage.measure.mesh_surface_area`. Regarding visualization of algorithm output, the ``mayavi`` package is recommended. To contour a volume named `myvolume` about the level 0.0:: @@ -122,8 +126,18 @@ def marching_cubes(volume, level, spacing=(1., 1., 1.)): # Returns a true mesh with no degenerate faces. verts, faces = _marching_cubes_cy.unpack_unique_verts(raw_faces) + verts = np.asarray(verts) + faces = np.asarray(faces) + + # Calculate gradient of `volume`, then interpolate to vertices in `verts` + grad_x, grad_y, grad_z = np.gradient(volume) + + # Fancy indexing to define two vector arrays from triangle vertices + faces = _correct_mesh_orientation(volume, verts[faces], faces, spacing, + gradient_direction) + # Adjust for non-isotropic spacing in `verts` at time of return - return np.asarray(verts) * np.r_[spacing], np.asarray(faces) + return verts * np.r_[spacing], faces def mesh_surface_area(verts, faces): @@ -187,7 +201,7 @@ def correct_mesh_orientation(volume, verts, faces, spacing=(1., 1., 1.), indexing dimensions (M, N, P) as in `volume`. gradient_direction : string Controls if the mesh was generated from an isosurface with gradient - ascent toward objects of interest (the default), or the opposite. + descent toward objects of interest (the default), or the opposite. The two options are: * descent : Object was greater than exterior * ascent : Exterior was greater than object @@ -225,17 +239,86 @@ def correct_mesh_orientation(volume, verts, faces, spacing=(1., 1., 1.), skimage.measure.mesh_surface_area """ - import scipy.ndimage as ndi + import warnings + warnings.warn( + DeprecationWarning("`correct_mesh_orientation` is deprecated for " + "removal as `marching_cubes` now guarantess " + "correct mesh orientation.")) - # Calculate gradient of `volume`, then interpolate to vertices in `verts` - grad_x, grad_y, grad_z = np.gradient(volume) + verts = verts.copy() + verts[:, 0] /= spacing[0] + verts[:, 1] /= spacing[1] + verts[:, 2] /= spacing[2] # Fancy indexing to define two vector arrays from triangle vertices actual_verts = verts[faces] - actual_verts[:, 0] /= spacing[0] - actual_verts[:, 1] /= spacing[1] - actual_verts[:, 2] /= spacing[2] - + + return _correct_mesh_orientation(volume, actual_verts, faces, spacing, + gradient_direction) + + +def _correct_mesh_orientation(volume, actual_verts, faces, + spacing=(1., 1., 1.), + gradient_direction='descent'): + """ + Correct orientations of mesh faces. + + Parameters + ---------- + volume : (M, N, P) array of doubles + Input data volume to find isosurfaces. Will be cast to `np.float64`. + actual_verts : (F, 3, 3) array of floats + Array with (face, vertex, coords) index coordinates. + faces : (F, 3) array of ints + List of length-3 lists of integers, referencing vertex coordinates as + provided in `verts`. + spacing : length-3 tuple of floats + Voxel spacing in spatial dimensions corresponding to numpy array + indexing dimensions (M, N, P) as in `volume`. + gradient_direction : string + Controls if the mesh was generated from an isosurface with gradient + descent toward objects of interest (the default), or the opposite. + The two options are: + * descent : Object was greater than exterior + * ascent : Exterior was greater than object + + Returns + ------- + faces_corrected (F, 3) array of ints + Corrected list of faces referencing vertex coordinates in `verts`. + + Notes + ----- + Certain applications and mesh processing algorithms require all faces + to be oriented in a consistent way. Generally, this means a normal vector + points "out" of the meshed shapes. This algorithm corrects the output from + `skimage.measure.marching_cubes` by flipping the orientation of + mis-oriented faces. + + Because marching cubes could be used to find isosurfaces either on + gradient descent (where the desired object has greater values than the + exterior) or ascent (where the desired object has lower values than the + exterior), the ``gradient_direction`` kwarg allows the user to inform this + algorithm which is correct. If the resulting mesh appears to be oriented + completely incorrectly, try changing this option. + + The arguments expected by this function are the exact outputs from + `skimage.measure.marching_cubes` except `actual_verts`, which is an + uncorrected version of the fancy indexing operation `verts[faces]`. + Only `faces` is corrected and returned as the vertices do not change, + only the order in which they are referenced. + + This algorithm assumes ``faces`` provided are exclusively triangles. + + See Also + -------- + skimage.measure.marching_cubes + skimage.measure.mesh_surface_area + + """ + # Calculate gradient of `volume`, then interpolate to vertices in `verts` + grad_x, grad_y, grad_z = np.gradient(volume) + a = actual_verts[:, 0, :] - actual_verts[:, 1, :] b = actual_verts[:, 0, :] - actual_verts[:, 2, :] diff --git a/skimage/measure/tests/test_marching_cubes.py b/skimage/measure/tests/test_marching_cubes.py index 60adbb0d..0e3d77ec 100644 --- a/skimage/measure/tests/test_marching_cubes.py +++ b/skimage/measure/tests/test_marching_cubes.py @@ -39,7 +39,24 @@ def test_invalid_input(): def test_correct_mesh_orientation(): sphere_small = ellipsoid(1, 1, 1, levelset=True) - verts, faces = marching_cubes(sphere_small, 0.) + + # Mesh with incorrectly oriented faces which was previously returned from + # `marching_cubes`, before it guaranteed correct mesh orientation + verts = np.array([[1., 2., 2.], + [2., 2., 1.], + [2., 1., 2.], + [2., 2., 3.], + [2., 3., 2.], + [3., 2., 2.]]) + + faces = np.array([[0, 1, 2], + [2, 0, 3], + [1, 0, 4], + [4, 0, 3], + [1, 2, 5], + [2, 3, 5], + [1, 4, 5], + [5, 4, 3]]) # Correct mesh orientation - descent corrected_faces1 = correct_mesh_orientation(sphere_small, verts, faces,