mirror of
https://github.com/wassname/scikit-image.git
synced 2026-07-02 11:37:39 +08:00
FEAT: Add function to correct mesh face orientations
This commit is contained in:
Josh Warner (Mac)
@@ -75,11 +75,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 +90,7 @@ def marching_cubes(volume, level, spacing=(1., 1., 1.)):
|
||||
|
||||
See Also
|
||||
--------
|
||||
skimage.measure.correct_mesh_orientation
|
||||
skimage.measure.mesh_surface_area
|
||||
|
||||
"""
|
||||
@@ -106,17 +107,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 +146,122 @@ 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.linalg.norm(grad_centroids, axis=1)[:, np.newaxis])
|
||||
|
||||
# Calculate and normalize cross products for each face
|
||||
crosses = np.cross(a, b)
|
||||
crosses = crosses / np.linalg.norm(crosses, axis=1)[:, 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))
|
||||
|
||||
# Swap orientation and returh, without modifying original data
|
||||
faces_corrected = faces.copy()
|
||||
faces_corrected[indices] = faces_corrected[indices, ::-1]
|
||||
|
||||
return faces_corrected
|
||||
|
||||
Reference in New Issue
Block a user