Morfometría del embrión de Tribolium a lo largo del tiempo en Napari

Autores: Robert Haase, Daniela Vorkel, 2020

Esta es la versión en pyclesperanto de un flujo de trabajo publicado anteriormente para clij2. Original en ImageJ Macro

Este script es un ejemplo de procesamiento intensivo acelerado por GPU. Se recomienda utilizar una tarjeta gráfica dedicada con al menos 8 GB de memoria GDDR6. De lo contrario, puede ser bastante lento.

Comencemos verificando que pyclesperanto esté instalado y qué GPU utiliza.

import pyclesperanto_prototype as cle
import numpy as np

cle.select_device("RTX")

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

Cargar un conjunto de datos

El conjunto de datos muestra un embrión de Tribolium castaneum, capturado por un microscopio de hoja de luz personalizado, a una longitud de onda de 488nm (Créditos de imagen: Daniela Vorkel, laboratorio Myers, MPI CBG). Se puede descargar desde zenodo. El conjunto de datos ha sido remuestreado a un tamaño de vóxel de 1x1x1 micras. El embrión expresa núcleos-GFP. Utilizaremos el conjunto de datos para detectar núcleos y generar una segmentación celular estimada.

Todos los pasos de procesamiento se realizan en el espacio 3D.

from skimage.io import imread

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

# imprimir las dimensiones espaciales de la imagen
print(timelapse.shape)

(25, 140, 532, 266)

def process_image(image):
    import time

    start_time = time.time()
    
    # enviar imagen a la memoria de la GPU y mostrarla
    gpu_input = cle.push(image)
    # print(gpu_input)
    
    # desenfoque gaussiano
    sigma = 2.0
    gpu_blurred = cle.gaussian_blur(gpu_input, sigma_x=sigma, sigma_y=sigma, sigma_z=sigma)

    # detectar máximos
    gpu_detected_maxima = cle.detect_maxima_box(gpu_blurred)
    
    # umbral
    threshold = 300.0
    gpu_thresholded = cle.greater_constant(gpu_blurred, constant=threshold)

    # máscara
    gpu_masked_spots = cle.mask(gpu_detected_maxima, gpu_thresholded)

    # etiquetar puntos
    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("Número de puntos detectados: " + str(number_of_spots))
    
    # cierre del mapa de etiquetas
    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)

    # erosionar etiquetas
    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)
    
    # obtener el resultado de la GPU como un array numpy
    result = cle.pull(gpu_labels).astype(np.uint16)
        
    print("El procesamiento tomó " + 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

# tamaño de vóxel z,y,x
calibration = [1, 1, 1]

# función de conveniencia para visualización
def show(image):
    viewer.add_image(image, scale=calibration)
    
def show_labels(labels):
    viewer.add_labels(labels, scale=calibration)

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

# crear pila dask de lectores de imágenes perezosos
lazy_process_image = dask.delayed(process_image)  # lector perezoso
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
]
# Apilar en un gran 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)