From ed1f05427e741175c17467fabbc20f0def6b18aa Mon Sep 17 00:00:00 2001 From: emmanuelle Date: Mon, 24 Oct 2011 21:31:24 +0200 Subject: [PATCH] ENH address comments on PR #65 made by Stefan and Neil --- doc/examples/plot_medial_transform.py | 2 +- .../{_cpmorphology2.pyx => _skeletonize.pyx} | 34 +++++------ skimage/morphology/setup.py | 4 ++ skimage/morphology/skeletonize.py | 58 +++++++++++-------- 4 files changed, 52 insertions(+), 46 deletions(-) rename scikits/image/morphology/{_cpmorphology2.pyx => _skeletonize.pyx} (92%) diff --git a/doc/examples/plot_medial_transform.py b/doc/examples/plot_medial_transform.py index 47fcd5bd..af0b4c38 100644 --- a/doc/examples/plot_medial_transform.py +++ b/doc/examples/plot_medial_transform.py @@ -15,7 +15,7 @@ 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 less branches, there exists another skeletonization +For a skeleton with fewer branches, there exists another skeletonization algorithm in ``skimage``: ``skimage.morphology.skeletonize``, that computes a skeleton by iterative morphological thinnings. """ diff --git a/scikits/image/morphology/_cpmorphology2.pyx b/scikits/image/morphology/_skeletonize.pyx similarity index 92% rename from scikits/image/morphology/_cpmorphology2.pyx rename to scikits/image/morphology/_skeletonize.pyx index 23f9736f..0b2f8a8c 100644 --- a/scikits/image/morphology/_cpmorphology2.pyx +++ b/scikits/image/morphology/_skeletonize.pyx @@ -14,21 +14,9 @@ import numpy as np cimport numpy as np cimport cython -cdef extern from "Python.h": - ctypedef int Py_intptr_t - -cdef extern from "numpy/arrayobject.h": - ctypedef class numpy.ndarray [object PyArrayObject]: - cdef char *data - cdef Py_intptr_t *dimensions - cdef Py_intptr_t *strides - cdef void import_array() - cdef int PyArray_ITEMSIZE(np.ndarray) - -import_array() @cython.boundscheck(False) -def skeletonize_loop(np.ndarray[dtype=np.uint8_t, ndim=2, +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, @@ -44,15 +32,18 @@ def skeletonize_loop(np.ndarray[dtype=np.uint8_t, ndim=2, Parameters ---------- - result: ndarray of uint8 + result : ndarray of uint8 On input, the image to be skeletonized, on output the skeletonized image. - i, j: ndarrays + + i, j : ndarrays The coordinates of each foreground pixel in the image - order: ndarray + + order : ndarray The index of each pixel, in the order of processing (order[0] is the first pixel to process, etc.) - table: ndarray + + table : ndarray The 512-element lookup table of values after transformation (whether to keep or not each configuration in a binary 3x3 array) @@ -103,21 +94,24 @@ def skeletonize_loop(np.ndarray[dtype=np.uint8_t, ndim=2, result[ii, jj] = table[accumulator] @cython.boundscheck(False) -def table_lookup_index(np.ndarray[dtype=np.uint8_t, ndim=2, +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: + looks like this:: + 1 2 4 8 16 32 64 128 256 - This code could be replaced by a convolution with the kernel: + 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. """ 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 748d1b80..b7c54c9e 100644 --- a/skimage/morphology/skeletonize.py +++ b/skimage/morphology/skeletonize.py @@ -5,7 +5,7 @@ Algorithms for computing the skeleton of a binary image import numpy as np from scipy import ndimage -from _cpmorphology2 import skeletonize_loop, table_lookup_index +from _skeletonize import _skeletonize_loop, _table_lookup_index # --------- Skeletonization by morphological thinning --------- @@ -149,7 +149,7 @@ def skeletonize(image): # --------- Skeletonization by medial axis transform -------- -eight_connect = ndimage.generate_binary_structure(2, 2) +_eight_connect = ndimage.generate_binary_structure(2, 2) def medial_axis(image, mask=None, return_distance=False): @@ -159,22 +159,22 @@ def medial_axis(image, mask=None, return_distance=False): Parameters ---------- - image: binary ndarray + image : binary ndarray - mask: binary ndarray, optional + 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 + return_distance : bool, optional If true, the distance transform is returned as well as the skeleton. Returns ------- - out: ndarray of bools + out : ndarray of bools Medial axis transform of the image - dist: ndarray of ints + dist : ndarray of ints Distance transform of the image (only returned if `return_distance` is True) @@ -201,8 +201,9 @@ def medial_axis(image, mask=None, return_distance=False): * 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 not a pixel according to the lookup table. Because of the - ordering, it is possible to process all pixels in only one pass. + 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 -------- @@ -226,7 +227,7 @@ def medial_axis(image, mask=None, return_distance=False): [0, 0, 0, 0, 0, 0, 0]], dtype=uint8) """ - global eight_connect + global _eight_connect if mask is None: masked_image = image.astype(np.bool) else: @@ -243,12 +244,17 @@ def medial_axis(image, mask=None, return_distance=False): # 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 & - (np.array([ndimage.label(_pattern_of(index), eight_connect)[1] != + 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)]) | - np.array([np.sum(_pattern_of(index)) < 3 for index in range(512)]))) + _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) @@ -289,7 +295,7 @@ def medial_axis(image, mask=None, return_distance=False): table = np.ascontiguousarray(table, np.uint8) # Remove pixels not belonging to the medial axis - skeletonize_loop(result, i, j, order, table) + _skeletonize_loop(result, i, j, order, table) result = result.astype(bool) if not mask is None: @@ -316,20 +322,22 @@ def _table_lookup(image, table): Parameters ---------- - image - a binary image - table - a 512-element table giving the transform of each pixel given - the values of that pixel and its 8-connected neighbors. - border_value - the value of pixels beyond the border of the image. - This should test as True or False. + 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` + result : ndarray of same shape as `image` Transformed image Notes ----- - The pixels are numbered like this: + The pixels are numbered like this:: 0 1 2 3 4 5 @@ -343,7 +351,7 @@ def _table_lookup(image, table): # if image.shape[0] < 3 or image.shape[1] < 3: image = image.astype(bool) - indexer = np.zeros(image.shape,int) + 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 @@ -356,7 +364,7 @@ def _table_lookup(image, table): indexer[:-1, :] += image[1:, :] * 2**7 indexer[:-1, :-1] += image[1:, 1:] * 2**8 else: - indexer = table_lookup_index(np.ascontiguousarray(image, np.uint8)) + indexer = _table_lookup_index(np.ascontiguousarray(image, np.uint8)) image = table[indexer] return image