#cython: cdivision=True #cython: boundscheck=False #cython: nonecheck=False #cython: wraparound=False import numpy as np cimport numpy as cnp cnp.import_array() 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.uint8_t [:, ::1] 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) with nogil: 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 # nogil can't iterate through `(True, False)` because it is a Python # tuple. Use the fact that 0 is Falsy, and 1 is truthy in C # for the iteration instead. # for first_pass in (True, False): for pass_num in range(2): first_pass = (pass_num == 0) 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[:, :] return _skeleton[1:nrows-1, 1:ncols-1].astype(bool) """ Originally part of CellProfiler, code licensed under both GPL and BSD licenses. Website: http://www.cellprofiler.org Copyright (c) 2003-2009 Massachusetts Institute of Technology Copyright (c) 2009-2011 Broad Institute All rights reserved. Original author: Lee Kamentsky """ def _skeletonize_loop(cnp.uint8_t[:, ::1] result, Py_ssize_t[::1] i, Py_ssize_t[::1] j, cnp.int32_t[::1] order, cnp.uint8_t[::1] table): """ Inner loop of skeletonize function Parameters ---------- result : ndarray of uint8 On input, the image to be skeletonized, on output the skeletonized image. i, j : ndarrays The coordinates of each foreground pixel in the image order : ndarray The index of each pixel, in the order of processing (order[0] is the first pixel to process, etc.) table : ndarray The 512-element lookup table of values after transformation (whether to keep or not each configuration in a binary 3x3 array) Notes ----- The loop determines whether each pixel in the image can be removed without changing the Euler number of the image. The pixels are ordered by increasing distance from the background which means a point nearer to the quench-line of the brushfire will be evaluated later than a point closer to the edge. Note that the neighbourhood of a pixel may evolve before the loop arrives at this pixel. This is why it is possible to compute the skeleton in only one pass, thanks to an adapted ordering of the pixels. """ cdef: Py_ssize_t accumulator Py_ssize_t index, order_index Py_ssize_t ii, jj Py_ssize_t rows = result.shape[0] Py_ssize_t cols = result.shape[1] with nogil: for index in range(order.shape[0]): accumulator = 16 order_index = order[index] ii = i[order_index] jj = j[order_index] # Compute the configuration around the pixel if ii > 0: if jj > 0 and result[ii - 1, jj - 1]: accumulator += 1 if result[ii - 1, jj]: accumulator += 2 if jj < cols - 1 and result[ii - 1, jj + 1]: accumulator += 4 if jj > 0 and result[ii, jj - 1]: accumulator += 8 if jj < cols - 1 and result[ii, jj + 1]: accumulator += 32 if ii < rows - 1: if jj > 0 and result[ii + 1, jj - 1]: accumulator += 64 if result[ii + 1, jj]: accumulator += 128 if jj < cols - 1 and result[ii + 1, jj + 1]: accumulator += 256 # Assign the value of table corresponding to the configuration result[ii, jj] = table[accumulator] def _table_lookup_index(cnp.uint8_t[:, ::1] image): """ Return an index into a table per pixel of a binary image Take the sum of true neighborhood pixel values where the neighborhood looks like this:: 1 2 4 8 16 32 64 128 256 This code could be replaced by a convolution with the kernel:: 256 128 64 32 16 8 4 2 1 but this runs about twice as fast because of inlining and the hardwired kernel. """ cdef: Py_ssize_t[:, ::1] indexer Py_ssize_t *p_indexer cnp.uint8_t *p_image Py_ssize_t i_stride Py_ssize_t i_shape Py_ssize_t j_shape Py_ssize_t i Py_ssize_t j Py_ssize_t offset i_shape = image.shape[0] j_shape = image.shape[1] indexer = np.zeros((i_shape, j_shape), dtype=np.intp) p_indexer = &indexer[0, 0] p_image = &image[0, 0] i_stride = image.strides[0] assert i_shape >= 3 and j_shape >= 3, \ "Please use the slow method for arrays < 3x3" with nogil: for i in range(1, i_shape-1): offset = i_stride* i + 1 for j in range(1, j_shape - 1): if p_image[offset]: p_indexer[offset + i_stride + 1] += 1 p_indexer[offset + i_stride] += 2 p_indexer[offset + i_stride - 1] += 4 p_indexer[offset + 1] += 8 p_indexer[offset] += 16 p_indexer[offset - 1] += 32 p_indexer[offset - i_stride + 1] += 64 p_indexer[offset - i_stride] += 128 p_indexer[offset - i_stride - 1] += 256 offset += 1 # # Do the corner cases (literally) # if image[0, 0]: indexer[0, 0] += 16 indexer[0, 1] += 8 indexer[1, 0] += 2 indexer[1, 1] += 1 if image[0, j_shape - 1]: indexer[0, j_shape - 2] += 32 indexer[0, j_shape - 1] += 16 indexer[1, j_shape - 2] += 4 indexer[1, j_shape - 1] += 2 if image[i_shape - 1, 0]: indexer[i_shape - 2, 0] += 128 indexer[i_shape - 2, 1] += 64 indexer[i_shape - 1, 0] += 16 indexer[i_shape - 1, 1] += 8 if image[i_shape - 1, j_shape - 1]: indexer[i_shape - 2, j_shape - 2] += 256 indexer[i_shape - 2, j_shape - 1] += 128 indexer[i_shape - 1, j_shape - 2] += 32 indexer[i_shape - 1, j_shape - 1] += 16 # # Do the edges # for j in range(1, j_shape - 1): if image[0, j]: indexer[0, j - 1] += 32 indexer[0, j] += 16 indexer[0, j + 1] += 8 indexer[1, j - 1] += 4 indexer[1, j] += 2 indexer[1, j + 1] += 1 if image[i_shape - 1, j]: indexer[i_shape - 2, j - 1] += 256 indexer[i_shape - 2, j] += 128 indexer[i_shape - 2, j + 1] += 64 indexer[i_shape - 1, j - 1] += 32 indexer[i_shape - 1, j] += 16 indexer[i_shape - 1, j + 1] += 8 for i in range(1, i_shape - 1): if image[i, 0]: indexer[i - 1, 0] += 128 indexer[i, 0] += 16 indexer[i + 1, 0] += 2 indexer[i - 1, 1] += 64 indexer[i, 1] += 8 indexer[i + 1, 1] += 1 if image[i, j_shape - 1]: indexer[i - 1, j_shape - 2] += 256 indexer[i, j_shape - 2] += 32 indexer[i + 1, j_shape - 2] += 4 indexer[i - 1, j_shape - 1] += 128 indexer[i, j_shape - 1] += 16 indexer[i + 1, j_shape - 1] += 2 return np.asarray(indexer)