mirror of
https://github.com/wassname/scikit-image.git
synced 2026-07-20 12:40:31 +08:00
Josh Warner (Mac)
5eba5a42a7
Also includes surface area calculation algorithm from generated mesh. Convenient output to visualize with `mayavi.mlab`. Efficient Cython implementation.
201 lines
7.3 KiB
Python
201 lines
7.3 KiB
Python
import numpy as np
|
|
from . import _marching_cubes
|
|
|
|
|
|
def marching_cubes(volume, level, sampling=(1., 1., 1.)):
|
|
"""
|
|
Marching cubes algorithm to find iso-valued surfaces in 3d volumetric data
|
|
|
|
Parameters
|
|
----------
|
|
volume : (M, N, P) array of doubles
|
|
Input data volume to find isosurfaces. Will be cast to `np.float64` if
|
|
not provided in this format.
|
|
level : float
|
|
Contour value to search for isosurfaces in `volume`.
|
|
sampling : length-3 tuple of floats
|
|
Voxel spacing in spatial dimensions corresponding to numpy array
|
|
indexing dimensions (M, N, P) as in `volume`.
|
|
|
|
Returns
|
|
-------
|
|
vert_list : list
|
|
Every entry in this list is a unique vertex on the isosurface.
|
|
tri_list : list
|
|
Every entry in this list is a length-3 list of integers. These
|
|
represent triangular faces; the integers in each sub-list correspond
|
|
to vertices held in `vert_list`.
|
|
|
|
Notes
|
|
-----
|
|
The marching cubes algorithm is implemented as described in [1]_.
|
|
A simple explanation is available here::
|
|
|
|
http://www.essi.fr/~lingrand/MarchingCubes/algo.html
|
|
|
|
There are several known ambiguous cases in the marching cubes algorithm.
|
|
Using point labeling as in [1]_, Figure 4, as shown:
|
|
|
|
v8 ------ v7
|
|
/ | / | y
|
|
/ | / | ^ z
|
|
v4 ------ v3 | | /
|
|
| v5 ----|- v6 |/ (note: NOT right handed!)
|
|
| / | / ----> x
|
|
|/ | /
|
|
v1 ------ v2
|
|
|
|
Most notably, if v4, v8, v2, and v6 are all >= `level` (or any
|
|
generalization of this case) two parallel planes are generated by this
|
|
algorithm, separating v4 and v8 from v2 and v6. An equally valid
|
|
interpretation would be a single connected thin surface enclosing all
|
|
four points. This is the best known ambiguity, though there are others.
|
|
|
|
This algorithm does not attempt to resolve such ambiguities; it is a naive
|
|
implementation of marching cubes as in [1]_, but may be a good beginning
|
|
for work with more recent techniques (Dual Marching Cubes, Extended
|
|
Marching Cubes, Cubic Marching Squares, etc.).
|
|
|
|
Because of interactions between neighboring cubes, the isosurface(s)
|
|
generated by this algorithm are NOT guaranteed to be closed, particularly
|
|
for complicated contours. Furthermore, this algorithm does not guarantee
|
|
a single contour will be returned. Indeed, ALL isosurfaces which cross
|
|
`level` will be found, regardless of connectivity.
|
|
|
|
The output is a triangular mesh consisting of a set of unique vertices and
|
|
connecting triangles. The order of these vertices and triangles in the
|
|
output list is determined by the position of the smallest ``x,y,z`` (in
|
|
lexicographical order) coordinate in the contour. This is a side-effect
|
|
of how the input array is traversed, but can be relied upon.
|
|
|
|
To quantify the area of an isosurface generated by this algorithm, pass
|
|
the output 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:
|
|
|
|
>>> from mayavi import mlab
|
|
>>> verts, tris = marching_cubes(myvolume, 0.0, (1., 1., 2.))
|
|
>>> mlab.triangular_mesh([vert[0] for vert in verts],
|
|
[vert[1] for vert in verts],
|
|
[vert[2] for vert in verts],
|
|
tris)
|
|
>>> mlab.show()
|
|
|
|
References
|
|
----------
|
|
.. [1] Lorensen, William and Harvey E. Cline. Marching Cubes: A High
|
|
Resolution 3D Surface Construction Algorithm. Computer Graphics
|
|
(SIGGRAPH 87 Proceedings) 21(4) July 1987, p. 163-170).
|
|
|
|
See Also
|
|
--------
|
|
skimage.measure.mesh_surface_area
|
|
|
|
"""
|
|
# Check inputs
|
|
if volume.ndim != 3:
|
|
raise ValueError("Input volume must be 3d.")
|
|
if volume.dtype.kind == 'f':
|
|
volume = volume.astype(np.float)
|
|
else:
|
|
from skimage.util import img_as_float
|
|
# If incorrect type provided, convert BOTH contour value and input
|
|
# volume using same method
|
|
level = img_as_float(np.array(level, dtype=volume.dtype))[0]
|
|
volume = img_as_float(volume)
|
|
|
|
# Extract raw triangles using marching cubes in Cython
|
|
# Returns a list of length-3 lists, each sub-list containing three
|
|
# tuples. The tuples hold (x, y, z) coordinates for triangle vertices.
|
|
# 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.iterate_and_store_3d(volume, float(level),
|
|
sampling)
|
|
|
|
# Find and collect unique vertices, storing triangle verts as indices.
|
|
# Removes much redundancy and eliminates degenerate "triangles".
|
|
vert_list, tri_list = _unpack_unique_verts(raw_tris)
|
|
|
|
return vert_list, tri_list
|
|
|
|
|
|
def _unpack_unique_verts(trilist):
|
|
"""
|
|
Converts a list of lists of tuples corresponding to triangle vertices into
|
|
a unique vertex list, and a list of triangle faces w/indices corresponding
|
|
to entries of the vertex list.
|
|
|
|
"""
|
|
idx = 0
|
|
vert_index = {}
|
|
vert_list = []
|
|
tri_list = []
|
|
|
|
# Iterate over triangles
|
|
for i in range(len(trilist)):
|
|
templist = []
|
|
|
|
# Only parse vertices from non-degenerate triangles
|
|
if not ((trilist[i][0] == trilist[i][1]) or
|
|
(trilist[i][0] == trilist[i][2]) or
|
|
(trilist[i][1] == trilist[i][2])):
|
|
|
|
# Iterate over vertices within each triangle
|
|
for j in range(3):
|
|
vert = trilist[i][j]
|
|
|
|
# Check if a new unique vertex found
|
|
if vert not in vert_index:
|
|
vert_index[vert] = idx
|
|
templist.append(idx)
|
|
vert_list.append(vert)
|
|
idx += 1
|
|
else:
|
|
templist.append(vert_index[vert])
|
|
|
|
tri_list.append(templist)
|
|
|
|
return vert_list, tri_list
|
|
|
|
|
|
def mesh_surface_area(verts, tris):
|
|
"""
|
|
Compute surface area, given vertices & triangular faces
|
|
|
|
Parameters
|
|
----------
|
|
verts : list
|
|
List of length-3 NumPy arrays containing vertex coordinates.
|
|
Units in each dimension should be consistent.
|
|
|
|
tris : list
|
|
List of length-3 lists of integers, referencing vertex coordinates as
|
|
provided in `verts`
|
|
|
|
Returns
|
|
-------
|
|
area : float
|
|
Surface area of mesh. Units in coordinates maintained, but squared.
|
|
|
|
Notes
|
|
-----
|
|
The arguments expected by this function are the exact outputs from
|
|
`skimage.measure.marching_cubes`. For unit correct output, ensure correct
|
|
`spacing` was passed to `skimage.measure.marching_cubes`.
|
|
|
|
See Also
|
|
--------
|
|
skimage.measure.marching_cubes
|
|
|
|
"""
|
|
# Define two vector arrays `a` and `b` from triangle vertices
|
|
actual_verts = np.array([[verts[i] for i in tri] for tri in tris])
|
|
a = actual_verts[:, 0, :] - actual_verts[:, 1, :]
|
|
b = actual_verts[:, 0, :] - actual_verts[:, 2, :]
|
|
del actual_verts
|
|
|
|
# Area of triangle = 1/2 * Euclidean norm of cross product
|
|
return ((np.cross(a, b) ** 2).sum(axis=1) ** 0.5).sum() / 2.
|