Morphométrie de l’embryon de Tribolium au fil du temps dans Napari

Auteurs : Robert Haase, Daniela Vorkel, 2020

Ceci est la version pyclesperanto d’un flux de travail précédemment publié pour clij2. Original en macro ImageJ

Ce script est un exemple de traitement lourd accéléré par GPU. Il est recommandé d’utiliser une carte graphique dédiée avec au moins 8 Go de mémoire GDDR6. Sinon, cela peut être assez lent.

Commençons par vérifier que pyclesperanto est installé et quel GPU il utilise.

import pyclesperanto_prototype as cle
import numpy as np

cle.select_device("RTX")

<Apple M1 Max on Platform: Apple (2 refs)>

Charger un ensemble de données

Le jeu de données montre un embryon de Tribolium castaneum, imagé par un microscope à feuille de lumière personnalisé, à une longueur d’onde de 488nm (Crédits d’imagerie : Daniela Vorkel, laboratoire Myers, MPI CBG). Il peut être téléchargé depuis zenodo. L’ensemble de données a été rééchantillonné à une taille de voxel de 1x1x1 microns. L’embryon exprime des noyaux-GFP. Nous utiliserons le jeu de données pour détecter les noyaux et générer une segmentation cellulaire estimée.

Toutes les étapes de traitement sont effectuées dans l’espace 3D.

from skimage.io import imread

timelapse = imread('/Users/haase/data/Lund_18.0_22.0_Hours-resampled.tif')

# afficher les dimensions spatiales de l'image
print(timelapse.shape)

(25, 140, 532, 266)

def process_image(image):
    import time

    start_time = time.time()
    
    # pousser l'image dans la mémoire GPU et l'afficher
    gpu_input = cle.push(image)
    # print(gpu_input)
    
    # flou gaussien
    sigma = 2.0
    gpu_blurred = cle.gaussian_blur(gpu_input, sigma_x=sigma, sigma_y=sigma, sigma_z=sigma)

    # détecter les maxima
    gpu_detected_maxima = cle.detect_maxima_box(gpu_blurred)
    
    # seuil
    threshold = 300.0
    gpu_thresholded = cle.greater_constant(gpu_blurred, constant=threshold)

    # masque
    gpu_masked_spots = cle.mask(gpu_detected_maxima, gpu_thresholded)

    # étiqueter les spots
    gpu_labelled_spots = cle.connected_components_labeling_box(gpu_masked_spots)
    # show_labels(gpu_labelled_spots)
    
    number_of_spots = int(cle.maximum_of_all_pixels(gpu_labelled_spots))
    # print("Nombre de spots détectés : " + str(number_of_spots))
    
    # fermeture de la carte d'étiquettes
    number_of_dilations = 10
    flip = cle.create_labels_like(gpu_labelled_spots)
    flop = cle.create_labels_like(gpu_labelled_spots)
    flag = cle.create([1,1,1])
    cle.copy(gpu_labelled_spots, flip)

    for i in range (0, number_of_dilations) :
        cle.onlyzero_overwrite_maximum_box(flip, flag, flop)
        cle.onlyzero_overwrite_maximum_diamond(flop, flag, flip)

    # éroder les étiquettes
    flap = cle.greater_constant(flip, constant=1)
    number_of_erosions = 4
    for i in range(0, number_of_erosions):
        cle.erode_box(flap, flop)
        cle.erode_box(flop, flap)

    gpu_labels = cle.mask(flip, flap)
    
    # récupérer le résultat du GPU sous forme de tableau numpy
    result = cle.pull(gpu_labels).astype(np.uint16)
        
    print("Le traitement a pris " + str(time.time() - start_time) + " s")

    return result

from skimage import data
import napari
viewer = napari.Viewer()

Warning: Could not find scikit-tensor which is needed for separable approximations...
If you want to compute separable approximations, please install it with
pip install scikit-tensor-py3

# taille de voxel z,y,x
calibration = [1, 1, 1]

# fonction pratique pour la visualisation
def show(image):
    viewer.add_image(image, scale=calibration)
    
def show_labels(labels):
    viewer.add_labels(labels, scale=calibration)

# adapté de : https://github.com/tlambert03/napari-dask-example/blob/master/dask_napari.ipynb
import dask
import dask.array as da

# créer une pile dask de lecteurs d'images paresseux
lazy_process_image = dask.delayed(process_image)  # lecteur paresseux
lazy_arrays = [lazy_process_image(timelapse[n]) for n in range(0, timelapse.shape[0])]
dask_arrays = [
    da.from_delayed(lazy_array, shape=timelapse[0].shape, dtype=np.uint16)
    for lazy_array in lazy_arrays
]
# Empiler en un grand dask.array
dask_stack = da.stack(dask_arrays, axis=0)
dask_stack

Array Chunk Bytes 0.92 GiB 37.79 MiB Shape (25, 140, 532, 266) (1, 140, 532, 266) Count 75 Tasks 25 Chunks Type uint16 numpy.ndarray

print(timelapse.shape[0])

25

show(timelapse)
show_labels(dask_stack)

Processing took 0.3502938747406006 s
Processing took 0.268435001373291 s
Processing took 0.2793867588043213 s

from napari.utils import nbscreenshot
nbscreenshot(viewer)