Une introduction à la déconvolution d’image#

En microscopie, les images sont couramment déformées en raison de la nature du système optique, le microscope. Une structure biologique émettant des photons n’apparaîtra pas visuellement dans une image de la scène reflétant 100% la réalité. L’image sera floue, convoluée avec la fonction d’étalement du point (PSF) du montage optique. Si nous connaissons la PSF, il est techniquement possible d’annuler cette convolution dans une certaine mesure.

import numpy as np
from skimage.io import imread, imsave
from skimage.data import cells3d
from skimage.filters import gaussian
from scipy.ndimage import convolve
from pyclesperanto_prototype import imshow
import pyclesperanto_prototype as cle
from napari_simpleitk_image_processing import richardson_lucy_deconvolution
import matplotlib.pyplot as plt

Pour démontrer la convolution et la déconvolution, nous supposerons pour un moment que cette image montrant un noyau reflète la réalité dans l’échantillon biologique.

image = cells3d()[30,1,120:190,80:150]
imshow(image)
../_images/2021282d89929f35e9a6759815d04191be7840f800b96f2b20237b281bf9f60d.png

De plus, nous construisons une image PSF artificielle. Nous commençons par une image PSF parfaite, qui est noire (0) dans tous les pixels sauf un, typiquement mais pas nécessairement au centre.

perfect_psf = np.zeros( (25,25) )
perfect_psf[12, 12] = 1

imshow(perfect_psf, colorbar=True)
../_images/2222ae3d5e3ae13dee9ac16a7b3ef7d5c95be17456d88e68154ef90890e9ddff.png

Notez que dans l’exemple ci-dessus, l’intensité totale de l’image est de 1. L’image floue montrée ensuite aura également une intensité totale d’environ 1. Ceci est crucial lors de la déconvolution ultérieure, car cela permet que l’intensité totale d’une image (ou de toute structure biologique imagée) ne change pas à travers la déconvolution.

À partir de cette PSF parfaite, nous dérivons un exemple de PSF plus réaliste en la floutant avec un flou gaussien.

psf = gaussian(perfect_psf, sigma=2)

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

Convolution#

Comme introduit dans la section sur le filtrage d’image, la convolution est le processus de combinaison de chaque valeur de pixel d’une image d’entrée donnée avec les pixels voisins pondérés selon la PSF.

convolved = convolve(image, psf)

imshow(convolved)
../_images/deb5a3a9284ca00bc99c8b73fa7e8f59a4e6c9e7137847025b1b724589820187.png

Déconvolution#

Un algorithme courant pour restaurer l’image originale à partir de l’image convoluée est la déconvolution de Richardson-Lucy qui est implémentée comme RichardsonLucyDeconvolutionImageFilter de SimpleITK. Ici, nous utilisons une couche de commodité qui fait partie de [napari-simpleitk-image-processing].

number_of_iterations = 10
deconvolved = richardson_lucy_deconvolution(convolved, psf, number_of_iterations)

imshow(deconvolved)
../_images/90e5ef549e09a09032046555638f0ec604a90a25d0e18967370b81f19e728102.png

Comme vous pouvez le voir, cette image déconvoluée n’est pas égale à l’image originale montrée ci-dessus. La déconvolution n’est pas magique. C’est un filtre de traitement d’image et tous les filtres ont des limites. De plus, dans le cas de cet algorithme de déconvolution, le résultat dépend du nombre d’itérations. Nous pouvons empiriquement trouver un bon nombre pour ce paramètre en testant différentes valeurs.

fig, axs = plt.subplots(3, 2, figsize=(20,20))

imshow(image, plot=axs[0, 0])
axs[0,0].set_title("original")

imshow(convolved, plot=axs[0, 1])
axs[0,1].set_title("convoluée")

for i, number_of_iterations in enumerate([10, 20, 40, 80]):

    deconvolved = richardson_lucy_deconvolution(convolved, psf, number_of_iterations)

    axis = axs[(i) % 2 + 1, int((i) / 2)]
    imshow(deconvolved, plot=axis)
    
    axis.set_title("déconvoluée pour " + str(number_of_iterations) + " itérations")
../_images/b54827e4a67830a1c2ca869cd67e5505404ca66247e10e8495278bbe73162576.png

Exercice#

Programmez une boucle for qui déconvolue l’image convolved en utilisant différents number_of_iterations et en déterminant l’erreur quadratique moyenne entre l’image originale et l’image déconvoluée. Indice : Vous pouvez utiliser sklearn.metrics.mean_squared_error.