From 54f219400b9b312bb8e3b91dc7c750cc1439ce09 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <juan.n@unimelb.edu.au>
Date: Thu, 21 Jul 2016 11:20:58 +1000
Subject: [PATCH 01/33] Replace <> with more common != in heap_watershed

---
 skimage/morphology/heap_watershed.pxi | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/skimage/morphology/heap_watershed.pxi b/skimage/morphology/heap_watershed.pxi
index 07b29f5cf10..caa32503256 100644
--- a/skimage/morphology/heap_watershed.pxi
+++ b/skimage/morphology/heap_watershed.pxi
@@ -19,7 +19,7 @@ cdef struct Heapitem:
 
 
 cdef inline int smaller(Heapitem *a, Heapitem *b):
-    if a.value <> b.value:
+    if a.value != b.value:
         return a.value < b.value
     return a.age < b.age
 

From 08e038bcbcaf6df565b4801504a2fdcc857aad08 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <juan.n@unimelb.edu.au>
Date: Thu, 21 Jul 2016 11:31:52 +1000
Subject: [PATCH 02/33] Use np.pad instead of bespoke one in watershed.py

---
 skimage/morphology/watershed.py | 22 ++++++++--------------
 1 file changed, 8 insertions(+), 14 deletions(-)

diff --git a/skimage/morphology/watershed.py b/skimage/morphology/watershed.py
index 8a39e19e784..34c41429c4b 100644
--- a/skimage/morphology/watershed.py
+++ b/skimage/morphology/watershed.py
@@ -141,22 +141,16 @@ def watershed(image, markers, connectivity=None, offset=None, mask=None):
         #
         offset = np.array(c_connectivity.shape) // 2
 
+    if mask is None:
+        # Use a complete `True` mask if none is provided
+        mask = np.ones(image.shape, bool)
+
     # pad the image, markers, and mask so that we can use the mask to
     # keep from running off the edges
-    pads = offset
-
-    def pad(im):
-        new_im = np.zeros(
-            [i + 2 * p for i, p in zip(im.shape, pads)], im.dtype)
-        new_im[[slice(p, -p, None) for p in pads]] = im
-        return new_im
-
-    if mask is not None:
-        mask = pad(mask)
-    else:
-        mask = pad(np.ones(image.shape, bool))
-    image = pad(image)
-    markers = pad(markers)
+    pad_width = [(p, p) for p in offset]
+    image = np.pad(image, pad_width, mode='constant')
+    mask = np.pad(mask, pad_width, mode='constant')
+    markers = np.pad(markers, pad_width, mode='constant')
 
     c_image = rank_order(image)[0].astype(np.int32)
     c_markers = np.ascontiguousarray(markers, dtype=np.int32)

From b8db4b2ab9c62198c8a011ebeb9369913f34e5bf Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <juan.n@unimelb.edu.au>
Date: Thu, 21 Jul 2016 12:53:41 +1000
Subject: [PATCH 03/33] Move input validation code out of watershed function

---
 skimage/morphology/watershed.py | 24 ++++++++++--------------
 1 file changed, 10 insertions(+), 14 deletions(-)

diff --git a/skimage/morphology/watershed.py b/skimage/morphology/watershed.py
index 34c41429c4b..5cb60542abe 100644
--- a/skimage/morphology/watershed.py
+++ b/skimage/morphology/watershed.py
@@ -32,6 +32,14 @@
 from . import _watershed
 
 
+def _validate_inputs(image, mask, markers):
+    if markers.shape != image.shape:
+        raise ValueError("Markers (shape %s) must have same shape "
+                         "as image (shape %s)" % (markers.ndim, image.ndim))
+    if mask is not None and mask.shape != image.shape:
+        raise ValueError("mask must have same shape as image")
+
+
 def watershed(image, markers, connectivity=None, offset=None, mask=None):
     """
     Return a matrix labeled using the watershed segmentation algorithm
@@ -125,6 +133,7 @@ def watershed(image, markers, connectivity=None, offset=None, mask=None):
     The algorithm works also for 3-D images, and can be used for example to
     separate overlapping spheres.
     """
+    _validate_inputs(image, markers, mask)
 
     if connectivity is None:
         c_connectivity = ndi.generate_binary_structure(image.ndim, 1)
@@ -154,20 +163,7 @@ def watershed(image, markers, connectivity=None, offset=None, mask=None):
 
     c_image = rank_order(image)[0].astype(np.int32)
     c_markers = np.ascontiguousarray(markers, dtype=np.int32)
-    if c_markers.ndim != c_image.ndim:
-        raise ValueError("markers (ndim=%d) must have same # of dimensions "
-                         "as image (ndim=%d)" % (c_markers.ndim, c_image.ndim))
-    if c_markers.shape != c_image.shape:
-        raise ValueError("image and markers must have the same shape")
-    if mask is not None:
-        c_mask = np.ascontiguousarray(mask, dtype=bool)
-        if c_mask.ndim != c_markers.ndim:
-            raise ValueError("mask must have same # of dimensions as image")
-        if c_markers.shape != c_mask.shape:
-            raise ValueError("mask must have same shape as image")
-        c_markers[np.logical_not(mask)] = 0
-    else:
-        c_mask = None
+    c_mask = np.ascontiguousarray(mask, dtype=bool)
     c_output = c_markers.copy()
 
     #

From 4c1707aa0a14643e569c544be211ce719f8baa81 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <juan.n@unimelb.edu.au>
Date: Thu, 21 Jul 2016 16:11:03 +1000
Subject: [PATCH 04/33] Move connectivity validation out of watershed main
 function

---
 skimage/morphology/watershed.py | 82 ++++++++++++++++++++++++++-------
 1 file changed, 65 insertions(+), 17 deletions(-)

diff --git a/skimage/morphology/watershed.py b/skimage/morphology/watershed.py
index 5cb60542abe..331b267fbf3 100644
--- a/skimage/morphology/watershed.py
+++ b/skimage/morphology/watershed.py
@@ -32,7 +32,23 @@
 from . import _watershed
 
 
-def _validate_inputs(image, mask, markers):
+def _validate_inputs(image, markers, mask):
+    """Ensure that all inputs to watershed have matching shapes.
+
+    Parameters
+    ----------
+    image : array
+        The input image.
+    markers : array
+        The marker image.
+    mask : array, or None
+        A boolean mask, True where we want to compute the watershed.
+
+    Raises
+    ------
+    ValueError
+        If the shapes of the given arrays don't match.
+    """
     if markers.shape != image.shape:
         raise ValueError("Markers (shape %s) must have same shape "
                          "as image (shape %s)" % (markers.ndim, image.ndim))
@@ -40,7 +56,52 @@ def _validate_inputs(image, mask, markers):
         raise ValueError("mask must have same shape as image")
 
 
-def watershed(image, markers, connectivity=None, offset=None, mask=None):
+def _validate_connectivity(image_dim, connectivity, offset):
+    """Convert any valid connectivity to a structuring element and offset.
+
+    Parameters
+    ----------
+    image_dim : int
+        The number of dimensions of the input image.
+    connectivity : int, array, or None
+        The neighborhood connectivity. An integer is interpreted as in
+        ``scipy.ndimage.generate_binary_structure``, as the maximum number
+        of orthogonal steps to reach a neighbor. An array is directly
+        interpreted as a structuring element and its shape is validated against
+        the input image shape. ``None`` is interpreted as a connectivity of 1.
+    offset : tuple of int, or None
+        The coordinates of the center of the structuring element.
+
+    Returns
+    -------
+    c_connectivity : array of bool
+        The structuring element corresponding to the input `connectivity`.
+    offset : array of int
+        The offset corresponding to the center of the structuring element.
+
+    Raises
+    ------
+    ValueError:
+        If the image dimension and the connectivity or offset dimensions don't
+        match.
+    """
+    if connectivity is None:
+        connectivity = 1
+    if np.isscalar(connectivity):
+        c_connectivity = ndi.generate_binary_structure(image_dim, connectivity)
+    else:
+        c_connectivity = np.array(connectivity, bool)
+        if c_connectivity.ndim != image_dim:
+            raise ValueError("Connectivity dimension must be same as image")
+    if offset is None:
+        if any([x % 2 == 0 for x in c_connectivity.shape]):
+            raise ValueError("Connectivity array must have an unambiguous "
+                             "center")
+        offset = np.array(c_connectivity.shape) // 2
+    return c_connectivity, offset
+
+
+def watershed(image, markers, connectivity=1, offset=None, mask=None):
     """
     Return a matrix labeled using the watershed segmentation algorithm
 
@@ -134,21 +195,8 @@ def watershed(image, markers, connectivity=None, offset=None, mask=None):
     separate overlapping spheres.
     """
     _validate_inputs(image, markers, mask)
-
-    if connectivity is None:
-        c_connectivity = ndi.generate_binary_structure(image.ndim, 1)
-    else:
-        c_connectivity = np.array(connectivity, bool)
-        if c_connectivity.ndim != image.ndim:
-            raise ValueError("Connectivity dimension must be same as image")
-    if offset is None:
-        if any([x % 2 == 0 for x in c_connectivity.shape]):
-            raise ValueError("Connectivity array must have an unambiguous "
-                             "center")
-        #
-        # offset to center of connectivity array
-        #
-        offset = np.array(c_connectivity.shape) // 2
+    c_connectivity, offset = _validate_connectivity(image.ndim, connectivity,
+                                                    offset)
 
     if mask is None:
         # Use a complete `True` mask if none is provided

From 1d5830374a56bee58877997757742639ae5cef51 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <juan.n@unimelb.edu.au>
Date: Thu, 21 Jul 2016 16:27:39 +1000
Subject: [PATCH 05/33] Use float values instead of int32s in watershed

---
 skimage/morphology/_watershed.pyx     | 8 ++++----
 skimage/morphology/heap_watershed.pxi | 2 +-
 skimage/morphology/watershed.py       | 2 +-
 3 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/skimage/morphology/_watershed.pyx b/skimage/morphology/_watershed.pyx
