marching cubes tests

skimage/measure/tests/test_marching_cubes.py

@@ -2,16 +2,23 @@ import numpy as np
 from numpy.testing import assert_raises
 
 from skimage.draw import ellipsoid, ellipsoid_stats
-from skimage.measure import (marching_cubes, mesh_surface_area,
-                             correct_mesh_orientation)
+from skimage.measure import (marching_cubes, marching_cubes_lewiner, 
+                             mesh_surface_area, correct_mesh_orientation)
 
 
 def test_marching_cubes_isotropic():
     ellipsoid_isotropic = ellipsoid(6, 10, 16, levelset=True)
     _, surf = ellipsoid_stats(6, 10, 16)
+    
+    # Classic
     verts, faces = marching_cubes(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
+    
+    # Lewiner
+    verts, faces, *_ = marching_cubes_lewiner(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
 
@@ -21,20 +28,36 @@ def test_marching_cubes_anisotropic():
     ellipsoid_anisotropic = ellipsoid(6, 10, 16, spacing=spacing,
                                       levelset=True)
     _, surf = ellipsoid_stats(6, 10, 16)
+    
+    # Classic
     verts, faces = marching_cubes(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
+    
+    # Lewiner
+    verts, faces, *_ = marching_cubes_lewiner(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
 
 
 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,
                   spacing=(1, 2))
     assert_raises(ValueError, marching_cubes, np.zeros((20, 20)), 0)
+    
+    # Lewiner
+    assert_raises(ValueError, marching_cubes_lewiner, np.zeros((2, 2, 1)), 0)
+    assert_raises(ValueError, marching_cubes_lewiner, np.zeros((2, 2, 1)), 1)
+    assert_raises(ValueError, marching_cubes_lewiner, np.ones((3, 3, 3)), 1,
+                  spacing=(1, 2))
+    assert_raises(ValueError, marching_cubes_lewiner, np.zeros((20, 20)), 0)
 
 
 def test_correct_mesh_orientation():
@@ -74,5 +97,58 @@ def test_correct_mesh_orientation():
     np.testing.assert_array_equal(expected, corrected_faces1)
 
 
+def test_both_algs_same_result_ellipse():
+    # Performing this test on data that does not have ambiguities
+    
+    sphere_small = ellipsoid(1, 1, 1, levelset=True)
+    
+    vertices1, faces1, *_ = marching_cubes(sphere_small, 0)
+    vertices2, faces2, *_ = marching_cubes_lewiner(sphere_small, 0, allow_degenerate=False)
+    vertices3, faces3, *_ = marching_cubes_lewiner(sphere_small, 0, allow_degenerate=False, use_classic=True)
+    
+    # Order id different, best we can do is test equal shape and same vertices present
+    assert _same_mesh(vertices1, faces1, vertices2, faces2)
+    assert _same_mesh(vertices1, faces1, vertices3, faces3)
+
+
+def _same_mesh(vertices1, faces1, vertices2, faces2):
+    faces1 = [sorted(f) for f in faces1]
+    faces2 = [sorted(f) for f in faces2]
+    triangles1 = vertices1[np.array(faces1)]
+    triangles2 = vertices2[np.array(faces2)]
+    triang1 = set([tuple(t.flat) for t in triangles1])
+    triang2 = set([tuple(t.flat) for t in triangles2])
+    return triang1 == triang2
+
+
+def test_both_algs_same_result_donut():
+    # Performing this test on data that does not have ambiguities
+    
+    n = 48
+    a, b = 2.5/n, -1.25
+    isovalue = 0.0
+    #
+    vol = np.empty((n,n,n), 'float32')
+    for iz in range(vol.shape[0]):
+        for iy in range(vol.shape[1]):
+            for ix in range(vol.shape[2]):
+                z, y, x = float(iz)*a+b, float(iy)*a+b, float(ix)*a+b
+                vol[iz,iy,ix] = ( ( 
+                    (8*x)**2 + (8*y-2)**2 + (8*z)**2 + 16 - 1.85*1.85 ) * ( (8*x)**2 +
+                    (8*y-2)**2 + (8*z)**2 + 16 - 1.85*1.85 ) - 64 * ( (8*x)**2 + (8*y-2)**2 )
+                    ) * ( ( (8*x)**2 + ((8*y-2)+4)*((8*y-2)+4) + (8*z)**2 + 16 - 1.85*1.85 )
+                    * ( (8*x)**2 + ((8*y-2)+4)*((8*y-2)+4) + (8*z)**2 + 16 - 1.85*1.85 ) -
+                    64 * ( ((8*y-2)+4)*((8*y-2)+4) + (8*z)**2 
+                    ) ) + 1025
+    
+    vertices1, faces1, *_ = marching_cubes(vol, 0)
+    vertices2, faces2, *_ = marching_cubes_lewiner(vol, 0)
+    vertices3, faces3, *_ = marching_cubes_lewiner(vol, 0, use_classic=True)
+    
+    assert not _same_mesh(vertices1, faces1, vertices2, faces2)
+    #assert _same_mesh(vertices1, faces1, vertices3, faces3)  # would have been nice
+    assert not _same_mesh(vertices2, faces2, vertices3, faces3)
+
+
 if __name__ == '__main__':
     np.testing.run_module_suite()