Segmentación de imágenes en 3D

La segmentación de imágenes en 3D es desafiante por varias razones: En muchas técnicas de imagen por microscopía, la calidad de la imagen varía en el espacio: Por ejemplo, la intensidad y/o el contraste se degradan cuanto más profundo se imagen dentro de una muestra. Además, los núcleos que se tocan son difíciles de diferenciar de manera automatizada. Por último, pero no menos importante, la anisotropía es difícil de manejar dependiendo de los algoritmos aplicados y los respectivos parámetros dados. Algunos algoritmos, como el enfoque de Etiquetado Voronoi-Otsu demostrado aquí, solo funcionan con datos isotrópicos.

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)>

Para demostrar el flujo de trabajo, estamos utilizando datos de imagen recortados y remuestreados del 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 en 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

Para fines de visualización, mostramos proyecciones de intensidad a lo largo de X, Y y Z.

def show(image_to_show, labels=False):
    """
    Esta función genera tres proyecciones: en dirección X, Y y Z y las muestra.
    """
    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)

Obviamente, el tamaño del vóxel no es isotrópico. Por lo tanto, escalamos la imagen con el tamaño del vóxel utilizado como factor de escala para obtener una pila de imágenes con vóxeles isotrópicos.

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)

Corrección de intensidad y fondo

Como podemos ver, la intensidad está disminuyendo en la dirección Z (de corte a corte) y el contraste también. Al menos la disminución de intensidad se puede corregir. En CLIJx, este método se conoce como 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):
    # obtener un solo corte de la pila
    cle.copy_slice(resampled, a_slice, z)
    # medir su intensidad
    mean_intensity_slice = cle.mean_of_all_pixels(a_slice)
    # corregir la intensidad
    correction_factor = mean_intensity_slice/mean_intensity_stack
    corrected_slice = cle.multiply_image_and_scalar(a_slice, corrected_slice, correction_factor)
    # copiar el corte de vuelta en una pila
    cle.copy_slice(corrected_slice, equalized_intensities_stack, z)

show(equalized_intensities_stack)

Además, la intensidad de fondo parece aumentar, potencialmente como resultado de una mayor dispersión en profundidad en la muestra. Podemos compensar esto utilizando una técnica de sustracción de fondo:

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

Segmentación

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

Como los resultados de la segmentación son difíciles de inspeccionar en 3D, generamos una pila de imágenes con las intensidades originales + contornos de la segmentación. Mostramos esta pila para un par de cortes.

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

    # obtener un solo corte de la imagen de intensidad y la imagen de etiquetas segmentada
    cle.copy_slice(resampled, a_slice, z)
    cle.copy_slice(segmented, segmented_slice, z)

    # determinar contornos alrededor de objetos etiquetados
    label_outlines = cle.detect_label_edges(segmented_slice, label_outlines)

    # combinar ambas imágenes
    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)

    # visualización
    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])

Visualización en 3D

Para la visualización real en 3D también puedes usar napari.

# iniciar napari
viewer = napari.Viewer()

# mostrar imágenes
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'

Podemos cambiar a una vista 3D haciendo clic en el botón 3D en la esquina inferior izquierda.

nbscreenshot(viewer)

Luego también podemos inclinar y girar la vista.

nbscreenshot(viewer)