Tribolium-Embryo-Morphometrie über die Zeit in Napari#

Autoren: Robert Haase, Daniela Vorkel, 2020

Dies ist die pyclesperanto-Version eines Workflows, der zuvor für clij2 veröffentlicht wurde. ImageJ Macro Original

Dieses Skript ist ein Beispiel für eine stark GPU-beschleunigte Verarbeitung. Es wird empfohlen, eine dedizierte Grafikkarte mit mindestens 8 GB GDDR6-Speicher zu verwenden. Andernfalls kann es recht langsam sein.

Beginnen wir damit zu überprüfen, ob pyclesperanto installiert ist und welche GPU es verwendet.

import pyclesperanto_prototype as cle
import numpy as np

cle.select_device("RTX")
<Apple M1 Max on Platform: Apple (2 refs)>

Laden eines Datensatzes#

Der Datensatz zeigt einen Tribolium castaneum Embryo, aufgenommen mit einem speziellen Lichtblattmikroskop, bei einer Wellenlänge von 488nm (Bildnachweis: Daniela Vorkel, Myers Lab, MPI CBG). Er kann von Zenodo heruntergeladen werden. Der Datensatz wurde auf eine Voxelgröße von 1x1x1 Mikrometer neu abgetastet. Der Embryo exprimiert Kern-GFP. Wir werden den Datensatz verwenden, um Zellkerne zu erkennen und eine geschätzte Zellsegmentierung zu generieren.

Alle Verarbeitungsschritte werden im 3D-Raum durchgeführt.

from skimage.io import imread

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

# Ausgabe der räumlichen Dimensionen des Bildes
print(timelapse.shape)
(25, 140, 532, 266)
def process_image(image):
    import time

    start_time = time.time()
    
    # Bild in den GPU-Speicher laden und anzeigen
    gpu_input = cle.push(image)
    # print(gpu_input)
    
    # Gaußsche Unschärfe
    sigma = 2.0
    gpu_blurred = cle.gaussian_blur(gpu_input, sigma_x=sigma, sigma_y=sigma, sigma_z=sigma)

    # Maxima erkennen
    gpu_detected_maxima = cle.detect_maxima_box(gpu_blurred)
    
    # Schwellenwert
    threshold = 300.0
    gpu_thresholded = cle.greater_constant(gpu_blurred, constant=threshold)

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

    # Spots markieren
    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("Anzahl der erkannten Spots: " + str(number_of_spots))
    
    # Markierungskarte schließen
    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)

    # Markierungen erodieren
    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)
    
    # Ergebnis von der GPU als numpy-Array zurückholen
    result = cle.pull(gpu_labels).astype(np.uint16)
        
    print("Verarbeitung dauerte " + 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
# Voxelgröße z,y,x
calibration = [1, 1, 1]

# Hilfsfunktion zur Visualisierung
def show(image):
    viewer.add_image(image, scale=calibration)
    
def show_labels(labels):
    viewer.add_labels(labels, scale=calibration)
# angepasst von: https://github.com/tlambert03/napari-dask-example/blob/master/dask_napari.ipynb
import dask
import dask.array as da

# Dask-Stack von trägen Bild-Lesern erstellen
lazy_process_image = dask.delayed(process_image)  # träger Leser
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
]
# In ein großes dask.array stapeln
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
25 1 266 532 140
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)