From a9bdaaf2747a7c786db1d494aea9312563b0c2ae Mon Sep 17 00:00:00 2001 From: kpk09 Date: Wed, 20 May 2015 16:41:32 +0100 Subject: [PATCH 1/2] proposed fix and test for broken watershed on flat areas (https://github.com/scikit-image/scikit-image/issues/803) --- skimage/morphology/_watershed.pyx | 7 ++++--- skimage/morphology/tests/test_watershed.py | 16 ++++++++++++++++ 2 files changed, 20 insertions(+), 3 deletions(-) diff --git a/skimage/morphology/_watershed.pyx b/skimage/morphology/_watershed.pyx index a0857707..812c4310 100644 --- a/skimage/morphology/_watershed.pyx +++ b/skimage/morphology/_watershed.pyx @@ -83,10 +83,11 @@ def watershed(DTYPE_INT32_t[::1] image, 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 = age + new_elem.index = index + output[index] = output[old_index] # # Push the neighbor onto the heap to work on it later diff --git a/skimage/morphology/tests/test_watershed.py b/skimage/morphology/tests/test_watershed.py index 788d999e..1bba6836 100644 --- a/skimage/morphology/tests/test_watershed.py +++ b/skimage/morphology/tests/test_watershed.py @@ -385,6 +385,22 @@ class TestWatershed(unittest.TestCase): scipy.ndimage.watershed_ift(image.astype(np.uint16), markers, self.eight) + def test_watershed10(self): + "watershed 10" + data = np.array([[1, 1, 1, 1], + [1, 1, 1, 1], + [1, 1, 1, 1], + [1, 1, 1, 1]], np.uint8) + markers = np.array([[1, 0, 0, 2], + [0, 0, 0, 0], + [0, 0, 0, 0], + [3, 0, 0, 4]], np.int8) + out = watershed(data, markers, self.eight) + error = diff([[1, 1, 2, 2], + [1, 1, 2, 2], + [3, 3, 4, 4], + [3, 3, 4, 4]], out) + self.assertTrue(error < eps) if __name__ == "__main__": np.testing.run_module_suite() From 54f8265603c1ec4f873f7aeda266e4bba157214e Mon Sep 17 00:00:00 2001 From: Lee Kamentsky Date: Tue, 26 May 2015 15:24:32 -0400 Subject: [PATCH 2/2] Issue #803 order the structuring element by distance to get diagonals processed last --- skimage/morphology/tests/test_watershed.py | 25 ++++++++++++++++++++++ skimage/morphology/watershed.py | 7 +++++- 2 files changed, 31 insertions(+), 1 deletion(-) diff --git a/skimage/morphology/tests/test_watershed.py b/skimage/morphology/tests/test_watershed.py index 1bba6836..39b464cd 100644 --- a/skimage/morphology/tests/test_watershed.py +++ b/skimage/morphology/tests/test_watershed.py @@ -402,5 +402,30 @@ class TestWatershed(unittest.TestCase): [3, 3, 4, 4]], out) self.assertTrue(error < eps) + def test_watershed11(self): + '''Make sure that all points on this plateau are assigned to closest seed''' + # https://github.com/scikit-image/scikit-image/issues/803 + # + # Make sure that no point in a level image is farther away + # from its seed than any other + # + image = np.zeros((21, 21)) + markers = np.zeros((21, 21), int) + markers[5, 5] = 1 + markers[5, 10] = 2 + markers[10, 5] = 3 + markers[10, 10] = 4 + + structure = np.array([[False, True, False], + [True, True, True], + [False, True, False]]) + out = watershed(image, markers, structure) + i, j = np.mgrid[0:21, 0:21] + d = np.dstack( + [np.sqrt((i.astype(float)-i0)**2, (j.astype(float)-j0)**2) + for i0, j0 in ((5, 5), (5, 10), (10, 5), (10, 10))]) + dmin = np.min(d, 2) + self.assertTrue(np.all(d[i, j, out[i, j]-1] == dmin)) + if __name__ == "__main__": np.testing.run_module_suite() diff --git a/skimage/morphology/watershed.py b/skimage/morphology/watershed.py index 5e4a9156..d788d0e2 100644 --- a/skimage/morphology/watershed.py +++ b/skimage/morphology/watershed.py @@ -186,6 +186,7 @@ def watershed(image, markers, connectivity=None, offset=None, mask=None): # and the second through last are the x,y...whatever offsets # (to do bounds checking). c = [] + distances = [] image_stride = np.array(image.strides) // image.itemsize for i in range(np.product(c_connectivity.shape)): multiplier = 1 @@ -202,10 +203,14 @@ def watershed(image, markers, connectivity=None, offset=None, mask=None): multiplier *= c_connectivity.shape[j] if (not ignore) and c_connectivity.__getitem__(tuple(indexes)): stride = np.dot(image_stride, np.array(offs)) + d = np.sum(np.abs(offs)) - 1 offs.insert(0, stride) c.append(offs) + distances.append(d) + c = np.array(c, dtype=np.int32) - + c = c[np.argsort(distances)] + pq, age = __heapify_markers(c_markers, c_image) pq = np.ascontiguousarray(pq, dtype=np.int32) if np.product(pq.shape) > 0: