Richardson-Lucy-Deconvolution on OpenCL-compatible GPUs#

Richardson-Lucy-Deconvolution is a common and yet basic algorithm for image deconvolution in microscopy. In this notebook we will use a GPU-accelerated version of it that is implemented in the napari-plugin RedLionFish. Hence, you can use the same algorithm from the graphical user interface in napari.

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

We will load an image showing fluorescent intensity along lines. This 3D image was taken with a confocal microscope.

image = imread('../../data/DeconvolutionSampleVerticalGrid1AU-crop.tif')
image.shape
(21, 150, 150)
imshow(image, colorbar=True)
../_images/dcb7253afa0708bc1466100c38ca84c604606d441c48a94631fc9f7b7d6212ac.png

The following PSF image was extracted from images taken with the same microscope using the procedure explained before.

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

imshow(psf, colorbar=True)
../_images/a9ef73e14ebfd8a086446101509d41b76b27e3fc686670b9e4d076e15758902b.png

We can now deconvolve the image using RedLionFish’s Richardson-Lucy-Deconvolution algorithm. We should specify that the algorith shall be executed on the 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'
../_images/6196b1264f6a72457f0b4c418c842f03d15c187f3d5a0e3b47dcb093596f0c33.png

To visualize more precisely how the original image and the deconvolved version differ, we can plot the intensity along a line from the left to the right. We retrieve these numbers from a maximum intensity projection along 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()
../_images/784be31367ebbb4a870b07ac894f215ce8a9376033cce6227646848b6b76325a.png

As you can see, the intensity range has change through deconvolution. This depends on the algorithm and implementation. Whenever applying deconvolution, consider checking if the total intensity in the original image and the deconvolved image are within the same range:

image.min(), image.max()
(1, 8)
deconvolved.min(), deconvolved.max()
(0.0, 28.122286)