mirror of
https://github.com/wassname/scikit-image.git
synced 2026-07-23 13:10:18 +08:00
Merge pull request #1572 from csachs/cython_skeletonize
fast Cython based skeletonize implementation
This commit is contained in:
Juan Nunez-Iglesias
@@ -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 --------
|
||||
|
||||
|
||||
@@ -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,
|
||||
|
||||
Reference in New Issue
Block a user