diff --git a/doc/examples/plot_medial_transform.py b/doc/examples/plot_medial_transform.py new file mode 100644 index 00000000..1421b4ce --- /dev/null +++ b/doc/examples/plot_medial_transform.py @@ -0,0 +1,70 @@ +""" +=========================== +Medial axis skeletonization +=========================== + +The medial axis of an object is the set of all points having more than one +closest point on the object's boundary. It is often called the **topological +skeleton**, because it is a 1-pixel wide skeleton of the object, with the same +connectivity as the original object. + +Here, we use the medial axis transform to compute the width of the foreground +objects. As the function ``medial_axis`` (``skimage.morphology.medial_axis``) +returns the distance transform in addition to the medial axis (with the keyword +argument ``return_distance=True``), it is possible to compute the distance to +the background for all points of the medial axis with this function. This gives +an estimate of the local width of the objects. + +For a skeleton with fewer branches, there exists another skeletonization +algorithm in ``skimage``: ``skimage.morphology.skeletonize``, that computes +a skeleton by iterative morphological thinnings. +""" + +import numpy as np +from scipy import ndimage +from skimage.morphology import medial_axis +import matplotlib.pyplot as plt + + +def microstructure(l=256): + """ + Synthetic binary data: binary microstructure with blobs. + + Parameters + ---------- + + l: int, optional + linear size of the returned image + + """ + n = 5 + x, y = np.ogrid[0:l, 0:l] + mask_outer = (x - l/2)**2 + (y - l/2)**2 < (l/2)**2 + mask = np.zeros((l, l)) + generator = np.random.RandomState(1) + points = l * generator.rand(2, n**2) + mask[(points[0]).astype(np.int), (points[1]).astype(np.int)] = 1 + mask = ndimage.gaussian_filter(mask, sigma=l/(4.*n)) + return mask > mask.mean() + +data = microstructure(l=64) + +# Compute the medial axis (skeleton) and the distance transform +skel, distance = medial_axis(data, return_distance=True) + +# Distance to the background for pixels of the skeleton +dist_on_skel = distance * skel + +plt.figure(figsize=(8, 4)) +plt.subplot(121) +plt.imshow(data, cmap=plt.cm.gray, interpolation='nearest') +plt.axis('off') +plt.subplot(122) +plt.imshow(dist_on_skel, cmap=plt.cm.spectral, interpolation='nearest') +plt.contour(data, [0.5], colors='w') +plt.axis('off') + +plt.subplots_adjust(hspace=0.01, wspace=0.01, top=1, bottom=0, left=0, + right=1) +plt.show() + diff --git a/doc/examples/plot_skeleton.py b/doc/examples/plot_skeleton.py index a61bbbdd..7cb33d95 100644 --- a/doc/examples/plot_skeleton.py +++ b/doc/examples/plot_skeleton.py @@ -15,8 +15,8 @@ results. The input is a 2D ndarray, with either boolean or integer elements. In the case of boolean, 'True' indicates foreground, and for integer arrays, the foreground is 1's. """ -from scikits.image.morphology import skeletonize -from scikits.image.draw import draw +from skimage.morphology import skeletonize +from skimage.draw import draw import numpy as np import matplotlib.pyplot as plt diff --git a/skimage/morphology/__init__.py b/skimage/morphology/__init__.py index d57cc859..efb92507 100644 --- a/skimage/morphology/__init__.py +++ b/skimage/morphology/__init__.py @@ -2,4 +2,4 @@ from grey import * from selem import * from .ccomp import label from watershed import watershed, is_local_maximum -from skeletonize import skeletonize +from skeletonize import skeletonize, medial_axis diff --git a/skimage/morphology/_skeletonize.pyx b/skimage/morphology/_skeletonize.pyx new file mode 100644 index 00000000..ff5fcdf2 --- /dev/null +++ b/skimage/morphology/_skeletonize.pyx @@ -0,0 +1,212 @@ +''' +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 +''' + +import numpy as np +cimport numpy as np +cimport cython + + +@cython.boundscheck(False) +def _skeletonize_loop(np.ndarray[dtype=np.uint8_t, ndim=2, + negative_indices=False, mode='c'] result, + np.ndarray[dtype=np.int32_t, ndim=1, + negative_indices=False, mode='c'] i, + np.ndarray[dtype=np.int32_t, ndim=1, + negative_indices=False, mode='c'] j, + np.ndarray[dtype=np.int32_t, ndim=1, + negative_indices=False, mode='c'] order, + np.ndarray[dtype=np.uint8_t, ndim=1, + negative_indices=False, mode='c'] 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: + np.int32_t accumulator + np.int32_t index, order_index + np.int32_t ii, jj + + 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 < result.shape[1] - 1 and result[ii - 1, jj + 1]: + accumulator += 4 + if jj > 0 and result[ii, jj - 1]: + accumulator += 8 + if jj < result.shape[1] - 1 and result[ii, jj + 1]: + accumulator += 32 + if ii < result.shape[0]-1: + if jj > 0 and result[ii + 1, jj - 1]: + accumulator += 64 + if result[ii + 1, jj]: + accumulator += 128 + if jj < result.shape[1] - 1 and result[ii + 1, jj + 1]: + accumulator += 256 + # Assign the value of table corresponding to the configuration + result[ii, jj] = table[accumulator] + +@cython.boundscheck(False) +def _table_lookup_index(np.ndarray[dtype=np.uint8_t, ndim=2, + negative_indices=False, mode='c'] 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: + np.ndarray[dtype=np.int32_t, ndim=2, + negative_indices=False, mode='c'] indexer + np.int32_t *p_indexer + np.uint8_t *p_image + np.int32_t i_stride + np.int32_t i_shape + np.int32_t j_shape + np.int32_t i + np.int32_t j + np.int32_t offset + + i_shape = image.shape[0] + j_shape = image.shape[1] + indexer = np.zeros((i_shape, j_shape), np.int32) + p_indexer = indexer.data + p_image = image.data + i_stride = image.strides[0] + assert i_shape >= 3 and j_shape >= 3, \ + "Please use the slow method for arrays < 3x3" + + 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 indexer diff --git a/skimage/morphology/setup.py b/skimage/morphology/setup.py index 9d8b5527..46010989 100644 --- a/skimage/morphology/setup.py +++ b/skimage/morphology/setup.py @@ -14,6 +14,7 @@ def configuration(parent_package='', top_path=None): cython(['ccomp.pyx'], working_path=base_path) cython(['cmorph.pyx'], working_path=base_path) cython(['_watershed.pyx'], working_path=base_path) + cython(['_skeletonize.pyx'], working_path=base_path) config.add_extension('ccomp', sources=['ccomp.c'], include_dirs=[get_numpy_include_dirs()]) @@ -21,6 +22,9 @@ def configuration(parent_package='', top_path=None): include_dirs=[get_numpy_include_dirs()]) config.add_extension('_watershed', sources=['_watershed.c'], include_dirs=[get_numpy_include_dirs()]) + config.add_extension('_skeletonize', sources=['_skeletonize.c'], + include_dirs=[get_numpy_include_dirs()]) + return config diff --git a/skimage/morphology/skeletonize.py b/skimage/morphology/skeletonize.py index ea478a9a..5ddf2c4c 100644 --- a/skimage/morphology/skeletonize.py +++ b/skimage/morphology/skeletonize.py @@ -1,10 +1,13 @@ -"""Use an iterative thinning algorithm to find the skeletons of binary -objects in an image. - +""" +Algorithms for computing the skeleton of a binary image """ import numpy as np -from scipy.ndimage import correlate +from scipy import ndimage + +from _skeletonize import _skeletonize_loop, _table_lookup_index + +# --------- Skeletonization by morphological thinning --------- def skeletonize(image): """Return the skeleton of a binary image. @@ -24,6 +27,10 @@ def skeletonize(image): skeleton : ndarray A matrix containing the thinned image. + See also + -------- + medial_axis + Notes ----- The algorithm [1] works by making successive passes of the image, @@ -107,7 +114,7 @@ def skeletonize(image): pixelRemoved = False; # assign each pixel a unique value based on its foreground neighbours - neighbours = correlate(skeleton, mask, mode='constant') + neighbours = ndimage.correlate(skeleton, mask, mode='constant') # ignore background neighbours *= skeleton @@ -126,7 +133,7 @@ def skeletonize(image): skeleton[code_mask] = 0 # pass 2 - remove the 2's and 3's - neighbours = correlate(skeleton, mask, mode='constant') + neighbours = ndimage.correlate(skeleton, mask, mode='constant') neighbours *= skeleton codes = np.take(lut, neighbours) code_mask = (codes == 2) @@ -139,3 +146,226 @@ def skeletonize(image): skeleton[code_mask] = 0 return skeleton + +# --------- Skeletonization by medial axis transform -------- + +_eight_connect = ndimage.generate_binary_structure(2, 2) + + +def medial_axis(image, mask=None, return_distance=False): + """ + Compute the medial axis transform of a binary image + + Parameters + ---------- + + image : binary ndarray + + mask : binary ndarray, optional + If a mask is given, only those elements 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. + + Returns + ------- + + out : ndarray of bools + Medial axis transform of the image + + dist : ndarray of ints + Distance transform of the image (only returned if `return_distance` + is True) + + See also + -------- + skeletonize + + 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=np.uint8) + >>> square[1:-1, 2:-2] = 1 + >>> square + 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) + >>> morphology.medial_axis(square).astype(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(np.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 neighbourhood 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([ndimage.label(_pattern_of(index), _eight_connect)[1] != + ndimage.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 = ndimage.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 neighbourhood + 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], np.int32) + j = np.ascontiguousarray(j[result], np.int32) + 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.RandomState(0) + tiebreaker = generator.permutation(np.arange(masked_image.sum())) + order = np.lexsort((tiebreaker, + corner_score[masked_image], + distance)) + order = np.ascontiguousarray(order, np.int32) + + table = np.ascontiguousarray(table, np.uint8) + # Remove pixels not belonging to the medial axis + _skeletonize_loop(result, i, j, order, table) + + result = result.astype(bool) + if not mask is 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. + border_value : bool + The value of pixels beyond the border of the image. + + 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 + diff --git a/skimage/morphology/tests/test_skeletonize.py b/skimage/morphology/tests/test_skeletonize.py index 8d900cd7..539e0ad8 100644 --- a/skimage/morphology/tests/test_skeletonize.py +++ b/skimage/morphology/tests/test_skeletonize.py @@ -1,5 +1,5 @@ import numpy as np -from skimage.morphology import skeletonize +from skimage.morphology import skeletonize, medial_axis import numpy.testing from skimage.draw import draw from scipy.ndimage import correlate @@ -90,6 +90,64 @@ class TestSkeletonize(): blocks = correlate(result, mask, mode='constant') assert not numpy.any(blocks == 4) +class TestMedialAxis(): + def test_00_00_zeros(self): + '''Test skeletonize on an array of all zeros''' + result = medial_axis(np.zeros((10, 10), bool)) + assert np.all(result == False) + + def test_00_01_zeros_masked(self): + '''Test skeletonize on an array that is completely masked''' + result = medial_axis(np.zeros((10, 10), bool), + np.zeros((10, 10), bool)) + assert np.all(result == False) + + def test_01_01_rectangle(self): + '''Test skeletonize on a rectangle''' + image = np.zeros((9, 15), bool) + image[1:-1, 1:-1] = True + # + # The result should be four diagonals from the + # corners, meeting in a horizontal line + # + expected = np.array([[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,1,0], + [0,0,1,0,0,0,0,0,0,0,0,0,1,0,0], + [0,0,0,1,0,0,0,0,0,0,0,1,0,0,0], + [0,0,0,0,1,1,1,1,1,1,1,0,0,0,0], + [0,0,0,1,0,0,0,0,0,0,0,1,0,0,0], + [0,0,1,0,0,0,0,0,0,0,0,0,1,0,0], + [0,1,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,0,0]], bool) + result = medial_axis(image) + assert np.all(result == expected) + result, distance = medial_axis(image, return_distance=True) + assert distance.max() == 4 + + def test_01_02_hole(self): + '''Test skeletonize on a rectangle with a hole in the middle''' + image = np.zeros((9, 15), bool) + image[1:-1, 1:-1] = True + image[4, 4:-4] = False + expected = np.array([[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,1,0], + [0,0,1,1,1,1,1,1,1,1,1,1,1,0,0], + [0,0,1,0,0,0,0,0,0,0,0,0,1,0,0], + [0,0,1,0,0,0,0,0,0,0,0,0,1,0,0], + [0,0,1,0,0,0,0,0,0,0,0,0,1,0,0], + [0,0,1,1,1,1,1,1,1,1,1,1,1,0,0], + [0,1,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,0,0]],bool) + result = medial_axis(image) + assert np.all(result == expected) + + def test_narrow_image(self): + """Test skeletonize on a 1-pixel thin strip""" + image = np.zeros((1, 5), bool) + image[:, 1:-1] = True + result = medial_axis(image) + assert np.all(result == image) + if __name__ == '__main__': np.testing.run_module_suite() diff --git a/skimage/transform/setup.py b/skimage/transform/setup.py index 58cfaf32..34b76832 100644 --- a/skimage/transform/setup.py +++ b/skimage/transform/setup.py @@ -23,6 +23,7 @@ def configuration(parent_package='', top_path=None): config.add_extension('_project', sources=['_project.c'], include_dirs=[get_numpy_include_dirs()]) + return config if __name__ == '__main__':