index 812c4310c97..7034f3aeb5d 100644
--- a/skimage/morphology/_watershed.pyx
+++ b/skimage/morphology/_watershed.pyx
@@ -10,20 +10,20 @@ All rights reserved.
 Original author: Lee Kamentsky
 """
 import numpy as np
-cimport numpy as np
+cimport numpy as cnp
 cimport cython
 
 
-ctypedef np.int32_t DTYPE_INT32_t
+ctypedef cnp.int32_t DTYPE_INT32_t
 DTYPE_BOOL = np.bool
-ctypedef np.int8_t DTYPE_BOOL_t
+ctypedef cnp.int8_t DTYPE_BOOL_t
 
 
 include "heap_watershed.pxi"
 
 
 @cython.boundscheck(False)
-def watershed(DTYPE_INT32_t[::1] image,
+def watershed(cnp.float64_t[::1] image,
               DTYPE_INT32_t[:, ::1] pq,
               Py_ssize_t age,
               DTYPE_INT32_t[:, ::1] structure,
diff --git a/skimage/morphology/heap_watershed.pxi b/skimage/morphology/heap_watershed.pxi
index caa32503256..aff627fd14c 100644
--- a/skimage/morphology/heap_watershed.pxi
+++ b/skimage/morphology/heap_watershed.pxi
@@ -13,7 +13,7 @@ cimport numpy as cnp
 
 
 cdef struct Heapitem:
-    cnp.int32_t value
+    cnp.float64_t value
     cnp.int32_t age
     Py_ssize_t index
 
diff --git a/skimage/morphology/watershed.py b/skimage/morphology/watershed.py
index 331b267fbf3..71c686f4e61 100644
--- a/skimage/morphology/watershed.py
+++ b/skimage/morphology/watershed.py
@@ -209,7 +209,7 @@ def watershed(image, markers, connectivity=1, offset=None, mask=None):
     mask = np.pad(mask, pad_width, mode='constant')
     markers = np.pad(markers, pad_width, mode='constant')
 
-    c_image = rank_order(image)[0].astype(np.int32)
+    c_image = image.astype(np.float64)
     c_markers = np.ascontiguousarray(markers, dtype=np.int32)
     c_mask = np.ascontiguousarray(mask, dtype=bool)
     c_output = c_markers.copy()

From afa2c6a92e5afc582aabfd962b3deb6a91fa798f Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <juan.n@unimelb.edu.au>
Date: Thu, 21 Jul 2016 16:46:37 +1000
Subject: [PATCH 06/33] Update the doctest example for watershed

---
 skimage/morphology/watershed.py | 15 ++++++++++-----
 1 file changed, 10 insertions(+), 5 deletions(-)

diff --git a/skimage/morphology/watershed.py b/skimage/morphology/watershed.py
index 71c686f4e61..4a1db3b1e2d 100644
--- a/skimage/morphology/watershed.py
+++ b/skimage/morphology/watershed.py
@@ -170,18 +170,20 @@ def watershed(image, markers, connectivity=1, offset=None, mask=None):
 
     Examples
     --------
-    The watershed algorithm is very useful to separate overlapping objects
+    The watershed algorithm is useful to separate overlapping objects.
+
+    We first generate an initial image with two overlapping circles:
 
-    >>> # Generate an initial image with two overlapping circles
     >>> x, y = np.indices((80, 80))
     >>> x1, y1, x2, y2 = 28, 28, 44, 52
     >>> r1, r2 = 16, 20
     >>> mask_circle1 = (x - x1)**2 + (y - y1)**2 < r1**2
     >>> mask_circle2 = (x - x2)**2 + (y - y2)**2 < r2**2
     >>> image = np.logical_or(mask_circle1, mask_circle2)
-    >>> # Now we want to separate the two objects in image
-    >>> # Generate the markers as local maxima of the distance
-    >>> # to the background
+
+    Next, we want to separate the two circles. We generate markers at the
+    maxima of the distance to the background:
+
     >>> from scipy import ndimage as ndi
     >>> distance = ndi.distance_transform_edt(image)
     >>> from skimage.feature import peak_local_max
@@ -189,6 +191,9 @@ def watershed(image, markers, connectivity=1, offset=None, mask=None):
     ...                             footprint=np.ones((3, 3)),
     ...                             indices=False)
     >>> markers = ndi.label(local_maxi)[0]
+
+    Finally, we run the watershed on the image and markers:
+
     >>> labels = watershed(-distance, markers, mask=image)
 
     The algorithm works also for 3-D images, and can be used for example to

From a773591630fa65688b7be3df205d8d2049205fdd Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <juan.n@unimelb.edu.au>
Date: Thu, 21 Jul 2016 16:47:00 +1000
Subject: [PATCH 07/33] Add source to the heap items

---
 skimage/morphology/_watershed.pyx     | 2 ++
 skimage/morphology/heap_watershed.pxi | 1 +
 2 files changed, 3 insertions(+)

diff --git a/skimage/morphology/_watershed.pyx b/skimage/morphology/_watershed.pyx
index 7034f3aeb5d..e2db0b2cd4f 100644
--- a/skimage/morphology/_watershed.pyx
+++ b/skimage/morphology/_watershed.pyx
@@ -65,6 +65,7 @@ def watershed(cnp.float64_t[::1] image,
         elem.value = pq[i, 0]
         elem.age = pq[i, 1]
         elem.index = pq[i, 2]
+        elem.source = pq[i, 2]
         heappush(hp, &elem)
 
     while hp.items > 0:
@@ -87,6 +88,7 @@ def watershed(cnp.float64_t[::1] image,
             new_elem.value = image[index]
             new_elem.age = age
             new_elem.index = index
+            new_elem.source = elem.source
             
             output[index] = output[old_index]
             #
diff --git a/skimage/morphology/heap_watershed.pxi b/skimage/morphology/heap_watershed.pxi
index aff627fd14c..830adabbdf7 100644
--- a/skimage/morphology/heap_watershed.pxi
+++ b/skimage/morphology/heap_watershed.pxi
@@ -16,6 +16,7 @@ cdef struct Heapitem:
     cnp.float64_t value
     cnp.int32_t age
     Py_ssize_t index
+    Py_ssize_t source
 
 
 cdef inline int smaller(Heapitem *a, Heapitem *b):

From a93a8698ddf48fa1e9a31d869980816178919e22 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <juan.n@unimelb.edu.au>
Date: Thu, 21 Jul 2016 18:09:18 +1000
Subject: [PATCH 08/33] Move empty mask logic to validation function

---
 skimage/morphology/watershed.py | 18 ++++++++++++------
 1 file changed, 12 insertions(+), 6 deletions(-)

diff --git a/skimage/morphology/watershed.py b/skimage/morphology/watershed.py
index 4a1db3b1e2d..72840d8ec91 100644
--- a/skimage/morphology/watershed.py
+++ b/skimage/morphology/watershed.py
@@ -33,7 +33,7 @@
 
 
 def _validate_inputs(image, markers, mask):
-    """Ensure that all inputs to watershed have matching shapes.
+    """Ensure that all inputs to watershed have matching shapes and types.
 
     Parameters
     ----------
@@ -44,6 +44,12 @@ def _validate_inputs(image, markers, mask):
     mask : array, or None
         A boolean mask, True where we want to compute the watershed.
 
+    Returns
+    -------
+    mask : array
+        The validated and formatted mask array. If ``None`` was given, it
+        is a volume of all ``True`` values.
+
     Raises
     ------
     ValueError
@@ -54,6 +60,10 @@ def _validate_inputs(image, markers, mask):
                          "as image (shape %s)" % (markers.ndim, image.ndim))
     if mask is not None and mask.shape != image.shape:
         raise ValueError("mask must have same shape as image")
+    if mask is None:
+        # Use a complete `True` mask if none is provided
+        mask = np.ones(image.shape, bool)
+    return mask
 
 
 def _validate_connectivity(image_dim, connectivity, offset):
@@ -199,14 +209,10 @@ def watershed(image, markers, connectivity=1, offset=None, mask=None):
     The algorithm works also for 3-D images, and can be used for example to
     separate overlapping spheres.
     """
-    _validate_inputs(image, markers, mask)
+    mask = _validate_inputs(image, markers, mask)
     c_connectivity, offset = _validate_connectivity(image.ndim, connectivity,
                                                     offset)
 
-    if mask is None:
-        # Use a complete `True` mask if none is provided
-        mask = np.ones(image.shape, bool)
-
     # pad the image, markers, and mask so that we can use the mask to
     # keep from running off the edges
     pad_width = [(p, p) for p in offset]

From d2a603d4626d47283ecad9334bb1c52b10198412 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <juan.n@unimelb.edu.au>
Date: Thu, 21 Jul 2016 18:10:58 +1000
Subject: [PATCH 09/33] Simplify logic computing watershed neighborhood

---
 skimage/morphology/_watershed.pyx |  4 +-
 skimage/morphology/watershed.py   | 75 ++++++++++++++++---------------
 2 files changed, 42 insertions(+), 37 deletions(-)

