Mapas de distorsión - Cómo hacer SVD

En esta charla técnica, me gustaría hablar de cómo (y probablemente de cómo no) desarrollar una corrección de la distorsión óptica. Se trata de una problema óptico comúny que, de hecho, puede compensarse mediante un algoritmo, ya que, aunque la imagen esté distorsionada, puede seguir estando libre de aberraciones del frente de onda. Supongo que el título ya lo dice todo, será una continuación de mi historia de amor con la descomposición de valor singular.

Este problema es interesante, no sólo porque es común, sino también porque demostrará muy claramente lo que es el SVD y sus puntos fuertes innatos. Ahora, antes de sumergirnos con algunos números, vamos a establecer un poco de contexto.

Supuestos sobre el problema a resolver

Para un sistema bien alineado con una gran distorsión diseñada, tenemos simetría rotacional, y el polinomio radial será la mejor manera de abordar el efecto. Sin embargo, este no es el problema que me interesa. El problema que me gustaría abordar es uno en el que la distorsión es pequeña y cuando la medimos, no podemos ignorar el ruido de medición. Sin embargo, aunque esta distorsión sea pequeña, no es lo suficientemente pequeña como para que la ignoremos. ¿Dónde encontraríamos una aplicación así? Creo que la litografía puede ser un ejemplo en el que podríamos aspirar a corregirla mejor que una parte entre cien mil. A este nivel de precisión, no podemos suponer nada sobre la naturaleza de la distorsión. La óptica es tan buena como puede serlo y sigue habiendo una distorsión que podemos medir.

El problema de corrección del problema de distorsión encaja bien en el marco de la SVD porque podemos recopilar datos que se ajusten al formalismo matricial. Por lo tanto, supondré que recogeremos la información sobre la distorsión en una cuadrícula que cubre uniformemente una imagen rectangular.

Así que vamos a meternos de lleno, con los dos pies en el suelo. Veamos qué podemos esperar de este enfoque. Para generar algunos datos de distorsión con los que podamos experimentar, he escrito un script Matlab/Octave. Los lectores que quieran experimentar con esto pueden descargarlo aquí. (Puede enviar un correo electrónico a póngase en contacto con para obtener acceso si aún no lo tiene)

N = 6;
[x,y] = meshgrid(-10:10, -13:13);
x = x/max(x(:));
y = y/max(y(:));
%a = randn(1,(N+1)*(N+2)/2).*(1:((N+1)*(N+2)/2))/N^2;
fx = a(1)*x + a(2)*y + ...
      a(3)*x.^2 + a(4)*x.*y + a(5)*y.^2 + ...
      a(6)*x.^3 + a(7)*x.^2.*y + a(8)*x.*y.^2 + a(9)*y.^3 + ...
      a(10)*x.^4 + a(11)*x.^3.*y + a(12)*x.^2.*y.^2 + a(13)*x.*y.^3 + a(14)*y.^4 + ...
      a(15)*x.^5 + a(16)*x.^4.*y + a(17)*x.^3.*y.^2 + a(18)*x.^2.*y.^3 + a(19)*x.*y.^4 + a(20)*y.^5 ;
% a(20)*x.^6 + a(21)*x.^5.*y + a(22)*x.^4.*y.^2 + a(23)*x.^3.*y.^3 + a(24)*x.^2.*y.^4 + a(25)*x.*y.^5 + a(26)*y.^6;

%b = randn(1,(N+1)*(N+2)/2).*(1:((N+1)*(N+2)/2))/N^2;
fy = b(1)*x + b(2)*y + ...
     b(3)*x.^2 + b(4)*x.*y + b(5)*y.^2 + ...
     b(6)*x.^3 + b(7)*x.^2.*y + b(8)*x.*y.^2 + b(9)*y.^3 + ...
     b(10)*x.^4 + b(11)*x.^3.*y + b(12)*x.^2.*y.^2 + b(13)*x.*y.^3 + b(14)*y.^4 + ...
     b(15)*x.^5 + b(16)*x.^4.*y + b(17)*x.^3.*y.^2 + b(18)*x.^2.*y.^3 + b(19)*x.*y.^4 + b(20)*y.^5;
% b(21)*x.^6 + b(22)*x.^5.*y + b(23)*x.^4.*y.^2 + b(24)*x.^3.*y.^3 + b(25)*x.^2.*y.^4 + b(26)*x.*y.^5 + a(27)*y.^6;

fx = fx + 1*randn(tamaño(fx))/400;
fy = fy + 1*randn(tamaño(fy))/400;
[ux,sx,vx] = svd(fx);
[uy,sy,vy] = svd(fy);

SVD y separación de variables: una analogía

