proposed fix and test for broken watershed on flat areas (https://github.com/scikit-image/scikit-image/issues/803)

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,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()