diff --git a/skimage/morphology/_watershed.pyx b/skimage/morphology/_watershed.pyx
index e2db0b2cd4f..7bfa28d69c4 100644
--- a/skimage/morphology/_watershed.pyx
+++ b/skimage/morphology/_watershed.pyx
@@ -26,7 +26,7 @@ include "heap_watershed.pxi"
 def watershed(cnp.float64_t[::1] image,
               DTYPE_INT32_t[:, ::1] pq,
               Py_ssize_t age,
-              DTYPE_INT32_t[:, ::1] structure,
+              DTYPE_INT32_t[::1] structure,
               DTYPE_BOOL_t[::1] mask,
               DTYPE_INT32_t[::1] output):
     """Do heavy lifting of watershed algorithm
@@ -79,7 +79,7 @@ def watershed(cnp.float64_t[::1] image,
         old_index = elem.index
         for i in range(nneighbors):
             # get the flattened address of the neighbor
-            index = structure[i, 0] + old_index
+            index = structure[i] + old_index
             if index < 0 or index >= max_index or output[index] or \
                     not mask[index]:
                 continue
diff --git a/skimage/morphology/watershed.py b/skimage/morphology/watershed.py
index 72840d8ec91..a9d13528140 100644
--- a/skimage/morphology/watershed.py
+++ b/skimage/morphology/watershed.py
@@ -111,6 +111,43 @@ def _validate_connectivity(image_dim, connectivity, offset):
     return c_connectivity, offset
 
 
+def _compute_neighbors(image, structure, offset):
+    #
+    # We pass a connectivity array that pre-calculates the stride for each
+    # neighbor.
+    #
+    # The result of this bit of code is an array with one row per
+    # point to be considered. The first column is the pre-computed stride
+    # 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(structure.shape)):
+        multiplier = 1
+        offs = []
+        indexes = []
+        ignore = True
+        for j in range(len(structure.shape)):
+            idx = (i // multiplier) % structure.shape[j]
+            off = idx - offset[j]
+            if off:
+                ignore = False
+            offs.append(off)
+            indexes.append(idx)
+            multiplier *= structure.shape[j]
+        if (not ignore) and structure.__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)
+    neighborhood = np.ascontiguousarray(c[np.argsort(distances), 0])
+    return neighborhood
+
+
 def watershed(image, markers, connectivity=1, offset=None, mask=None):
     """
     Return a matrix labeled using the watershed segmentation algorithm
@@ -225,40 +262,8 @@ def watershed(image, markers, connectivity=1, offset=None, mask=None):
     c_mask = np.ascontiguousarray(mask, dtype=bool)
     c_output = c_markers.copy()
 
-    #
-    # We pass a connectivity array that pre-calculates the stride for each
-    # neighbor.
-    #
-    # The result of this bit of code is an array with one row per
-    # point to be considered. The first column is the pre-computed stride
-    # 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
-        offs = []
-        indexes = []
-        ignore = True
-        for j in range(len(c_connectivity.shape)):
-            idx = (i // multiplier) % c_connectivity.shape[j]
-            off = idx - offset[j]
-            if off:
-                ignore = False
-            offs.append(off)
-            indexes.append(idx)
-            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)]
-    
+    flat_neighborhood = _compute_neighbors(image, c_connectivity, offset)
+
     pq, age = __heapify_markers(c_markers, c_image)
     pq = np.ascontiguousarray(pq, dtype=np.int32)
     if np.product(pq.shape) > 0:
