From fe2edd51491f960343af3f6f5fcb2925130cfa85 Mon Sep 17 00:00:00 2001 From: "Josh Warner (Mac)" Date: Sun, 9 Feb 2014 13:16:19 -0600 Subject: [PATCH] FEAT: Add function to correct mesh face orientations --- skimage/measure/_marching_cubes.py | 124 +++++++++++++++++++++++++++-- 1 file changed, 117 insertions(+), 7 deletions(-) diff --git a/skimage/measure/_marching_cubes.py b/skimage/measure/_marching_cubes.py index 7ea0667e..b42c8c9d 100644 --- a/skimage/measure/_marching_cubes.py +++ b/skimage/measure/_marching_cubes.py @@ -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