Una introducción a la deconvolución de imágenes

En microscopía, las imágenes comúnmente se distorsionan debido a la naturaleza del sistema óptico, el microscopio. Una estructura biológica que emite fotones no aparecerá visualmente en una imagen de la escena reflejando la realidad al 100%. La imagen se verá borrosa, convolucionada con la función de dispersión de punto (PSF) del sistema óptico. Si conocemos la PSF, técnicamente es posible deshacer esta convolución hasta cierto punto.

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

Para demostrar la convolución y deconvolución, asumiremos por un momento que esta imagen que muestra un núcleo refleja la realidad en la muestra biológica.

image = cells3d()[30,1,120:190,80:150]
imshow(image)

Además, construimos una imagen PSF artificial. Comenzamos con una imagen PSF perfecta, que es negra (0) en todos los píxeles excepto uno, típicamente pero no necesariamente en el centro.

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

imshow(perfect_psf, colorbar=True)

Ten en cuenta que en el ejemplo anterior la intensidad total en la imagen es 1. También la imagen borrosa que se muestra a continuación tendrá aproximadamente una intensidad total de 1. Esto es crucial al realizar la deconvolución más adelante, porque hace posible que la intensidad total en una imagen (o cualquier estructura biológica fotografiada) no cambie a través de la deconvolución.

A partir de esta PSF perfecta, derivamos un ejemplo de PSF más realista difuminándola con un desenfoque gaussiano.

psf = gaussian(perfect_psf, sigma=2)

imshow(psf, colorbar=True)

Convolución

Como se introdujo en la sección sobre filtrado de imágenes, la convolución es el proceso de combinar cada valor de píxel de una imagen de entrada dada con los píxeles vecinos ponderados según la PSF.

convolved = convolve(image, psf)

imshow(convolved)

Deconvolución

Un algoritmo común para restaurar la imagen original a partir de la imagen convolucionada es la deconvolución de Richardson-Lucy que está implementada como RichardsonLucyDeconvolutionImageFilter de SimpleITK. Aquí usamos una capa de conveniencia que es parte de [napari-simpleitk-image-processing].

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

imshow(deconvolved)

Como puedes ver, esta imagen deconvolucionada no es igual a la imagen original mostrada arriba. La deconvolución no es magia. Es un filtro de procesamiento de imágenes y todos los filtros tienen limitaciones. Además, en el caso de este algoritmo de deconvolución, el resultado depende del número de iteraciones. Podemos encontrar empíricamente un buen número para este parámetro probando diferentes valores.

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("convolucionada")

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("deconvolucionada para " + str(number_of_iterations) + " iteraciones")

Ejercicio

Programa un bucle for que deconvolucione la imagen convolved usando diferentes number_of_iterations y determine el error cuadrático medio entre la imagen original y la imagen deconvolucionada. Pista: Puedes usar sklearn.metrics.mean_squared_error.