diff --git a/doc/examples/applications/plot_human_mitosis.py b/doc/examples/applications/plot_human_mitosis.py
index 377180cbecc..e0abc2e1252 100644
--- a/doc/examples/applications/plot_human_mitosis.py
+++ b/doc/examples/applications/plot_human_mitosis.py
@@ -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')
@@ -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))
@@ -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)
 
@@ -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()
 
 #####################################################################
@@ -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
 
@@ -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())
 
 #####################################################################
@@ -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