From 2b4854a977fde4a4bd9c3f4747d5543ac919d00d Mon Sep 17 00:00:00 2001 From: emmanuelle Date: Thu, 6 Oct 2011 09:04:47 +0200 Subject: [PATCH] Added the is_local_maximum function, that is useful to determine markers for the watershed. --- scikits/image/morphology/_watershed.pyx | 38 +++--- .../image/morphology/tests/test_watershed.py | 101 ++++++++++++++- scikits/image/morphology/watershed.py | 119 +++++++++++++++++- 3 files changed, 241 insertions(+), 17 deletions(-) diff --git a/scikits/image/morphology/_watershed.pyx b/scikits/image/morphology/_watershed.pyx index f8afcc27..c86d8744 100644 --- a/scikits/image/morphology/_watershed.pyx +++ b/scikits/image/morphology/_watershed.pyx @@ -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 diff --git a/scikits/image/morphology/tests/test_watershed.py b/scikits/image/morphology/tests/test_watershed.py index 8499442c..98c8100a 100644 --- a/scikits/image/morphology/tests/test_watershed.py +++ b/scikits/image/morphology/tests/test_watershed.py @@ -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)) + + + diff --git a/scikits/image/morphology/watershed.py b/scikits/image/morphology/watershed.py index 95721e5f..a159a6ab 100644 --- a/scikits/image/morphology/watershed.py +++ b/scikits/image/morphology/watershed.py @@ -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