The watershed algorithm (see http://en.wikipedia.org/wiki/Watershed_(algorithm)) is used to split an image into distinct components.

Suppose that we have the following image, composed of three whites disks (pixels of value 1) and a black background (pixels of value 0). We want to obtain a new array where each pixel is labeled with the index of the component to which it belongs, that is a segmentation of the orginal array, as shown in the image below. We will use the watershed algorithm provided by scipy.ndimage, scipy.ndimage.watershed_ift.

```
1 >>> # Create the initial black and white image
2 >>> import numpy as np
3 >>> from scipy import ndimage
4 >>> a = np.zeros((512, 512)).astype(uint8) #unsigned integer type needed by watershed
5 >>> y, x = np.ogrid[0:512, 0:512]
6 >>> m1 = ((y-200)**2 + (x-100)**2 < 30**2)
7 >>> m2 = ((y-350)**2 + (x-400)**2 < 20**2)
8 >>> m3 = ((y-260)**2 + (x-200)**2 < 20**2)
9 >>> a[m1+m2+m3]=1
10 >>> imshow(a, cmap = cm.gray)# left plot in the image above
```

The watershed algorithm relies on the *flooding* of different basins, so we need to put markers in the image to initiate the flooding. If one knows approximately where the objects are, and there are only a few objects, it is possible to set the markers by hand

```
1 >>> markers = np.zeros_like(a).astype(int16)
2 >>> markers[0, 0] = 1
3 >>> markers[200, 100] = 2
4 >>> markers[350, 400] = 3
5 >>> markers[260, 200] = 4
6 >>> res1 = ndimage.watershed_ift(a.astype(uint8), markers)
7 >>> np.unique1d(res1) # pixels are tagged according to the object they belong to
8 array([1, 2, 3, 4], dtype=int16)
9 >>> imshow(res1, cmap=cm.jet) # central plot in the image above
```

If you don't know where to put the markers, and you know the minimal size of the objects, you can spread a lot of markers on a grid finer than the objects.

```
1 >>> xm, ym = np.ogrid[0:512:10, 0:512:10]
2 >>> markers = np.zeros_like(a).astype(int16)
3 >>> markers[xm, ym]= np.arange(xm.size*ym.size).reshape((xm.size,ym.size))
4 >>> res2 = ndimage.watershed_ift(a.astype(uint8), markers)
5 >>> res2[xm, ym] = res2[xm-1, ym-1] # remove the isolate seeds
6 >>> imshow(res2, cmap=cm.jet)
7 <matplotlib.image.AxesImage object at 0xf1fd1ac>
```