Merge pull request #1530 from LeeKamentsky/master

Watershed fixes for issue #803 (incorrect propagation in plateaus)

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

skimage/morphology/tests/test_watershed.py

@@ -385,6 +385,47 @@ 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)
+
+    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()

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: