Deconvolución de Richardson-Lucy en GPUs compatibles con OpenCL

La deconvolución de Richardson-Lucy es un algoritmo común y básico para la deconvolución de imágenes en microscopía. En este notebook usaremos una versión acelerada por GPU que está implementada en el plugin de napari RedLionFish. Por lo tanto, puedes usar el mismo algoritmo desde la interfaz gráfica de usuario en napari.

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

Cargaremos una imagen que muestra la intensidad fluorescente a lo largo de líneas. Esta imagen 3D fue tomada con un microscopio confocal.

image = imread('../../data/DeconvolutionSampleVerticalGrid1AU-crop.tif')
image.shape

(21, 150, 150)

imshow(image, colorbar=True)

La siguiente imagen PSF fue extraída de imágenes tomadas con el mismo microscopio utilizando el procedimiento explicado anteriormente.

psf = imread('../../data/psf.tif')

imshow(psf, colorbar=True)

Ahora podemos deconvolucionar la imagen usando el algoritmo de Deconvolución de Richardson-Lucy de RedLionFish. Debemos especificar que el algoritmo se ejecutará en la gpu.

iterations = 50

deconvolved = rl.doRLDeconvolutionFromNpArrays(image, 
                                               psf, 
                                               niter=iterations, 
                                               method='gpu', 
                                               resAsUint8=False )
imshow(deconvolved)

ERROR:root:Failed to setup Reikna with OpenCL.
ERROR:root:No module named 'reikna'

Para visualizar de manera más precisa cómo difieren la imagen original y la versión deconvolucionada, podemos graficar la intensidad a lo largo de una línea de izquierda a derecha. Obtenemos estos números de una proyección de intensidad máxima a lo largo de Z.

max_intensity_image = image.max(axis=0)
max_intensity_image.shape

(150, 150)

max_intensity_deconvolved = deconvolved.max(axis=0)
max_intensity_deconvolved.shape

(150, 150)

plt.plot(max_intensity_image[80])
plt.plot(max_intensity_deconvolved[80])
plt.show()

Como puedes ver, el rango de intensidad ha cambiado a través de la deconvolución. Esto depende del algoritmo y la implementación. Cuando apliques deconvolución, considera verificar si la intensidad total en la imagen original y la imagen deconvolucionada están dentro del mismo rango:

image.min(), image.max()

(1, 8)

deconvolved.min(), deconvolved.max()

(0.0, 28.122286)