From a90f875d15333230ad3bfbd7f526af9a889513d5 Mon Sep 17 00:00:00 2001 From: Almar Klein Date: Mon, 25 Jul 2016 12:25:57 +0200 Subject: [PATCH] Rename MC algs to make lewiner the default --- bento.info | 4 +- skimage/measure/__init__.py | 8 ++- ...ng_cubes.py => _marching_cubes_classic.py} | 40 +++++------ ..._cy.pyx => _marching_cubes_classic_cy.pyx} | 0 skimage/measure/_marching_cubes_lewiner.py | 67 +++++++++++++++++-- skimage/measure/mc_meta/visual_test.py | 4 +- skimage/measure/setup.py | 6 +- skimage/measure/tests/test_marching_cubes.py | 49 +++++++++++--- 8 files changed, 132 insertions(+), 46 deletions(-) rename skimage/measure/{_marching_cubes.py => _marching_cubes_classic.py} (91%) rename skimage/measure/{_marching_cubes_cy.pyx => _marching_cubes_classic_cy.pyx} (100%) diff --git a/bento.info b/bento.info index 5597045c..53ef3e08 100644 --- a/bento.info +++ b/bento.info @@ -50,9 +50,9 @@ Library: Extension: skimage.measure._moments_cy Sources: skimage/measure/_moments_cy.pyx - Extension: skimage.measure._marching_cubes_cy + Extension: skimage.measure._marching_cubes_classic_cy Sources: - skimage/measure/_marching_cubes_cy.pyx + skimage/measure/_marching_cubes_classic_cy.pyx Extension: skimage.measure._marching_cubes_lewiner_cy Sources: skimage/measure/_marching_cubes_lewiner_cy.pyx diff --git a/skimage/measure/__init__.py b/skimage/measure/__init__.py index c03587d4..015f67d5 100755 --- a/skimage/measure/__init__.py +++ b/skimage/measure/__init__.py @@ -1,7 +1,8 @@ from ._find_contours import find_contours -from ._marching_cubes import (marching_cubes, mesh_surface_area, - correct_mesh_orientation) -from ._marching_cubes_lewiner import marching_cubes_lewiner +from ._marching_cubes_lewiner import marching_cubes, marching_cubes_lewiner +from ._marching_cubes_classic import (marching_cubes_classic, + mesh_surface_area, + correct_mesh_orientation) from ._regionprops import regionprops, perimeter from .simple_metrics import compare_mse, compare_nrmse, compare_psnr from ._structural_similarity import compare_ssim, structural_similarity @@ -31,6 +32,7 @@ __all__ = ['find_contours', 'moments_hu', 'marching_cubes', 'marching_cubes_lewiner', + 'marching_cubes_classic', 'mesh_surface_area', 'correct_mesh_orientation', 'profile_line', diff --git a/skimage/measure/_marching_cubes.py b/skimage/measure/_marching_cubes_classic.py similarity index 91% rename from skimage/measure/_marching_cubes.py rename to skimage/measure/_marching_cubes_classic.py index fa91de07..b6e77c5d 100644 --- a/skimage/measure/_marching_cubes.py +++ b/skimage/measure/_marching_cubes_classic.py @@ -1,15 +1,15 @@ import numpy as np import scipy.ndimage as ndi from .._shared.utils import warn -from . import _marching_cubes_cy +from . import _marching_cubes_classic_cy -def marching_cubes(volume, level=None, spacing=(1., 1., 1.), - gradient_direction='descent'): +def marching_cubes_classic(volume, level=None, spacing=(1., 1., 1.), + gradient_direction='descent'): """ Classic marching cubes algorithm to find surfaces in 3d volumetric data. - Note that the ``marching_cubes_lewiner()`` algorithm is recommended over + Note that the ``marching_cubes()`` algorithm is recommended over this algorithm, because it's faster and produces better results. Parameters @@ -90,7 +90,7 @@ def marching_cubes(volume, level=None, spacing=(1., 1., 1.), named `myvolume` about the level 0.0, using the ``mayavi`` package:: >>> from mayavi import mlab # doctest: +SKIP - >>> verts, faces = marching_cubes(myvolume, 0.0) # doctest: +SKIP + >>> verts, faces = marching_cubes_classic(myvolume, 0.0) # doctest: +SKIP >>> mlab.triangular_mesh([vert[0] for vert in verts], ... [vert[1] for vert in verts], ... [vert[2] for vert in verts], @@ -100,7 +100,7 @@ def marching_cubes(volume, level=None, spacing=(1., 1., 1.), Similarly using the ``visvis`` package:: >>> import visvis as vv # doctest: +SKIP - >>> verts, faces = marching_cubes(myvolume, 0.0) # doctest: +SKIP + >>> verts, faces = marching_cubes_classic(myvolume, 0.0) # doctest: +SKIP >>> vv.mesh(np.fliplr(verts), faces) # doctest: +SKIP >>> vv.use().Run() # doctest: +SKIP @@ -113,7 +113,7 @@ def marching_cubes(volume, level=None, spacing=(1., 1., 1.), See Also -------- - skimage.measure.marching_cubes_lewiner + skimage.measure.marching_cubes skimage.measure.mesh_surface_area """ # Check inputs and ensure `volume` is C-contiguous for memoryviews @@ -136,11 +136,12 @@ def marching_cubes(volume, level=None, 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_faces = _marching_cubes_cy.iterate_and_store_3d(volume, float(level)) + raw_faces = _marching_cubes_classic_cy.iterate_and_store_3d(volume, + float(level)) # 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_faces) + verts, faces = _marching_cubes_classic_cy.unpack_unique_verts(raw_faces) verts = np.asarray(verts) faces = np.asarray(faces) @@ -172,7 +173,7 @@ def mesh_surface_area(verts, faces): Notes ----- - The arguments expected by this function are the exact outputs from + The arguments expected by this function are the first two outputs from `skimage.measure.marching_cubes`. For unit correct output, ensure correct `spacing` was passed to `skimage.measure.marching_cubes`. @@ -182,6 +183,7 @@ def mesh_surface_area(verts, faces): See Also -------- skimage.measure.marching_cubes + skimage.measure.marching_cubes_classic skimage.measure.correct_mesh_orientation """ @@ -229,7 +231,7 @@ def correct_mesh_orientation(volume, verts, faces, spacing=(1., 1., 1.), 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 + `skimage.measure.marching_cubes_classic` by flipping the orientation of mis-oriented faces. Because marching cubes could be used to find isosurfaces either on @@ -240,21 +242,21 @@ def correct_mesh_orientation(volume, verts, faces, spacing=(1., 1., 1.), 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 + `skimage.measure.marching_cubes_classic`. 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.marching_cubes_classic skimage.measure.mesh_surface_area """ warn(DeprecationWarning("`correct_mesh_orientation` is deprecated for " - "removal as `marching_cubes` now guarantess " - "correct mesh orientation.")) + "removal as `marching_cubes`_classic now " + "guarantees correct mesh orientation.")) verts = verts.copy() verts[:, 0] /= spacing[0] @@ -303,7 +305,7 @@ def _correct_mesh_orientation(volume, actual_verts, faces, 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 + `skimage.measure.marching_cubes_classic` by flipping the orientation of mis-oriented faces. Because marching cubes could be used to find isosurfaces either on @@ -314,7 +316,7 @@ def _correct_mesh_orientation(volume, actual_verts, faces, completely incorrectly, try changing this option. The arguments expected by this function are the exact outputs from - `skimage.measure.marching_cubes` except `actual_verts`, which is an + `skimage.measure.marching_cubes_classic` except `actual_verts`, which is an uncorrected version of the fancy indexing operation `verts[faces]`. Only `faces` is corrected and returned as the vertices do not change, only the order in which they are referenced. @@ -323,7 +325,7 @@ def _correct_mesh_orientation(volume, actual_verts, faces, See Also -------- - skimage.measure.marching_cubes + skimage.measure.marching_cubes_classic skimage.measure.mesh_surface_area """ diff --git a/skimage/measure/_marching_cubes_cy.pyx b/skimage/measure/_marching_cubes_classic_cy.pyx similarity index 100% rename from skimage/measure/_marching_cubes_cy.pyx rename to skimage/measure/_marching_cubes_classic_cy.pyx diff --git a/skimage/measure/_marching_cubes_lewiner.py b/skimage/measure/_marching_cubes_lewiner.py index 2c6f1fd6..067396db 100644 --- a/skimage/measure/_marching_cubes_lewiner.py +++ b/skimage/measure/_marching_cubes_lewiner.py @@ -1,24 +1,55 @@ 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 marching_cubes_lewiner(volume, level=None, spacing=(1., 1., 1.), - gradient_direction='descent', step_size=1, - allow_degenerate=True, use_classic=False): +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()``, this algorithm is faster, + 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. @@ -54,8 +85,8 @@ def marching_cubes_lewiner(volume, level=None, spacing=(1., 1., 1.), 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()`` - function. + this option are *not* generally the same as the + ``marching_cubes_classic()`` function. Returns ------- @@ -92,10 +123,32 @@ def marching_cubes_lewiner(volume, level=None, spacing=(1., 1., 1.), See Also -------- - skimage.measure.marching_cubes + 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.') diff --git a/skimage/measure/mc_meta/visual_test.py b/skimage/measure/mc_meta/visual_test.py index bc95758a..64504fd9 100644 --- a/skimage/measure/mc_meta/visual_test.py +++ b/skimage/measure/mc_meta/visual_test.py @@ -8,7 +8,7 @@ import time import numpy as np import visvis as vv -from skimage.measure import marching_cubes, marching_cubes_lewiner +from skimage.measure import marching_cubes_classic, marching_cubes_lewiner from skimage.draw import ellipsoid @@ -57,7 +57,7 @@ vertices1, faces1, *_ = marching_cubes_lewiner(vol, isovalue, gradient_direction print('finding surface lewiner took %1.0f ms' % (1000*(time.time()-t0)) ) t0 = time.time() -vertices2, faces2, *_ = marching_cubes(vol, isovalue, gradient_direction=gradient_dir) +vertices2, faces2, *_ = marching_cubes_classic(vol, isovalue, gradient_direction=gradient_dir) print('finding surface classic took %1.0f ms' % (1000*(time.time()-t0)) ) # Show diff --git a/skimage/measure/setup.py b/skimage/measure/setup.py index cea657c5..df69ff5d 100644 --- a/skimage/measure/setup.py +++ b/skimage/measure/setup.py @@ -15,7 +15,7 @@ def configuration(parent_package='', top_path=None): cython(['_ccomp.pyx'], working_path=base_path) cython(['_find_contours_cy.pyx'], working_path=base_path) cython(['_moments_cy.pyx'], working_path=base_path) - cython(['_marching_cubes_cy.pyx'], working_path=base_path) + cython(['_marching_cubes_classic_cy.pyx'], working_path=base_path) cython(['_marching_cubes_lewiner_cy.pyx'], working_path=base_path) cython(['_pnpoly.pyx'], working_path=base_path) @@ -25,8 +25,8 @@ def configuration(parent_package='', top_path=None): include_dirs=[get_numpy_include_dirs()]) config.add_extension('_moments_cy', sources=['_moments_cy.c'], include_dirs=[get_numpy_include_dirs()]) - config.add_extension('_marching_cubes_cy', - sources=['_marching_cubes_cy.c'], + config.add_extension('_marching_cubes_classic_cy', + sources=['_marching_cubes_classic_cy.c'], include_dirs=[get_numpy_include_dirs()]) config.add_extension('_marching_cubes_lewiner_cy', sources=['_marching_cubes_lewiner_cy.c'], diff --git a/skimage/measure/tests/test_marching_cubes.py b/skimage/measure/tests/test_marching_cubes.py index 2d0c766d..3fa85da3 100644 --- a/skimage/measure/tests/test_marching_cubes.py +++ b/skimage/measure/tests/test_marching_cubes.py @@ -1,9 +1,38 @@ +import sys import numpy as np from numpy.testing import assert_raises from skimage.draw import ellipsoid, ellipsoid_stats -from skimage.measure import (marching_cubes, marching_cubes_lewiner, +from skimage.measure import (marching_cubes, + marching_cubes_classic, marching_cubes_lewiner, mesh_surface_area, correct_mesh_orientation) +from skimage.measure._marching_cubes_lewiner import expected_output_args + + +def test_expected_output_args(): + + res = [] + + def foo(): + nout = expected_output_args() + print(nout) + res.append(nout) + return [nout] * int(nout) + + foo() + a = foo() + a, b = foo() + a, b, c = foo() + assert res == [0, 1, 2, 3] or res == [0, 0, 2, 3] + # ``a = foo()`` somehow yields 0 in test, which is ok for us; + # we only want to distinguish between > 2 args or not + + if sys.version_info >= (3, 3): + res = [] + exec('*a, b, c = foo()') + exec('a, b, c, *d = foo()') + exec('a, b, *c, d, e = foo()') + assert res == [2.1, 3.1, 4.1] def test_marching_cubes_isotropic(): @@ -11,7 +40,7 @@ def test_marching_cubes_isotropic(): _, surf = ellipsoid_stats(6, 10, 16) # Classic - verts, faces = marching_cubes(ellipsoid_isotropic, 0.) + verts, faces = marching_cubes_classic(ellipsoid_isotropic, 0.) surf_calc = mesh_surface_area(verts, faces) # Test within 1% tolerance for isotropic. Will always underestimate. assert surf > surf_calc and surf_calc > surf * 0.99 @@ -30,8 +59,8 @@ def test_marching_cubes_anisotropic(): _, surf = ellipsoid_stats(6, 10, 16) # Classic - verts, faces = marching_cubes(ellipsoid_anisotropic, 0., - spacing=spacing) + verts, faces = marching_cubes_classic(ellipsoid_anisotropic, 0., + spacing=spacing) surf_calc = mesh_surface_area(verts, faces) # Test within 1.5% tolerance for anisotropic. Will always underestimate. assert surf > surf_calc and surf_calc > surf * 0.985 @@ -45,11 +74,11 @@ def test_marching_cubes_anisotropic(): def test_invalid_input(): # Classic - assert_raises(ValueError, marching_cubes, np.zeros((2, 2, 1)), 0) - assert_raises(ValueError, marching_cubes, np.zeros((2, 2, 1)), 1) - assert_raises(ValueError, marching_cubes, np.ones((3, 3, 3)), 1, + assert_raises(ValueError, marching_cubes_classic, np.zeros((2, 2, 1)), 0) + assert_raises(ValueError, marching_cubes_classic, np.zeros((2, 2, 1)), 1) + assert_raises(ValueError, marching_cubes_classic, np.ones((3, 3, 3)), 1, spacing=(1, 2)) - assert_raises(ValueError, marching_cubes, np.zeros((20, 20)), 0) + assert_raises(ValueError, marching_cubes_classic, np.zeros((20, 20)), 0) # Lewiner assert_raises(ValueError, marching_cubes_lewiner, np.zeros((2, 2, 1)), 0) @@ -101,7 +130,7 @@ def test_both_algs_same_result_ellipse(): sphere_small = ellipsoid(1, 1, 1, levelset=True) - vertices1, faces1 = marching_cubes(sphere_small, 0)[:2] + vertices1, faces1 = marching_cubes_classic(sphere_small, 0)[:2] vertices2, faces2 = marching_cubes_lewiner(sphere_small, 0, allow_degenerate=False)[:2] vertices3, faces3 = marching_cubes_lewiner(sphere_small, 0, allow_degenerate=False, use_classic=True)[:2] @@ -147,7 +176,7 @@ def test_both_algs_same_result_donut(): 64 * ( ((8*y-2)+4)*((8*y-2)+4) + (8*z)**2 ) ) + 1025 - vertices1, faces1 = marching_cubes(vol, 0)[:2] + vertices1, faces1 = marching_cubes_classic(vol, 0)[:2] vertices2, faces2 = marching_cubes_lewiner(vol, 0)[:2] vertices3, faces3 = marching_cubes_lewiner(vol, 0, use_classic=True)[:2]