ENH address comments on PR #65 made by Stefan and Neil

This commit is contained in:
emmanuelle
2011-10-24 23:23:52 +02:00
parent 736802bad4
commit ed1f05427e
4 changed files with 52 additions and 46 deletions
+1 -1
View File
@@ -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.
"""
@@ -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.
"""
+4
View File
@@ -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
+33 -25
View File
@@ -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