diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt index 44214d85..e4c852e6 100644 --- a/CONTRIBUTORS.txt +++ b/CONTRIBUTORS.txt @@ -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. diff --git a/skimage/measure/_marching_cubes.py b/skimage/measure/_marching_cubes.py index 09332d1c..476f4e39 100644 --- a/skimage/measure/_marching_cubes.py +++ b/skimage/measure/_marching_cubes.py @@ -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: diff --git a/skimage/measure/_marching_cubes_lewiner.py b/skimage/measure/_marching_cubes_lewiner.py index f105b4b3..7f8339e5 100644 --- a/skimage/measure/_marching_cubes_lewiner.py +++ b/skimage/measure/_marching_cubes_lewiner.py @@ -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 diff --git a/skimage/measure/_marching_cubes_lewiner_cy.pyx b/skimage/measure/_marching_cubes_lewiner_cy.pyx index d1aa2280..b37fb718 100644 --- a/skimage/measure/_marching_cubes_lewiner_cy.pyx +++ b/skimage/measure/_marching_cubes_lewiner_cy.pyx @@ -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. """ diff --git a/skimage/measure/_marching_cubes_lewiner_luts.py b/skimage/measure/_marching_cubes_lewiner_luts.py index 59243662..07fb5728 100644 --- a/skimage/measure/_marching_cubes_lewiner_luts.py +++ b/skimage/measure/_marching_cubes_lewiner_luts.py @@ -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), """ diff --git a/skimage/measure/mc_meta/createluts.py b/skimage/measure/mc_meta/createluts.py index c6be213c..670a27c8 100644 --- a/skimage/measure/mc_meta/createluts.py +++ b/skimage/measure/mc_meta/createluts.py @@ -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 diff --git a/skimage/measure/mc_meta/visual_test.py b/skimage/measure/mc_meta/visual_test.py index 7a8ea699..e8b7e717 100644 --- a/skimage/measure/mc_meta/visual_test.py +++ b/skimage/measure/mc_meta/visual_test.py @@ -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: diff --git a/skimage/measure/tests/test_marching_cubes.py b/skimage/measure/tests/test_marching_cubes.py index 4d62abe0..25262682 100644 --- a/skimage/measure/tests/test_marching_cubes.py +++ b/skimage/measure/tests/test_marching_cubes.py @@ -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__':