Issue #803 order the structuring element by distance to get diagonals processed last

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

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: