wassname/scikit-image
1
0
Fork 0
mirror of https://github.com/wassname/scikit-image.git synced 2026-07-27 11:27:08 +08:00
Code Issues Packages Projects Releases Wiki Activity

Added the is_local_maximum function, that is useful to determine

Browse Source 
markers for the watershed.
This commit is contained in:
emmanuelle
2011-10-06 09:04:47 +02:00
parent f0702fb938
commit 2b4854a977
3 changed files with 241 additions and 17 deletions
Download Patch File Download Diff File Expand all files Collapse all files
scikits/image/morphology/_watershed.pyx
+24 -14
View File
@@ -26,16 +26,25 @@ ctypedef np.int8_t DTYPE_BOOL_t
include "heap_watershed.pxi"
@cython.boundscheck(False)
def watershed(np.ndarray[DTYPE_INT32_t,ndim=1,negative_indices=False, mode='c'] image,
np.ndarray[DTYPE_INT32_t,ndim=2,negative_indices=False, mode='c'] pq,
def watershed(np.ndarray[DTYPE_INT32_t, ndim=1, negative_indices=False,
mode='c'] image,
np.ndarray[DTYPE_INT32_t, ndim=2, negative_indices=False,
mode='c'] pq,
DTYPE_INT32_t age,
np.ndarray[DTYPE_INT32_t,ndim=2,negative_indices=False, mode='c'] structure,
np.ndarray[DTYPE_INT32_t, ndim=2, negative_indices=False,
mode='c'] structure,
DTYPE_INT32_t ndim,
np.ndarray[DTYPE_BOOL_t,ndim=1,negative_indices=False, mode='c'] mask,
np.ndarray[DTYPE_INT32_t,ndim=1,negative_indices=False, mode='c'] image_shape,
np.ndarray[DTYPE_INT32_t,ndim=1,negative_indices=False, mode='c'] output):
np.ndarray[DTYPE_BOOL_t, ndim=1, negative_indices=False,
mode='c'] mask,
np.ndarray[DTYPE_INT32_t, ndim=1, negative_indices=False,
mode='c'] image_shape,
np.ndarray[DTYPE_INT32_t, ndim=1, negative_indices=False,
mode='c'] output):
"""Do heavy lifting of watershed algorithm
Parameters
----------
image - the flattened image pixels, converted to rank-order
pq - the priority queue, starts with the marked pixels
the first element in each row is the image intensity
@@ -49,7 +58,7 @@ def watershed(np.ndarray[DTYPE_INT32_t,ndim=1,negative_indices=False, mode='c']
in a flattened array. The remaining elements are the
offsets from the point to its neighbor in the various
dimensions
ndim - # of dimensions in the image
ndim - # of dimensions in the image
mask - numpy boolean (char) array indicating which pixels to consider
and which to ignore. Also flattened.
image_shape - the dimensions of the image, for boundary checking,
@@ -83,14 +92,15 @@ def watershed(np.ndarray[DTYPE_INT32_t,ndim=1,negative_indices=False, mode='c']
old_index = elem.index
for i in range(nneighbors):
# get the flattened address of the neighbor
index = structure[i,0]+old_index
if index < 0 or index >= max_index or output[index] or not mask[index]:
index = structure[i, 0] + old_index
if index < 0 or index >= max_index or output[index] or \
not mask[index]:
continue
new_elem.value = image[index]
new_elem.age = elem.age + 1
new_elem.index = index
age += 1
new_elem.value = image[index]
new_elem.age = elem.age + 1
new_elem.index = index
age += 1
output[index] = output[old_index]
#
# Push the neighbor onto the heap to work on it later
scikits/image/morphology/tests/test_watershed.py
+100 -1
View File
@@ -50,7 +50,7 @@ import numpy as np
import scipy.ndimage
from scikits.image.morphology.watershed import watershed,fast_watershed, \
_slow_watershed
_slow_watershed, is_local_maximum
eps = 1e-12
@@ -394,3 +394,102 @@ class TestFastWatershed(unittest.TestCase):
elapsed = time.clock() - before
print "Scipy watershed ran a megapixel image in %f seconds"%(elapsed)
class TestIsLocalMaximum(unittest.TestCase):
def test_00_00_empty(self):
image = np.zeros((10, 20))
labels = np.zeros((10, 20), int)
result = is_local_maximum(image, labels, np.ones((3, 3), bool))
self.assertTrue(np.all(~ result))
def test_01_01_one_point(self):
image = np.zeros((10, 20))
labels = np.zeros((10, 20), int)
image[5, 5] = 1
labels[5, 5] = 1
result = is_local_maximum(image, labels, np.ones((3, 3), bool))
self.assertTrue(np.all(result == (labels == 1)))
def test_01_02_adjacent_and_same(self):
image = np.zeros((10, 20))
labels = np.zeros((10, 20), int)
image[5, 5:6] = 1
labels[5, 5:6] = 1
result = is_local_maximum(image, labels, np.ones((3, 3), bool))
self.assertTrue(np.all(result == (labels == 1)))
def test_01_03_adjacent_and_different(self):
image = np.zeros((10, 20))
labels = np.zeros((10, 20), int)
image[5, 5] = 1
image[5, 6] = .5
labels[5, 5:6] = 1
expected = (image == 1)
result = is_local_maximum(image, labels, np.ones((3, 3), bool))
self.assertTrue(np.all(result == expected))
result = is_local_maximum(image, labels)
self.assertTrue(np.all(result == expected))
def test_01_04_not_adjacent_and_different(self):
image = np.zeros((10,20))
labels = np.zeros((10,20), int)
image[5,5] = 1
image[5,8] = .5
labels[image > 0] = 1
expected = (labels == 1)
result = is_local_maximum(image, labels, np.ones((3,3), bool))
self.assertTrue(np.all(result == expected))
def test_01_05_two_objects(self):
image = np.zeros((10,20))
labels = np.zeros((10,20), int)
image[5,5] = 1
image[5,15] = .5
labels[5,5] = 1
labels[5,15] = 2
expected = (labels > 0)
result = is_local_maximum(image, labels, np.ones((3,3), bool))
self.assertTrue(np.all(result == expected))
def test_01_06_adjacent_different_objects(self):
image = np.zeros((10,20))
labels = np.zeros((10,20), int)
image[5,5] = 1
image[5,6] = .5
labels[5,5] = 1
labels[5,6] = 2
expected = (labels > 0)
result = is_local_maximum(image, labels, np.ones((3,3), bool))
self.assertTrue(np.all(result == expected))
def test_02_01_four_quadrants(self):
np.random.seed(21)
image = np.random.uniform(size=(40,60))
i,j = np.mgrid[0:40,0:60]
labels = 1 + (i >= 20) + (j >= 30) * 2
i,j = np.mgrid[-3:4,-3:4]
footprint = (i*i + j*j <=9)
expected = np.zeros(image.shape, float)
for imin, imax in ((0, 20), (20, 40)):
for jmin, jmax in ((0, 30), (30, 60)):
expected[imin:imax,jmin:jmax] = scipy.ndimage.maximum_filter(
image[imin:imax, jmin:jmax], footprint = footprint)
expected = (expected == image)
result = is_local_maximum(image, labels, footprint)
self.assertTrue(np.all(result == expected))
def test_03_01_disk_1(self):
'''regression test of img-1194, footprint = [1]
Test is_local_maximum when every point is a local maximum
'''
np.random.seed(31)
image = np.random.uniform(size=(10,20))
footprint = np.array([[1]])
result = is_local_maximum(image, np.ones((10,20)), footprint)
self.assertTrue(np.all(result))
result = is_local_maximum(image, footprint=footprint)
self.assertTrue(np.all(result))
scikits/image/morphology/watershed.py
+117 -2
View File
@@ -93,7 +93,7 @@ def fast_watershed(image, markers, connectivity=None, offset=None, mask=None):
See also
--------
`regional_maximum`, `is_local_maximum`, `scipy.ndimage.label`
`is_local_maximum`, `scipy.ndimage.label`
Examples
--------
@@ -111,7 +111,7 @@ def fast_watershed(image, markers, connectivity=None, offset=None, mask=None):
>>> # to the background
>>> from scipy import ndimage
>>> distance = ndimage.distance_transform_edt(image)
>>> local_maxi = regional_maximum(distance)
>>> local_maxi = is_local_maximum(distance, image, np.ones((3, 3)))
>>> markers = ndimage.label(local_maxi)[0]
>>> labels = watershed(-distance, markers, mask=image)
@@ -224,6 +224,121 @@ def fast_watershed(image, markers, connectivity=None, offset=None, mask=None):
watershed = fast_watershed
def is_local_maximum(image, labels=None, footprint=None):
"""
Return a boolean array of points that are local maxima
Parameters
----------
image: ndarray (2-D, 3-D, ...)
intensity image
labels: ndarray, optional
find maxima only within labels. Zero is reserved for background.
footprint: ndarray of bools, optional
binary mask indicating the neighborhood to be examined
`footprint` must be a matrix with odd dimensions, the center is taken
to be the point in question.
Returns
-------
result: ndarray of bools
mask that is True for pixels that are local maxima of `image`
Examples
--------
>>> image = np.zeros((4, 4))
>>> image[1, 2] = 2
>>> image[3, 3] = 1
>>> image
array([[ 0., 0., 0., 0.],
[ 0., 0., 2., 0.],
[ 0., 0., 0., 0.],
[ 0., 0., 0., 1.]])
>>> is_local_maximum(image)
array([[ True, False, False, False],
[ True, False, True, False],
[ True, False, False, False],
[ True, True, False, True]], dtype=bool)
>>> image = np.arange(16).reshape((4, 4))
>>> labels = np.array([[1, 2], [3, 4]])
>>> labels = np.repeat(np.repeat(labels, 2, axis=0), 2, axis=1)
>>> labels
array([[1, 1, 2, 2],
[1, 1, 2, 2],
[3, 3, 4, 4],
[3, 3, 4, 4]])
>>> image
array([[ 0, 1, 2, 3],
[ 4, 5, 6, 7],
[ 8, 9, 10, 11],
[12, 13, 14, 15]])
>>> is_local_maximum(image, labels=labels)
array([[False, False, False, False],
[False, True, False, True],
[False, False, False, False],
[False, True, False, True]], dtype=bool)
"""
if labels is None:
labels = np.ones(image.shape, dtype=np.uint8)
if footprint is None:
footprint = np.ones([3] * image.ndim, dtype=np.uint8)
assert((np.all(footprint.shape) & 1) == 1)
footprint = (footprint != 0)
footprint_extent = (np.array(footprint.shape)-1) / 2
if np.all(footprint_extent == 0):
return labels > 0
result = (labels > 0).copy()
#
# Create a labels matrix with zeros at the borders that might be
# hit by the footprint.
#
big_labels = np.zeros(np.array(labels.shape) + footprint_extent*2,
labels.dtype)
big_labels[[slice(fe,-fe) for fe in footprint_extent]] = labels
#
# Find the relative indexes of each footprint element
#
image_strides = np.array(image.strides) / image.dtype.itemsize
big_strides = np.array(big_labels.strides) / big_labels.dtype.itemsize
result_strides = np.array(result.strides) / result.dtype.itemsize
footprint_offsets = np.mgrid[[slice(-fe,fe+1) for fe in footprint_extent]]
fp_image_offsets = np.sum(image_strides[:, np.newaxis] *
footprint_offsets[:, footprint], 0)
fp_big_offsets = np.sum(big_strides[:, np.newaxis] *
footprint_offsets[:, footprint], 0)
#
# Get the index of each labeled pixel in the image and big_labels arrays
#
indexes = np.mgrid[[slice(0,x) for x in labels.shape]][:, labels > 0]
image_indexes = np.sum(image_strides[:, np.newaxis] * indexes, 0)
big_indexes = np.sum(big_strides[:, np.newaxis] *
(indexes + footprint_extent[:, np.newaxis]), 0)
result_indexes = np.sum(result_strides[:, np.newaxis] * indexes, 0)
#
# Now operate on the raveled images
#
big_labels_raveled = big_labels.ravel()
image_raveled = image.ravel()
result_raveled = result.ravel()
#
# A hit is a hit if the label at the offset matches the label at the pixel
# and if the intensity at the pixel is greater or equal to the intensity
# at the offset.
#
for fp_image_offset, fp_big_offset in zip(fp_image_offsets, fp_big_offsets):
same_label = (big_labels_raveled[big_indexes + fp_big_offset] ==
big_labels_raveled[big_indexes])
less_than = (image_raveled[image_indexes[same_label]] <
image_raveled[image_indexes[same_label]+ fp_image_offset])
result_raveled[result_indexes[same_label][less_than]] = False
return result
# ---------------------- deprecated ------------------------------
# Deprecate slower pure-Python code, that we keep only for
# pedagogical purposes