address Stefans comments

This commit is contained in:
Almar Klein
2016-05-20 13:45:12 +02:00
parent f44e55b36d
commit 41857a0736
8 changed files with 62 additions and 48 deletions
+2 -1
View File
@@ -35,7 +35,8 @@
and more.
- Almar Klein
Binary heap class for graph algorithms
Binary heap class and other improvements for graph algorithms
Lewiner variant of marching cubes algorithm
- Lee Kamentsky and Thouis Jones of the CellProfiler team, Broad Institute, MIT
Constant time per pixel median filter, edge detectors, and more.
+2 -2
View File
@@ -109,9 +109,9 @@ def marching_cubes(volume, level=None, spacing=(1., 1., 1.),
See Also
--------
skimage.measure.correct_mesh_orientation
skimage.measure.marching_cubes_lewiner
skimage.measure.mesh_surface_area
"""
# Check inputs and ensure `volume` is C-contiguous for memoryviews
if volume.ndim != 3:
+36 -29
View File
@@ -24,8 +24,9 @@ def marching_cubes_lewiner(volume, level=None, spacing=(1., 1., 1.),
Parameters
----------
volume : (M, N, P) array (the data is internally converted to
float32 if necessary)
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.
@@ -73,18 +74,24 @@ def marching_cubes_lewiner(volume, level=None, spacing=(1., 1., 1.),
Notes about the algorithm
-------------------------
This is an implementation of:
Efficient implementation of Marching Cubes' cases with
topological guarantees. Thomas Lewiner, Helio Lopes, Antonio
Wilson Vieira and Geovan Tavares. Journal of Graphics Tools
8(2): pp. 1-15 (december 2003)
The algorithm is an improved version of Chernyaev's Marching Cubes 33
The algorithm [1] is an improved version of Chernyaev's Marching Cubes 33
algorithm, originally written in C++. It is an efficient algorithm
that relies on heavy use of lookup tables to handle the many different
cases. This keeps the algorithm relatively easy. The current algorithm
is a port of Lewiner's algorithm and written in Cython.
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).
See Also
--------
skimage.measure.marching_cubes
skimage.measure.mesh_surface_area
"""
# Check volume and ensure its in the format that the alg needs
@@ -113,7 +120,7 @@ def marching_cubes_lewiner(volume, level=None, spacing=(1., 1., 1.),
use_classic = bool(use_classic)
# Get LutProvider class (reuse if possible)
L = _getMCLuts()
L = _get_mc_luts()
# Apply algorithm
func = _marching_cubes_lewiner_cy.marching_cubes
@@ -144,7 +151,7 @@ def marching_cubes_lewiner(volume, level=None, spacing=(1., 1., 1.),
return fun(vertices, faces, normals, values)
def _toArray(args):
def _to_array(args):
shape, text = args
byts = base64decode(text.encode('utf-8'))
ar = np.frombuffer(byts, dtype='int8')
@@ -163,7 +170,7 @@ EDGETORELATIVEPOSY = np.array([ [0,0],[0,1],[1,1],[1,0], [0,0],[0,1],[1,1],[1,0]
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 _getMCLuts():
def _get_mc_luts():
""" Kind of lazy obtaining of the luts.
"""
if not hasattr(mcluts, 'THE_LUTS'):
@@ -171,24 +178,24 @@ def _getMCLuts():
mcluts.THE_LUTS = _marching_cubes_lewiner_cy.LutProvider(
EDGETORELATIVEPOSX, EDGETORELATIVEPOSY, EDGETORELATIVEPOSZ,
_toArray(mcluts.CASESCLASSIC), _toArray(mcluts.CASES),
_to_array(mcluts.CASESCLASSIC), _to_array(mcluts.CASES),
_toArray(mcluts.TILING1), _toArray(mcluts.TILING2), _toArray(mcluts.TILING3_1), _toArray(mcluts.TILING3_2),
_toArray(mcluts.TILING4_1), _toArray(mcluts.TILING4_2), _toArray(mcluts.TILING5), _toArray(mcluts.TILING6_1_1),
_toArray(mcluts.TILING6_1_2), _toArray(mcluts.TILING6_2), _toArray(mcluts.TILING7_1),
_toArray(mcluts.TILING7_2), _toArray(mcluts.TILING7_3), _toArray(mcluts.TILING7_4_1),
_toArray(mcluts.TILING7_4_2), _toArray(mcluts.TILING8), _toArray(mcluts.TILING9),
_toArray(mcluts.TILING10_1_1), _toArray(mcluts.TILING10_1_1_), _toArray(mcluts.TILING10_1_2),
_toArray(mcluts.TILING10_2), _toArray(mcluts.TILING10_2_), _toArray(mcluts.TILING11),
_toArray(mcluts.TILING12_1_1), _toArray(mcluts.TILING12_1_1_), _toArray(mcluts.TILING12_1_2),
_toArray(mcluts.TILING12_2), _toArray(mcluts.TILING12_2_), _toArray(mcluts.TILING13_1),
_toArray(mcluts.TILING13_1_), _toArray(mcluts.TILING13_2), _toArray(mcluts.TILING13_2_),
_toArray(mcluts.TILING13_3), _toArray(mcluts.TILING13_3_), _toArray(mcluts.TILING13_4),
_toArray(mcluts.TILING13_5_1), _toArray(mcluts.TILING13_5_2), _toArray(mcluts.TILING14),
_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),
_toArray(mcluts.TEST3), _toArray(mcluts.TEST4), _toArray(mcluts.TEST6),
_toArray(mcluts.TEST7), _toArray(mcluts.TEST10), _toArray(mcluts.TEST12),
_toArray(mcluts.TEST13), _toArray(mcluts.SUBCONFIG13),
_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
@@ -1,6 +1,4 @@
# -*- coding: utf-8 -*-
# Copyright (C) 2012, Almar Klein
# Copyright (C) 2002, Thomas Lewiner
#
#cython: cdivision=True
#cython: boundscheck=False
@@ -14,11 +12,12 @@ Efficient implementation of Marching Cubes' cases with topological guarantees.
Thomas Lewiner, Helio Lopes, Antonio Wilson Vieira and Geovan Tavares.
Journal of Graphics Tools 8(2): pp. 1-15 (december 2003)
I selected this algorithm because it provides topologically correct results,
and because the algorithms implementation is relatively simple. Most of
the magic is in the lookup tables, which are provided as open source.
This algorithm has the advantage that it provides topologically correct
results, and the algorithms implementation is relatively simple. Most
of the magic is in the lookup tables, which are provided as open source.
This code is distributed under the terms of the (new) BSD License.
Originally implemented in C++ by Thomas Lewiner in 2002, ported to Cython
by Almar Klein in 2012. Adapted for scikit-image in 2016.
"""
@@ -1,8 +1,7 @@
# -*- coding: utf-8 -*-
# Copyright (C) 2012, Almar Klein
# Copyright (C) 2002, Thomas Lewiner
# This file was auto-generated from LookUpTable.h by createluts.py.
# The luts are Copyright (C) 2002 by Thomas Lewiner
#static const char casesClassic[256][16]
CASESCLASSIC = (256, 16), """
-1
View File
@@ -1,5 +1,4 @@
# -*- coding: utf-8 -*-
# Copyright (C) 2012, Almar Klein
""" Create lookup tables for the marching cubes algorithm, by parsing
the file "LookUpTable.h". This prints a text to the stdout wich
+1 -1
View File
@@ -13,7 +13,7 @@ from skimage.draw import ellipsoid
# Create test volume
SELECT = 4
SELECT = 3
gradient_dir = 'descent' # ascent or descent
if SELECT == 1:
+15 -6
View File
@@ -111,15 +111,21 @@ def test_both_algs_same_result_ellipse():
assert _same_mesh(vertices1, faces1, vertices3, faces3)
def _same_mesh(vertices1, faces1, vertices2, faces2):
rounder = lambda x: int(x*1000)/1000 # to take into account small variations
def _same_mesh(vertices1, faces1, vertices2, faces2, tol=1e-10):
""" Compare two meshes, using a certain tolerance and invariant to
the order of the faces.
"""
# Unwind vertices
triangles1 = vertices1[np.array(faces1)]
triangles2 = vertices2[np.array(faces2)]
# Sort vertices within each triangle
triang1 = [np.concatenate(sorted(t, key=lambda x:tuple(x))) for t in triangles1]
triang1 = set([tuple([rounder(i) for i in t]) for t in triang1])
triang2 = [np.concatenate(sorted(t, key=lambda x:tuple(x))) for t in triangles2]
triang2 = set([tuple([rounder(i) for i in t]) for t in triang2])
return triang1 == triang2
# Sort the resulting 9-element "tuples"
triang1 = np.array(sorted([tuple(x) for x in triang1]))
triang2 = np.array(sorted([tuple(x) for x in triang2]))
return triang1.shape == triang2.shape and np.allclose(triang1, triang2, 0, tol)
def test_both_algs_same_result_donut():
# Performing this test on data that does not have ambiguities
@@ -145,9 +151,12 @@ def test_both_algs_same_result_donut():
vertices2, faces2, *_ = marching_cubes_lewiner(vol, 0)
vertices3, faces3, *_ = marching_cubes_lewiner(vol, 0, use_classic=True)
# Old and new alg are different
assert not _same_mesh(vertices1, faces1, vertices2, faces2)
#assert _same_mesh(vertices1, faces1, vertices3, faces3) # would have been nice
# New classic and new Lewiner are different
assert not _same_mesh(vertices2, faces2, vertices3, faces3)
# Would have been nice if old and new classic would have been the same
# assert _same_mesh(vertices1, faces1, vertices3, faces3, 5)
if __name__ == '__main__':