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

Merge pull request #882 from JDWarner/correct_mesh_orientation

Browse Source 
Correct mesh orientation
This commit is contained in:
Johannes Schönberger
2014-03-14 18:10:28 -04:00
parent 9501c2dc25 5ac910b0b8
commit 944c9c6bde
3 changed files with 149 additions and 9 deletions
Download Patch File Download Diff File Expand all files Collapse all files
skimage/measure/__init__.py
+3 -1
View File
@@ -1,5 +1,6 @@
from ._find_contours import find_contours
from ._marching_cubes import marching_cubes, mesh_surface_area
from ._marching_cubes import (marching_cubes, mesh_surface_area,
correct_mesh_orientation)
from ._regionprops import regionprops, perimeter
from ._structural_similarity import structural_similarity
from ._polygon import approximate_polygon, subdivide_polygon
@@ -26,4 +27,5 @@ __all__ = ['find_contours',
'moments_hu',
'marching_cubes',
'mesh_surface_area',
'correct_mesh_orientation',
'profile_line']
skimage/measure/_marching_cubes.py
+124 -7
View File
@@ -68,6 +68,12 @@ 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.
To quantify the area of an isosurface generated by this algorithm, pass
the outputs directly into `skimage.measure.mesh_surface_area`.
@@ -75,11 +81,11 @@ def marching_cubes(volume, level, spacing=(1., 1., 1.)):
is recommended. To contour a volume named `myvolume` about the level 0.0::
>>> from mayavi import mlab # doctest: +SKIP
>>> verts, tris = marching_cubes(myvolume, 0.0, (1., 1., 2.)) # doctest: +SKIP
>>> verts, faces = marching_cubes(myvolume, 0.0, (1., 1., 2.)) # doctest: +SKIP
>>> mlab.triangular_mesh([vert[0] for vert in verts],
... [vert[1] for vert in verts],
... [vert[2] for vert in verts],
... tris) # doctest: +SKIP
... faces) # doctest: +SKIP
>>> mlab.show() # doctest: +SKIP
References
@@ -90,6 +96,7 @@ def marching_cubes(volume, level, spacing=(1., 1., 1.)):
See Also
--------
skimage.measure.correct_mesh_orientation
skimage.measure.mesh_surface_area
"""
@@ -106,17 +113,17 @@ 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_tris = _marching_cubes_cy.iterate_and_store_3d(volume, float(level),
spacing)
raw_faces = _marching_cubes_cy.iterate_and_store_3d(volume, float(level),
spacing)
# 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_tris)
verts, faces = _marching_cubes_cy.unpack_unique_verts(raw_faces)
return np.asarray(verts), np.asarray(faces)
def mesh_surface_area(verts, tris):
def mesh_surface_area(verts, faces):
"""
Compute surface area, given vertices & triangular faces
@@ -145,13 +152,123 @@ def mesh_surface_area(verts, tris):
See Also
--------
skimage.measure.marching_cubes
skimage.measure.correct_mesh_orientation
"""
# Fancy indexing to define two vector arrays from triangle vertices
actual_verts = verts[tris]
actual_verts = verts[faces]
a = actual_verts[:, 0, :] - actual_verts[:, 1, :]
b = actual_verts[:, 0, :] - actual_verts[:, 2, :]
del actual_verts
# Area of triangle in 3D = 1/2 * Euclidean norm of cross product
return ((np.cross(a, b) ** 2).sum(axis=1) ** 0.5).sum() / 2.
def correct_mesh_orientation(volume, 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`.
verts : (V, 3) array of floats
Array containing (x, y, z) coordinates for V unique mesh vertices.
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
ascent 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`. 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 all triangles.
See Also
--------
skimage.measure.marching_cubes
skimage.measure.mesh_surface_area
"""
import scipy.ndimage as ndi
# Calculate gradient of `volume`, then interpolate to vertices in `verts`
grad_x, grad_y, grad_z = np.gradient(volume, *spacing)
# Fancy indexing to define two vector arrays from triangle vertices
actual_verts = verts[faces]
a = actual_verts[:, 0, :] - actual_verts[:, 1, :]
b = actual_verts[:, 0, :] - actual_verts[:, 2, :]
# Find triangle centroids
centroids = (actual_verts.sum(axis=1) / 3.).T
del actual_verts
# Interpolate face centroids into each gradient axis
grad_centroids_x = ndi.map_coordinates(grad_x, centroids)
grad_centroids_y = ndi.map_coordinates(grad_y, centroids)
grad_centroids_z = ndi.map_coordinates(grad_z, centroids)
# Combine and normalize interpolated gradients
grad_centroids = np.c_[grad_centroids_x, grad_centroids_y,
grad_centroids_z]
grad_centroids = (grad_centroids /
(np.sum(grad_centroids ** 2,
axis=1) ** 0.5)[:, np.newaxis])
# Find normal vectors for each face via cross product
crosses = np.cross(a, b)
crosses = crosses / (np.sum(crosses ** 2, axis=1) ** (0.5))[:, np.newaxis]
# Take dot product
dotproducts = (grad_centroids * crosses).sum(axis=1)
# Find mis-oriented faces
if 'descent' in gradient_direction:
# Faces with incorrect orientations have dot product < 0
indices = (dotproducts < 0).nonzero()[0]
elif 'ascent' in gradient_direction:
# Faces with incorrection orientation have dot product > 0
indices = (dotproducts > 0).nonzero()[0]
else:
raise ValueError("Incorrect input %s in `gradient_direction`, see "
"docstring." % (gradient_direction))
# Correct orientation and return, without modifying original data
faces_corrected = faces.copy()
faces_corrected[indices] = faces_corrected[indices, ::-1]
return faces_corrected
skimage/measure/tests/test_marching_cubes.py
+22 -1
View File
@@ -2,7 +2,8 @@ import numpy as np
from numpy.testing import assert_raises
from skimage.draw import ellipsoid, ellipsoid_stats
from skimage.measure import marching_cubes, mesh_surface_area
from skimage.measure import (marching_cubes, mesh_surface_area,
correct_mesh_orientation)
def test_marching_cubes_isotropic():
@@ -36,5 +37,25 @@ def test_invalid_input():
assert_raises(ValueError, marching_cubes, np.zeros((20, 20)), 0)
def test_correct_mesh_orientation():
sphere_small = ellipsoid(1, 1, 1, levelset=True)
verts, faces = marching_cubes(sphere_small, 0.)
# Correct mesh orientation - descent
corrected_faces1 = correct_mesh_orientation(sphere_small, verts, faces,
gradient_direction='descent')
corrected_faces2 = correct_mesh_orientation(sphere_small, verts, faces,
gradient_direction='ascent')
# Ensure ascent is opposite of descent for all faces
np.testing.assert_array_equal(corrected_faces1, corrected_faces2[:, ::-1])
# Ensure correct faces have been reversed: 1, 4, and 5
idx = [1, 4, 5]
expected = faces.copy()
expected[idx] = expected[idx, ::-1]
np.testing.assert_array_equal(expected, corrected_faces1)
if __name__ == '__main__':
np.testing.run_module_suite()