Segmentation d’image 3D

La segmentation d’image en 3D est difficile pour plusieurs raisons : Dans de nombreuses techniques d’imagerie microscopique, la qualité de l’image varie dans l’espace : Par exemple, l’intensité et/ou le contraste se dégradent plus on image profondément à l’intérieur d’un échantillon. De plus, les noyaux en contact sont difficiles à différencier de manière automatisée. Enfin, l’anisotropie est difficile à gérer en fonction des algorithmes appliqués et des paramètres respectifs donnés. Certains algorithmes, comme l’approche Voronoi-Otsu-Labeling démontrée ici, ne fonctionnent que pour des données isotropes.

from skimage.io import imread
from pyclesperanto_prototype import imshow
import pyclesperanto_prototype as cle
import matplotlib.pyplot as plt

import napari
from napari.utils import nbscreenshot

# For 3D processing, powerful graphics
# processing units might be necessary
cle.select_device('TX')

<NVIDIA GeForce GTX 1650 with Max-Q Design on Platform: NVIDIA CUDA (1 refs)>

Pour démontrer le flux de travail, nous utilisons des données d’image recadrées et rééchantillonnées du Broad Bio Image Challenge : Ljosa V, Sokolnicki KL, Carpenter AE (2012). Annotated high-throughput microscopy image sets for validation. Nature Methods 9(7):637 / doi. PMID: 22743765 PMCID: PMC3627348. Disponible sur http://dx.doi.org/10.1038/nmeth.2083

input_image = imread("../../data/BMP4blastocystC3-cropped_resampled_8bit.tif")

voxel_size_x = 0.202
voxel_size_y = 0.202
voxel_size_z = 1

À des fins de visualisation, nous montrons les projections d’intensité le long de X, Y et Z.

def show(image_to_show, labels=False):
    """
    Cette fonction génère trois projections : dans les directions X, Y et Z et les affiche.
    """
    projection_x = cle.maximum_x_projection(image_to_show)
    projection_y = cle.maximum_y_projection(image_to_show)
    projection_z = cle.maximum_z_projection(image_to_show)

    fig, axs = plt.subplots(1, 3, figsize=(15, 15))
    cle.imshow(projection_x, plot=axs[0], labels=labels)
    cle.imshow(projection_y, plot=axs[1], labels=labels)
    cle.imshow(projection_z, plot=axs[2], labels=labels)
    plt.show()

show(input_image)
print(input_image.shape)

(86, 396, 393)

Évidemment, la taille des voxels n’est pas isotrope. Ainsi, nous mettons à l’échelle l’image en utilisant la taille des voxels comme facteur d’échelle pour obtenir une pile d’images avec des voxels isotropes.

resampled = cle.scale(input_image, factor_x=voxel_size_x, factor_y=voxel_size_y, factor_z=voxel_size_z, auto_size=True)

show(resampled)
print(resampled.shape)

(86, 79, 79)

Correction d’intensité et de fond

Comme nous pouvons le voir, l’intensité diminue dans la direction Z (de coupe en coupe) et le contraste également. Au moins la décroissance d’intensité peut être corrigée. Dans CLIJx, cette méthode est connue sous le nom de equalize_mean_intensities_of_slices.

equalized_intensities_stack = cle.create_like(resampled)
a_slice = cle.create([resampled.shape[1], resampled.shape[0]])

num_slices = resampled.shape[0]
mean_intensity_stack = cle.mean_of_all_pixels(resampled)

corrected_slice = None
for z in range(0, num_slices):
    # extraire une seule coupe de la pile
    cle.copy_slice(resampled, a_slice, z)
    # mesurer son intensité
    mean_intensity_slice = cle.mean_of_all_pixels(a_slice)
    # corriger l'intensité
    correction_factor = mean_intensity_slice/mean_intensity_stack
    corrected_slice = cle.multiply_image_and_scalar(a_slice, corrected_slice, correction_factor)
    # copier la coupe dans une pile
    cle.copy_slice(corrected_slice, equalized_intensities_stack, z)

show(equalized_intensities_stack)

De plus, l’intensité de fond semble augmenter, potentiellement à cause d’une diffusion plus importante en profondeur dans l’échantillon. Nous pouvons compenser cela en utilisant une technique de soustraction de fond :

backgrund_subtracted = cle.top_hat_box(equalized_intensities_stack, radius_x=5, radius_y=5, radius_z=5)
show(backgrund_subtracted)

Segmentation

segmented = cle.voronoi_otsu_labeling(backgrund_subtracted, spot_sigma=3, outline_sigma=1)
show(segmented, labels=True)

Comme les résultats de segmentation sont difficiles à inspecter en 3D, nous générons une pile d’images avec les intensités originales + les contours de la segmentation. Nous montrons cette pile pour quelques coupes.

a_slice = cle.create([resampled.shape[1], resampled.shape[0]])
segmented_slice = cle.create([resampled.shape[1], resampled.shape[0]])

for z in range(0, resampled.shape[2], 20):
    label_outlines = None
    combined = None

    # obtenir une seule coupe de l'image d'intensité et de l'image d'étiquettes segmentée
    cle.copy_slice(resampled, a_slice, z)
    cle.copy_slice(segmented, segmented_slice, z)

    # déterminer les contours autour des objets étiquetés
    label_outlines = cle.detect_label_edges(segmented_slice, label_outlines)

    # combiner les deux images
    outline_intensity_factor = cle.maximum_of_all_pixels(a_slice)
    combined = cle.add_images_weighted(a_slice, label_outlines, combined, 1.0, outline_intensity_factor)

    # visualisation
    fig, axs = plt.subplots(1, 3, figsize=(15, 15))
    cle.imshow(a_slice, plot=axs[0])
    cle.imshow(segmented_slice, plot=axs[1], labels=True)
    cle.imshow(combined, plot=axs[2])

Visualisation en 3D

Pour une visualisation réelle en 3D, vous pouvez également utiliser napari.

# démarrer napari
viewer = napari.Viewer()

# afficher les images
viewer.add_image(cle.pull(resampled))
viewer.add_image(cle.pull(equalized_intensities_stack))
viewer.add_labels(cle.pull(segmented))

INFO:xmlschema:Resource 'XMLSchema.xsd' is already loaded

<Labels layer 'Labels' at 0x1eba6d51dc0>

viewer.dims.current_step = (40, 0, 0)
nbscreenshot(viewer)

INFO:OpenGL.acceleratesupport:No OpenGL_accelerate module loaded: No module named 'OpenGL_accelerate'

Nous pouvons passer à une vue 3D en cliquant sur le bouton 3D dans le coin inférieur gauche.

nbscreenshot(viewer)

Nous pouvons ensuite aussi incliner et faire pivoter la vue.

nbscreenshot(viewer)