Convolution

Lorsque nous appliquons un filtre dit linéaire à une image, nous calculons chaque nouveau pixel comme une somme pondérée de ses voisins. Ce processus s’appelle convolution et la matrice définissant les poids est appelée noyau de convolution. Dans le contexte de la microscopie, on parle souvent de la fonction d’étalement du point (PSF) des microscopes. Cette PSF décrit techniquement comment une image est convoluée par le microscope avant que nous ne l’enregistrions sur le disque.

import numpy as np
import pyclesperanto_prototype as cle
from skimage.io import imread
from pyclesperanto_prototype import imshow
from skimage import filters
from skimage.morphology import ball
from scipy.ndimage import convolve
import matplotlib.pyplot as plt

cle.select_device('RTX')

<gfx90c on Platform: AMD Accelerated Parallel Processing (2 refs)>

Pour démontrer le principe de la convolution, nous définissons d’abord une image d’exemple plutôt simple.

image = np.asarray([
  [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
  [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
  [0, 0, 1, 0, 0, 0, 0, 0, 0, 0],
  [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
  [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
  [0, 0, 0, 0, 0, 0, 0, 2, 0, 0],
  [0, 0, 0, 0, 0, 0, 1, 0, 0, 0],
  [0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
  [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
]).astype(float)

imshow(image)

Ensuite, nous définissons un noyau de convolution simple, qui est représenté par une petite image.

kernel = np.asarray([
  [0, 1, 0],
  [1, 1, 1],
  [0, 1, 0],
])

Ensuite, nous convoluons l’image avec le noyau en utilisant scipy.ndimage.convolve. Lorsque nous affichons le résultat, nous pouvons voir comment un 1 dans l’image originale se propage, car pour chaque pixel directement voisin, le noyau somme les intensités des voisins. S’il y a plusieurs pixels avec une intensité > 0 dans l’image originale, l’image résultante calculera leur somme dans leur voisinage. On pourrait appeler ce noyau un noyau de somme locale.

convolved = convolve(image, kernel)

imshow(convolved, colorbar=True)

convolved

array([[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 1., 0., 0., 0., 0., 0., 0., 0.],
       [0., 1., 1., 1., 0., 0., 0., 0., 0., 0.],
       [0., 0., 1., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 2., 0., 0.],
       [0., 0., 0., 0., 0., 0., 3., 2., 2., 0.],
       [0., 0., 0., 0., 0., 1., 1., 3., 0., 0.],
       [0., 0., 0., 0., 0., 0., 1., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]])

Autres noyaux

Selon le noyau utilisé pour la convolution, les images peuvent avoir des aspects très différents. Un noyau de moyenne, par exemple, calcule l’intensité moyenne des pixels localement :

mean_kernel = np.asarray([
  [0, 0.2, 0],
  [0.2, 0.2, 0.2],
  [0, 0.2, 0],
])

mean_convolved = convolve(image, mean_kernel)

imshow(mean_convolved, colorbar=True)

mean_convolved

array([[0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ],
       [0. , 0. , 0.2, 0. , 0. , 0. , 0. , 0. , 0. , 0. ],
       [0. , 0.2, 0.2, 0.2, 0. , 0. , 0. , 0. , 0. , 0. ],
       [0. , 0. , 0.2, 0. , 0. , 0. , 0. , 0. , 0. , 0. ],
       [0. , 0. , 0. , 0. , 0. , 0. , 0. , 0.4, 0. , 0. ],
       [0. , 0. , 0. , 0. , 0. , 0. , 0.6, 0.4, 0.4, 0. ],
       [0. , 0. , 0. , 0. , 0. , 0.2, 0.2, 0.6, 0. , 0. ],
       [0. , 0. , 0. , 0. , 0. , 0. , 0.2, 0. , 0. , 0. ],
       [0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ]])

Le noyau suivant est une forme simple d’un opérateur de Laplace.

laplace_operator = np.asarray([
  [0, 1, 0],
  [1, -4, 1],
  [0, 1, 0],
])

laplacian = convolve(image, laplace_operator)

imshow(laplacian, colorbar=True)

laplacian

array([[ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  1., -4.,  1.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  2.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  3., -8.,  2.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  1., -4.,  3.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.]])

Pour démontrer ce que font ces différents noyaux, nous les appliquons à l’image IRM montrée précédemment.

# ouvrir le jeu de données et extraire un seul plan
noisy_mri = imread('../../data/Haase_MRT_tfl3d1.tif')[90].astype(float)

# zoomer en recadrant une partie
noisy_mri_zoom = noisy_mri[50:100, 50:100]

convolved_mri = convolve(noisy_mri_zoom, kernel)
mean_mri = convolve(noisy_mri_zoom, mean_kernel)
laplacian_mri = convolve(noisy_mri_zoom, laplace_operator)

fig, axes = plt.subplots(2, 2, figsize=(10,10))

imshow(noisy_mri_zoom, plot=axes[0,0])
axes[0,0].set_title("original")
imshow(laplacian_mri, plot=axes[0,1])
axes[0,1].set_title("Laplacien")
imshow(convolved_mri, plot=axes[1,0])
axes[1,0].set_title("noyau de somme")
imshow(mean_mri, plot=axes[1,1])
axes[1,1].set_title("Noyau de moyenne")

Text(0.5, 1.0, 'Mean-kernel')