mirror of
https://github.com/wassname/scikit-image.git
synced 2026-06-27 17:31:50 +08:00
Almar Klein
258 lines
11 KiB
Python
258 lines
11 KiB
Python
import sys
|
|
import base64
|
|
import dis
|
|
import inspect
|
|
|
|
import numpy as np
|
|
|
|
if sys.version_info >= (3, ):
|
|
base64decode = base64.decodebytes
|
|
ordornot = lambda x:x
|
|
else:
|
|
base64decode = base64.decodestring
|
|
ordornot = ord
|
|
|
|
from . import _marching_cubes_lewiner_luts as mcluts
|
|
from . import _marching_cubes_lewiner_cy
|
|
from .._shared.utils import skimage_deprecation, warn
|
|
|
|
|
|
def expected_output_args():
|
|
""" Get number of expected output args.
|
|
|
|
Please don't use this to influence the algorithmic bahaviour of a function.
|
|
For ``a, b, rest*, c = ...`` syntax, returns n + 0.1 (3.1 in this example).
|
|
"""
|
|
f = inspect.currentframe().f_back.f_back
|
|
i = f.f_lasti + 3
|
|
bytecode = f.f_code.co_code
|
|
instruction = ordornot(bytecode[i])
|
|
while True:
|
|
if instruction == dis.opmap['DUP_TOP']:
|
|
if ordornot(bytecode[i + 1]) == dis.opmap['UNPACK_SEQUENCE']:
|
|
return ordornot(bytecode[i + 2])
|
|
i += 4
|
|
instruction = ordornot(bytecode[i])
|
|
continue
|
|
if instruction == dis.opmap['STORE_NAME']:
|
|
return 1
|
|
if instruction == dis.opmap['UNPACK_SEQUENCE']:
|
|
return ordornot(bytecode[i + 1])
|
|
if instruction == dis.opmap.get('UNPACK_EX', -1): # py3k
|
|
return ordornot(bytecode[i + 1]) + ordornot(bytecode[i + 2]) + 0.1
|
|
return 0
|
|
|
|
|
|
def marching_cubes(volume, level=None, spacing=(1., 1., 1.),
|
|
gradient_direction='descent', step_size=1,
|
|
allow_degenerate=True, use_classic=False):
|
|
"""
|
|
Lewiner marching cubes algorithm to find surfaces in 3d volumetric data.
|
|
|
|
In contrast to ``marching_cubes_classic()``, this algorithm is faster,
|
|
resolves ambiguities, and guarantees topologically correct results.
|
|
Therefore, this algorithm generally a better choice, unless there
|
|
is a specific need for the classic algorithm.
|
|
|
|
Parameters
|
|
----------
|
|
volume : (M, N, P) array
|
|
Input data volume to find isosurfaces. Will internally be
|
|
converted to float32 if necessary.
|
|
level : float
|
|
Contour value to search for isosurfaces in `volume`. If not
|
|
given or None, the average of the min and max of vol is used.
|
|
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,
|
|
considering the *left-hand* rule.
|
|
The two options are:
|
|
* descent : Object was greater than exterior
|
|
* ascent : Exterior was greater than object
|
|
step_size : int
|
|
Step size in voxels. Default 1. Larger steps yield faster but
|
|
coarser results. The result will always be topologically correct
|
|
though.
|
|
allow_degenerate : bool
|
|
Whether to allow degenerate (i.e. zero-area) triangles in the
|
|
end-result. Default True. If False, degenerate triangles are
|
|
removed, at the cost of making the algorithm slower.
|
|
use_classic : bool
|
|
If given and True, the classic marching cubes by Lorensen (1987)
|
|
is used. This option is included for reference purposes. Note
|
|
that this algorithm has ambiguities and is not guaranteed to
|
|
produce a topologically correct result. The results with using
|
|
this option are *not* generally the same as the
|
|
``marching_cubes_classic()`` function.
|
|
|
|
Returns
|
|
-------
|
|
verts : (V, 3) array
|
|
Spatial coordinates for V unique mesh vertices. Coordinate order
|
|
matches input `volume` (M, N, P).
|
|
faces : (F, 3) array
|
|
Define triangular faces via referencing vertex indices from ``verts``.
|
|
This algorithm specifically outputs triangles, so each face has
|
|
exactly three indices.
|
|
normals : (V, 3) array
|
|
The normal direction at each vertex, as calculated from the
|
|
data.
|
|
values : (V, ) array
|
|
Gives a measure for the maximum value of the data in the local region
|
|
near each vertex. This can be used by visualization tools to apply
|
|
a colormap to the mesh.
|
|
|
|
Notes
|
|
-----
|
|
The algorithm [1] is an improved version of Chernyaev's Marching
|
|
Cubes 33 algorithm. It is an efficient algorithm that relies on
|
|
heavy use of lookup tables to handle the many different cases,
|
|
keeping the algorithm relatively easy. This implementation is
|
|
written in Cython, ported from Lewiner's C++ implementation.
|
|
|
|
References
|
|
----------
|
|
.. [1] Thomas Lewiner, Helio Lopes, Antonio Wilson Vieira and Geovan
|
|
Tavares. Efficient implementation of Marching Cubes' cases with
|
|
topological guarantees. Journal of Graphics Tools 8(2)
|
|
pp. 1-15 (december 2003).
|
|
DOI: 10.1080/10867651.2003.10487582
|
|
|
|
See Also
|
|
--------
|
|
skimage.measure.marching_cubes_classic
|
|
skimage.measure.mesh_surface_area
|
|
"""
|
|
|
|
# This signature (output args) of this func changed after 0.12
|
|
try:
|
|
nout = expected_output_args()
|
|
except Exception:
|
|
nout = 0
|
|
if nout <= 2:
|
|
warn(skimage_deprecation('`marching_cubes` now uses a better and '
|
|
'faster algorithm, and returns four instead '
|
|
'of two outputs (see docstring for details). '
|
|
'Backwards compatibility with 0.12 and prior '
|
|
'is available with `marching_cubes_classic`.'))
|
|
|
|
return marching_cubes_lewiner(volume, level, spacing, gradient_direction,
|
|
step_size, allow_degenerate, use_classic)
|
|
|
|
|
|
def marching_cubes_lewiner(volume, level=None, spacing=(1., 1., 1.),
|
|
gradient_direction='descent', step_size=1,
|
|
allow_degenerate=True, use_classic=False):
|
|
""" Alias for ``marching_cubes()``.
|
|
"""
|
|
|
|
# Check volume and ensure its in the format that the alg needs
|
|
if not isinstance(volume, np.ndarray) or (volume.ndim != 3):
|
|
raise ValueError('Input volume should be a 3D numpy array.')
|
|
if volume.shape[0] < 2 or volume.shape[1] < 2 or volume.shape[2] < 2:
|
|
raise ValueError("Input array must be at least 2x2x2.")
|
|
volume = np.ascontiguousarray(volume, np.float32) # no copy if not necessary
|
|
|
|
# Check/convert other inputs:
|
|
# level
|
|
if level is None:
|
|
level = 0.5 * (volume.min() + volume.max())
|
|
else:
|
|
level = float(level)
|
|
if level < volume.min() or level > volume.max():
|
|
raise ValueError("Surface level must be within volume data range.")
|
|
# spacing
|
|
if len(spacing) != 3:
|
|
raise ValueError("`spacing` must consist of three floats.")
|
|
# step_size
|
|
step_size = int(step_size)
|
|
if step_size < 1:
|
|
raise ValueError('step_size must be at least one.')
|
|
# use_classic
|
|
use_classic = bool(use_classic)
|
|
|
|
# Get LutProvider class (reuse if possible)
|
|
L = _get_mc_luts()
|
|
|
|
# Apply algorithm
|
|
func = _marching_cubes_lewiner_cy.marching_cubes
|
|
vertices, faces , normals, values = func(volume, level, L, step_size, use_classic)
|
|
|
|
if not len(vertices):
|
|
raise RuntimeError('No surface found at the given iso value.')
|
|
|
|
# Output in z-y-x order, as is common in skimage
|
|
vertices = np.fliplr(vertices)
|
|
normals = np.fliplr(normals)
|
|
|
|
# Finishing touches to output
|
|
faces.shape = -1, 3
|
|
if gradient_direction == 'descent':
|
|
# MC implementation is right-handed, but gradient_direction is left-handed
|
|
faces = np.fliplr(faces)
|
|
elif not gradient_direction == 'ascent':
|
|
raise ValueError("Incorrect input %s in `gradient_direction`, see "
|
|
"docstring." % (gradient_direction))
|
|
if spacing != (1, 1, 1):
|
|
vertices = vertices * np.r_[spacing]
|
|
|
|
if allow_degenerate:
|
|
return vertices, faces, normals, values
|
|
else:
|
|
fun = _marching_cubes_lewiner_cy.remove_degenerate_faces
|
|
return fun(vertices, faces, normals, values)
|
|
|
|
|
|
def _to_array(args):
|
|
shape, text = args
|
|
byts = base64decode(text.encode('utf-8'))
|
|
ar = np.frombuffer(byts, dtype='int8')
|
|
ar.shape = shape
|
|
return ar
|
|
|
|
|
|
# Map an edge-index to two relative pixel positions. The ege index
|
|
# represents a point that lies somewhere in between these pixels.
|
|
# Linear interpolation should be used to determine where it is exactly.
|
|
# 0
|
|
# 3 1 -> 0x
|
|
# 2 xx
|
|
EDGETORELATIVEPOSX = np.array([ [0,1],[1,1],[1,0],[0,0], [0,1],[1,1],[1,0],[0,0], [0,0],[1,1],[1,1],[0,0] ], 'int8')
|
|
EDGETORELATIVEPOSY = np.array([ [0,0],[0,1],[1,1],[1,0], [0,0],[0,1],[1,1],[1,0], [0,0],[0,0],[1,1],[1,1] ], 'int8')
|
|
EDGETORELATIVEPOSZ = np.array([ [0,0],[0,0],[0,0],[0,0], [1,1],[1,1],[1,1],[1,1], [0,1],[0,1],[0,1],[0,1] ], 'int8')
|
|
|
|
|
|
def _get_mc_luts():
|
|
""" Kind of lazy obtaining of the luts.
|
|
"""
|
|
if not hasattr(mcluts, 'THE_LUTS'):
|
|
|
|
mcluts.THE_LUTS = _marching_cubes_lewiner_cy.LutProvider(
|
|
EDGETORELATIVEPOSX, EDGETORELATIVEPOSY, EDGETORELATIVEPOSZ,
|
|
|
|
_to_array(mcluts.CASESCLASSIC), _to_array(mcluts.CASES),
|
|
|
|
_to_array(mcluts.TILING1), _to_array(mcluts.TILING2), _to_array(mcluts.TILING3_1), _to_array(mcluts.TILING3_2),
|
|
_to_array(mcluts.TILING4_1), _to_array(mcluts.TILING4_2), _to_array(mcluts.TILING5), _to_array(mcluts.TILING6_1_1),
|
|
_to_array(mcluts.TILING6_1_2), _to_array(mcluts.TILING6_2), _to_array(mcluts.TILING7_1),
|
|
_to_array(mcluts.TILING7_2), _to_array(mcluts.TILING7_3), _to_array(mcluts.TILING7_4_1),
|
|
_to_array(mcluts.TILING7_4_2), _to_array(mcluts.TILING8), _to_array(mcluts.TILING9),
|
|
_to_array(mcluts.TILING10_1_1), _to_array(mcluts.TILING10_1_1_), _to_array(mcluts.TILING10_1_2),
|
|
_to_array(mcluts.TILING10_2), _to_array(mcluts.TILING10_2_), _to_array(mcluts.TILING11),
|
|
_to_array(mcluts.TILING12_1_1), _to_array(mcluts.TILING12_1_1_), _to_array(mcluts.TILING12_1_2),
|
|
_to_array(mcluts.TILING12_2), _to_array(mcluts.TILING12_2_), _to_array(mcluts.TILING13_1),
|
|
_to_array(mcluts.TILING13_1_), _to_array(mcluts.TILING13_2), _to_array(mcluts.TILING13_2_),
|
|
_to_array(mcluts.TILING13_3), _to_array(mcluts.TILING13_3_), _to_array(mcluts.TILING13_4),
|
|
_to_array(mcluts.TILING13_5_1), _to_array(mcluts.TILING13_5_2), _to_array(mcluts.TILING14),
|
|
|
|
_to_array(mcluts.TEST3), _to_array(mcluts.TEST4), _to_array(mcluts.TEST6),
|
|
_to_array(mcluts.TEST7), _to_array(mcluts.TEST10), _to_array(mcluts.TEST12),
|
|
_to_array(mcluts.TEST13), _to_array(mcluts.SUBCONFIG13),
|
|
)
|
|
|
|
return mcluts.THE_LUTS
|
|
|