diff --git a/skimage/morphology/_skeletonize.py b/skimage/morphology/_skeletonize.py index e892c185..0ea9cf7d 100644 --- a/skimage/morphology/_skeletonize.py +++ b/skimage/morphology/_skeletonize.py @@ -5,11 +5,10 @@ Algorithms for computing the skeleton of a binary image import numpy as np from scipy import ndimage as ndi -from ._skeletonize_cy import _skeletonize_loop, _table_lookup_index +from ._skeletonize_cy import _fast_skeletonize, _skeletonize_loop, _table_lookup_index # --------- Skeletonization by morphological thinning --------- - def skeletonize(image): """Return the skeleton of a binary image. @@ -34,7 +33,7 @@ def skeletonize(image): Notes ----- - The algorithm [1] works by making successive passes of the image, + The algorithm [1]_ works by making successive passes of the image, removing pixels on object borders. This continues until no more pixels can be removed. The image is correlated with a mask that assigns each pixel a number in the range [0...255] @@ -50,8 +49,8 @@ def skeletonize(image): References ---------- .. [1] A fast parallel algorithm for thinning digital patterns, - T. Y. ZHANG and C. Y. SUEN, Communications of the ACM, - March 1984, Volume 27, Number 3 + T. Y. Zhang and C. Y. Suen, Communications of the ACM, + March 1984, Volume 27, Number 3. Examples @@ -81,76 +80,20 @@ def skeletonize(image): [0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=uint8) """ - # look up table - there is one entry for each of the 2^8=256 possible - # combinations of 8 binary neighbours. 1's, 2's and 3's are candidates - # for removal at each iteration of the algorithm. - lut = [0, 0, 0, 1, 0, 0, 1, 3, 0, 0, 3, 1, 1, 0, 1, 3, 0, 0, 0, 0, 0, 0, - 0, 0, 2, 0, 2, 0, 3, 0, 3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0, - 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 3, 0, 2, 2, 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, 2, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 2, 0, - 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0, 3, 0, 2, 0, 0, 0, 3, 1, - 0, 0, 1, 3, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, - 0, 0, 0, 0, 0, 1, 3, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, - 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 3, 1, 3, 0, 0, - 1, 3, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, - 0, 0, 0, 0, 2, 3, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 3, 3, - 0, 1, 0, 0, 0, 0, 2, 2, 0, 0, 2, 0, 0, 0] # convert to unsigned int (this should work for boolean values) - skeleton = image.astype(np.uint8) + image = image.astype(np.uint8) # check some properties of the input image: # - 2D # - binary image with only 0's and 1's - if skeleton.ndim != 2: + if image.ndim != 2: raise ValueError('Skeletonize requires a 2D array') - if not np.all(np.in1d(skeleton.flat, (0, 1))): + if not np.all(np.in1d(image.flat, (0, 1))): raise ValueError('Image contains values other than 0 and 1') - # create the mask that will assign a unique value based on the - # arrangement of neighbouring pixels - mask = np.array([[ 1, 2, 4], - [128, 0, 8], - [ 64, 32, 16]], np.uint8) + return _fast_skeletonize(image) - pixel_removed = True - while pixel_removed: - pixel_removed = False - - # assign each pixel a unique value based on its foreground neighbours - neighbours = ndi.correlate(skeleton, mask, mode='constant') - - # ignore background - neighbours *= skeleton - - # use LUT to categorize each foreground pixel as a 0, 1, 2 or 3 - codes = np.take(lut, neighbours) - - # pass 1 - remove the 1's and 3's - code_mask = (codes == 1) - if np.any(code_mask): - pixel_removed = True - skeleton[code_mask] = 0 - code_mask = (codes == 3) - if np.any(code_mask): - pixel_removed = True - skeleton[code_mask] = 0 - - # pass 2 - remove the 2's and 3's - neighbours = ndi.correlate(skeleton, mask, mode='constant') - neighbours *= skeleton - codes = np.take(lut, neighbours) - code_mask = (codes == 2) - if np.any(code_mask): - pixel_removed = True - skeleton[code_mask] = 0 - code_mask = (codes == 3) - if np.any(code_mask): - pixel_removed = True - skeleton[code_mask] = 0 - - return skeleton.astype(bool) # --------- Skeletonization by medial axis transform -------- diff --git a/skimage/morphology/_skeletonize_cy.pyx b/skimage/morphology/_skeletonize_cy.pyx index 8dc8eb03..54e90266 100644 --- a/skimage/morphology/_skeletonize_cy.pyx +++ b/skimage/morphology/_skeletonize_cy.pyx @@ -3,7 +3,115 @@ #cython: nonecheck=False #cython: wraparound=False -''' +import numpy as np +cimport numpy as cnp + +def _fast_skeletonize(image): + """Optimized parts of the Zhang-Suen [1]_ skeletonization. + Iteratively, pixels meeting removal criteria are removed, + till only the skeleton remains (that is, no further removable pixel + was found). + + Performs a hard-coded correlation to assign every neighborhood of 8 a + unique number, which in turn is used in conjunction with a look up + table to select the appropriate thinning criteria. + + Parameters + ---------- + image : numpy.ndarray + A binary image containing the objects to be skeletonized. '1' + represents foreground, and '0' represents background. + + Returns + ------- + skeleton : ndarray + A matrix containing the thinned image. + + References + ---------- + .. [1] A fast parallel algorithm for thinning digital patterns, + T. Y. Zhang and C. Y. Suen, Communications of the ACM, + March 1984, Volume 27, Number 3. + + """ + + # look up table - there is one entry for each of the 2^8=256 possible + # combinations of 8 binary neighbours. 1's, 2's and 3's are candidates + # for removal at each iteration of the algorithm. + cdef int *lut = \ + [0, 0, 0, 1, 0, 0, 1, 3, 0, 0, 3, 1, 1, 0, 1, 3, 0, 0, 0, 0, 0, 0, + 0, 0, 2, 0, 2, 0, 3, 0, 3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 3, 0, 2, 2, 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, 2, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 2, 0, + 0, 0, 3, 0, 0, 0, 0, 0, 0, 0, 3, 0, 0, 0, 3, 0, 2, 0, 0, 0, 3, 1, + 0, 0, 1, 3, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 1, 3, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 3, 1, 3, 0, 0, + 1, 3, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 2, 3, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 3, 3, + 0, 1, 0, 0, 0, 0, 2, 2, 0, 0, 2, 0, 0, 0] + + cdef int pixel_removed, first_pass, neighbors + + # indices for fast iteration + cdef Py_ssize_t row, col, nrows = image.shape[0]+2, ncols = image.shape[1]+2 + + # we copy over the image into a larger version with a single pixel border + # this removes the need to handle border cases below + _skeleton = np.zeros((nrows, ncols), dtype=np.uint8) + _skeleton[1:nrows-1, 1:ncols-1] = image > 0 + + _cleaned_skeleton = _skeleton.copy() + + # cdef'd numpy-arrays for fast, typed access + cdef cnp.ndarray[cnp.uint8_t, ndim=2] skeleton, cleaned_skeleton + + skeleton = _skeleton + cleaned_skeleton = _cleaned_skeleton + + pixel_removed = True + + # the algorithm reiterates the thinning till + # no further thinning occurred (variable pixel_removed set) + + while pixel_removed: + pixel_removed = False + + # there are two phases, in the first phase, pixels labeled (see below) + # 1 and 3 are removed, in the second 2 and 3 + + for first_pass in (True, False): + for row in range(1, nrows-1): + for col in range(1, ncols-1): + # all set pixels ... + if skeleton[row, col]: + # are correlated with a kernel (coefficients spread around here ...) + # to apply a unique number to every possible neighborhood ... + + # which is used with the lut to find the "connectivity type" + + neighbors = lut[ 1*skeleton[row - 1, col - 1] + 2*skeleton[row - 1, col] +\ + 4*skeleton[row - 1, col + 1] + 8*skeleton[row, col + 1] +\ + 16*skeleton[row + 1, col + 1] + 32*skeleton[row + 1, col] +\ + 64*skeleton[row + 1, col - 1] + 128*skeleton[row, col - 1]] + + # if the condition is met, the pixel is removed (unset) + if ((neighbors == 1 and first_pass) or + (neighbors == 2 and not first_pass) or + (neighbors == 3)): + cleaned_skeleton[row, col] = 0 + pixel_removed = True + + # once a step has been processed, the original skeleton + # is overwritten with the cleaned version + _skeleton = _cleaned_skeleton.copy() + skeleton = _skeleton + + return _skeleton[1:nrows-1, 1:ncols-1].astype(np.bool) + + +""" Originally part of CellProfiler, code licensed under both GPL and BSD licenses. Website: http://www.cellprofiler.org @@ -12,11 +120,7 @@ Copyright (c) 2009-2011 Broad Institute All rights reserved. Original author: Lee Kamentsky -''' - -import numpy as np -cimport numpy as cnp - +""" def _skeletonize_loop(cnp.uint8_t[:, ::1] result, Py_ssize_t[::1] i, Py_ssize_t[::1] j,