Convolución

Cuando aplicamos un filtro llamado lineal a una imagen, calculamos cada nuevo píxel como una suma ponderada de sus vecinos. El proceso se llama convolución y la matriz que define los pesos se denomina kernel de convolución. En el contexto de la microscopía, a menudo hablamos de la función de dispersión de punto (PSF) de los microscopios. Esta PSF describe técnicamente cómo el microscopio convoluciona una imagen antes de que la guardemos en el disco.

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

Para demostrar el principio de la convolución, primero definimos una imagen de ejemplo que es bastante 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)

A continuación, definimos un kernel de convolución simple, que está representado por una pequeña imagen.

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

Luego, convolucionamos la imagen con el kernel usando scipy.ndimage.convolve. Cuando imprimimos el resultado, podemos ver cómo un 1 en la imagen original se extiende, porque para cada píxel vecino directo, el kernel suma las intensidades de los vecinos. Si hay múltiples píxeles con intensidad > 0 en la imagen original, la imagen resultante calculará la suma en su vecindad. Podrías llamar a este kernel un kernel de suma local.

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

Otros kernels

Dependiendo del kernel que se use para la convolución, las imágenes pueden verse bastante diferentes. Un kernel de media, por ejemplo, calcula la intensidad promedio de los píxeles localmente:

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

El siguiente kernel es una forma simple del operador 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.]])

Para demostrar lo que hacen estos diferentes kernels, los aplicamos a la imagen de MRI mostrada anteriormente.

# abrir el conjunto de datos y extraer un solo plano
noisy_mri = imread('../../data/Haase_MRT_tfl3d1.tif')[90].astype(float)

# hacer zoom recortando una parte
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("Laplaciano")
imshow(convolved_mri, plot=axes[1,0])
axes[1,0].set_title("kernel de suma")
imshow(mean_mri, plot=axes[1,1])
axes[1,1].set_title("Kernel de media")

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