""" Algorithms for computing the skeleton of a binary image """ import numpy as np from scipy import ndimage as ndi from .._shared.utils import check_nD from ..util import crop from ._skeletonize_lee_cy import _compute_thin_image from ._skeletonize_various_cy import ( _fast_skeletonize, _skeletonize_loop, _table_lookup_index, ) def skeletonize(image, *, method=None): """Compute the skeleton of the input image via thinning. Parameters ---------- image : (M, N[, P]) ndarray of bool or int The image containing the objects to be skeletonized. Each connected component in the image is reduced to a single-pixel wide skeleton. The image is binarized prior to thinning; thus, adjacent objects of different intensities are considered as one. Zero or ``False`` values represent the background, nonzero or ``True`` values -- foreground. method : {'zhang', 'lee'}, optional Which algorithm to use. Zhang's algorithm [Zha84]_ only works for 2D images, and is the default for 2D. Lee's algorithm [Lee94]_ works for 2D or 3D images and is the default for 3D. Returns ------- skeleton : (M, N[, P]) ndarray of bool The thinned image. See Also -------- medial_axis 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. .. [Zha84] 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. Examples -------- >>> X, Y = np.ogrid[0:9, 0:9] >>> ellipse = (1./3 * (X - 4)**2 + (Y - 4)**2 < 3**2).astype(bool) >>> ellipse.view(np.uint8) array([[0, 0, 0, 1, 1, 1, 0, 0, 0], [0, 0, 1, 1, 1, 1, 1, 0, 0], [0, 0, 1, 1, 1, 1, 1, 0, 0], [0, 0, 1, 1, 1, 1, 1, 0, 0], [0, 0, 1, 1, 1, 1, 1, 0, 0], [0, 0, 1, 1, 1, 1, 1, 0, 0], [0, 0, 1, 1, 1, 1, 1, 0, 0], [0, 0, 1, 1, 1, 1, 1, 0, 0], [0, 0, 0, 1, 1, 1, 0, 0, 0]], dtype=uint8) >>> skel = skeletonize(ellipse) >>> skel.view(np.uint8) 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, 1, 0, 0, 0, 0], [0, 0, 0, 0, 1, 0, 0, 0, 0], [0, 0, 0, 0, 1, 0, 0, 0, 0], [0, 0, 0, 0, 1, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=uint8) """ image = image.astype(bool, order="C", copy=False) if method not in {'zhang', 'lee', None}: raise ValueError( f'skeletonize method should be either "lee" or "zhang", ' f'got {method}.' ) if image.ndim == 2 and (method is None or method == 'zhang'): skeleton = _skeletonize_zhang(image) elif image.ndim == 3 and method == 'zhang': raise ValueError('skeletonize method "zhang" only works for 2D ' 'images.') elif image.ndim == 3 or (image.ndim == 2 and method == 'lee'): skeleton = _skeletonize_lee(image) else: raise ValueError( f'skeletonize requires a 2D or 3D image as input, ' f'got {image.ndim}D.' ) return skeleton def _skeletonize_zhang(image): """Return the skeleton of a 2D binary image. Thinning is used to reduce each connected component in a binary image to a single-pixel wide skeleton. Parameters ---------- image : numpy.ndarray An image containing the objects to be skeletonized. Zeros or ``False`` represent background, nonzero values or ``True`` are foreground. Returns ------- skeleton : ndarray A matrix containing the thinned image. See Also -------- medial_axis, skeletonize, thin Notes ----- The algorithm [Zha84]_ 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] corresponding to each possible pattern of its 8 neighboring pixels. A look up table is then used to assign the pixels a value of 0, 1, 2 or 3, which are selectively removed during the iterations. Note that this algorithm will give different results than a medial axis transform, which is also often referred to as "skeletonization". References ---------- .. [Zha84] 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. Examples -------- >>> X, Y = np.ogrid[0:9, 0:9] >>> ellipse = (1./3 * (X - 4)**2 + (Y - 4)**2 < 3**2).astype(bool) >>> ellipse.view(np.uint8) array([[0, 0, 0, 1, 1, 1, 0, 0, 0], [0, 0, 1, 1, 1, 1, 1, 0, 0], [0, 0, 1, 1, 1, 1, 1, 0, 0], [0, 0, 1, 1, 1, 1, 1, 0, 0], [0, 0, 1, 1, 1, 1, 1, 0, 0], [0, 0, 1, 1, 1, 1, 1, 0, 0], [0, 0, 1, 1, 1, 1, 1, 0, 0], [0, 0, 1, 1, 1, 1, 1, 0, 0], [0, 0, 0, 1, 1, 1, 0, 0, 0]], dtype=uint8) >>> skel = skeletonize(ellipse) >>> skel.view(np.uint8) 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, 1, 0, 0, 0, 0], [0, 0, 0, 0, 1, 0, 0, 0, 0], [0, 0, 0, 0, 1, 0, 0, 0, 0], [0, 0, 0, 0, 1, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=uint8) """ if image.ndim != 2: raise ValueError("Zhang's skeletonize method requires a 2D array") return _fast_skeletonize(image) # --------- Skeletonization and thinning based on Guo and Hall 1989 --------- def _generate_thin_luts(): """generate LUTs for thinning algorithm (for reference)""" def nabe(n): return np.array([n >> i & 1 for i in range(0, 9)]).astype(bool) def G1(n): s = 0 bits = nabe(n) for i in (0, 2, 4, 6): if not (bits[i]) and (bits[i + 1] or bits[(i + 2) % 8]): s += 1 return s == 1 g1_lut = np.array([G1(n) for n in range(256)]) def G2(n): n1, n2 = 0, 0 bits = nabe(n) for k in (1, 3, 5, 7): if bits[k] or bits[k - 1]: n1 += 1 if bits[k] or bits[(k + 1) % 8]: n2 += 1 return min(n1, n2) in [2, 3] g2_lut = np.array([G2(n) for n in range(256)]) g12_lut = g1_lut & g2_lut def G3(n): bits = nabe(n) return not ((bits[1] or bits[2] or not (bits[7])) and bits[0]) def G3p(n): bits = nabe(n) return not ((bits[5] or bits[6] or not (bits[3])) and bits[4]) g3_lut = np.array([G3(n) for n in range(256)]) g3p_lut = np.array([G3p(n) for n in range(256)]) g123_lut = g12_lut & g3_lut g123p_lut = g12_lut & g3p_lut return g123_lut, g123p_lut # fmt: off G123_LUT = np.array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 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, 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, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0], dtype=bool) G123P_LUT = np.array([0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 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, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 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, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 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, 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 1, 0, 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], dtype=bool) # fmt: on def thin(image, max_num_iter=None): """ Perform morphological thinning of a binary image. Parameters ---------- image : binary (M, N) ndarray The image to thin. If this input isn't already a binary image, it gets converted into one: In this case, zero values are considered background (False), nonzero values are considered foreground (True). max_num_iter : int, number of iterations, optional Regardless of the value of this parameter, the thinned image is returned immediately if an iteration produces no change. If this parameter is specified it thus sets an upper bound on the number of iterations performed. Returns ------- out : ndarray of bool Thinned image. See Also -------- skeletonize, medial_axis Notes ----- This algorithm [1]_ works by making multiple passes over the image, removing pixels matching a set of criteria designed to thin connected regions while preserving eight-connected components and 2 x 2 squares [2]_. In each of the two sub-iterations the algorithm correlates the intermediate skeleton image with a neighborhood mask, then looks up each neighborhood in a lookup table indicating whether the central pixel should be deleted in that sub-iteration. References ---------- .. [1] Z. Guo and R. W. Hall, "Parallel thinning with two-subiteration algorithms," Comm. ACM, vol. 32, no. 3, pp. 359-373, 1989. :DOI:`10.1145/62065.62074` .. [2] Lam, L., Seong-Whan Lee, and Ching Y. Suen, "Thinning Methodologies-A Comprehensive Survey," IEEE Transactions on Pattern Analysis and Machine Intelligence, Vol 14, No. 9, p. 879, 1992. :DOI:`10.1109/34.161346` Examples -------- >>> square = np.zeros((7, 7), dtype=bool) >>> square[1:-1, 2:-2] = 1 >>> square[0, 1] = 1 >>> square.view(np.uint8) array([[0, 1, 0, 0, 0, 0, 0], [0, 0, 1, 1, 1, 0, 0], [0, 0, 1, 1, 1, 0, 0], [0, 0, 1, 1, 1, 0, 0], [0, 0, 1, 1, 1, 0, 0], [0, 0, 1, 1, 1, 0, 0], [0, 0, 0, 0, 0, 0, 0]], dtype=uint8) >>> skel = thin(square) >>> skel.view(np.uint8) array([[0, 1, 0, 0, 0, 0, 0], [0, 0, 1, 0, 0, 0, 0], [0, 0, 0, 1, 0, 0, 0], [0, 0, 0, 1, 0, 0, 0], [0, 0, 0, 1, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0, 0]], dtype=uint8) """ # check that image is 2d check_nD(image, 2) # convert image to uint8 with values in {0, 1} skel = np.asanyarray(image, dtype=bool).copy().view(np.uint8) # neighborhood mask mask = np.array([[8, 4, 2], [16, 0, 1], [32, 64, 128]], dtype=np.uint8) # iterate until convergence, up to the iteration limit max_num_iter = max_num_iter or np.inf num_iter = 0 n_pts_old, n_pts_new = np.inf, np.sum(skel) while n_pts_old != n_pts_new and num_iter < max_num_iter: n_pts_old = n_pts_new # perform the two "subiterations" described in the paper for lut in [G123_LUT, G123P_LUT]: # correlate image with neighborhood mask N = ndi.correlate(skel, mask, mode='constant') # take deletion decision from this subiteration's LUT D = np.take(lut, N) # perform deletion skel[D] = 0 n_pts_new = np.sum(skel) # count points after thinning num_iter += 1 return skel.astype(bool) # --------- Skeletonization by medial axis transform -------- _eight_connect = ndi.generate_binary_structure(2, 2) def medial_axis(image, mask=None, return_distance=False, *, rng=None): """Compute the medial axis transform of a binary image. Parameters ---------- image : binary ndarray, shape (M, N) The image of the shape to skeletonize. If this input isn't already a binary image, it gets converted into one: In this case, zero values are considered background (False), nonzero values are considered foreground (True). mask : binary ndarray, shape (M, N), optional If a mask is given, only those elements in `image` with a true value in `mask` are used for computing the medial axis. return_distance : bool, optional If true, the distance transform is returned as well as the skeleton. rng : {`numpy.random.Generator`, int}, optional Pseudo-random number generator. By default, a PCG64 generator is used (see :func:`numpy.random.default_rng`). If `rng` is an int, it is used to seed the generator. The PRNG determines the order in which pixels are processed for tiebreaking. .. versionadded:: 0.19 Returns ------- out : ndarray of bools Medial axis transform of the image dist : ndarray of ints, optional Distance transform of the image (only returned if `return_distance` is True) See Also -------- skeletonize, thin Notes ----- This algorithm computes the medial axis transform of an image as the ridges of its distance transform. The different steps of the algorithm are as follows * A lookup table is used, that assigns 0 or 1 to each configuration of the 3x3 binary square, whether the central pixel should be removed or kept. We want a point to be removed if it has more than one neighbor and if removing it does not change the number of connected components. * The distance transform to the background is computed, as well as the cornerness of the pixel. * The foreground (value of 1) points are ordered by the distance transform, then the cornerness. * A cython function is called to reduce the image to its skeleton. It processes pixels in the order determined at the previous step, and removes or maintains a pixel according to the lookup table. Because of the ordering, it is possible to process all pixels in only one pass. Examples -------- >>> square = np.zeros((7, 7), dtype=bool) >>> square[1:-1, 2:-2] = 1 >>> square.view(np.uint8) array([[0, 0, 0, 0, 0, 0, 0], [0, 0, 1, 1, 1, 0, 0], [0, 0, 1, 1, 1, 0, 0], [0, 0, 1, 1, 1, 0, 0], [0, 0, 1, 1, 1, 0, 0], [0, 0, 1, 1, 1, 0, 0], [0, 0, 0, 0, 0, 0, 0]], dtype=uint8) >>> medial_axis(square).view(np.uint8) array([[0, 0, 0, 0, 0, 0, 0], [0, 0, 1, 0, 1, 0, 0], [0, 0, 0, 1, 0, 0, 0], [0, 0, 0, 1, 0, 0, 0], [0, 0, 0, 1, 0, 0, 0], [0, 0, 1, 0, 1, 0, 0], [0, 0, 0, 0, 0, 0, 0]], dtype=uint8) """ global _eight_connect if mask is None: masked_image = image.astype(bool) else: masked_image = image.astype(bool).copy() masked_image[~mask] = False # # Build lookup table - three conditions # 1. Keep only positive pixels (center_is_foreground array). # AND # 2. Keep if removing the pixel results in a different connectivity # (if the number of connected components is different with and # without the central pixel) # OR # 3. Keep if # pixels in neighborhood is 2 or less # Note that table is independent of image center_is_foreground = (np.arange(512) & 2**4).astype(bool) table = ( center_is_foreground # condition 1. & ( np.array( [ ndi.label(_pattern_of(index), _eight_connect)[1] != ndi.label(_pattern_of(index & ~(2**4)), _eight_connect)[1] for index in range(512) ] ) # condition 2 | np.array([np.sum(_pattern_of(index)) < 3 for index in range(512)]) ) # condition 3 ) # Build distance transform distance = ndi.distance_transform_edt(masked_image) if return_distance: store_distance = distance.copy() # Corners # The processing order along the edge is critical to the shape of the # resulting skeleton: if you process a corner first, that corner will # be eroded and the skeleton will miss the arm from that corner. Pixels # with fewer neighbors are more "cornery" and should be processed last. # We use a cornerness_table lookup table where the score of a # configuration is the number of background (0-value) pixels in the # 3x3 neighborhood cornerness_table = np.array( [9 - np.sum(_pattern_of(index)) for index in range(512)] ) corner_score = _table_lookup(masked_image, cornerness_table) # Define arrays for inner loop i, j = np.mgrid[0 : image.shape[0], 0 : image.shape[1]] result = masked_image.copy() distance = distance[result] i = np.ascontiguousarray(i[result], dtype=np.intp) j = np.ascontiguousarray(j[result], dtype=np.intp) result = np.ascontiguousarray(result, np.uint8) # Determine the order in which pixels are processed. # We use a random # for tiebreaking. Assign each pixel in the image a # predictable, random # so that masking doesn't affect arbitrary choices # of skeletons # generator = np.random.default_rng(rng) tiebreaker = generator.permutation(np.arange(masked_image.sum())) order = np.lexsort((tiebreaker, corner_score[masked_image], distance)) order = np.ascontiguousarray(order, dtype=np.int32) table = np.ascontiguousarray(table, dtype=np.uint8) # Remove pixels not belonging to the medial axis _skeletonize_loop(result, i, j, order, table) result = result.astype(bool) if mask is not None: result[~mask] = image[~mask] if return_distance: return result, store_distance else: return result def _pattern_of(index): """ Return the pattern represented by an index value Byte decomposition of index """ return np.array( [ [index & 2**0, index & 2**1, index & 2**2], [index & 2**3, index & 2**4, index & 2**5], [index & 2**6, index & 2**7, index & 2**8], ], bool, ) def _table_lookup(image, table): """ Perform a morphological transform on an image, directed by its neighbors Parameters ---------- image : ndarray A binary image table : ndarray A 512-element table giving the transform of each pixel given the values of that pixel and its 8-connected neighbors. Returns ------- result : ndarray of same shape as `image` Transformed image Notes ----- The pixels are numbered like this:: 0 1 2 3 4 5 6 7 8 The index at a pixel is the sum of 2** for pixels that evaluate to true. """ # # We accumulate into the indexer to get the index into the table # at each point in the image # if image.shape[0] < 3 or image.shape[1] < 3: image = image.astype(bool) indexer = np.zeros(image.shape, int) indexer[1:, 1:] += image[:-1, :-1] * 2**0 indexer[1:, :] += image[:-1, :] * 2**1 indexer[1:, :-1] += image[:-1, 1:] * 2**2 indexer[:, 1:] += image[:, :-1] * 2**3 indexer[:, :] += image[:, :] * 2**4 indexer[:, :-1] += image[:, 1:] * 2**5 indexer[:-1, 1:] += image[1:, :-1] * 2**6 indexer[:-1, :] += image[1:, :] * 2**7 indexer[:-1, :-1] += image[1:, 1:] * 2**8 else: indexer = _table_lookup_index(np.ascontiguousarray(image, np.uint8)) image = table[indexer] return image def _skeletonize_lee(image): """Compute the skeleton of a binary image. Thinning is used to reduce each connected component in a binary image to a single-pixel wide skeleton. Parameters ---------- image : ndarray, 2D or 3D An image containing the objects to be skeletonized. Zeros or ``False`` represent background, nonzero values or ``True`` are foreground. Returns ------- skeleton : ndarray of bool The thinned image. See Also -------- skeletonize, medial_axis Notes ----- The method of [Lee94]_ uses an octree data structure to examine a 3x3x3 neighborhood of a pixel. The algorithm proceeds by iteratively sweeping over the image, and removing pixels at each iteration until the image stops changing. Each iteration consists of two steps: first, a list of candidates for removal is assembled; then pixels from this list are rechecked sequentially, to better preserve connectivity of the image. The algorithm this function implements is different from the algorithms used by either `skeletonize` or `medial_axis`, thus for 2D images the results produced by this function are generally different. 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. """ # make sure the image is 3D or 2D if image.ndim < 2 or image.ndim > 3: raise ValueError( "skeletonize can only handle 2D or 3D images; " f"got image.ndim = {image.ndim} instead." ) image_o = image.astype(bool, order="C", copy=False) # make a 2D input image 3D and pad it w/ zeros to simplify dealing w/ boundaries # NB: careful here to not clobber the original *and* minimize copying if image.ndim == 2: image_o = image_o[np.newaxis, ...] image_o = np.pad(image_o, pad_width=1, mode='constant') # copies # do the computation image_o = _compute_thin_image(image_o) # crop it back and restore the original intensity range image_o = crop(image_o, crop_width=1) if image.ndim == 2: image_o = image_o[0] return image_o