@@ -270,7 +275,7 @@ def watershed(image, markers, connectivity=1, offset=None, mask=None):
         else:
             c_mask = c_mask.astype(np.int8).flatten()
         _watershed.watershed(c_image.flatten(),
-                             pq, age, c,
+                             pq, age, flat_neighborhood,
                              c_mask,
                              c_output)
     c_output = c_output.reshape(c_image.shape)[[slice(1, -1, None)] *

From b0b629c6be895dd6b301643daf5a3473ca28a714 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <juan.n@unimelb.edu.au>
Date: Thu, 21 Jul 2016 19:04:19 +1000
Subject: [PATCH 10/33] Simplify passing of markers to cython watershed

---
 skimage/morphology/_watershed.pyx | 15 ++++++++-------
 skimage/morphology/watershed.py   | 21 ++++++---------------
 2 files changed, 14 insertions(+), 22 deletions(-)

diff --git a/skimage/morphology/_watershed.pyx b/skimage/morphology/_watershed.pyx
index 7bfa28d69c4..49559061632 100644
--- a/skimage/morphology/_watershed.pyx
+++ b/skimage/morphology/_watershed.pyx
@@ -24,8 +24,7 @@ include "heap_watershed.pxi"
 
 @cython.boundscheck(False)
 def watershed(cnp.float64_t[::1] image,
-              DTYPE_INT32_t[:, ::1] pq,
-              Py_ssize_t age,
+              DTYPE_INT32_t[::1] marker_locations,
               DTYPE_INT32_t[::1] structure,
               DTYPE_BOOL_t[::1] mask,
               DTYPE_INT32_t[::1] output):
@@ -55,17 +54,19 @@ def watershed(cnp.float64_t[::1] image,
     cdef Heapitem new_elem
     cdef Py_ssize_t nneighbors = structure.shape[0]
     cdef Py_ssize_t i = 0
+    cdef Py_ssize_t age = 1
     cdef Py_ssize_t index = 0
     cdef Py_ssize_t old_index = 0
     cdef Py_ssize_t max_index = image.shape[0]
 
     cdef Heap *hp = <Heap *> heap_from_numpy2()
 
-    for i in range(pq.shape[0]):
-        elem.value = pq[i, 0]
-        elem.age = pq[i, 1]
-        elem.index = pq[i, 2]
-        elem.source = pq[i, 2]
+    for i in range(marker_locations.shape[0]):
+        index = marker_locations[i]
+        elem.value = image[index]
+        elem.age = 0
+        elem.index = index
+        elem.source = index
         heappush(hp, &elem)
 
     while hp.items > 0:
diff --git a/skimage/morphology/watershed.py b/skimage/morphology/watershed.py
index a9d13528140..6f8660bb946 100644
--- a/skimage/morphology/watershed.py
+++ b/skimage/morphology/watershed.py
@@ -258,24 +258,15 @@ def watershed(image, markers, connectivity=1, offset=None, mask=None):
     markers = np.pad(markers, pad_width, mode='constant')
 
     c_image = image.astype(np.float64)
-    c_markers = np.ascontiguousarray(markers, dtype=np.int32)
-    c_mask = np.ascontiguousarray(mask, dtype=bool)
-    c_output = c_markers.copy()
+    c_mask = np.ascontiguousarray(mask, dtype=np.int8).ravel()
+    c_output = np.array(markers, dtype=np.int32).ravel()
 
     flat_neighborhood = _compute_neighbors(image, c_connectivity, offset)
 
-    pq, age = __heapify_markers(c_markers, c_image)
-    pq = np.ascontiguousarray(pq, dtype=np.int32)
-    if np.product(pq.shape) > 0:
-        # If nothing is labeled, the output is empty and we don't have to
-        # do anything
-        c_output = c_output.flatten()
-        if c_mask is None:
-            c_mask = np.ones(c_image.shape, np.int8).flatten()
-        else:
-            c_mask = c_mask.astype(np.int8).flatten()
-        _watershed.watershed(c_image.flatten(),
-                             pq, age, flat_neighborhood,
+    marker_locations = np.flatnonzero(markers).astype(np.int32)
+    if len(marker_locations) > 0:
+        _watershed.watershed(c_image.ravel(),
+                             marker_locations, flat_neighborhood,
                              c_mask,
                              c_output)
     c_output = c_output.reshape(c_image.shape)[[slice(1, -1, None)] *

From e6020f747e01b2867f4b31c34a0088951dc654c7 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <juan.n@unimelb.edu.au>
Date: Thu, 21 Jul 2016 19:55:44 +1000
Subject: [PATCH 11/33] Add compactness parameter to watershed

---
 skimage/morphology/_watershed.pyx | 26 ++++++++++++++++++++++++++
 skimage/morphology/watershed.py   | 13 +++++++++++--
 2 files changed, 37 insertions(+), 2 deletions(-)

diff --git a/skimage/morphology/_watershed.pyx b/skimage/morphology/_watershed.pyx
index 49559061632..e0dcb5cc2b6 100644
--- a/skimage/morphology/_watershed.pyx
+++ b/skimage/morphology/_watershed.pyx
@@ -22,11 +22,34 @@ ctypedef cnp.int8_t DTYPE_BOOL_t
 include "heap_watershed.pxi"
 
 
+cdef extern from "math.h":
+    double sqrt(double x)
+
+
+@cython.boundscheck(False)
+@cython.cdivision(True)
+@cython.overflowcheck(False)
+@cython.unraisable_tracebacks(False)
+cdef inline double _euclid_dist(cnp.int32_t pt0, cnp.int32_t pt1,
+                                cnp.int32_t[::1] strides):
+    """Return the Euclidean distance between raveled points pt0 and pt1."""
+    cdef float result = 0
+    cdef float curr = 0
+    for i in range(strides.shape[0]):
+        curr = (pt0 // strides[i]) - (pt1 // strides[i])
+        result += curr * curr
+        pt0 = pt0 % strides[i]
+        pt1 = pt1 % strides[i]
+    return sqrt(result)
+
+
 @cython.boundscheck(False)
 def watershed(cnp.float64_t[::1] image,
               DTYPE_INT32_t[::1] marker_locations,
               DTYPE_INT32_t[::1] structure,
               DTYPE_BOOL_t[::1] mask,
+              cnp.int32_t[::1] strides,
+              cnp.float32_t compactness,
               DTYPE_INT32_t[::1] output):
     """Do heavy lifting of watershed algorithm
 
@@ -87,6 +110,9 @@ def watershed(cnp.float64_t[::1] image,
 
             age += 1
             new_elem.value = image[index]
+            if compactness > 0:
+                new_elem.value += (compactness *
+                                   _euclid_dist(index, elem.source, strides))
             new_elem.age = age
             new_elem.index = index
             new_elem.source = elem.source
diff --git a/skimage/morphology/watershed.py b/skimage/morphology/watershed.py
index 6f8660bb946..48d727b1082 100644
--- a/skimage/morphology/watershed.py
+++ b/skimage/morphology/watershed.py
@@ -148,7 +148,8 @@ def _compute_neighbors(image, structure, offset):
     return neighborhood
 
 
-def watershed(image, markers, connectivity=1, offset=None, mask=None):
+def watershed(image, markers, connectivity=1, offset=None, mask=None,
+              compactness=0):
     """
     Return a matrix labeled using the watershed segmentation algorithm
 
@@ -171,6 +172,8 @@ def watershed(image, markers, connectivity=1, offset=None, mask=None):
     mask: ndarray of bools or 0s and 1s, optional
         Array of same shape as `image`. Only points at which mask == True
         will be labeled.
+    compactness : float, optional
+        Use compact watershed [3]_ with given compactness parameter.
 
     Returns
     -------
@@ -215,6 +218,11 @@ def watershed(image, markers, connectivity=1, offset=None, mask=None):
 
     .. [2] http://cmm.ensmp.fr/~beucher/wtshed.html
 
+    .. [3] Peer Neubert & Peter Protzel (2014). Compact Watershed and
+           Preemptive SLIC: On Improving Trade-offs of Superpixel Segmentation
+           Algorithms. ICPR 2014, pp 996-1001. DOI:10.1109/ICPR.2014.181
+           https://www.tu-chemnitz.de/etit/proaut/forschung/rsrc/cws_pSLIC_ICPR.pdf
+
     Examples
     --------
     The watershed algorithm is useful to separate overlapping objects.
@@ -264,10 +272,11 @@ def watershed(image, markers, connectivity=1, offset=None, mask=None):
     flat_neighborhood = _compute_neighbors(image, c_connectivity, offset)
 
     marker_locations = np.flatnonzero(markers).astype(np.int32)
+    image_strides = np.array(image.strides, dtype=np.int32) // image.itemsize
     if len(marker_locations) > 0:
         _watershed.watershed(c_image.ravel(),
                              marker_locations, flat_neighborhood,
-                             c_mask,
+                             c_mask, image_strides, compactness,
                              c_output)
     c_output = c_output.reshape(c_image.shape)[[slice(1, -1, None)] *
                                                image.ndim]

From 8c8623cf1272b7d19b27994fe594846501a24e63 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <juan.n@unimelb.edu.au>
Date: Fri, 22 Jul 2016 15:20:11 +1000
Subject: [PATCH 12/33] Add pixels when processed, not when pushed

This ensures that the lowest-cost path to a pixel is always used.
---
 skimage/morphology/_watershed.pyx | 18 +++++-------------
 1 file changed, 5 insertions(+), 13 deletions(-)

diff --git a/skimage/morphology/_watershed.pyx b/skimage/morphology/_watershed.pyx
index e0dcb5cc2b6..418fa216fbf 100644
--- a/skimage/morphology/_watershed.pyx
+++ b/skimage/morphology/_watershed.pyx
@@ -79,7 +79,6 @@ def watershed(cnp.float64_t[::1] image,
     cdef Py_ssize_t i = 0
     cdef Py_ssize_t age = 1
     cdef Py_ssize_t index = 0
-    cdef Py_ssize_t old_index = 0
     cdef Py_ssize_t max_index = image.shape[0]
 
     cdef Heap *hp = <Heap *> heap_from_numpy2()
@@ -93,17 +92,14 @@ def watershed(cnp.float64_t[::1] image,
         heappush(hp, &elem)
 
     while hp.items > 0:
-        #
-        # Pop off an item to work on
-        #
         heappop(hp, &elem)
-        ####################################################
-        # loop through each of the structuring elements
-        #
-        old_index = elem.index
+        if output[elem.index] and elem.index != elem.source:
+            # non-marker, already visited from another neighbor
+            continue
+        output[elem.index] = output[elem.source]
         for i in range(nneighbors):
             # get the flattened address of the neighbor
-            index = structure[i] + old_index
+            index = structure[i] + elem.index
             if index < 0 or index >= max_index or output[index] or \
                     not mask[index]:
                 continue
@@ -117,9 +113,5 @@ def watershed(cnp.float64_t[::1] image,
             new_elem.index = index
             new_elem.source = elem.source
             
-            output[index] = output[old_index]
-            #
-            # Push the neighbor onto the heap to work on it later
-            #
             heappush(hp, &new_elem)
     heap_done(hp)

From a27fa83ca4109e67425587060cffed9fba7ff552 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <juan.n@unimelb.edu.au>
Date: Fri, 22 Jul 2016 15:27:59 +1000
Subject: [PATCH 13/33] Add test for compact watershed

---
 skimage/morphology/tests/test_watershed.py | 20 ++++++++++++++++++++
 1 file changed, 20 insertions(+)

diff --git a/skimage/morphology/tests/test_watershed.py b/skimage/morphology/tests/test_watershed.py
index cdaed77d98e..023c1df3c08 100644
--- a/skimage/morphology/tests/test_watershed.py
+++ b/skimage/morphology/tests/test_watershed.py
@@ -426,5 +426,25 @@ def test_watershed11(self):
         dmin = np.min(d, 2)
         self.assertTrue(np.all(d[i, j, out[i, j]-1] == dmin))
 
+
+def test_compact_watershed():
+    image = np.zeros((5, 6))
+    image[:, 3:] = 1
+    seeds = np.zeros((5, 6), dtype=int)
+    seeds[2, 0] = 1
+    seeds[2, 3] = 2
+    compact = watershed(image, seeds, compactness=0.01)
+    expected = np.array([[1, 1, 1, 2, 2, 2],
+                         [1, 1, 1, 2, 2, 2],
+                         [1, 1, 1, 2, 2, 2],
+                         [1, 1, 1, 2, 2, 2],
+                         [1, 1, 1, 2, 2, 2]], dtype=int)
+    np.testing.assert_equal(compact, expected)
+    normal = watershed(image, seeds)
+    expected = np.ones(image.shape, dtype=int)
+    expected[2, 3:] = 2
+    np.testing.assert_equal(normal, expected)
+
+
 if __name__ == "__main__":
     np.testing.run_module_suite()

From f53a5fefb82629f865985c7d7f43fde9148a6b93 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <juan.n@unimelb.edu.au>
Date: Fri, 22 Jul 2016 15:37:37 +1000
Subject: [PATCH 14/33] Remove deprecated slow watershed code

Remove test of _slow_watershed
---
 skimage/morphology/tests/test_watershed.py |   5 +-
 skimage/morphology/watershed.py            | 111 ---------------------
 2 files changed, 1 insertion(+), 115 deletions(-)

diff --git a/skimage/morphology/tests/test_watershed.py b/skimage/morphology/tests/test_watershed.py
index 023c1df3c08..d057aae701f 100644
--- a/skimage/morphology/tests/test_watershed.py
+++ b/skimage/morphology/tests/test_watershed.py
@@ -48,7 +48,7 @@
 import numpy as np
 from scipy import ndimage as ndi
 
-from skimage.morphology.watershed import watershed, _slow_watershed
+from skimage.morphology.watershed import watershed
 
 eps = 1e-12
 
@@ -114,9 +114,6 @@ def test_watershed01(self):
                       [-1, -1, -1, -1, -1, -1, -1]])
         error = diff(expected, out)
         assert error < eps
-        out = _slow_watershed(data, markers, 8)
-        error = diff(expected, out)
-        assert error < eps
 
     def test_watershed02(self):
         "watershed 2"
diff --git a/skimage/morphology/watershed.py b/skimage/morphology/watershed.py
index 48d727b1082..a2eed291ae6 100644
--- a/skimage/morphology/watershed.py
+++ b/skimage/morphology/watershed.py
@@ -284,114 +284,3 @@ def watershed(image, markers, connectivity=1, offset=None, mask=None,
         return c_output.astype(markers.dtype)
     except:
         return c_output
-
-
-# ---------------------- deprecated ------------------------------
-# Deprecate slower pure-Python code, that we keep only for
-# pedagogical purposes
-def __heapify_markers(markers, image):
-    """Create a priority queue heap with the markers on it"""
-    stride = np.array(image.strides) // image.itemsize
-    coords = np.argwhere(markers != 0)
-    ncoords = coords.shape[0]
-    if ncoords > 0:
-        pixels = image[markers != 0]
-        age = np.arange(ncoords)
-        offset = np.zeros(coords.shape[0], int)
-        for i in range(image.ndim):
-            offset = offset + stride[i] * coords[:, i]
-        pq = np.column_stack((pixels, age, offset, coords))
-        # pixels = top priority, age=second
-        ordering = np.lexsort((age, pixels))
-        pq = pq[ordering, :]
-    else:
-        pq = np.zeros((0, markers.ndim + 3), int)
-    return (pq, ncoords)
-
-
-def _slow_watershed(image, markers, connectivity=8, mask=None):
-    """Return a matrix labeled using the watershed algorithm
-
-    Use the `watershed` function for a faster execution.
-    This pure Python function is solely for pedagogical purposes.
-
-    Parameters
-    ----------
-    image: 2-d ndarray of integers
-        a two-dimensional matrix where the lowest value points are
-        labeled first.
-    markers: 2-d ndarray of integers
-        a two-dimensional matrix marking the basins with the values
-        to be assigned in the label matrix. Zero means not a marker.
-    connectivity: {4, 8}, optional
-        either 4 for four-connected or 8 (default) for eight-connected
-    mask: 2-d ndarray of bools, optional
-        don't label points in the mask
-
-    Returns
-    -------
-    out: ndarray
-        A labeled matrix of the same type and shape as markers
-
-
-    Notes
-    -----
-    This function implements a watershed algorithm [1]_that apportions pixels
-    into marked basins. The algorithm uses a priority queue to hold the pixels
-    with the metric for the priority queue being pixel value, then the time of
-    entry into the queue - this settles ties in favor of the closest marker.
-
-    Some ideas taken from
-    Soille, "Automated Basin Delineation from Digital Elevation Models Using
-    Mathematical Morphology", Signal Processing 20 (1990) 171-182
-
-    The most important insight in the paper is that entry time onto the queue
-    solves two problems: a pixel should be assigned to the neighbor with the
-    largest gradient or, if there is no gradient, pixels on a plateau should
-    be split between markers on opposite sides.
-
-    This implementation converts all arguments to specific, lowest common
-    denominator types, then passes these to a C algorithm.
-
-    Markers can be determined manually, or automatically using for example
-    the local minima of the gradient of the image, or the local maxima of the
-    distance function to the background for separating overlapping objects.
-    """
-    if connectivity not in (4, 8):
-        raise ValueError("Connectivity was %d: it should be either \
-        four or eight" % (connectivity))
-
-    image = np.array(image)
-    markers = np.array(markers)
-    labels = markers.copy()
-    max_x = markers.shape[0]
-    max_y = markers.shape[1]
-    if connectivity == 4:
-        connect_increments = ((1, 0), (0, 1), (-1, 0), (0, -1))
-    else:
-        connect_increments = ((1, 0), (1, 1), (0, 1), (-1, 1),
-                              (-1, 0), (-1, -1), (0, -1), (1, -1))
-    pq, age = __heapify_markers(markers, image)
-    pq = pq.tolist()
-    #
-    # The second step pops a value off of the queue, then labels and pushes
-    # the neighbors
-    #
-    while len(pq):
-        pix_value, pix_age, ignore, pix_x, pix_y = heappop(pq)
-        pix_label = labels[pix_x, pix_y]
-        for xi, yi in connect_increments:
-            x = pix_x + xi
-            y = pix_y + yi
-            if x < 0 or y < 0 or x >= max_x or y >= max_y:
-                continue
-            if labels[x, y]:
-                continue
-            if mask is not None and not mask[x, y]:
-                continue
-            # label the pixel
-            labels[x, y] = pix_label
-            # put the pixel onto the queue
-            heappush(pq, [image[x, y], age, 0, x, y])
-            age += 1
-    return labels

From b4d0236cca0bc82db612fae18ffe18715b85cb29 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <juan.n@unimelb.edu.au>
Date: Fri, 22 Jul 2016 16:56:24 +1000
Subject: [PATCH 15/33] Remove redundant checks from Cython watershed

---
 skimage/morphology/_watershed.pyx | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/skimage/morphology/_watershed.pyx b/skimage/morphology/_watershed.pyx
index 418fa216fbf..bf7b89037e2 100644
--- a/skimage/morphology/_watershed.pyx
+++ b/skimage/morphology/_watershed.pyx
@@ -79,7 +79,6 @@ def watershed(cnp.float64_t[::1] image,
     cdef Py_ssize_t i = 0
     cdef Py_ssize_t age = 1
     cdef Py_ssize_t index = 0
-    cdef Py_ssize_t max_index = image.shape[0]
 
     cdef Heap *hp = <Heap *> heap_from_numpy2()
 
@@ -100,8 +99,8 @@ def watershed(cnp.float64_t[::1] image,
         for i in range(nneighbors):
             # get the flattened address of the neighbor
             index = structure[i] + elem.index
-            if index < 0 or index >= max_index or output[index] or \
-                    not mask[index]:
+            if output[index] or not mask[index]:
+                # previously visited, masked, or border pixel
                 continue
 
             age += 1

From 259fdbc7c57871fcfa51497a366bbd14476e085a Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <juan.n@unimelb.edu.au>
Date: Fri, 22 Jul 2016 16:57:40 +1000
Subject: [PATCH 16/33] Remove no longer used imports in watershed.py

---
 skimage/morphology/watershed.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/skimage/morphology/watershed.py b/skimage/morphology/watershed.py
index a2eed291ae6..61d388cdd1a 100644
--- a/skimage/morphology/watershed.py
+++ b/skimage/morphology/watershed.py
@@ -24,10 +24,8 @@
 Original author: Lee Kamentsky
 """
 
-from _heapq import heappush, heappop
 import numpy as np
 from scipy import ndimage as ndi
-from ..filters import rank_order
 
 from . import _watershed
 

From b4bfb0de9c7237c9d404d5565f6d77c21bec5b48 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <juan.n@unimelb.edu.au>
Date: Fri, 22 Jul 2016 16:58:11 +1000
Subject: [PATCH 17/33] General cleanup of watershed function

---
 skimage/morphology/watershed.py | 49 ++++++++++++++++-----------------
 1 file changed, 23 insertions(+), 26 deletions(-)

diff --git a/skimage/morphology/watershed.py b/skimage/morphology/watershed.py
index 61d388cdd1a..b4e65099c5b 100644
--- a/skimage/morphology/watershed.py
+++ b/skimage/morphology/watershed.py
@@ -28,6 +28,7 @@
 from scipy import ndimage as ndi
 
 from . import _watershed
+from ..util import crop
 
 
 def _validate_inputs(image, markers, mask):
@@ -44,9 +45,10 @@ def _validate_inputs(image, markers, mask):
 
     Returns
     -------
-    mask : array
-        The validated and formatted mask array. If ``None`` was given, it
-        is a volume of all ``True`` values.
+    image, markers, mask : arrays
+        The validated and formatted arrays. Image will have dtype float64,
+        markers int32, and mask int8. If ``None`` was given for the mask,
+        it is a volume of all 1s.
 
     Raises
     ------
@@ -61,7 +63,9 @@ def _validate_inputs(image, markers, mask):
     if mask is None:
         # Use a complete `True` mask if none is provided
         mask = np.ones(image.shape, bool)
-    return mask
+    return (image.astype(np.float64),
+            markers.astype(np.int32),
+            mask.astype(np.int8))
 
 
 def _validate_connectivity(image_dim, connectivity, offset):
@@ -252,33 +256,26 @@ def watershed(image, markers, connectivity=1, offset=None, mask=None,
     The algorithm works also for 3-D images, and can be used for example to
     separate overlapping spheres.
     """
-    mask = _validate_inputs(image, markers, mask)
-    c_connectivity, offset = _validate_connectivity(image.ndim, connectivity,
-                                                    offset)
+    image, markers, mask = _validate_inputs(image, markers, mask)
+    connectivity, offset = _validate_connectivity(image.ndim, connectivity,
+                                                  offset)
 
     # pad the image, markers, and mask so that we can use the mask to
     # keep from running off the edges
     pad_width = [(p, p) for p in offset]
     image = np.pad(image, pad_width, mode='constant')
-    mask = np.pad(mask, pad_width, mode='constant')
-    markers = np.pad(markers, pad_width, mode='constant')
+    mask = np.pad(mask, pad_width, mode='constant').ravel()
+    output = np.pad(markers, pad_width, mode='constant')
 
-    c_image = image.astype(np.float64)
-    c_mask = np.ascontiguousarray(mask, dtype=np.int8).ravel()
-    c_output = np.array(markers, dtype=np.int32).ravel()
+    flat_neighborhood = _compute_neighbors(image, connectivity, offset)
+    marker_locations = np.flatnonzero(output).astype(np.int32)
+    image_strides = np.array(image.strides, dtype=np.int32) // image.itemsize
 
-    flat_neighborhood = _compute_neighbors(image, c_connectivity, offset)
+    _watershed.watershed(image.ravel(),
+                         marker_locations, flat_neighborhood,
+                         mask, image_strides, compactness,
+                         output.ravel())
 
-    marker_locations = np.flatnonzero(markers).astype(np.int32)
-    image_strides = np.array(image.strides, dtype=np.int32) // image.itemsize
-    if len(marker_locations) > 0:
-        _watershed.watershed(c_image.ravel(),
-                             marker_locations, flat_neighborhood,
-                             c_mask, image_strides, compactness,
-                             c_output)
-    c_output = c_output.reshape(c_image.shape)[[slice(1, -1, None)] *
-                                               image.ndim]
-    try:
-        return c_output.astype(markers.dtype)
-    except:
-        return c_output
+    output = crop(output, pad_width, copy=True)
+
+    return output

From df3536e3ef06c1d7460fcbb089f1fa40a1827d56 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <juan.n@unimelb.edu.au>
Date: Fri, 22 Jul 2016 17:55:11 +1000
Subject: [PATCH 18/33] Rename Cython watershed and update docstring

---
 skimage/morphology/_watershed.pyx | 56 ++++++++++++++++++-------------
 skimage/morphology/watershed.py   |  8 ++---
 2 files changed, 36 insertions(+), 28 deletions(-)

diff --git a/skimage/morphology/_watershed.pyx b/skimage/morphology/_watershed.pyx
index bf7b89037e2..8702ea44818 100644
--- a/skimage/morphology/_watershed.pyx
+++ b/skimage/morphology/_watershed.pyx
@@ -44,34 +44,42 @@ cdef inline double _euclid_dist(cnp.int32_t pt0, cnp.int32_t pt1,
 
 
 @cython.boundscheck(False)
-def watershed(cnp.float64_t[::1] image,
-              DTYPE_INT32_t[::1] marker_locations,
-              DTYPE_INT32_t[::1] structure,
-              DTYPE_BOOL_t[::1] mask,
-              cnp.int32_t[::1] strides,
-              cnp.float32_t compactness,
-              DTYPE_INT32_t[::1] output):
-    """Do heavy lifting of watershed algorithm
+def watershed_raveled(cnp.float64_t[::1] image,
+                      DTYPE_INT32_t[::1] marker_locations,
+                      DTYPE_INT32_t[::1] structure,
+                      DTYPE_BOOL_t[::1] mask,
+                      cnp.int32_t[::1] strides,
+                      cnp.float32_t compactness,
+                      DTYPE_INT32_t[::1] output):
+    """Perform watershed algorithm using a raveled image and neighborhood.
 
     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
-            the second element is the age at entry into the queue
-            the third element is the index into the flattened image or labels
-            the remaining elements are the coordinates of the point
-    age   - the next age to assign to a pixel
-    structure - a numpy int32 array containing the structuring elements
-                that define nearest neighbors. For each row, the first
-                element is the stride from the point to its neighbor
-                in a flattened array. The remaining elements are the
-                offsets from the point to its neighbor in the various
-                dimensions
-    mask  - numpy boolean (char) array indicating which pixels to consider
-            and which to ignore. Also flattened.
-    output - put the image labels in here
+    image : array of float
+        The flattened image pixels.
+    marker_locations : array of int
+        The raveled coordinates of the initial markers (aka seeds) for the
+        watershed. NOTE: these should *all* point to nonzero entries in the
+        output, or the algorithm will never terminate and blow up your memory!
+    structure : array of int
+        A list of coordinate offsets to compute the raveled coordinates of each
+        neighbor from the raveled coordinates of the current pixel.
+    mask : array of int
+        An array of the same shape as `image` where each pixel contains a
+        nonzero value if it is to be considered for flooding with watershed,
+        zero otherwise. NOTE: it is *essential* that the border pixels (those
+        with neighbors falling outside the volume) are all set to zero, or
+        segfaults could occur.
+    strides : array of int
+        An array representing the number of steps to move along each dimension.
+        This is used in computing the Euclidean distance between raveled
+        indices.
+    compactness : float
+        A value >0 implements the compact watershed algorithm (see .py file).
+    output : array of int
+        The output array, which must already contain nonzero entries at all the
+        seed locations.
     """
     cdef Heapitem elem
     cdef Heapitem new_elem
diff --git a/skimage/morphology/watershed.py b/skimage/morphology/watershed.py
index b4e65099c5b..6082fed90f2 100644
--- a/skimage/morphology/watershed.py
+++ b/skimage/morphology/watershed.py
@@ -271,10 +271,10 @@ def watershed(image, markers, connectivity=1, offset=None, mask=None,
     marker_locations = np.flatnonzero(output).astype(np.int32)
     image_strides = np.array(image.strides, dtype=np.int32) // image.itemsize
 
-    _watershed.watershed(image.ravel(),
-                         marker_locations, flat_neighborhood,
-                         mask, image_strides, compactness,
-                         output.ravel())
+    _watershed.watershed_raveled(image.ravel(),
+                                 marker_locations, flat_neighborhood,
+                                 mask, image_strides, compactness,
+                                 output.ravel())
 
     output = crop(output, pad_width, copy=True)
 

From 4816a7e7c24474e37d9dfe4ecde8079c92eb66e1 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <juan.n@unimelb.edu.au>
Date: Fri, 22 Jul 2016 18:43:53 +1000
Subject: [PATCH 19/33] Add compact watershed gallery example

---
 .../segmentation/plot_compact_watershed.py    | 42 +++++++++++++++++++
 1 file changed, 42 insertions(+)
 create mode 100644 doc/examples/segmentation/plot_compact_watershed.py

diff --git a/doc/examples/segmentation/plot_compact_watershed.py b/doc/examples/segmentation/plot_compact_watershed.py
new file mode 100644
index 00000000000..96618d6b677
--- /dev/null
+++ b/doc/examples/segmentation/plot_compact_watershed.py
@@ -0,0 +1,42 @@
+"""
+=============================================
+Find Regular Segments Using Compact Watershed
+=============================================
+
+The watershed transform is commonly used as a starting point for many
+segmentation algorithms. However, without a judicious choice of seeds, it
+can produce very uneven fragment sizes, which can be difficult to deal with
+in downstream analyses.
+
+The *compact* watershed transform remedies this by favoring seeds that are
+close to the pixel being considered.
+
+Both are implemented in the :py:func:`skimage.morphology.watershed` function.
+To use the compact form, simply pass a ``compactness`` value greater than 0.
+"""
+
+import numpy as np
+from skimage import data, util, filters, color
+from skimage.morphology import watershed
+import matplotlib.pyplot as plt
+
+coins = data.coins()
+edges = filters.sobel(coins)
+
+grid = util.regular_grid(coins.shape, n_points=468)
+
+seeds = np.zeros(coins.shape, dtype=int)
+seeds[grid] = np.arange(seeds[grid].size).reshape(seeds[grid].shape) + 1
+
+w0 = watershed(edges, seeds)
+w1 = watershed(edges, seeds, compactness=0.01)
+
+fig, (ax0, ax1) = plt.subplots(1, 2)
+
+ax0.imshow(color.label2rgb(w0, coins))
+ax0.set_title('Classical watershed')
+
+ax1.imshow(color.label2rgb(w1, coins))
+ax1.set_title('Compact watershed')
+
+plt.show()

From 278745d8f81a976f8720685942028b8e38973b1b Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <juan.n@unimelb.edu.au>
Date: Thu, 28 Jul 2016 11:26:31 +1000
Subject: [PATCH 20/33] Clarify watershed docstring about compactness

---
 skimage/morphology/watershed.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/skimage/morphology/watershed.py b/skimage/morphology/watershed.py
index 6082fed90f2..769a9abce40 100644
--- a/skimage/morphology/watershed.py
+++ b/skimage/morphology/watershed.py
@@ -152,8 +152,7 @@ def _compute_neighbors(image, structure, offset):
 
 def watershed(image, markers, connectivity=1, offset=None, mask=None,
               compactness=0):
-    """
-    Return a matrix labeled using the watershed segmentation algorithm
+    """Find watershed basins [1]_ in `image` flooded from given `markers`.
 
     Parameters
     ----------
@@ -176,6 +175,7 @@ def watershed(image, markers, connectivity=1, offset=None, mask=None,
         will be labeled.
     compactness : float, optional
         Use compact watershed [3]_ with given compactness parameter.
+        This results in more regularly-shaped watershed basins.
 
     Returns
     -------

From c41cd0af63be199f74a50d3e7ea4607ed8b92265 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <juan.n@unimelb.edu.au>
Date: Thu, 28 Jul 2016 11:33:17 +1000
Subject: [PATCH 21/33] Clarify intro to compact watershed example

---
 doc/examples/segmentation/plot_compact_watershed.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/doc/examples/segmentation/plot_compact_watershed.py b/doc/examples/segmentation/plot_compact_watershed.py
index 96618d6b677..bb472b24150 100644
--- a/doc/examples/segmentation/plot_compact_watershed.py
+++ b/doc/examples/segmentation/plot_compact_watershed.py
@@ -11,8 +11,9 @@
 The *compact* watershed transform remedies this by favoring seeds that are
 close to the pixel being considered.
 
-Both are implemented in the :py:func:`skimage.morphology.watershed` function.
-To use the compact form, simply pass a ``compactness`` value greater than 0.
+Both algorithms are implemented in the :py:func:`skimage.morphology.watershed`
+function. To use the compact form, simply pass a ``compactness`` value greater
+than 0.
 """
 
 import numpy as np

From e85e5d8792061bafb141f8bd8b5bf98a7673432d Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <juan.n@unimelb.edu.au>
Date: Thu, 28 Jul 2016 14:05:06 +1000
Subject: [PATCH 22/33] Clean up neighbor index computation

---
 skimage/morphology/watershed.py | 46 +++++++++------------------------
 1 file changed, 12 insertions(+), 34 deletions(-)

diff --git a/skimage/morphology/watershed.py b/skimage/morphology/watershed.py
index 769a9abce40..78e42d6a96f 100644
--- a/skimage/morphology/watershed.py
+++ b/skimage/morphology/watershed.py
@@ -114,40 +114,18 @@ def _validate_connectivity(image_dim, connectivity, offset):
 
 
 def _compute_neighbors(image, structure, offset):
-    #
-    # We pass a connectivity array that pre-calculates the stride for each
-    # neighbor.
-    #
-    # The result of this bit of code is an array with one row per
-    # point to be considered. The first column is the pre-computed stride
-    # 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(structure.shape)):
-        multiplier = 1
-        offs = []
-        indexes = []
-        ignore = True
-        for j in range(len(structure.shape)):
-            idx = (i // multiplier) % structure.shape[j]
-            off = idx - offset[j]
-            if off:
-                ignore = False
-            offs.append(off)
-            indexes.append(idx)
-            multiplier *= structure.shape[j]
-        if (not ignore) and structure.__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)
-    neighborhood = np.ascontiguousarray(c[np.argsort(distances), 0])
-    return neighborhood
+    """Compute neighborhood as an array of linear offsets into the image.
+
+    These are sorted according to Euclidean distance from the center (given
+    by `offset`), ensuring that immediate neighbors are visited first.
+    """
+    structure[tuple(offset)] = 0  # ignore the center; it's not a neighbor
+    locations = np.transpose(np.nonzero(structure))
+    sqdistances = np.sum((locations - offset)**2, axis=1)
+    neighborhood = (np.ravel_multi_index(locations.T, image.shape) -
+                    np.ravel_multi_index(offset, image.shape)).astype(np.int32)
+    sorted_neighborhood = neighborhood[np.argsort(sqdistances)]
+    return sorted_neighborhood
 
 
 def watershed(image, markers, connectivity=1, offset=None, mask=None,

From c4f39ec603958c32c41ced4954875f943a64c3c5 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <juan.n@unimelb.edu.au>
Date: Thu, 28 Jul 2016 14:09:59 +1000
Subject: [PATCH 23/33] Make some types double to maybe save on casting

---
 skimage/morphology/_watershed.pyx | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/skimage/morphology/_watershed.pyx b/skimage/morphology/_watershed.pyx
index 8702ea44818..2e3dc3b66e5 100644
--- a/skimage/morphology/_watershed.pyx
+++ b/skimage/morphology/_watershed.pyx
@@ -33,8 +33,8 @@ cdef extern from "math.h":
 cdef inline double _euclid_dist(cnp.int32_t pt0, cnp.int32_t pt1,
                                 cnp.int32_t[::1] strides):
     """Return the Euclidean distance between raveled points pt0 and pt1."""
-    cdef float result = 0
-    cdef float curr = 0
+    cdef double result = 0
+    cdef double curr = 0
     for i in range(strides.shape[0]):
         curr = (pt0 // strides[i]) - (pt1 // strides[i])
         result += curr * curr
@@ -49,7 +49,7 @@ def watershed_raveled(cnp.float64_t[::1] image,
                       DTYPE_INT32_t[::1] structure,
                       DTYPE_BOOL_t[::1] mask,
                       cnp.int32_t[::1] strides,
-                      cnp.float32_t compactness,
+                      cnp.double_t compactness,
                       DTYPE_INT32_t[::1] output):
     """Perform watershed algorithm using a raveled image and neighborhood.
 

From 343609ce55ae1b06f7e22f1d7b38884fcff4d06b Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <juan.n@unimelb.edu.au>
Date: Thu, 28 Jul 2016 14:17:36 +1000
Subject: [PATCH 24/33] Minor Cython cleanups

---
 skimage/morphology/_watershed.pyx | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/skimage/morphology/_watershed.pyx b/skimage/morphology/_watershed.pyx
index 2e3dc3b66e5..38de44a4cdb 100644
--- a/skimage/morphology/_watershed.pyx
+++ b/skimage/morphology/_watershed.pyx
@@ -15,15 +15,13 @@ cimport cython
 
 
 ctypedef cnp.int32_t DTYPE_INT32_t
-DTYPE_BOOL = np.bool
 ctypedef cnp.int8_t DTYPE_BOOL_t
 
 
 include "heap_watershed.pxi"
 
 
-cdef extern from "math.h":
-    double sqrt(double x)
+from libc.math cimport sqrt
 
 
 @cython.boundscheck(False)

From 2d1fabf5ee06716365c946e7750f51c705edd929 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <juan.n@unimelb.edu.au>
Date: Fri, 29 Jul 2016 14:20:17 +1000
Subject: [PATCH 25/33] Add regular_seeds function to util

---
 skimage/util/__init__.py      |  3 ++-
 skimage/util/_regular_grid.py | 33 +++++++++++++++++++++++++++++++++
 2 files changed, 35 insertions(+), 1 deletion(-)

diff --git a/skimage/util/__init__.py b/skimage/util/__init__.py
index b64836f1ca4..5a16f084eff 100644
--- a/skimage/util/__init__.py
+++ b/skimage/util/__init__.py
@@ -5,7 +5,7 @@
 from .apply_parallel import apply_parallel
 
 from .arraypad import pad, crop
-from ._regular_grid import regular_grid
+from ._regular_grid import regular_grid, regular_seeds
 from .unique import unique_rows
 from ._invert import invert
 
@@ -22,6 +22,7 @@
            'crop',
            'random_noise',
            'regular_grid',
+           'regular_seeds',
            'apply_parallel',
            'invert',
            'unique_rows']
diff --git a/skimage/util/_regular_grid.py b/skimage/util/_regular_grid.py
index f4b77a4a490..cf7acaf2ba5 100644
--- a/skimage/util/_regular_grid.py
+++ b/skimage/util/_regular_grid.py
@@ -70,3 +70,36 @@ def regular_grid(ar_shape, n_points):
               start, step in zip(starts, stepsizes)]
     slices = [slices[i] for i in unsort_dim_idxs]
     return slices
+
+
+def regular_seeds(ar_shape, n_points, dtype=int):
+    """Return an image with ~`n_points` regularly-spaced nonzero pixels.
+
+    Parameters
+    ----------
+    ar_shape : tuple of int
+        The shape of the desired output image.
+    n_points : int
+        The desired number of nonzero points.
+    dtype : numpy data type, optional
+        The desired data type of the output.
+
+    Returns
+    -------
+    seed_img : array of int or bool
+        The desired image.
+
+    Examples
+    --------
+    >>> regular_seeds((5, 5), 4)
+    array([[0, 0, 0, 0, 0],
+           [0, 1, 0, 2, 0],
+           [0, 0, 0, 0, 0],
+           [0, 3, 0, 4, 0],
+           [0, 0, 0, 0, 0]])
+    """
+    grid = regular_grid(ar_shape, n_points)
+    seed_img = np.zeros(ar_shape, dtype=dtype)
+    seed_img[grid] = 1 + np.reshape(np.arange(seed_img[grid].size),
+                                    seed_img[grid].shape)
+    return seed_img

From 5c490ecdbf7acff8a64de43efa5c5d94f7f59b14 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <juan.n@unimelb.edu.au>
Date: Fri, 29 Jul 2016 14:45:35 +1000
Subject: [PATCH 26/33] Allow specifying watershed seeds with just int

This makes compact watershed substantially simpler to use.
---
 skimage/morphology/tests/test_watershed.py | 13 +++++++++++++
 skimage/morphology/watershed.py            | 16 +++++++++-------
 2 files changed, 22 insertions(+), 7 deletions(-)

diff --git a/skimage/morphology/tests/test_watershed.py b/skimage/morphology/tests/test_watershed.py
index d057aae701f..d0c7ca8d7fb 100644
--- a/skimage/morphology/tests/test_watershed.py
+++ b/skimage/morphology/tests/test_watershed.py
@@ -443,5 +443,18 @@ def test_compact_watershed():
     np.testing.assert_equal(normal, expected)
 
 
+def test_numeric_seed_watershed():
+    """Test that passing just the number of seeds to watershed works."""
+    image = np.zeros((5, 6))
+    image[:, 3:] = 1
+    compact = watershed(image, 2, compactness=0.01)
+    expected = np.array([[1, 1, 1, 1, 2, 2],
+                         [1, 1, 1, 1, 2, 2],
+                         [1, 1, 1, 1, 2, 2],
+                         [1, 1, 1, 1, 2, 2],
+                         [1, 1, 1, 1, 2, 2]], dtype=np.int32)
+    np.testing.assert_equal(compact, expected)
+
+
 if __name__ == "__main__":
     np.testing.run_module_suite()
diff --git a/skimage/morphology/watershed.py b/skimage/morphology/watershed.py
index 78e42d6a96f..3fb7df09a25 100644
--- a/skimage/morphology/watershed.py
+++ b/skimage/morphology/watershed.py
@@ -28,7 +28,7 @@
 from scipy import ndimage as ndi
 
 from . import _watershed
-from ..util import crop
+from ..util import crop, regular_seeds
 
 
 def _validate_inputs(image, markers, mask):
@@ -38,7 +38,7 @@ def _validate_inputs(image, markers, mask):
     ----------
     image : array
         The input image.
-    markers : array
+    markers : int or array of int
         The marker image.
     mask : array, or None
         A boolean mask, True where we want to compute the watershed.
@@ -55,7 +55,10 @@ def _validate_inputs(image, markers, mask):
     ValueError
         If the shapes of the given arrays don't match.
     """
-    if markers.shape != image.shape:
+    if not isinstance(markers, (np.ndarray, list, tuple)):
+        # not array-like, assume int
+        markers = regular_seeds(image.shape, markers)
+    elif markers.shape != image.shape:
         raise ValueError("Markers (shape %s) must have same shape "
                          "as image (shape %s)" % (markers.ndim, image.ndim))
     if mask is not None and mask.shape != image.shape:
@@ -137,10 +140,9 @@ def watershed(image, markers, connectivity=1, offset=None, mask=None,
 
     image: ndarray (2-D, 3-D, ...) of integers
         Data array where the lowest value points are labeled first.
-    markers: ndarray of the same shape as `image`
-        An array marking the basins with the values to be assigned in the
-        label matrix. Zero means not a marker. This array should be of an
-        integer type.
+    markers: int, or ndarray of int, same shape as `image`
+        The desired number of markers, or an array marking the basins with the
+        values to be assigned in the label matrix. Zero means not a marker.
     connectivity: ndarray, optional
         An array with the same number of dimensions as `image` whose
         non-zero elements indicate neighbors for connection.

From 74b00f4538a6d44f8e8fcb6abe9081733273e16c Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <juan.n@unimelb.edu.au>
Date: Fri, 29 Jul 2016 14:47:13 +1000
Subject: [PATCH 27/33] Alias watershed within the segmentation module

---
 skimage/segmentation/__init__.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/skimage/segmentation/__init__.py b/skimage/segmentation/__init__.py
index 63f7d66d8ba..dfcecb28cd1 100644
--- a/skimage/segmentation/__init__.py
+++ b/skimage/segmentation/__init__.py
@@ -6,6 +6,7 @@
 from .boundaries import find_boundaries, mark_boundaries
 from ._clear_border import clear_border
 from ._join import join_segmentations, relabel_from_one, relabel_sequential
+from ..morphology import watershed
 
 
 __all__ = ['random_walker',
@@ -18,4 +19,5 @@
            'clear_border',
            'join_segmentations',
            'relabel_from_one',
-           'relabel_sequential']
+           'relabel_sequential',
+           'watershed']

From 4f467271d7c72aa515bee2769aa54e54daefd560 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <juan.n@unimelb.edu.au>
Date: Fri, 29 Jul 2016 14:47:31 +1000
Subject: [PATCH 28/33] Add watershed to plot_segmentations

---
 .../segmentation/plot_segmentations.py        | 51 +++++++++++++++----
 1 file changed, 40 insertions(+), 11 deletions(-)

diff --git a/doc/examples/segmentation/plot_segmentations.py b/doc/examples/segmentation/plot_segmentations.py
index 33008069669..e06538ceb7d 100644
--- a/doc/examples/segmentation/plot_segmentations.py
+++ b/doc/examples/segmentation/plot_segmentations.py
@@ -59,6 +59,29 @@
 .. [3] Radhakrishna Achanta, Appu Shaji, Kevin Smith, Aurelien Lucchi,
     Pascal Fua, and Sabine Suesstrunk, SLIC Superpixels Compared to
     State-of-the-art Superpixel Methods, TPAMI, May 2012.
+
+
+Compact watershed segmentation of gradient images
+-------------------------------------------------
+
+Instead of taking a color image as input, watershed requires a grayscale
+*gradient* image, where bright pixels denote a boundary between regions.
+The algorithm views the image as a landscape, with bright pixels forming high
+peaks. This landscape is then flooded from the given *markers*, until separate
+flood basins meet at the peaks. Each distinct basin then forms a different
+image segment. [4]_
+
+As with SLIC, there is an additional *compactness* argument that makes it
+harder for markers to flood faraway pixels. This makes the watershed regions
+more regularly shaped. [5]_
+
+
+.. [4] http://en.wikipedia.org/wiki/Watershed_%28image_processing%29
+
+.. [5] Peer Neubert & Peter Protzel (2014). Compact Watershed and
+       Preemptive SLIC: On Improving Trade-offs of Superpixel Segmentation
+       Algorithms. ICPR 2014, pp 996-1001. DOI:10.1109/ICPR.2014.181
+       https://www.tu-chemnitz.de/etit/proaut/forschung/rsrc/cws_pSLIC_ICPR.pdf
 """
 from __future__ import print_function
 
@@ -66,7 +89,9 @@
 import numpy as np
 
 from skimage.data import astronaut
-from skimage.segmentation import felzenszwalb, slic, quickshift
+from skimage.color import rgb2gray
+from skimage.filters import sobel
+from skimage.segmentation import felzenszwalb, slic, quickshift, watershed
 from skimage.segmentation import mark_boundaries
 from skimage.util import img_as_float
 
@@ -74,23 +99,27 @@
 segments_fz = felzenszwalb(img, scale=100, sigma=0.5, min_size=50)
 segments_slic = slic(img, n_segments=250, compactness=10, sigma=1)
 segments_quick = quickshift(img, kernel_size=3, max_dist=6, ratio=0.5)
+gradient = sobel(rgb2gray(img))
+segments_watershed = watershed(gradient, markers=250, compactness=0.001)
 
 print("Felzenszwalb's number of segments: %d" % len(np.unique(segments_fz)))
-print("Slic number of segments: %d" % len(np.unique(segments_slic)))
-print("Quickshift number of segments: %d" % len(np.unique(segments_quick)))
+print('SLIC number of segments: %d' % len(np.unique(segments_slic)))
+print('Quickshift number of segments: %d' % len(np.unique(segments_quick)))
 
-fig, ax = plt.subplots(1, 3, sharex=True, sharey=True,
+fig, ax = plt.subplots(2, 2, sharex=True, sharey=True,
                        subplot_kw={'adjustable': 'box-forced'})
 fig.set_size_inches(8, 3, forward=True)
 fig.tight_layout()
 
-ax[0].imshow(mark_boundaries(img, segments_fz))
-ax[0].set_title("Felzenszwalbs's method")
-ax[1].imshow(mark_boundaries(img, segments_slic))
-ax[1].set_title("SLIC")
-ax[2].imshow(mark_boundaries(img, segments_quick))
-ax[2].set_title("Quickshift")
-for a in ax:
+ax[0, 0].imshow(mark_boundaries(img, segments_fz))
+ax[0, 0].set_title("Felzenszwalbs's method")
+ax[0, 1].imshow(mark_boundaries(img, segments_slic))
+ax[0, 1].set_title('SLIC')
+ax[1, 0].imshow(mark_boundaries(img, segments_quick))
+ax[1, 0].set_title('Quickshift')
+ax[1, 1].imshow(mark_boundaries(img, segments_watershed))
+ax[1, 1].set_title('Compact watershed')
+for a in ax.ravel():
     a.set_xticks(())
     a.set_yticks(())
 plt.show()

From 1f011d340ee97c8306ba93100f268599dc77d97d Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <juan.n@unimelb.edu.au>
Date: Sat, 30 Jul 2016 12:47:16 +1000
Subject: [PATCH 29/33] Clarify compactness arg description in docstring

---
 skimage/morphology/watershed.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/skimage/morphology/watershed.py b/skimage/morphology/watershed.py
index 3fb7df09a25..568bf553183 100644
--- a/skimage/morphology/watershed.py
+++ b/skimage/morphology/watershed.py
@@ -155,7 +155,7 @@ def watershed(image, markers, connectivity=1, offset=None, mask=None,
         will be labeled.
     compactness : float, optional
         Use compact watershed [3]_ with given compactness parameter.
-        This results in more regularly-shaped watershed basins.
+        Higher values result in more regularly-shaped watershed basins.
 
     Returns
     -------

From 18552c5001fa405a377f3532828831ca662eeefb Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <juan.n@unimelb.edu.au>
Date: Mon, 1 Aug 2016 14:35:39 +1000
Subject: [PATCH 30/33] Fix incorrect in-text references

---
 doc/examples/segmentation/plot_segmentations.py | 2 +-
 skimage/morphology/watershed.py                 | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/doc/examples/segmentation/plot_segmentations.py b/doc/examples/segmentation/plot_segmentations.py
index e06538ceb7d..0847a9fb29f 100644
--- a/doc/examples/segmentation/plot_segmentations.py
+++ b/doc/examples/segmentation/plot_segmentations.py
@@ -51,7 +51,7 @@
 image location and is therefore closely related to quickshift. As the
 clustering method is simpler, it is very efficient. It is essential for this
 algorithm to work in Lab color space to obtain good results.  The algorithm
-quickly gained momentum and is now widely used. See [3] for details.  The
+quickly gained momentum and is now widely used. See [3]_ for details.  The
 ``compactness`` parameter trades off color-similarity and proximity, as in the
 case of Quickshift, while ``n_segments`` chooses the number of centers for
 kmeans.
diff --git a/skimage/morphology/watershed.py b/skimage/morphology/watershed.py
index 568bf553183..639182f9484 100644
--- a/skimage/morphology/watershed.py
+++ b/skimage/morphology/watershed.py
@@ -172,7 +172,7 @@ def watershed(image, markers, connectivity=1, offset=None, mask=None,
 
     Notes
     -----
-    This function implements a watershed algorithm [1]_that apportions pixels
+    This function implements a watershed algorithm [1]_ that apportions pixels
     into marked basins. The algorithm uses a priority queue to hold the pixels
     with the metric for the priority queue being pixel value, then the time of
     entry into the queue - this settles ties in favor of the closest marker.

From 47af70a5909fbec11edb0aae26db9f8bdf6bd750 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <juan.n@unimelb.edu.au>
Date: Tue, 2 Aug 2016 09:49:23 +1000
Subject: [PATCH 31/33] Tweak text spacing

---
 doc/examples/segmentation/plot_segmentations.py | 1 -
 skimage/morphology/_watershed.pyx               | 3 ++-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/doc/examples/segmentation/plot_segmentations.py b/doc/examples/segmentation/plot_segmentations.py
index 0847a9fb29f..e179b30a474 100644
--- a/doc/examples/segmentation/plot_segmentations.py
+++ b/doc/examples/segmentation/plot_segmentations.py
@@ -75,7 +75,6 @@
 harder for markers to flood faraway pixels. This makes the watershed regions
 more regularly shaped. [5]_
 
-
 .. [4] http://en.wikipedia.org/wiki/Watershed_%28image_processing%29
 
 .. [5] Peer Neubert & Peter Protzel (2014). Compact Watershed and
diff --git a/skimage/morphology/_watershed.pyx b/skimage/morphology/_watershed.pyx
index 38de44a4cdb..cac8894f9ae 100644
--- a/skimage/morphology/_watershed.pyx
+++ b/skimage/morphology/_watershed.pyx
@@ -74,7 +74,8 @@ def watershed_raveled(cnp.float64_t[::1] image,
         This is used in computing the Euclidean distance between raveled
         indices.
     compactness : float
-        A value >0 implements the compact watershed algorithm (see .py file).
+        A value greater than 0 implements the compact watershed algorithm
+        (see .py file).
     output : array of int
         The output array, which must already contain nonzero entries at all the
         seed locations.

From 5b0ba4b16bf0a73912650d29befef00246678f2f Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <juan.n@unimelb.edu.au>
Date: Wed, 3 Aug 2016 17:01:06 +1000
Subject: [PATCH 32/33] Add reference to ref 2 in text

---
 skimage/morphology/watershed.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/skimage/morphology/watershed.py b/skimage/morphology/watershed.py
index 639182f9484..7e12c77ff24 100644
--- a/skimage/morphology/watershed.py
+++ b/skimage/morphology/watershed.py
@@ -133,7 +133,7 @@ def _compute_neighbors(image, structure, offset):
 
 def watershed(image, markers, connectivity=1, offset=None, mask=None,
               compactness=0):
-    """Find watershed basins [1]_ in `image` flooded from given `markers`.
+    """Find watershed basins [1]_ [2]_ in `image` flooded from given `markers`.
 
     Parameters
     ----------

From d34f7ab51bf6f39f0c67c6ad65ccafdba2d33754 Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <juan.n@unimelb.edu.au>
Date: Thu, 25 Aug 2016 12:27:48 +1000
Subject: [PATCH 33/33] Reorganise refs in watershed docstring

---
 skimage/morphology/watershed.py | 11 ++++++-----
 1 file changed, 6 insertions(+), 5 deletions(-)

diff --git a/skimage/morphology/watershed.py b/skimage/morphology/watershed.py
index 7e12c77ff24..6e6120db471 100644
--- a/skimage/morphology/watershed.py
+++ b/skimage/morphology/watershed.py
@@ -133,7 +133,7 @@ def _compute_neighbors(image, structure, offset):
 
 def watershed(image, markers, connectivity=1, offset=None, mask=None,
               compactness=0):
-    """Find watershed basins [1]_ [2]_ in `image` flooded from given `markers`.
+    """Find watershed basins in `image` flooded from given `markers`.
 
     Parameters
     ----------
@@ -172,10 +172,11 @@ def watershed(image, markers, connectivity=1, offset=None, mask=None,
 
     Notes
     -----
-    This function implements a watershed algorithm [1]_ that apportions pixels
-    into marked basins. The algorithm uses a priority queue to hold the pixels
-    with the metric for the priority queue being pixel value, then the time of
-    entry into the queue - this settles ties in favor of the closest marker.
+    This function implements a watershed algorithm [1]_ [2]_ that apportions
+    pixels into marked basins. The algorithm uses a priority queue to hold
+    the pixels with the metric for the priority queue being pixel value, then
+    the time of entry into the queue - this settles ties in favor of the
+    closest marker.
 
     Some ideas taken from
     Soille, "Automated Basin Delineation from Digital Elevation Models Using