Skip to content

Commit

Permalink
Include border when generating basin markers for watershed. (#7362)
Browse files Browse the repository at this point in the history
Updated the Segment human cells in mitosis (link to stable):

- removed the very last figure, which doesn't add much value;
- added a few words about the watershed segmentation step;
- fixed the way markers are computed, so that basins touching the image border are included. You can see the difference especially if you look at the RHS vertical border of the image.

Co-authored-by: Lars Grüter <lagru+github@mailbox.org>
  • Loading branch information
mkcor and lagru committed Apr 4, 2024
1 parent 8cc2c24 commit eab8481
Showing 1 changed file with 32 additions and 26 deletions.
58 changes: 32 additions & 26 deletions doc/examples/applications/plot_human_mitosis.py
Expand Up @@ -21,10 +21,9 @@
import numpy as np
from scipy import ndimage as ndi

from skimage import color, feature, filters, measure, morphology, segmentation, util
from skimage.data import human_mitosis
import skimage as ski

image = human_mitosis()
image = ski.data.human_mitosis()

fig, ax = plt.subplots()
ax.imshow(image, cmap='gray')
Expand Down Expand Up @@ -70,7 +69,7 @@
# To separate these three different classes of pixels, we
# resort to :ref:`sphx_glr_auto_examples_segmentation_plot_multiotsu.py`.

thresholds = filters.threshold_multiotsu(image, classes=3)
thresholds = ski.filters.threshold_multiotsu(image, classes=3)
regions = np.digitize(image, bins=thresholds)

fig, ax = plt.subplots(ncols=2, figsize=(10, 5))
Expand All @@ -90,8 +89,8 @@

cells = image > thresholds[0]
dividing = image > thresholds[1]
labeled_cells = measure.label(cells)
labeled_dividing = measure.label(dividing)
labeled_cells = ski.measure.label(cells)
labeled_dividing = ski.measure.label(dividing)
naive_mi = labeled_dividing.max() / labeled_cells.max()
print(naive_mi)

Expand All @@ -112,12 +111,12 @@
ax[0].imshow(image)
ax[0].set_title('Original')
ax[0].set_axis_off()
ax[2].imshow(cells)
ax[2].set_title('All nuclei?')
ax[2].set_axis_off()
ax[1].imshow(dividing)
ax[1].set_title('Dividing nuclei?')
ax[1].set_axis_off()
ax[2].imshow(cells)
ax[2].set_title('All nuclei?')
ax[2].set_axis_off()
plt.show()

#####################################################################
Expand All @@ -140,7 +139,9 @@
higher_threshold = 125
dividing = image > higher_threshold

smoother_dividing = filters.rank.mean(util.img_as_ubyte(dividing), morphology.disk(4))
smoother_dividing = ski.filters.rank.mean(
ski.util.img_as_ubyte(dividing), ski.morphology.disk(4)
)

binary_smoother_dividing = smoother_dividing > 20

Expand All @@ -152,7 +153,7 @@

#####################################################################
# We are left with
cleaned_dividing = measure.label(binary_smoother_dividing)
cleaned_dividing = ski.measure.label(binary_smoother_dividing)
print(cleaned_dividing.max())

#####################################################################
Expand All @@ -163,38 +164,43 @@
# ==============
# To separate overlapping nuclei, we resort to
# :ref:`sphx_glr_auto_examples_segmentation_plot_watershed.py`.
# To visualize the segmentation conveniently, we colour-code the labelled
# regions using the `color.label2rgb` function, specifying the background
# label with argument `bg_label=0`.
# The idea of the algorithm is to find watershed basins as if flooded from
# a set of `markers`. We generate these markers as the local maxima of the
# distance function to the background. Given the typical size of the nuclei,
# we pass ``min_distance=7`` so that local maxima and, hence, markers will lie
# at least 7 pixels away from one another.
# We also use ``exclude_border=False`` so that all nuclei touching the image
# border will be included.

distance = ndi.distance_transform_edt(cells)

local_max_coords = feature.peak_local_max(distance, min_distance=7)
local_max_coords = ski.feature.peak_local_max(
distance, min_distance=7, exclude_border=False
)
local_max_mask = np.zeros(distance.shape, dtype=bool)
local_max_mask[tuple(local_max_coords.T)] = True
markers = measure.label(local_max_mask)
markers = ski.measure.label(local_max_mask)

segmented_cells = segmentation.watershed(-distance, markers, mask=cells)
segmented_cells = ski.segmentation.watershed(-distance, markers, mask=cells)

#####################################################################
# To visualize the segmentation conveniently, we colour-code the labelled
# regions using the `color.label2rgb` function, specifying the background
# label with argument `bg_label=0`.

fig, ax = plt.subplots(ncols=2, figsize=(10, 5))
ax[0].imshow(cells, cmap='gray')
ax[0].set_title('Overlapping nuclei')
ax[0].set_axis_off()
ax[1].imshow(color.label2rgb(segmented_cells, bg_label=0))
ax[1].imshow(ski.color.label2rgb(segmented_cells, bg_label=0))
ax[1].set_title('Segmented nuclei')
ax[1].set_axis_off()
plt.show()

#####################################################################
# Additionally, we may use function `color.label2rgb` to overlay the original
# image with the segmentation result, using transparency (alpha parameter).

color_labels = color.label2rgb(segmented_cells, image, alpha=0.4, bg_label=0)
# Make sure that the watershed algorithm has led to identifying more nuclei:

fig, ax = plt.subplots(figsize=(5, 5))
ax.imshow(color_labels)
ax.set_title('Segmentation result over raw image')
plt.show()
assert segmented_cells.max() > labeled_cells.max()

#####################################################################
# Finally, we find a total number of
Expand Down

0 comments on commit eab8481

Please sign in to comment.