Recientemente hemos estado hablando de la deconvolución, pero he notado que quizá los conceptos básicos no están suficientemente claros en general. Voy a intentar dar una descripción sencilla de este problema, que es esencial en muchos campos del procesamiento de señales, incluyendo muy especialmente el de las imágenes astronómicas.

Aquí tenemos una ecuación muy sencilla:

O = I*P + N

que conocemos como ecuación de formación de la imagen. En esta ecuación:

O es la imagen observada.
I es la imagen real, que obtendríamos en condiciones ideales de observación.
N representa el ruido.
P es la función de dispersión de punto (point-spread function, o PSF).

El símbolo * representa aquí a la convolución. La convolución es una operación matemática que tiende a conferir a una función la forma de otra función. Por ejemplo, al aplicar una convolución con una función Gaussiana a una imagen, lo que hacemos es forzar a que las variaciones de brillo existentes en la imagen se parezcan más a la forma de una Gaussiana, que es muy suave. Esto produce un efecto parecido a un desenfoque, de ahí el nombre de desenfoque Gaussiano. Al contrario, si convolucionamos con una función puntiaguda que tiene una salto brusco de valores positivos a negativos, vamos a hacer que las variaciones de brillo se intensifiquen localmente, lo cual produce un efecto de realce de bordes (sharpening).

La ecuación de formación de la imagen nos dice que lo que obtenemos cuando fotografiamos cualquier objeto es básicamente el resultado de una convolución. La PSF (P en la ecuación) es una función que caracteriza las imperfecciones del telescopio y sus limitaciones, pero sobre todo, tratándose de imágenes obtenidas desde la superficie terrestre, los efectos de las perturbaciones atmosféricas (seeing). Podemos pensar en la PSF como en la forma que tiene una estrella en la imagen observada. Como una estrella es puntual, su representación en la imagen observada es una buena medida de lo que ha pasado para que las cosas se hayan estropeado y no hayamos obtenido I, en vez de nuestra patética O

En condiciones ideales de observación, es decir con un telescopio increíblemente grande y desde el espacio exterior, la PSF sería una delta de Dirac, es decir una función consistente en un pico único cuya anchura tiende a cero. La convolución con una delta de Dirac es una identidad, es decir que no modifica la función convolucionada.

Entonces, ¿qué es la deconvolución? Muy sencillo: la deconvolución es la inversa de la convolución. En ausencia de ruido (N=0), conocido el resultado de una convolución:

O = I*P,

y sabiendo qué forma tiene la PSF, es posible revertir el proceso para obtener la función original:

I = O/P

donde el símbolo / representa a la deconvolución.

Esto es cierto en teoría, pero no es tan secillo en el mundo real, ni mucho menos (¿o acaso esperábais otra cosa?). En la práctica, trabajando con imágenes reales, no podemos llegar a un resultado razonable aplicando una deconvolución de forma tan inocente. Hay varias razones para ello. Una, muy importante, es que nuestras imágenes son representaciones (funciones) discretas, lo cual es realmente necesario para construir imágenes digitales. Todo proceso de muestreo necesariamente se realiza con un rango finito (y de hecho bastante limitado, en nuestro caso) de frecuencias, con lo cual una parte de la información se pierde por el camino de forma irreversible. Otra razón, igualmente importante, es que en realidad no conocemos la PSF. Las perturbaciones atmosféricas que intervienen durante la captura de una imagen astronómica se deben a procesos que entran de lleno en la matemática del caos; sólo podemos modelizarlos de forma aproximada, y una pequeña variación en la PSF tiene efectos muy importantes en el resultado de la deconvolución. Y por último, cómo no, tenemos a nuestro compañero inseparable: el ruido, la N en la ecuación. El ruido es una componente aditiva de naturaleza caótica que complica la solución del problema terriblemente, primero porque es muy difícil caracterizarlo correctamente, y además porque es aún más difícil separarlo de la señal (o sea, saber qué es O y qué es N con exactitud).

Por todas estas razones, la deconvolución de imágenes reales es un problema bastante mal condicionado por naturaleza, que no acepta una solución sencilla. Y por todo esto, los algoritmos de deconvolución que diseñamos e implementamos son bastante complicados y refinados. Hay principalmente tres formas de enfocar el problema: los métodos basados en alcanzar una solución de máxima probabilidad (maximum likelihood), los basados en la probabilidad Bayesiana, y los métodos de deconvolución a ciegas (blind deconvolution). No quiero enrollarme demasiado, y además no soy la persona adecuada para entrar en explicaciones más técnicas al respecto. Entre los métodos Bayesianos destaca el algoritmo Pixon, bastante famoso por sus aplicaciones militares. Pixon es un enfoque interesante, pero para aplicaciones astronómicas está superado, sin mencionar el hecho de que hay que pagar royalties para poder implementarlo. Los algoritmos iterativos de Richardson-Lucy y Van Cittert son muy conocidos y se basan en alcanzar una solución de máxima probabilidad, pero en sus formas clásicas tienen grandes limitaciones (ringing e intensificación del ruido). Richardson-Lucy se hizo famoso porque en NASA/JPL lo implementaron para tratar de corregir las imágenes de Hubble, cuando éste estaba miope. Los métodos de deconvolución a ciegas pretenden extraer una PSF plausible a partir de la información existente en la imagen y de cierto conocimiento a priori del contenido de la misma.

En nuestro caso (PixInsight) hemos implementado las que creemos que son las técnicas más evolucionadas y refinadas disponibles actualmente para la deconvolución de imágenes astronómicas, con aplicaciones específicas en la astrofotografía: los algoritmos de Richardson-Lucy y Van Cittert regularizados mediante técnicas de análisis multiescala, basadas en wavelets. La regularización consiste en aplicar una reducción de ruido mientras progresa la deconvolución, iteración a iteración. Durante este proceso, los wavelets nos ayudan a separar el ruido de las estructuras significativas de la imagen, para poder eliminar o atenuar el ruido y potenciar sólo lo que es información relevante. Estas técnicas han sido creadas por Jean-Luc Starck, uno de los gurús del procesamiento multiescala de imágenes (PI está repleto de ideas suyas). Finalmente, nuestro algoritmo de deringing, también basado en gran medida en wavelets, implementa nuestras propias ideas (de Vicent Peris, de Carlos Milovic y mías), y resuelve el problema del efecto Gibbs (ringing) de una forma práctica y muy eficiente. El resultado es una herramienta que realmente funciona y que permite aplicar la deconvolución a imágenes astronómicas que con otras soluciones más simples han estado hasta ahora fuera del alcance de estas técnicas.

Espero que este ladrillo sirva para aclarar un poco de qué va todo esto, y no al revés