SeededWatershed and 'uphill only'

[acycliq]

Component
diplib 3.3.0

Describe the bug
I dont know is this a bug but I am looking into watershed and I might be missing something. I am experimenting with a toy example, two circles, one white the other gray on a black background. When I am using the uphill only flag then my segmentation is just two crosses. Is this something to be expected? Does the algorithm walk along the vertical and horizontal dimension first, then it hits the edges of the basin and it stays there because all other pixels under the mask have lower values?

See code below for a reproducible example

import skimage.io
from scipy import ndimage
import diplib as dip
import numpy as np
from skimage.draw import disk
import skimage


img = np.zeros(shape=[512, 512], dtype=np.uint8)
rr, cc = disk((200, 200), 100, shape=img.shape)
img[rr, cc] = 255

rr, cc = disk((340, 340), 100, shape=img.shape)
img[rr, cc] = 127

img = img.astype(np.uint8)
skimage.io.imshow(img, cmap="gray")
skimage.io.show()

distance = ndimage.distance_transform_edt(img)

seeds = np.zeros(img.shape)
seeds[200, 200] = 1
seeds[340, 340] = 2
seeds = seeds.astype(np.uint8)

labels = dip.SeededWatershed(-distance.astype(np.uint8), dip.Maxima(seeds), np.array(img)>0,
                         maxDepth=86, # 86: doesnt merge, 87: merges. Need to remove 'uphill only'
                         flags={'labels', 'uphill only'}
                         )
labels = 255 * np.array(labels)/np.array(labels).max()
skimage.io.imshow(labels.astype(np.uint8), cmap="gray")
skimage.io.show()

[crisluengo]

The "uphill only" flag literally does that: it only propagates to higher-valued pixels. If a neighbor has the same value, it won't be propagated into.

Because you discretize the distance transform (-distance.astype(np.uint8)), you create many neighbors that have the same value value, and hence the propagation is blocked.

If you remove the .astype(np.uint8) conversion, it works as expected. You still get an ugly unlabeled square where the two circles meet, because the distance transform there is not strictly increasing.

labels = dip.SeededWatershed(-distance, seeds, img>0, maxDepth=86, flags={'labels', 'uphill only'})

---

Note that in DIPlib, most functionality is implemented for images of arbitrary types, and floating-point types tend to be preferred. There is hardly ever a need to explicitly cast and image to an integer type. Of course there are functions that only work on binary images, for example, and labelled images must always be an unsigned integer type.

---

Here's your example code using only DIPlib, in the hopes of helping you find functionality in the library:

import diplib as dip

img = dip.Image((512, 512), dt="UINT8")
img.Fill(0)
dip.DrawEllipsoid(img, (200, 200), (200, 200), 255)
dip.DrawEllipsoid(img, (200, 200), (340, 340), 127)
img.Show()

distance = dip.EuclideanDistanceTransform(img > 0)

seeds = img.Similar()
seeds.Fill(0)
seeds[200, 200] = 1
seeds[340, 340] = 2

labels = dip.SeededWatershed(distance, seeds, img > 0,
                         maxDepth=86, # 86: doesnt merge, 87: merges. Need to remove 'uphill only'
                         flags={'labels', 'uphill only', 'high first'}
                         )
labels.Show('labels')

[acycliq]

Thanks a lot Cris. Yes that makes sense. I had also noticed that too; When the int conversion is removed it works better, I was a bit puzzled about that square in the middle though.
Thanks for the example code and for the library!