Voronoi-Otsu-labeling

This workflow for image segmentation is a rather simple and yet powerful approach, e.g. for detecting and segmenting nuclei in fluorescence micropscopy images. A nuclei marker such as nuclei-GFP, DAPI or histone-RFP in combination with various microscopy techniques can be used to generate images of suitable kind.

To demonstrate the workflow, we're using image data from the Broad Bio Image Challenge: We used image set BBBC022v1 Gustafsdottir et al., PLOS ONE, 2013, available from the Broad Bioimage Benchmark Collection Ljosa et al., Nature Methods, 2012.

We start by opening an example image.

In [1]:
from skimage.io import imread, imshow
import matplotlib.pyplot as plt

input_image = imread("BBBC022/IXMtest_A02_s9.tif")[:,:,0]

# define an interesting sub-region
bb_x=200
bb_y=0
bb_width=200
bb_height=200

Next, we initialize the GPU

In [2]:
import pyclesperanto_prototype as cle

# select a specific OpenCL / GPU device and see which one was chosen
cle.select_device('RTX')
Out[2]:
<NVIDIA GeForce RTX 3050 Ti Laptop GPU on Platform: NVIDIA CUDA (1 refs)>

Before segmenting the image, need to push it to GPU memory. For visualisation purposes we crop out a sub-region:

In [3]:
input_gpu = cle.push(input_image)

input_crop = cle.crop(input_gpu, start_x=bb_x, start_y=bb_y, width=bb_width, height=bb_height)

fig, axs = plt.subplots(1, 2, figsize=(15, 15))
cle.imshow(input_gpu, plot=axs[0])
cle.imshow(input_crop, plot=axs[1])

Applying the algorithm

Voronoi-Otsu-labeling is a command in clesperanto, which asks for two sigma parameters. The first sigma controls how close detected cells can be (spot_sigma) and second controls how precise segmented objects are outlined (outline_sigma).

In [4]:
sigma_spot_detection = 5
sigma_outline = 1

segmented = cle.voronoi_otsu_labeling(input_gpu, spot_sigma=sigma_spot_detection, outline_sigma=sigma_outline)
segmented_crop = cle.crop(segmented, start_x=bb_x, start_y=bb_y, width=bb_width, height=bb_height)

fig, axs = plt.subplots(1, 2, figsize=(15, 15))
cle.imshow(segmented, labels=True, plot=axs[0])
cle.imshow(segmented_crop, labels=True, plot=axs[1])

How does it work?

The Voronoi-Otsu-Labeling workflow is a combination of Gaussian blur, spot detection, thresholding and binary watershed. The interested reader might want to see the open source code. The approach is similar to applying a seeded watershed to a binary image, e.g. in MorphoLibJ or scikit-image. However, the seeds are computed automatically and cannot be passed.

For demonstration purposes we do that only on the 2D cropped image shown above. If this algorithm is applied to 3D data, it is recommended to make it isotropic first.

In [5]:
image_to_segment = input_crop
print(image_to_segment.shape)
(200, 200)

As a first step, we blur the image with a given sigma and detect maxima in the resulting image.

In [6]:
blurred = cle.gaussian_blur(image_to_segment, sigma_x=sigma_spot_detection, sigma_y=sigma_spot_detection, sigma_z=sigma_spot_detection)

detected_spots = cle.detect_maxima_box(blurred, radius_x=0, radius_y=0, radius_z=0)

number_of_spots = cle.sum_of_all_pixels(detected_spots)
print("number of detected spots", number_of_spots)

fig, axs = plt.subplots(1, 2, figsize=(15, 15))
cle.imshow(blurred, plot=axs[0])
cle.imshow(detected_spots, plot=axs[1])
number of detected spots 29.0

Furthermore, we start again from the cropped image and blur it again, with a different sigma. Afterwards, we threshold the image using Otsu's thresholding method (Otsu et al 1979).

In [7]:
blurred = cle.gaussian_blur(image_to_segment, sigma_x=sigma_outline, sigma_y=sigma_outline, sigma_z=sigma_outline)

binary = cle.threshold_otsu(blurred)

fig, axs = plt.subplots(1, 2, figsize=(15, 15))
cle.imshow(blurred, plot=axs[0])
cle.imshow(binary, plot=axs[1])

Afterwards, we take the binary spots image and the binary segmentation image and apply a binary_and operation to exclude spots which were detected in the background area. Those likely corresponded to noise.

In [8]:
selected_spots = cle.binary_and(binary, detected_spots)

number_of_spots = cle.sum_of_all_pixels(selected_spots)
print("number of selected spots", number_of_spots)

fig, axs = plt.subplots(1, 3, figsize=(15, 15))
cle.imshow(detected_spots, plot=axs[0])
cle.imshow(binary, plot=axs[1])
cle.imshow(selected_spots, plot=axs[2])
number of selected spots 11.0

Next, we separate the image space between the selected spots using a Voronoi diagram which is limited to the positive pixels in the binary image.

In [9]:
voronoi_diagram = cle.masked_voronoi_labeling(selected_spots, binary)

fig, axs = plt.subplots(1, 3, figsize=(15, 15))
cle.imshow(selected_spots, plot=axs[0])
cle.imshow(binary, plot=axs[1])
cle.imshow(voronoi_diagram, labels=True, plot=axs[2])
In [ ]: