TST: skel3d: add more tests

skimage/morphology/__init__.py

@@ -6,7 +6,7 @@ from .selem import (square, rectangle, diamond, disk, cube, octahedron, ball,
                     octagon, star)
 from .watershed import watershed
 from ._skeletonize import skeletonize, medial_axis
-from ._skeletonize_3d import skeletonize_3d as compute_thin_image
+from ._skeletonize_3d import skeletonize_3d
 from .convex_hull import convex_hull_image, convex_hull_object
 from .greyreconstruct import reconstruction
 from .misc import remove_small_objects, remove_small_holes

skimage/morphology/_skeletonize_3d.py

@@ -1,10 +1,11 @@
 from __future__ import division, print_function, absolute_import
 
 import numpy as np
+from ..util import img_as_ubyte
 from ._skeletonize_3d_cy import _compute_thin_image
 
 
-def skeletonize_3d(img_in):
+def skeletonize_3d(img):
     """Compute the skeleton of a binary image.
 
     Thinning is used to reduce each connected component in a binary image
@@ -12,7 +13,7 @@ def skeletonize_3d(img_in):
 
     Parameters
     ----------
-    image : ndarray, 2D or 3D
+    img : ndarray, 2D or 3D
         A binary image containing the objects to be skeletonized. Zeros
         represent background, nonzero values are foreground.
 
@@ -33,17 +34,19 @@ def skeletonize_3d(img_in):
 
     """
     # make sure the image is 3D or 2D (if it is, temporarily upcast to 3D)
-    if img_in.ndim < 2 or img_in.ndim > 3:
-        raise ValueError('expect 2D, got ndim = %s' % img_in.ndim)
+    if img.ndim < 2 or img.ndim > 3:
+        raise ValueError('expect 2D, got ndim = %s' % img.ndim)
 
-    img = img_in.copy()
+    img = img_as_ubyte(img)
+    img = np.ascontiguousarray(img)
+
+    img = img.copy()
     if img.ndim == 2:
         img = img[None, ...]
 
     # normalize to binary
     maxval = img.max()
     img[img != 0] = 1
-    img = img.astype(np.uint8)
 
     # pad w/ zeros to simplify dealing w/ boundaries
     img_o = np.pad(img, pad_width=1, mode='constant')

skimage/morphology/tests/test_skeletonize_3d.py

@@ -1,14 +1,111 @@
 from __future__ import division, print_function, absolute_import
 
 import os
+import warnings
 
 import numpy as np
-from numpy.testing import assert_equal, run_module_suite
+from numpy.testing import (assert_equal, run_module_suite, assert_raises,
+                           assert_)
+
+import scipy.ndimage as ndi
 
 import skimage
-from skimage import io
+from skimage import io, draw, data_dir
+#from skimage import draw
+from skimage.util import img_as_ubyte
 
-from skimage.morphology import compute_thin_image
+from skimage.morphology import skeletonize_3d
+
+
+# basic behavior tests (mostly copied over from 2D skeletonize)
+
+def test_skeletonize_wrong_dim():
+    im = np.zeros(5, dtype=np.uint8)
+    assert_raises(ValueError, skeletonize_3d, im)
+
+    im = np.zeros((5, 5, 5, 5), dtype=np.uint8)
+    assert_raises(ValueError, skeletonize_3d, im)
+
+
+def test_skeletonize_no_foreground():
+    im = np.zeros((5, 5), dtype=np.uint8)
+    result = skeletonize_3d(im)
+    assert_equal(result, im)
+
+
+def test_skeletonize_all_foreground():
+    im = np.ones((3, 4), dtype=np.uint8)
+    assert_equal(skeletonize_3d(im),
+                 np.array([[0, 0, 0, 0],
+                           [1, 1, 1, 1],
+                           [0, 0, 0, 0]], dtype=np.uint8))
+
+
+def test_skeletonize_single_point():
+    im = np.zeros((5, 5), dtype=np.uint8)
+    im[3, 3] = 1
+    result = skeletonize_3d(im)
+    assert_equal(result, im)
+
+
+def test_skeletonize_already_thinned():
+    im = np.zeros((5, 5), dtype=np.uint8)
+    im[3, 1:-1] = 1
+    im[2, -1] = 1
+    im[4, 0] = 1
+    result = skeletonize_3d(im)
+    assert_equal(result, im)
+
+
+def test_dtype_conv():
+    # check that the operation does the right thing with floats etc
+    # also check non-contiguous input
+    img = np.random.random((16, 16))[::2, ::2]
+    img[img < 0.5] = 0
+
+    orig = img.copy()
+
+    with warnings.catch_warnings():
+        # UserWarning for possible precision loss, expected
+        warnings.simplefilter('ignore', UserWarning)
+        res = skeletonize_3d(img)
+
+    assert_equal(res.dtype, np.uint8)
+    assert_equal(img, orig)  # operation does not clobber the original 
+    assert_equal(res.max(),
+                 img_as_ubyte(img).max())    # the intensity range is preserved
+
+
+def test_skeletonize_num_neighbours():
+    # an empty image
+    image = np.zeros((300, 300))
+
+    # foreground object 1
+    image[10:-10, 10:100] = 1
+    image[-100:-10, 10:-10] = 1
+    image[10:-10, -100:-10] = 1
+
+    # foreground object 2
+    rs, cs = draw.line(250, 150, 10, 280)
+    for i in range(10):
+        image[rs + i, cs] = 1
+    rs, cs = draw.line(10, 150, 250, 280)
+    for i in range(20):
+        image[rs + i, cs] = 1
+
+    # foreground object 3
+    ir, ic = np.indices(image.shape)
+    circle1 = (ic - 135)**2 + (ir - 150)**2 < 30**2
+    circle2 = (ic - 135)**2 + (ir - 150)**2 < 20**2
+    image[circle1] = 1
+    image[circle2] = 0
+    result = skeletonize_3d(image)
+
+    # there should never be a 2x2 block of foreground pixels in a skeleton
+    mask = np.array([[1,  1],
+                     [1,  1]], np.uint8)
+    blocks = ndi.correlate(result, mask, mode='constant')
+    assert_(not np.any(blocks == 4))
 
 
 # nose test generators:
@@ -44,7 +141,7 @@ def check_skel(fname):
                      dtype=np.uint8)
 
     # compute
-    img1_2d = compute_thin_image(img)
+    img1_2d = skeletonize_3d(img)
 
     # and compare to FIJI
     img_f = np.loadtxt(os.path.join(get_data_path(), fname + '_fiji.txt'),
@@ -57,7 +154,7 @@ def check_skel_3d(fname):
     img = io.imread(os.path.join(get_data_path(), fname + '.tif'))
     img_f = io.imread(os.path.join(get_data_path(), fname + '_fiji.tif'))
 
-    img_s = compute_thin_image(img)
+    img_s = skeletonize_3d(img)
     assert_equal(img_s, img_f)