+ * More information at Skeletonize3D homepage: + * http://fiji.sc/Skeletonize3D + * + * @version 1.0 11/13/2015 (unique BSD licensed version for scikit-image) + * @author Ignacio Arganda-Carreras (iargandacarreras at gmail.com) + +References +---------- +.. [Lee94] T.-C. Lee, R.L. Kashyap and C.-N. Chu, Building skeleton models + via 3-D medial surface/axis thinning algorithms. + Computer Vision, Graphics, and Image Processing, 56(6):462-478, 1994. + +.. [IAC15] Ignacio Arganda-Carreras, 2015. Skeletonize3D plugin for ImageJ(C). + http://fiji.sc/Skeletonize3D + +""" +from __future__ import division, print_function, absolute_import + +from libc.string cimport memcpy + +import numpy as np +from numpy cimport npy_intp, npy_uint8, ndarray +cimport cython + +ctypedef npy_uint8 pixel_type + + +@cython.boundscheck(False) +@cython.wraparound(False) +def _compute_thin_image(pixel_type[:, :, ::1] img not None): + """Compute a thin image. + + Loop through the image multiple times, removing "simple" points, i.e. + those point which can be removed without changing local connectivity in the + 3x3x3 neighborhood of a point. + + This routine implements the two-pass algorithm of [Lee94]_. Namely, + for each of the six border types (positive and negative x-, y- and z-), + the algorithm first collects all possibly deletable points, and then + performs a sequential rechecking. + + The input, `img`, is assumed to be a 3D binary image in the + (p, r, c) format [i.e., C ordered array], filled by zeros (background) and + ones. Furthermore, `img` is assumed to be padded by zeros from all + directions --- this way the zero boundary conditions are automatic + and there is need to guard against out-of-bounds access. + + """ + cdef: + int unchanged_borders = 0, curr_border, num_borders + int borders[6] + npy_intp p, r, c + bint no_change + + list simple_border_points + Py_ssize_t num_border_points, i, j + (npy_intp, npy_intp, npy_intp) point + + pixel_type neighb[27] + # loop over the six directions in this order (for consistency with ImageJ) + borders[:] = [4, 3, 2, 1, 5, 6] + + # no need to worry about the z direction if the original image is 2D. + if img.shape[0] == 3: + num_borders = 4 + else: + num_borders = 6 + + # loop through the image several times until there is no change for all + # the six border types + while unchanged_borders < num_borders: + unchanged_borders = 0 + for j in range(num_borders): + curr_border = borders[j] + + simple_border_points = find_simple_point_candidates(img, curr_border) + + # sequential re-checking to preserve connectivity when deleting + # in a parallel way + no_change = True + num_border_points = len(simple_border_points) + for i in range(num_border_points): + point = simple_border_points[i] + p, r, c = point + get_neighborhood(img, p, r, c, neighb) + if is_simple_point(neighb): + img[p, r, c] = 0 + no_change = False + + if no_change: + unchanged_borders += 1 + + return np.asarray(img) + + +@cython.boundscheck(False) +@cython.wraparound(False) +cdef list find_simple_point_candidates(pixel_type[:, :, ::1] img, + int curr_border): + """Inner loop of compute_thin_image. + + The algorithm of [Lee94]_ proceeds in two steps: (1) six directions are + checked for simple border points to remove, and (2) these candidates are + sequentially rechecked, see Sec 3 of [Lee94]_ for rationale and discussion. + + This routine implements the first step above: it loops over the image + for a given direction and assembles candidates for removal. + + """ + # This routine looks like it could be nogil, but actually it cannot be, + # because of `simple_border_points` being a python list which is being + # mutated. + cdef: + list simple_border_points = [] + pixel_type neighborhood[27] + npy_intp p, r, c + bint is_border_pt + (npy_intp, npy_intp, npy_intp) point + + # rebind a global name to avoid lookup. The table is filled in + # at import time. + int[::1] Euler_LUT = LUT + + # loop through the image + # NB: each loop is from 1 to size-1: img is padded from all sides + for p in range(1, img.shape[0] - 1): + for r in range(1, img.shape[1] - 1): + for c in range(1, img.shape[2] - 1): + + # check if pixel is foreground + if img[p, r, c] != 1: + continue + + is_border_pt = (curr_border == 1 and img[p, r, c-1] == 0 or #N + curr_border == 2 and img[p, r, c+1] == 0 or #S + curr_border == 3 and img[p, r+1, c] == 0 or #E + curr_border == 4 and img[p, r-1, c] == 0 or #W + curr_border == 5 and img[p+1, r, c] == 0 or #U + curr_border == 6 and img[p-1, r, c] == 0) #B + if not is_border_pt: + # current point is not deletable + continue + + get_neighborhood(img, p, r, c, neighborhood) + + # check if (p, r, c) can be deleted: + # * it must not be an endpoint; + # * it must be Euler invariant (condition 1 in [Lee94]_); and + # * it must be simple (i.e., its deletion does not change + # connectivity in the 3x3x3 neighborhood) + # this is conditions 2 and 3 in [Lee94]_ + if (is_endpoint(neighborhood) or + not is_Euler_invariant(neighborhood, Euler_LUT) or + not is_simple_point(neighborhood)): + continue + + # ok, add (p, r, c) to the list of simple border points + point = p, r, c + simple_border_points.append(point) + return simple_border_points + + +@cython.boundscheck(False) +@cython.wraparound(False) +cdef void get_neighborhood(pixel_type[:, :, ::1] img, + npy_intp p, npy_intp r, npy_intp c, + pixel_type neighborhood[]): + """Get the neighborhood of a pixel. + + Assume zero boundary conditions. + Image is already padded, so no out-of-bounds checking. + + For the numbering of points see Fig. 1a. of [Lee94]_, where the numbers + do *not* include the center point itself. OTOH, this numbering below + includes it as number 13. The latter is consistent with [IAC15]_. + """ + neighborhood[0] = img[p-1, r-1, c-1] + neighborhood[1] = img[p-1, r, c-1] + neighborhood[2] = img[p-1, r+1, c-1] + + neighborhood[ 3] = img[p-1, r-1, c] + neighborhood[ 4] = img[p-1, r, c] + neighborhood[ 5] = img[p-1, r+1, c] + + neighborhood[ 6] = img[p-1, r-1, c+1] + neighborhood[ 7] = img[p-1, r, c+1] + neighborhood[ 8] = img[p-1, r+1, c+1] + + neighborhood[ 9] = img[p, r-1, c-1] + neighborhood[10] = img[p, r, c-1] + neighborhood[11] = img[p, r+1, c-1] + + neighborhood[12] = img[p, r-1, c] + neighborhood[13] = img[p, r, c] + neighborhood[14] = img[p, r+1, c] + + neighborhood[15] = img[p, r-1, c+1] + neighborhood[16] = img[p, r, c+1] + neighborhood[17] = img[p, r+1, c+1] + + neighborhood[18] = img[p+1, r-1, c-1] + neighborhood[19] = img[p+1, r, c-1] + neighborhood[20] = img[p+1, r+1, c-1] + + neighborhood[21] = img[p+1, r-1, c] + neighborhood[22] = img[p+1, r, c] + neighborhood[23] = img[p+1, r+1, c] + + neighborhood[24] = img[p+1, r-1, c+1] + neighborhood[25] = img[p+1, r, c+1] + neighborhood[26] = img[p+1, r+1, c+1] + + +###### look-up tables +def fill_Euler_LUT(): + """ Look-up table for preserving Euler characteristic. + + This is column $\delta G_{26}$ of Table 2 of [Lee94]_. + """ + cdef int arr[128] + arr[:] = [1, -1, -1, 1, -3, -1, -1, 1, -1, 1, 1, -1, 3, 1, 1, -1, -3, -1, + 3, 1, 1, -1, 3, 1, -1, 1, 1, -1, 3, 1, 1, -1, -3, 3, -1, 1, 1, + 3, -1, 1, -1, 1, 1, -1, 3, 1, 1, -1, 1, 3, 3, 1, 5, 3, 3, 1, + -1, 1, 1, -1, 3, 1, 1, -1, -7, -1, -1, 1, -3, -1, -1, 1, -1, + 1, 1, -1, 3, 1, 1, -1, -3, -1, 3, 1, 1, -1, 3, 1, -1, 1, 1, + -1, 3, 1, 1, -1, -3, 3, -1, 1, 1, 3, -1, 1, -1, 1, 1, -1, 3, + 1, 1, -1, 1, 3, 3, 1, 5, 3, 3, 1, -1, 1, 1, -1, 3, 1, 1, -1] + cdef ndarray LUT = np.zeros(256, dtype=np.intc) + LUT[1::2] = arr + return LUT +cdef int[::1] LUT = fill_Euler_LUT() + + +# Fill the look-up table for indexing octants for computing the Euler +# characteristic. See is_Euler_invariant routine below. +{{py: +_neighb_idx = [[2, 1, 11, 10, 5, 4, 14], # NEB + [0, 9, 3, 12, 1, 10, 4], # NWB + [8, 7, 17, 16, 5, 4, 14], # SEB + [6, 15, 7, 16, 3, 12, 4], # SWB + [20, 23, 19, 22, 11, 14, 10], # NEU + [18, 21, 9, 12, 19, 22, 10], # NWU + [26, 23, 17, 14, 25, 22, 16], # SEU + [24, 25, 15, 16, 21, 22, 12], # SWU + ] +}} + + +@cython.boundscheck(False) +@cython.wraparound(False) +cdef bint is_Euler_invariant(pixel_type neighbors[], + int[::1] lut): + """Check if a point is Euler invariant. + + Calculate Euler characteristic for each octant and sum up. + + Parameters + ---------- + neighbors + neighbors of a point + lut + The look-up table for preserving the Euler characteristic. + + Returns + ------- + bool (C bool, that is) + + """ + cdef int n, euler_char = 0 + {{for _octant in range(8)}} + + # octant {{_octant}}: + n = 1 + {{for _j in range(7):}} + {{py: _idx = _neighb_idx[_octant][_j]}} + if neighbors[{{_idx}}] == 1: + n |= {{1 << (7 - _j)}} + + {{endfor}} + euler_char += lut[n] + {{endfor}} + return euler_char == 0 + + +cdef inline bint is_endpoint(pixel_type neighbors[]): + """An endpoint has exactly one neighbor in the 26-neighborhood. + """ + # The center pixel is counted, thus r.h.s. is 2 + cdef int s = 0, j + for j in range(27): + s += neighbors[j] + return s == 2 + + +cdef bint is_simple_point(pixel_type neighbors[]): + """Check is a point is a Simple Point. + + A point is simple iff its deletion does not change connectivity in + the 3x3x3 neighborhood. (cf conditions 2 and 3 in [Lee94]_). + + This method is named "N(v)_labeling" in [Lee94]_. + + Parameters + ---------- + neighbors : uint8 C array, shape(27,) + neighbors of the point + + Returns + ------- + bool + Whether the point is simple or not. + + """ + # copy neighbors for labeling + # ignore center pixel (i=13) when counting (see [Lee94]_) + cdef pixel_type cube[26] + memcpy(cube, neighbors, 13*sizeof(pixel_type)) + memcpy(cube+13, neighbors+14, 13*sizeof(pixel_type)) + + # set initial label + cdef int label = 2, i + + # for all point in the neighborhood + for i in range(26): + if cube[i] == 1: + # voxel has not been labeled yet + # start recursion with any octant that contains the point i + if i in (0, 1, 3, 4, 9, 10, 12): + octree_labeling(1, label, cube) + elif i in (2, 5, 11, 13): + octree_labeling(2, label, cube) + elif i in (6, 7, 14, 15): + octree_labeling(3, label, cube) + elif i in (8, 16): + octree_labeling(4, label, cube) + elif i in (17, 18, 20, 21): + octree_labeling(5, label, cube) + elif i in (19, 22): + octree_labeling(6, label, cube) + elif i in (23, 24): + octree_labeling(7, label, cube) + elif i == 25: + octree_labeling(8, label, cube) + label += 1 + if label - 2 >= 2: + return False + return True + + +# Octree structure for labeling in `octree_labeling` routine below. +# NB: this is only available at build time, and is used by Tempita templating. +{{py: +_octree = [ + # octant 1 + ([0, 1, 3, 4, 9, 10, 12], + [[], [2], [3], [2, 3, 4], [5], [2, 5, 6], [3, 5, 7]]), + # octant 2 + ([1, 4, 10, 2, 5, 11, 13], + [[1], [1, 3, 4], [1, 5, 6], [], [4], [6], [4, 6, 8]]), + # octant 3 + ([3, 4, 12, 6, 7, 14, 15], + [[1], [1, 2, 4], [1, 5, 7], [], [4], [7], [4, 7, 8]]), + # octant 4 + ([4, 5, 13, 7, 15, 8, 16], + [[1, 2, 3], [2], [2, 6, 8], [3], [3, 7, 8], [], [8]]), + # octant 5 + ([9, 10, 12, 17, 18, 20, 21], + [[1], [1, 2, 6], [1, 3, 7], [], [6], [7], [6, 7, 8]]), + # octant 6 + ([10, 11, 13, 18, 21, 19, 22], + [[1, 2, 5], [2], [2, 4, 8], [5], [5, 7, 8], [], [8]]), + # octant 7 + ([12, 14, 15, 20, 21, 23, 24], + [[1, 3, 5], [3], [3, 4, 8], [5], [5, 6, 8], [], [8]]), + # octant 8 + ([13, 15, 16, 21, 22, 24, 25], + [[2, 4, 6], [3, 4, 7], [4], [5, 6, 7], [6], [7], []]) +] +}} + +@cython.boundscheck(False) +@cython.wraparound(False) +cdef void octree_labeling(int octant, int label, pixel_type cube[]): + """This is a recursive method that calculates the number of connected + components in the 3D neighborhood after the center pixel would + have been removed. + + See Figs. 6 and 7 of [Lee94]_ for the values of indices. + + Parameters + ---------- + octant : int + octant index + label : int + the current label of the center point + cube : uint8 C array, shape(26,) + local neighborhood of the point + + """ + # This routine checks if there are points in the octant with value 1 + # Then sets points in this octant to current label + # and recursive labeling of adjacent octants. + # + # Below, leading underscore means build-time variables. + {{for _oct in range(1, 9)}} + + if octant == {{_oct}}: + {{py: _indices, _list_octants = _octree[_oct-1]}} + {{for _idx, _new_octants in zip(_indices, _list_octants)}} + if cube[{{_idx}}] == 1: + cube[{{_idx}}] = label + {{for _new_octant in _new_octants}} + octree_labeling({{_new_octant}}, label, cube) + {{endfor}} + {{endfor}} + {{endfor}} diff --git a/skimage/morphology/setup.py b/skimage/morphology/setup.py index dbbcad8b..a1945d91 100644 --- a/skimage/morphology/setup.py +++ b/skimage/morphology/setup.py @@ -16,6 +16,7 @@ def configuration(parent_package='', top_path=None): cython(['_skeletonize_cy.pyx'], working_path=base_path) cython(['_convex_hull.pyx'], working_path=base_path) cython(['_greyreconstruct.pyx'], working_path=base_path) + cython(['_skeletonize_3d_cy.pyx.in'], working_path=base_path) config.add_extension('_watershed', sources=['_watershed.c'], include_dirs=[get_numpy_include_dirs()]) @@ -25,6 +26,8 @@ def configuration(parent_package='', top_path=None): include_dirs=[get_numpy_include_dirs()]) config.add_extension('_greyreconstruct', sources=['_greyreconstruct.c'], include_dirs=[get_numpy_include_dirs()]) + config.add_extension('_skeletonize_3d_cy', sources=['_skeletonize_3d_cy.c'], + include_dirs=[get_numpy_include_dirs()]) return config diff --git a/skimage/morphology/tests/test_skeletonize_3d.py b/skimage/morphology/tests/test_skeletonize_3d.py new file mode 100644 index 00000000..d5859a5d --- /dev/null +++ b/skimage/morphology/tests/test_skeletonize_3d.py @@ -0,0 +1,186 @@ +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, assert_raises, + assert_) + +import scipy.ndimage as ndi + +import skimage +from skimage import io, draw, data_dir +from skimage.data import binary_blobs +from skimage.util import img_as_ubyte + +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_1D(): + # a corner case of an image of a shape(1, N) + im = np.ones((5, 1), dtype=np.uint8) + res = skeletonize_3d(im) + assert_equal(res, 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_input(): + # check that the input is not clobbered + # for 2D and 3D images of varying dtypes + imgs = [np.ones((8, 8), dtype=float), np.ones((4, 8, 8), dtype=float), + np.ones((8, 8), dtype=np.uint8), np.ones((4, 8, 8), dtype=np.uint8), + np.ones((8, 8), dtype=bool), np.ones((4, 8, 8), dtype=bool)] + for img in imgs: + yield check_input, img + + +def check_input(img): + orig = img.copy() + with warnings.catch_warnings(): + # UserWarning for possible precision loss, expected + warnings.simplefilter('ignore', UserWarning) + res = skeletonize_3d(img) + assert_equal(img, orig) + + +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)) + + +def test_two_hole_image(): + # test a simple 2D image against FIJI + img_o = np.array([[0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0], + [0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0], + [0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0], + [0, 1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0], + [0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0], + [0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0], + [0, 0, 0, 1, 1, 0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0], + [0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0], + [0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0], + [0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0], + [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], + dtype=np.uint8) + img_f = np.array([[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0], + [0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0], + [0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0], + [0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0], + [0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0], + [0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0], + [0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0], + [0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], + dtype=np.uint8) + res = skeletonize_3d(img_o) + assert_equal(res, img_f) + + +def test_3d_vs_fiji(): + # generate an image with blobs and compate its skeleton to + # the skeleton generated by FIJI + img = binary_blobs(32, 0.05, n_dim=3, seed=1234) + img = img[:-2, ...] + img = img.astype(np.uint8)*255 + + img_s = skeletonize_3d(img) + img_f = io.imread(os.path.join(data_dir, "_blobs_3d_fiji_skeleton.tif")) + assert_equal(img_s, img_f) + + +if __name__ == '__main__': + run_module_suite() diff --git a/tools/travis_before_install.sh b/tools/travis_before_install.sh index b4baeaf4..98eef785 100755 --- a/tools/travis_before_install.sh +++ b/tools/travis_before_install.sh @@ -30,7 +30,7 @@ retry () { } # add build dependencies -echo "cython>=0.21" >> requirements.txt +echo "cython>=0.23.4" >> requirements.txt # require networkx 1.9.1 on 2.6, as 2.6 support was dropped in 1.10 # require matplotlib 1.4.3 on 2.6, as 2.6 support was dropped in 1.5