El código anterior no es más que una expansión MacLaurin hasta 5º orden con coeficientes aleatorios. La expansión define desplazamientos a lo largo del eje de coordenadas sobre una malla que analizamos mediante la descomposición en valores singulares. El resultado de un análisis SVD es una descomposición matricial que divide la matriz en tres términos,

La matriz "S" es diagonal y (convencionalmente) está ordenada de forma descendente. Las columnas de la matriz U y V son ortogonales y (normalmente) normalizadas,

Antes de sumergirnos en los detalles del análisis de la distorsión de nuestro modelo, me gustaría relacionar la SVD con algo con lo que la mayoría de los lectores pueden estar ya más familiarizados: el método de separación de variables. Cuando queremos resolver una ecuación diferencial parcial analíticamente, suponemos que la solución se puede escribir como,

Esto nos lleva posteriormente a una solución que escribimos como,

Este es exactamente el propósito de la descomposición de valor singular. Ofrece una forma de escribir nuestra matriz M como un producto de funciones que abarcan independientemente el espacio de filas y columnas de la matriz. Cuando preparamos nuestras entradas, hicimos algo parecido: nuestras "funciones base" eran todo tipo de combinaciones de potencias de x y poderes de y multiplicado por algún coeficiente (aleatorio). El análisis SVD redujo nuestro conjunto de 20 funciones de entrada a dos conjuntos (U y V) de 6 funciones. Podemos comprobarlo ejecutando el script (puede descargarse una versión completa del script aquí)

Los resultados anteriores se generaron sin añadir ruido sólo para separar limpiamente lo que es una solución y lo que es (numéricamente) cero.

Análisis de la distorsión con ruido

Cuando nos enfrentamos a datos reales, al menos para los supuestos que he planteado más arriba, no podemos dar por sentado que conozcamos el tamaño de nuestra expansión. Esto es algo que aún tenemos que averiguar. Por lo tanto, vamos a añadir un poco de ruido y ver cómo se ve,

Con el ruido añadido, la estructura de los valores singulares parece diferente. Mientras que en el primer caso los valores singulares caían 12 órdenes de magnitud después de los 6 iniciales, ahora hay lo que parece un cruce entre dos líneas rectas. En realidad, esto es bastante habitual cuando añadimos ruido blanco. Para nuestro modelo de distorsión, el ruido blanco es probablemente una buena suposición porque las mediciones deberían ser independientes.

El script ahora también nos muestra que nuestra recuperación de la distorsión (up-sampled 32×32 veces) no es perfecta. El ruido no nos permite recuperar de forma fiable más de 5 componentes. También podríamos discutir si el 5º componente es ruido o son datos. Al ejecutar el script suficientes veces para recopilar algunas estadísticas, mi propia conclusión fue que el error era menor si lo consideramos datos y lo incluimos. Aún más valioso es el hecho de que, al incluir ese 5º punto, el mapeo suprimió los valores atípicos que, de otro modo, aparecerían 1,5% o 2% de las veces que volvimos a ejecutar la simulación.

Incluir el 5º término tiene un valor muy específico. Hace que la distribución de errores parezca bastante gaussiana. A ojos de nuestros ingenieros (o corazones), debería ser muy deseable. Los valores atípicos pueden ser muy costosos en un contexto de ingeniería, provocando retrasos y frustración. Para cerrar el tema, por valor atípico entiendo un suceso cuya frecuencia no predice la estadística gaussiana, como un error de 6 desviaciones estándar que se produce una vez en una muestra de mil. Tales sucesos deberían requerir dos mil millones de muestras.

Comentarios sobre el guión

Hay algunas cosas en la script sobre los que merece la pena reflexionar un poco. Uno de los que me parece más inquietante es que cuando interpolamos y volvemos a muestrear los vectores singulares izquierdo y derecho (columnas de las matrices U y V), dejan de ser ortogonales. Esto significa que no podemos utilizar los vectores interpolados para proyectar otras distribuciones sobre este espacio. Sin embargo, la falta de ortogonalidad no parece añadir ningún sesgo en la reconstrucción.

Además, las interpolaciones utilizan el número de valores singulares extraídos como orden del polinomio utilizado para la interpolación. Es fácil argumentar a favor de esta elección para un modelo como éste. Sin embargo, esta elección sería una de las primeras que intentaría confirmar (o rechazar) cuando me enfrentara a datos reales.

Si te ha parecido interesante, deja un comentario. No me importaría reescribirlo en un programa en C de código abierto e incluir el pre-warping de los datos de entrada. También hay algunas consideraciones numéricas bastante interesantes sobre cómo estructurar los datos para hacer una aplicación de esta teoría capaz de manejar tasas de mapeo de giga píxeles por segundo, lo que parece ser alcanzable con CPUs genéricas.

Gracias por su lectura.