Registro de imágenes mediante transformaciones lineales por trozos

Registro de Imágenes mediante Transformaciones Lineales por Trozos Tesis Doctoral para optar al grado de Doctor en Informática presentada por: Vicente M. Arévalo Espejo Dirigida por: Javier González Jiménez Málaga, Diciembre 2007 AUTOR: Vicente Manuel Arévalo Espejo http://orcid.org/0000-0003-0622-207X EDITA: Publicaciones y Divulgación Científica. Universidad de Málaga Esta obra está sujeta a una licencia Creative Commons: Reconocimiento - No comercial - SinObraDerivada (cc-by-nc-nd): Http://creativecommons.org/licences/by-nc-nd/3.0/es Cualquier parte de esta obra se puede reproducir sin autorización pero con el reconocimiento y atribución de los autores. No se puede hacer uso comercial de la obra y no se puede alterar, transformar o hacer obras derivadas. Esta Tesis Doctoral está depositada en el Repositorio Institucional de la Universidad de Málaga (RIUMA): riuma.uma.es El Dr. Javier González Jiménez, director de la tesis titulada: “Registro de Imágenes mediante Transformaciones Lineales por Trozos” realizada por Vicente M. Arévalo Espejo certifica su idoneidad para la obtención del título de Doctor en Informática. Málaga, Diciembre 2007 Javier González Jiménez Lo que con mucho trabajo se adquiere, más se ama. — Aristóteles (384 AC-322 AC) Filósofo griego. A Carmen y Pablo. Agradecimientos Dedico estas líneas a todas aquellas personas, compañeros y amigos, que me han ayudado durante estos últimos años a llegar a donde estoy. A muchos de ellos, olvidaré citarlos, y no por ello, son menos importantes para mí. A todos ellos, mucha gracias. A Javier, por confiar en mí hace casi 5 años, por su inestimable y gran ayuda, por su apoyo y comprensión en los malos momentos, que los ha habido y muchos, y sobre todo, por su paciencia. Sólo él y yo sabemos lo que significa esta tesis y el trabajo que ha supuesto para ambos. Sin él, no sería una realidad, por todo ello, Javier, muchas gracias. A todos mis compañeros, los que están y los que se fueron, por los buenos momentos que hemos pasado juntos, por su apoyo y ayuda en estos años. A Juan Antonio, Ana, Goyo, Jose Luis, Paco, Tony, Cabello, etc. y en particular, a Cipri, amigo además de compañero, por sus consejos y saber escuchar, por aguantarme en las horas bajas. A él y a todos los demás. Muchas gracias. A mis padres y hermana, por su amor y apoyo incondicional, por su sacrificio cuando las cosas no iban bien, por enseñarme a ser humilde y responsable, por estos últimos años, tan duros para todos. Os quiero y muchas gracias. A Carmen y Pablo, las dos personas más importantes del mundo, mi universo, sin ellos, sin su amor, paciencia y sacrificio esta tesis no sería hoy una realidad. Os quiero mucho, perdón por estos últimos meses y muchas gracias. A todos vosotros, esta tesis también es vuestra :-). i ii AGRADECIMIENTOS Resumen Multitud de aplicaciones de la visión artificial necesitan comparar o integrar imágenes de un mismo objeto pero obtenidas en instantes de tiempo diferentes, con distintos dispositivos (cámaras), desde distintas posiciones, bajo distintas condiciones, etc. Estas diferencias en la captura dan lugar a imágenes con importantes diferencias geométricas relativas que impiden que éstas “encajen” con precisión unas sobre otras. El registro elimina estas diferencias geométricas de forma que píxeles situados en las mismas coordenadas se correspondan con el mismo punto del objeto y, por tanto, ambas imágenes se puedan comparar o integrar fácilmente. El registro de imágenes es esencial en disciplinas como la teledetección, radiología, visión robótica, etc.; campos, todos ellos, que superponen imágenes para estudiar fenómenos medio-ambientales, monitorizar tumores cancerígenos o para reconstruir la escena observada. En esta tesis se aborda la problemática asociada a este proceso, se analizan experimentalmente las técnicas de registro no-rígido más representativas. También se estudian diferentes medidas de similitud utilizadas para medir su consistencia y se propone un novedoso procedimiento para mejorar la precisión del registro lineal por trozos. Concretamente: •Se analizan experimentalmente los elementos que influyen en la estimación de distribuciones de probabilidad de los niveles de intensidad de las imágenes. Estas distribuciones son la base para el cálculo de medidas de similitud basadas en la entropía como la información mutua (MI) o el coeficiente de correlación de entropía (ECC). Por tanto, la efectividad de estas medidas depende críticamente de su correcta estimación. iii xÍNDICE DE FIGURAS 3.1 Familia de transformaciones rígidas . . . . . . . . . . . . . . . 41 3.2 Transformaciones polinomiales. Algunos ejemplos. . . . . . . . 44 3.3 Registro lineal por trozos . . . . . . . . . . . . . . . . . . . . . 45 3.4 Observación fuera del nadir . . . . . . . . . . . . . . . . . . . 48 3.5 Imágenes de prueba utilizadas en el estudio comparativo . . . 51 3.6 Procedimiento utilizado en la selección de los CPs e ICPs . . . 52 3.7 Resultados obtenidos fijando el número de CPs, el ángulo de observación y la geometría de la escena . . . . . . . . . . . . . 54 3.8 Resultados agrupados por ángulo de observación y geometría delaescena ............................ 55 4 Registro lineal por trozos 59 4.1 Ilustración de las transformaciones thin-plate-spline y lineal portrozos ............................. 60 4.2 Redes compatibles e incompatibles con la geometría de la escena 61 4.3 Acciones basadas en aristas empleadas para modificar la topología/geometría de una red dada. . . . . . . . . . . . . . . . 63 4.4 Proyección para-perspectiva de la cámara . . . . . . . . . . . . 65 4.5 Registro lineal por trozos de imágenes afectadas de distorsión proyectivayafín ......................... 66 4.6 Configuración escena-cámara sintética . . . . . . . . . . . . . . 69 4.7 Una red triangular en R3que aproxima la escena de la figura 4.6.................................. 70 4.8 Redes triangulares en R2que aproximan las proyecciones de la escena de la figura 4.6 sobre los planos imagen de las cámaras. 71 4.9 Registro lineal por trozos de dos realizaciones geométricas . . 73 4.10 Cambios topológicos propuestos para eliminar patch reversal en las redes iniciales . . . . . . . . . . . . . . . . . . . . . . . 75 4.11 Inconsistencias de la red producidas por el método de generación de redes propuesto . . . . . . . . . . . . . . . . . . . . . . 76 4.12 Proceso de selección robusta de correspondencias . . . . . . . 83 ÍNDICE DE FIGURAS xi 5 Optimización de redes triangulares 85 5.1 Conjuntos de símplices característicos: star, closure, bound, etc 87 5.2 Acción intercambiar arista .................... 88 5.3 Acción dividir arista ....................... 89 5.4 Acción eliminar arista ...................... 90 5.5 Patch reversal producido la eliminación de la arista . . . . . . 92 5.6 Inserción óptima de nuevos vértices en las redes triangulares . 93 5.7 Proceso de inserción de nuevos vértices en las redes triangulares 95 5.8 Comparación entre “intercambiar arista” y “dividir arista” . . . 97 5.9 Algoritmo greedy propuesto para mejorar el registro lineal por trozos de dos imágenes. . . . . . . . . . . . . . . . . . . . . . . 100 5.10 Ilustración del proceso de optimización propuesto . . . . . . . 102 5.11 Partición binaria de una red triangular. . . . . . . . . . . . . . 104 5.12 Imágenes reales de escenas poliédricas, sus correspondientes redes triangulares de Delaunay y las redes optimizadas . . . . 107 5.13 Evolución de la consistencia del registro global para las imágenes de la figura 5.12(a-d) . . . . . . . . . . . . . . . . . . . . 110 5.14 Red triangular optimizada con el método propuesto (todas las acciones)..............................112 5.15 Imágenes de prueba utilizadas en el estudio comparativo . . . 113 5.16 Evolución de la consistencia global del registro para las imágenes de la figura 5.15 . . . . . . . . . . . . . . . . . . . . . . 116 5.17 Reconstrucciones 3D generadas a partir de las redes: iniciales y optimizadas, de la figura 5.12 . . . . . . . . . . . . . . . . . 118 5.18 Reconstrucción 3D generada a partir de la red optimizada de la zona residencial de la figura 5.15 . . . . . . . . . . . . . . . 119 Apéndices 126 A Entropía, entropía relativa e información mutua 127 A.1 Comparativa H(X)versus p. ..................130 A.2 Relación entre entropía e información mutua. . . . . . . . . . . 136 xii ÍNDICE DE FIGURAS B Métodos de regresión robusta 139 B.1 Comparativa del ajuste mínimo cuadrático y un procedimiento de estimación robusta. . . . . . . . . . . . . . . . . . . . . . . 140 C Topología algebraica 145 C.1 n-símplices, con n= 0,1,2,3. ..................146 C.2 Transformación canónica de un 2-simplex estándar a un 2simplexcualquiera.........................147 C.3 Ilustración de un complejo simplicial abstracto de dimensión 2. 149 C.4 Ilustración de un complejo simplicial de dimensión 3. . . . . . 150 C.5 Realización topológica y geométrica de un complejo simplicial. 151 C.6 Simplex abierto |σ| ⊂ R3vs. Simplex estándar ∆n⊂R2. . . . 152 D Detectores de esquinas. Harris y Lowe 157 D.1 Proceso de detección de esquinas propuesto por Harris. . . . . 159 D.2 Proceso de construcción del espacio de escalas. . . . . . . . . . 162 D.3 Proceso de búsqueda de extremos en las octavas construidas. . 162 D.4 Descriptor de longitud 2×2×8 = 32 obtenido a partir de regiones de 8×8muestras.....................163 Índice de tablas 5 Optimización de redes triangulares 85 5.1 Comparación del método propuesto con otros métodos de similares. Datos de las redes y tiempos . . . . . . . . . . . . . . 108 5.2 Comparación del método propuesto con otros métodos de registrono-rígidos..........................114 5.3 Datos de la red y tiempos empleados en el proceso de optimización ...............................115 Apéndices 126 B Métodos de regresión robusta 139 B.1 Datos utilizados en la comparativa de la figura B.1. . . . . . . 140 C Topología algebraica 145 C.1 Matriz de incidencia de un complejo simplicial. . . . . . . . . . 155 C.2 Matriz de 1-conexión (γ1) de un complejo simplicial. . . . . . . 155 xiii xiv ÍNDICE DE TABLAS Acrónimos CP Punto de Control DDT Triangulación Dependiente de los Datos DEM Modelo Digital de Elevación DTM Modelo Digital del Terreno ECC Coeficiente de Correlación de Entropía GCP Punto de Control Terrestre GIS Sistema de Información Geográfica GPS Satélite de Posicionamiento Global GPU Unidad de Procesamiento Gráfico ICP Punto Independiente de Control MI Información Mutua MLE Estimación de Máxima Probabilidad MR Resonancia Magnética NCC Correlación Cruzada Normalizada PCA Análisis de Componentes Principales PET Tomografía por Emisión de Positrones PIU Partitioned Intensity Uniformity RIU Ratio Image Uniformity xv xvi ACRÓNIMOS RPC Coeficientes de Polinomios Racionales SAD Suma de Diferencias en Valor Absoluto SAR Radar de Apertura Sintética SIFT Scale Invariant Feature Transformation SSD Suma de Diferencias al Cuadrado SURF Speeded Up Robust Features TC Tomografía Computarizada TM Mapa Temático Capítulo 1 Introducción 1.1 Registro de imágenes El registro de imágenes consiste en superponer dos imágenes de la misma escena adquiridas en diferentes instantes de tiempo (análisis multi-temporal), desde distintos puntos de vista (análisis multi-vista) y/o usando distintos sensores (análisis multi-modal). En este proceso, una de las imágenes permanece sin modificar (imagen de referencia ofija), mientras que la otra (imagen de entrada omóvil) se transforma geométricamente hasta que se ajusta a la de referencia. Este proceso es crucial en todas aquellas aplicaciones que necesitan combinar, comparar o fusionar información visual. Uno de los campos de aplicación es la teledetección donde, para poder monitorizar fenómenos medioambientales, como los efectos de unas inundaciones, o humanos, como el impacto del desarrollo urbano sobre una determinada región es necesario superponer con precisión imágenes de fechas diferentes para estudiar el nivel de cambio antes y después del fenómeno a evaluar. Algo similar ocurre en medicina. Por ejemplo, para analizar diversos procesos biológicos, como la evolución de un tumor cancerígeno o una lesión muscular, se comparan resonancias magnéticas (MR) de la zona afectada adquiridas en distintas fechas (antes y después del tratamiento). También, para generar atlas médicos en los que se combinan imágenes de diferentes sujetos 1 21. INTRODUCCIÓN de un determina área del cerebro. En robótica es muy habitual utilizar cámaras para extraer información visual del entorno para detectar y/o manipular objetos, localizarse, desplazarse sin colisionar, reconstruir la escena, etc. Para realizar estas tareas se capturan imágenes desde distintas posiciones que se combinan o se comparan con otras adquiridas previamente. Numerosos sistemas de seguridad basados en visión comparan imágenes adquiridas en distintos instantes de tiempo y/o desde distintos puntos de vista para detectar posibles intrusos en un área restringida. Por otro lado, la reciente aparición de sistemas de bajo coste capaces de captar imágenes retinales1ha favorecido la aparición de sistemas (biométricos) de control de acceso. Estos sistemas comparan la imagen retinal de un sujeto con las almacenadas en el sistema antes de autorizar su acceso. Estas son sólo algunas de las muchas aplicaciones donde interviene el registro de imágenes. En las figuras 1.1 y 1.2 el lector puede encontrar otros ejemplos: generación de mosaicos, diagnóstico de enfermedades, compensación del desenfoque, reconstrucción 3D, etc. 1.2 La precisión en el registro La precisión del registro es vital en la mayoría de las aplicaciones. Sirva como ilustración el siguiente ejemplo. Supóngase que se pretende monitorizar el grado de desertización en una determinada región de interés mediante la comparación de mapas temáticos (TM) del satélite Landsat. Puesto que la resolución espacial de estas imágenes es de 25 m./píxel, cualquier error en el ajuste superior a un píxel desvirtuaría por completo los resultados obtenidos e induciría claramente a medidas incorrectas. Aunque conceptualmente la superposición de dos imágenes es un problema relativamente sencillo, la precisión, y por tanto el éxito, de este proceso se ve condicionado por las condiciones de adquisición de las imágenes (el 1Al igual que las huellas dactilares, las retinas (más concretamente, la disposición de los capilares retinales) de un individuo son únicos. 1.2. LA PRECISIÓN EN EL REGISTRO 3 c) a) b) d) Figura 1.1: Aplicaciones del registro de imágenes: (a) Mosaico de imágenes retinales para el diagnostico de enfermedades oculares; (b) Para ganar en realismo, los video juegos de reciente aparición proyectan vídeo sobre modelos 3D de la escena; (c) Mosaicos esféricos utilizados, típicamente, en realidad virtual y teleoperación; (d) Mosaicos de lechos submarinos para la realización de inventarios arqueológicos. 10 1. INTRODUCCIÓN Diversos autores utilizan el término elástico para referirse a las transformaciones no-rígidas [85, 33]. En esta tesis se utilizarán indistintamente ambos términos. Las funciones rígidas o elásticas, a su vez, se pueden clasificar atendiendo a su alcance o ámbito de influencia en: Globales: La transformación afecta a la totalidad de los píxeles de la imagen. Locales: La influencia de la transformación depende de la posición de cada píxel en la imagen. En la literatura también se pueden encontrar enfoques híbridos, como el propuesto en [41], donde se propone combinar una afinidad global y un conjunto de nfunciones locales de base radial. Euclídea Afín Proyectiva GlobalLocal Elástica Imagen Figura 1.6: Clasificación de las funciones de transformación. En el capítulo 3 se hace una revisión de las funciones más representativas, así como un análisis comparativo de diversas transformaciones no-rígidas (sección 3.3). Los parámetros de las funciones se estiman mediante puntos de control (CP) identificados en las imágenes. El procedimiento utilizado para su identificación robusta se describe en la sección 4.5.1. Por completitud, en 1.3. DEFINICIÓN DEL PROBLEMA 11 el apéndice D se describen las técnicas de detección de esquinas utilizadas en este proceso. Esta tesis se centra, en especial, en la función lineal por trozos, una transformación local que aborda el registro de las imágenes mediante su división en regiones triangulares conjugadas que se registran individualmente mediante transformaciones afines. En los capítulos 4 y 5 se formaliza el registro lineal por trozos y se propone un procedimiento basado en la optimización de redes triangulares para mejorar su precisión. Para formalizar el problema del registro lineal por trozos se proponen unas estructuras matemáticas denominadas complejos simpliciales. Estas estructuras, desarrolladas en el campo de la topología algebraica, se describen formalmente en el apéndice C. 1.3.3 La medida de consistencia La medida de consistencia del registro cuantifica cómo de bien se superponen la imagen registrada y la imagen de referencia. Para ello se puede utilizar una amplia variedad de métricas, las cuales también se pueden agrupar atendiendo a las bases utilizadas: Basadas en puntos: La precisión del registro se determina a partir de los errores geométricos de ajuste (distancia) de un conjunto de correspondencias identificados en las dos imágenes, denominados puntos independientes de control (ICP). Una transformación inadecuada conlleva importantes desajustes, y viceversa. Para medirlos se emplea típicamente el error cuadrático medio (RMSE) y/o el error circular con 90% de confianza (CE90). Ambos se comentan con más detalle en la sección 2.1. Basadas en intensidad: Este tipo de medidas, a diferencia de las anteriores, cuantifican la precisión del registro comparando el contenido de las dos imágenes. Existen tantas posibilidades para ello como enfoques para evaluar la similitud de dos series numéricas de datos como, por ejemplo, el coeficiente de correlación de Pearson, la información mutua (MI), el análisis de componentes principales (PCA), etc. Nótese que, en 12 1. INTRODUCCIÓN este enfoque, la normalización radiométrica de las imágenes es crucial, ya que cambios en el nivel de intensidad debidos, por ejemplo, a una distinta iluminación, serían considerados como errores de ajuste. La elección de la medida de consistencia idónea, crucial en cualquier procedimiento de registro, depende de una variedad de factores como la disponibilidad de ICPs fiables, la naturaleza de las imágenes (esto es, mono-modales o multi-modales), su sensibilidad a las posibles diferencias radiométricas e incluso su coste computacional. En el capítulo 2 se realiza una revisión de las medidas de similitud más representativas. El procedimiento de registro propuesto en el capítulo 5 de esta tesis emplea el coeficiente de correlación de entropía (ECC), una variante normalizada de la MI, para dirigir eficientemente el proceso de optimización. El ECC, al igual que la MI, se estima a partir de las funciones de distribución de probabilidad marginal y conjunta de las intensidades de ambas imágenes. En la sección 2.4 se analiza experimentalmente diversos procedimientos utilizados para la correcta estimación de estas distribuciones. Muchos de los conceptos y términos que se emplean en esta sección son descritos con más detalle en el apéndice A. 1.3.4 La función de interpolación En la mayoría de las ocasiones, el resultado de transformar geométricamente las coordenadas de un píxel no es un par de coordenadas discretas, es decir, las coordenadas transformadas no coinciden con un píxel. Se requiere, por tanto, una función de interpolación. El proceso de interpolación se aborda, típicamente, del siguiente modo (ver fig. 1.7): 1. para cada píxel x=(x, y)>de la imagen interpolada se obtienen las coordenadas origen en la imagen de entrada mediante f−1(x), 2. dependiendo de la función de interpolación, se determina la intensidad/color a transferir a la imagen interpolada, 3. finalmente, se “rellena” el píxel xcon dicho valor. 1.3. DEFINICIÓN DEL PROBLEMA 13 Imagen interpolada Imagen de entrada x   1 fx 1 f  x y b)a) c) Transformación inversa Vecino más próximo 4 vecinos 16 vecinos Figura 1.7: Funciones de interpolación: a) “vecino más próximo”, b) bilinear y c) bicúbica. En la literatura se proponen diferentes funciones de interpolación, algunas poco costosas computacionalmente, como “el vecino más próximo”, pero que produce imágenes con contornos fragmentados (efecto escalón). Otras más costosas, pero que producen imágenes de una mayor calidad visual, como la interpolación bilinear, bicúbica (que utilizan, respectivamente, los 4 o 16 píxeles más cercanos para estimar el valor del píxel transformado), splines, etc. Algunos paquetes de procesamiento de imágenes incorporan técnicas mucho más elaboradas, como la propuesta por Intel en [53], denominada super-sampling, especialmente apropiada para registrar imágenes que presentan cambios de escala muy importantes. 1.3.5 El método de estimación El proceso de registro consiste en determinar los parámetros óptimos de una función de transformación geométrica, esto es, resolver los problemas de minimización o maximización propuestos en (1.2) o (1.3), respectivamente. En este sentido se han propuesto una amplia variedad de procedimientos de optimización y estimación robusta. Así, se pueden encontrar trabajos en los que dada la simplicidad de las función objetivo, por ejemplo, Euclídea o polinomial de orden bajo, se resuelven con formulaciones cerradas [65]. Si la dimensión del problema crece, como ocurre normalmente cuando se emplean funciones de transformación no-rígidas y/o locales, ya no es posible emplear soluciones cerradas y se recurre a técnicas de optimización basadas en algo- 14 1. INTRODUCCIÓN ritmos iterativos (como el descenso del gradiente o la búsqueda greedy [77]), modelos bayesianos [105], el recocido simulado [91], etc. El lector puede dirigirse al trabajo de Maes [68] para un estudio completo sobre diferentes estrategias de optimización. 1.4 Contribuciones de la Tesis Esta tesis aborda el registro de imágenes que presentan diferencias geométricas locales. En este caso sólo las funciones de transformación elásticas y/o locales ofrecen soluciones aceptables en términos de precisión. Cómo medir esta precisión es justamente uno de los aspectos claves en el procedimiento de registro, ya que esta medida es la que guía el proceso de estimación de los parámetros de la función de transformación. Aunque algunas de las contribuciones de este trabajo ya se han apuntado con anterioridad, a continuación se detalla cada una de ellas así como las publicaciones que se han derivado de este trabajo. •Análisis experimental sobre la estimación robusta de funciones distribución de probabilidad de intensidades de la imagen y su aplicación en el cálculo del ECC, una variante normalizada de la MI. Este análisis ha sido publicado parcialmente en [6, 11]. •Evaluación de las técnicas de registro no-rígidas más representativas. Este trabajo ha sido publicado en [2, 4]. En [1, 7, 37] se pueden encontrar aplicaciones de este tipo de transformaciones en el campo de la teledetección. •Formalización del registro lineal por trozos mediante complejos simpliciales. Esta formalización ha sido parcialmente publicada en [3]. •Optimización topológica y geométrica de redes triangulares para mejorar la precisión del registro lineal por trozos. La optimización topológica ha sido publicada en [5, 6, 11] (el primero en proceso de revisión). En la implementación del método de registro lineal por trozos propuesto, así como en las diferentes pruebas experimentales, se ha utilizado una 1.4. CONTRIBUCIONES DE LA TESIS 15 librería open-source de visión por computador denominada OpenCV [100]. Se han realizado diversos trabajos en los que se describen las características de esta librería y sus posibles aplicaciones en el campo de la investigación y la docencia [8, 9, 10]. Por otro lado, para poder extraer robustamente pares de correspondencias en imágenes de satélite de alta-resolución en las que proliferan las sombras, se ha desarrollado una novedosa técnica que permite detectar y delimitar con precisión las zonas de la imagen afectadas [3, 12]. En los trabajos [38, 64, 76], aún no estando directamente relacionados con esta tesis, se aplican conceptos e ideas adquiridos a lo largo de ella para la detección de olivos en imágenes de satélite, la estimación del movimiento mediante imágenes adquiridas con un sistema estéreo y la detección de matrículas, respectivamente. 16 1. INTRODUCCIÓN Capítulo 2 Medidas de consistencia del registro La estimación de los parámetros de la función de transformación se realiza típicamente mediante un proceso iterativo dirigido por una medida de consistencia del registro. Estas medidas tratan de cuantificar la precisión con la que la imagen de entrada, una vez registrada, se ajusta a la de referencia. En la literatura se han propuesto diversas métricas para este fin, las cuales miden cómo de bien se ajustan ambas imágenes en base a: 1. las diferencias radiométricas de la totalidad sus píxeles (medidas basada en intensidad) o 2. los errores de ajuste de un conjunto de correspondencias identificadas en ellas (medidas basadas en puntos). La elección de la medida de similitud idónea es, por tanto, crucial en cualquier procedimiento de registro y depende de diversos factores: la disponibilidad de correspondencias fiables, la diferente naturaleza de las imágenes (esto es, monomodales o multimodales), las diferencias radiométricas, e incluso, su coste computacional. En este capítulo se revisan y comparan las más representativas, haciendo especial hincapié en la información mutua (MI), una medida de la entropía relativa de dos variables aleatorias y que fue propuesta originariamente para el registro de imágenes médicas [67, 104]. La 17 18 2. MEDIDAS DE CONSISTENCIA DEL REGISTRO MI se calcula a partir de las funciones de distribución de probabilidad marginales y conjunta de los niveles de intensidad de las imágenes. La medida MI de dos imágenes depende, por tanto, de la correcta estimación de estas distribuciones. Esta tesis contribuye con una evaluación experimental en la que se comparan diferentes procedimientos para la estimación eficaz de las distribuciones de probabilidad. También se analiza la influencia en esta estimación de parámetros tales como el tamaño de la muestra, el número de intervalos (o bins) del histograma conjunto, la función de interpolación utilizada, etc. En el apéndice A se describen conceptos utilizados en este capítulo tales como la entropía, entropía relativa y la MI. 2.1 Medidas basadas en puntos Las medidas basadas en puntos cuantifican la precisión del registro en base al error en el ajuste de un conjunto de correspondencias, denominados puntos independientes de control (ICP), seleccionadas en ambas imágenes. Estas medidas se caracterizan por su reducido coste computacional y su robustez a las diferencias radiométricas o diferente modalidad de las imágenes. Sin embargo, su efectividad depende críticamente del número de ICPs utilizados, así como de la precisión con la que han sido seleccionados y de lo representativos que sean de las diferencias geométricas de ambas imágenes (esto es, su distribución). Un claro ejemplo de utilización de estas medidas tiene lugar en el campo de la teledetección, y se debe fundamentalmente a lo siguiente: el tamaño de las imágenes, de varios cientos de megabytes en ocasiones, desaconsejan el empleo de cualquier técnica basada en intensidad, mucho más costosas computacionalmente. A continuación se definen las medidas basadas en puntos más utilizadas. Dadas dos imágenes IeI0del mismo tamaño; un conjunto de ncorrespondencias {(xi,x0 i), i = 1, . . . , n}identificadas en ellas y la función de transformación geométrica f, se definen la siguientes medidas de consistencia: 2.1. MEDIDAS BASADAS EN PUNTOS 19 •Error cuadrático medio: RMSE =v u u t 1 n n X i=1 d2 i(2.1) •Error circular con 90% de confianza: P(di≤CE90) = 90 % (2.2) donde di=kxi−f(x0 i)k. La principal diferencia entre ambas radica en lo siguiente: mientras el RMSE considera los errores de todos los ICPs, incluyendo los correspondientes a las posibles correspondencias espurias, el CE90 proporciona el valor que acota superiormente el 90% de los errores, esto es, no proporciona información sobre cómo de mal se ajusta es el 10 % restante. La figura 2.1 ilustra gráficamente el significado de ambas medidas. ICP estimado ICP  x' i f a) b) i d  xx' i df  0d 90% d 2 1 1n i i RMSE d n ¦  90 90% i Pd CEd xii d Figura 2.1: Medidas de similitud basadas en puntos: (a) RMSE mide el error de ajuste medio de todos los ICPs, incluidos los correspondientes a los espurios. (b) CE90 mide el error máximo que acota superiormente el 90 % de los errores de ajuste. 26 2. MEDIDAS DE CONSISTENCIA DEL REGISTRO (con distribuciones diferentes) o cero (con distribuciones iguales), no es una distancia en el sentido más estricto de la palabra ya que no es simétrica ni verifica la desigualdad del triángulo, circunstancias que dificultan su empleo como medida de consistencia del registro [26]. 2.3.1 Información mutua La información mutua (MI) (ver def. A.5) mide la dependencia estadística o redundancia de información de dos variables aleatorias. Se deriva a partir de la entropía relativa pero, a diferencia de ésta, es simétrica y verifica la desigualdad del triángulo. En 1997, Viola [104] y Maes [67] proponen, por primera vez y casi simultáneamente, el empleo de la MI como medida de consistencia del registro. Tradicionalmente, el parecido o similitud entre dos imágenes se ha medido con la SSD o su variante normalizada, la NCC que, cómo se ha comentado anteriormente, admite diferencias lineales entre las paletas de las dos imágenes [23]. Cuando esta relación no se verifica para toda la paleta, la NCC es incapaz de captar el parecido de éstas. La MI, sin embargo, no asume una relación funcional (normalmente lineal) entre los niveles de intensidad de las imágenes a comparar, sino su relación estadística. Es decir, no asume que la naturaleza radiométrica de las imágenes sea la misma o muy similar, por tanto, se puede utilizar con imágenes multimodales o monomodales con perfiles radiométricos muy dispares. Para ilustrar la ventaja de la MI frente a la NCC se presentan dos experimentos. En el primero, se alinean (mediante una rotación) los imágenes sintéticas de la figura 2.2, midiendo la consistencia del registro con ambas métricas. Como se observa en la figura 2.4(a), la MI proporciona un máximo en cero grados (esto es, cuando ambas imágenes encajan perfectamente), a pesar de que las intensidades de las imágenes no coinciden; este tipo de correlación pasa completamente desapercibida para la NCC. Lo mismo ocurre con el ejemplo que se muestra en la figura 2.4(b). En este caso, un trozo de una imagen de satélite QuickBird es desplazada sobre otra mayor pero adquirida seis meses después (esto es, diferentes condiciones de 2.3. MEDIDAS BASADAS EN LA ENTROPÍA 27 iluminación y cambios en el contenido debidos al cambio estacional, cambios temporales en la escena, etc.). Estos resultados ponen de relieve la robustez y efectividad de la MI cuando se utiliza con pares de imágenes con diferencias radiométricas no-funcionales. 50 100 150 200 250 50 100 150 200 250 50 100 150 200 250 50 100 150 200 250 b) Imagen 1 Correlación cruzada normalizada vs. Información mutua Píxeles Píxeles Imagen 2 Máximo global Máximo global a) Imagen 1 Imagen 2 Correlación cruzada normalizada vs. Información mutua Grados -30 -20 -10 010 20 30 0 0.2 0.4 0.6 0.8 1MI NCC Figura 2.4: Experimentos que ilustran la idoneidad de la MI frente a la NCC como medida de consistencia. (a) Un patrón sintético es rotado (desde −30◦hasta +30◦) sobre otro con la misma estructura (esto es, un cuadro de igual tamaño) pero con distintos niveles de gris; y (b) un pequeño trozo de una imagen de satélite QuickBird es desplazado sobre otra imagen de la misma escena adquirida en otra fecha. Obsérvese como, a diferencia de la MI, la NCC falla en ambos experimentos, esto es, no detecta un máximo global en (a) y proporciona un máximo global erróneo en (b). La MI se relaciona con la entropía mediante las siguientes expresiones (ver sec. A.4 para un desarrollo completo): MI I,ˆ I0=H(I) + Hˆ I0−HI,ˆ I0(2.13) donde H(I) yHˆ I0son las entropías de Iyˆ I0, respectivamente, y HI,ˆ I0 28 2. MEDIDAS DE CONSISTENCIA DEL REGISTRO la entropía conjunta. En la práctica son frecuentes los trabajos que utilizan derivaciones normalizadas de la MI, por ejemplo en Maes [67] y Chen [23]. A continuación se definen las dos más representativas: •Información mutua normalizada [97]: NMI I,ˆ I0= H(I) + Hˆ I0 HI,ˆ I0(2.14) •Coeficiente de correlación de entropía [13]: ECC I,ˆ I0= 2 − 2HI,ˆ I0 H(I) + Hˆ I0(2.15) que toman valores en los intervalos cerrados [0,2] y[0,1], respectivamente. Estas variantes normalizadas aportan ciertas ventajas respecto a sus correspondiente no normalizada. La más importante es que permite comparar cómo de bien se ajustan dos regiones conjugadas de la imagen con respecto a otras dos cuando el tamaño de ambos pares es diferente. Esta circunstancia es clave en el procedimiento de registro propuesto en el capítulo 5. El cálculo de la MI depende de forma crucial de la correcta estimación de las funciones de distribución marginales y conjunta de los niveles de intensidad de ambas imágenes, convirtiéndose éste en el principal escollo para la utilización de esta medida. A continuación se analiza esta problemática, así como diferentes procedimientos para la correcta estimación de las distribuciones. También se analiza la influencia del procedimiento de interpolación en la estimación de estas funciones. 2.4. ESTIMACIÓN DE LA DISTRIBUCIÓN DE PROBABILIDAD 29 2.4 Estimación de la distribución de probabilidad Una estimación imprecisa de las distribuciones de probabilidad introduce errores en la MI que derivan en una medida errónea de la similitud (expresiones (2.7) y (2.13)). Existen diversos factores que influyen en la estimación de estas distribuciones como, por ejemplo, la elección del número de intervalos del histograma conjunto (esto es, el número de bins) y el tamaño de las imágenes (número de píxeles considerados). Estos elementos determinan cómo de representativo es el histograma conjunto obtenido a partir de los niveles de intensidad de las imágenes. Así, por ejemplo, si se desea medir la similitud de una imagen de 32 ×32 píxeles, de 256 niveles de gris, superpuesta sobre otra de mayor tamaño, el histograma conjunto será una matriz de 256 ×256 donde se acumulan sólo 32 ×32 pares, lo que da lugar a un histograma poco representativo. En estos casos, si se reduce el número de bins del histograma conjunto (expresión 2.8) utilizando, por ejemplo, 64 ó 32 en lugar de 256 (en el caso de imágenes en escala de gris) se obtiene un histograma más representativo y, por tanto, mejores estimaciones de las distribuciones de probabilidad. Con esta actuación también se logra reducir el tiempo de computo (menos términos en el sumatorio de la expresión (2.7)) y dotar de robustez al proceso frente al ruido. Reducir el número de intervalos del histograma en función del número de muestras disponibles, como se propone arriba, no siempre garantiza una estimación apropiada de las funciones de distribución. Existen métodos de estimación no-paramétricos, más costosos computacionalmente, que tratan de minimizar esta limitación infiriendo funciones de distribución a partir de un conjunto limitado de muestras. Este es el caso de la ventana de Parzen [104], que propone acumular un kernel gaussiano K, centrado en cada par (i, j), en lugar de unidades discretas. Más concretamente, la ventana de Parzen consiste en una acumulación ponderada, función de σ, con aportación máxima en la celda (i, j)y que decrece a medida que se aleja de ella, de modo que, cuanto más pequeña sea la desviación más se asemejará a la dis- 30 2. MEDIDAS DE CONSISTENCIA DEL REGISTRO tribución obtenida con el método paramétrico. La estimación vendría dada por la siguiente expresión: p(i, j)≈p∗(i, j) = 1 NX (a,b)∈S K(i−a, j −b)(2.16) donde Ses el subconjunto de muestras utilizado. Naturalmente, este procedimiento asume que el subconjunto de muestras utilizado en la estimación es representativo de la distribución real, en cuyo caso, se obtiene una mejor estimación de ésta. La figura 2.5(a) ilustra gráficamente este procedimiento mediante la estimación de una distribución Gaussiana a partir de un conjunto reducido de muestras. 050 100 150 200 250 0 0.005 0.01 0.015 0.02 0.025 0.03 050 100 150 200 250 0 0.005 0.01 0.015 0.02 0.025 0.03 Ventana de Parzen Sigma Niveles de gris b) c) Método paramétrico vs. Ventana de Parzen a) Figura 2.5: Estimación de una distribución de probabilidad mediante la ventana de Parzen. (a) Estimación de una Gaussiana a partir de un conjunto reducido de muestras (ejemplo extraído de [104]). (b) Histograma normalizado de niveles de gris (expresión (2.11)). (c) Estimación obtenida con la ventana de Parzen. Las figuras 2.5(b,c) muestran las distribuciones de probabilidad de una imagen en escala de gris estimadas acumulando directamente los niveles de gris (histograma normalizado) y mediante la ventana de Parzen, respectivamente. Obsérvese cómo la estimación obtenida con la ventana de Parzen es mucho menos ruidosa que la obtenida a partir del histograma. 2.4. ESTIMACIÓN DE LA DISTRIBUCIÓN DE PROBABILIDAD 31 2.4.1 Evaluación experimental Independientemente del método utilizado, paramétrico o la ventana de Parzen, la estimación de la función de distribución depende, entre otros, del número de bins y del tamaño de la imagen y, en el caso de la ventana de Parzen, de la desviación típica del kernel gaussiano. En esta sección se analiza la influencia de estos parámetros en las estimaciones proporcionadas por ambos procedimientos. Para realizar este análisis se han efectuado diversas pruebas consistentes, básicamente, en desplazar horizontalmente y píxel a píxel (∆x= 1), trozos de imágenes de satélite IRS-1D (L= 256) sobre otra imagen de la misma zona y de mayor tamaño (imagen de referencia), calculando el ECC en cada posición. Recuérdese que, de acuerdo con la expresión (2.15), para calcular el ECC se requieren las distribuciones de probabilidad conjunta y marginales de los niveles de intensidad de ambas imágenes (expresión (2.8)). Los resultados de estas pruebas se recogen en las figuras 2.6 y 2.7. Del análisis de estas gráficas se extraen diversas conclusiones, algunas de ellas vienen a corroborar algunos aspectos comentados anteriormente: 1. El método paramétrico no produce distribuciones de probabilidad representativas cuando el tamaño de la imagen, en comparación con la dimensión del histograma, no es importante (fig. 2.6 abajo). Cuando el número de muestras utilizado en la estimación es grande (fig. 2.6 arriba), las curvas muestran un pico único y distintivo, independientemente del número de bins considerado. Obsérvese, también, como una reducción del número de intervalos del histograma (32 en lugar de 64, 128 ó 256) produce mejores resultados (fig. 2.6 abajo). 2. La ventana de Parzen produce mejores estimaciones de las distribuciones de probabilidad incluso cuando el número de muestras, en relación con el tamaño del histograma, es pequeño (fig. 2.7(a) abajo), mostrando un mejor comportamiento que el método de estimación tradicional en idénticas condiciones (fig. 2.7(a) arriba). Obsérvese en la figura 2.7(b) cómo, lógicamente, cuando σdecrece la respuesta de la ventana de Parzen se asemeja a la obtenida con el método paramétrico. 32 2. MEDIDAS DE CONSISTENCIA DEL REGISTRO Imagen referencia Imagen entrada Sentido del desplazamiento Tamaño muestra = 128 u 128) 128x128 64x64 32x32 ventana 0.4 0.6 C C 256 128 64 0 0.2 0.4 E C 32 20 40 60 80 100 120 0 Tamaño muestra = 64 u 64) ventana 0.4 0.6 CC 0 0.2 E 20 40 60 80 100 120 140 160 180 Tamaño muestra = 32 u 32) ventana 0.4 0.6 E CC 20 40 60 80 100 120 140 160 180 200 220 0 0.2 E 20 40 60 80 100 120 140 160 180 200 220 Figura 2.6: Valores de ECC obtenidos al desplazar imágenes cuadradas de 128, 64 y 32 píxeles de lado sobre la imagen de referencia. La estimación de las funciones de distribución se ha realizado con el método paramétrico modificando el número de bins y el tamaño de la ventana. 2.4. ESTIMACIÓN DE LA DISTRIBUCIÓN DE PROBABILIDAD 33 Imagen referencia Imagen de entrada Máximo global 010 20 30 40 50 60 70 80 0.1 0.2 0.3 0.4 0.5 0.6 Método paramétrico ECC (bins=32) 010 20 30 40 50 60 70 80 0 0.05 0.1 0.15 0.2 Ventana parzen ECC ( V =1) 128 64 32 10 20 30 40 50 60 70 0 0.05 0.1 0.15 0.2 0.25 Método paramétrico ECC (bins=32) 10 20 30 40 50 60 70 0 0.05 0.1 0.15 0.2 0.25 Ventana parzen ECC (bins=32) 0.25 0.50 1 a) b) Máximo global Sentido del desplazamiento Figura 2.7: Valores de ECC obtenidos al desplazar un fragmento cuadrado de la imagen de entrada (de 32 píxeles de lado) sobre la imagen de referencia. (a) Estimación de las funciones de distribución realizada con el método paramétrico (arriba) y la ventana de Parzen (abajo) modificando el número de bins L= 128,64,32. (b) Estimación de la función de distribución con la ventana de Parzen utilizando diferentes desviaciones típicas σ= 0,25,0,5,1. 34 2. MEDIDAS DE CONSISTENCIA DEL REGISTRO Como ilustración, obsérvese la figura 2.8 donde se muestra el aspecto de las distribuciones de probabilidad obtenidas con el método paramétrico y la ventana de Parzen en la posición denotada como “máximo global” en la imagen de la figura 2.7, y que representa el posición óptima de la imagen de entrada sobre la de referencia. b) Distribución de probabilidad conjunta a) Figura 2.8: Función densidad de probabilidad conjunta. (a) Estimación obtenida con el método paramétrico. (b) Estimación obtenida con la ventana de Parzen. 2.4.2 La función de interpolación Una vez estimados los parámetros de la función de transformación se ha de proceder a transferir los niveles de intensidad de la imagen de entrada a la imagen transformada. Típicamente, el resultado de transformar una coordenada discreta es un valor real, esto es, el valor que hay que transferir no se corresponde con un píxel concreto de la imagen de entrada sino con una posición cercana. Para determinar el nivel de intensidad a transferir se emplean funciones de interpolación. Las funciones de interpolación generan los niveles de gris a partir del conjunto de píxeles que rodean la coordenada transformada, de tal forma, que se generan distribuciones diferentes dependiendo del conjunto de píxeles que se consideran en este cálculo. Este aspecto es crítico en los procedimientos de registro basados en intensidad, ya que una estimación no válida de la función de distribución puede impedir que el proceso de optimización converja. 2.4. ESTIMACIÓN DE LA DISTRIBUCIÓN DE PROBABILIDAD 35 1 Tamaño ventana = 128 128, bins = 256, 0.5x' u 0.5 1 ecc a ) 050 100 150 200 250 0 1 c ) 0 50 100 150 200 250 0 0.5 ec c 0 50 100 150 200 250 1 b) 0.5 ecc 1 050 100 150 200 250 0 c) 0.5 ecc 050 100 150 200 250 0 d) Figura 2.9: Valores de ECC obtenidos al desplazar cada 0.5 píxeles una imagen cuadrada de 128 píxeles de lado sobre la imagen de referencia. La interpolación de los niveles de gris se ha realizado con diferentes funciones: (a) “el vecino más próximo” (NN), (b) bilineal (B), (c) volumen parcial (PV) y (d) interpolación parcial (PI). Las gráficas 2.9(a-b) recogen los valores de ECC obtenidos al repetir el experimento de la figura 2.6, pero en esta ocasión, en lugar de píxel a píxel, la imagen (de 128 de lado) se desplaza cada 0,5píxeles (∆x= 0,5) lo que obliga a interpolar. Como se observa en la figura, las curvas correspondientes a las funciones de interpolación clásicas, el “vecino más próximo” (NN) y bilineal (BI), aunque muestran un pico claro y distintivo, exhiben comportamientos bastantes ruidosos. Supóngase, por ejemplo, un procedimiento de registro basado en intensidad dirigido por el ECC y que utiliza NN o BI para interpolar los niveles de intensidad de la imagen registrada. El método de optimización actúa sobre los 42 3. REGISTRO BASADO EN PUNTOS (ver fig. 3.1(d)). Afín: La transformación afín viene dada por la expresión: x=x0a1+y0a2+a3 y=x0b1+y0b2+b3 (3.2) donde los términos de la forma aiybi,i= 1,2,3, son los coeficientes de la transformación en xey, respectivamente. Esta es una transformación con 6 grados de libertad que se caracteriza por preservar el paralelismo, el ratio de áreas, el ratio de longitudes de líneas colineales o paralelas y las combinaciones lineales de vectores (por ejemplo, centroides). Se requieren al menos 3 correspondencias no colineales para su estimación (ver fig. 3.1(b)). Las funciones afines se utilizan en el registro lineal por trozos para transformar las diferentes regiones triangulares en las que se dividen las imágenes conjugadas. El lector encontrará información más detallada sobre estas transformaciones en las secciones 3.2 y 4.2. Proyectiva: Transformación con 8 grados de libertad dada por la expresión: x=x0a1+y0a2+a3 1 + x0c1+y0c2 y=x0b1+y0b2+b3 1 + x0c1+y0c2 (3.3) donde los términos de la forma ai,bi,i= 1,2,3ycj,j= 1,2son los coeficientes de la transformación en xey, y el vector de distorsión perspectiva, respectivamente. La transformación proyectiva preserva la concurrencia, la colinealidad, el orden de contacto, etc., pero a diferencia de las transformaciones anteriores, no preserva el paralelismo. Para su estimación se necesitan como mínimo 4 correspondencias no colineales (ver fig. 3.1(a)). 3.2. TRANSFORMACIONES NO-RÍGIDAS 43 Normalmente, el número de correspondencias identificadas en ambas imágenes es superior al mínimo requerido por cada transformación. Esta circunstancia es explotada por diferentes procedimientos de estimación para minimizar el error de ajuste, como por ejemplo el ajuste de mínimo error cuadrático (RMSE) o los estimadores robustos LMedS, RANSAC o MAPSAC, que se describen en el apéndice B. Por otro lado, obsérvese como cada transformación extiende a la anterior, es decir, la similaridad extiende la Euclídea (similaridad = euclídea + cambio de escala), la afinidad extiende la similaridad, y así sucesivamente. El lector puede dirigirse al libro de Hartley y Zisserman [47] para un estudio completo de este tipo de transformaciones y sus invarianzas geométricas. 3.2 Transformaciones no-rígidas Salvo en casos muy concretos, tales como imágenes de escenas planas o imágenes de escenas no-planas tomadas desde la misma posición o después de un movimiento de rotación puro, las transformaciones rígidas rara vez son capaces de modelar las distorsiones, típicamente no lineales, que presentan un par de imágenes conjugadas. Son en estas situaciones donde las transformaciones no-rígidas o elásticas son requeridas. A continuación se detallan las transformaciones elásticas más representativas: Polinomial: Las transformaciones polinomiales se utilizan típicamente con imágenes que presentan diferencias geométricas que afectan a toda la imagen por igual (alcance global) (ver fig. 3.2). Este sería el caso en el que la imagen de referencia y la de entrada son adquiridas desde la misma posición o la escena es prácticamente plana. Este tipo se transformaciones se emplean con frecuencia en el campo de la teledetección para el registro de imágenes que son adquiridas desde una órbita fija y sin cabeceo del satélite [30, 37, 80]. Las funciones polinomiales tienen N= (m+ 1) (m+ 2) (3.4) grados de libertad, siendo mel grado de las funciones polinomiales. 44 3. REGISTRO BASADO EN PUNTOS Requieren, por tanto, de N/2correspondencias como mínimo para su estimación. Lógicamente, funciones polinomiales con grados elevados (esto es, más grados de libertad) modelan diferencias geométricas más severas y viceversa. Matemáticamente, una transformación polinomial se escribe como: x= m X i=0 i X j=0 aij(x0)i−j(y0)j y= m X i=0 i X j=0 bij(x0)i−j(y0)j (3.5) donde aij,bij son los coeficientes de los polinomios en xey, respectivamente. Este tipo de transformaciones no ofrecen resultados aceptables con imágenes que presentan diferencias locales, o sea, que no afectan a toda la imagen por igual, lo que limita significativamente su campo de aplicación. Este hecho queda patente en los resultados de la evaluación experimental que se presenta en la sección 3.3. Imagen original Algunas deformaciones polinomiales Figura 3.2: Transformaciones polinomiales. Algunos ejemplos. Lineal por trozos: Las funciones lineales por trozos abordan el registro dividiendo las imágenes en regiones triangulares conjugadas (por ejem- 3.2. TRANSFORMACIONES NO-RÍGIDAS 45 a) ' j t ' i t ˆ' j t ˆ' i t Líneas “quebradas” Imagen de entrada Imagen registrada Dominio de la función lineal por trozos Imagen Envolvente convexa de los CPs b) Figura 3.3: Registro lineal por trozos. (a) Efecto de líneas “quebradas” en las transiciones de los triángulos. (b) Dominio de la función lineal por trozos. Esta función sólo está definida dentro de la envolvente convexa del conjunto de puntos identificados en la imagen. plo, mediante el método de triangulación de Delaunay [93]) que son registradas individualmente mediante transformaciones lineales [6, 39]. Aunque este enfoque garantiza la continuidad en triángulos adyacentes, no produce transiciones suaves entre ellos, lo que puede causar efectos no deseados en la imagen transformada, tales como el “quebramiento” de las líneas, esto es, las rectas pueden no preservarse en las transiciones (ver fig. 3.3(a)). Matemáticamente se define como: x=fx(x0, y0) =        a11 +a12x0+a13y0si (x0, y0)∈t1 . . . an1+an2x0+an3y0si (x0, y0)∈tn y=fy(x0, y0) =        b11 +b12x0+b13y0si (x0, y0)∈t1 . . . bn1+bn2x0+bn3y0if (x0, y0)∈tn (3.6) donde ties un elemento triangular construido sobre tres corresponden- 46 3. REGISTRO BASADO EN PUNTOS cias; aij ybij, con j= 1,2,3, son los coeficientes de las transformación linear correspondiente a ti; y nes el número de elementos que configuran la malla triangular. Nótese que las funciones lineales por trozos sólo están definidas dentro de la envolvente convexa del conjunto de correspondencias (ver fig. 3.3(b)). Aunque es posible extrapolar fuera de esta región, en la práctica no es aconsejable ya que pueden introducir errores que distorsionan la imagen registrada. En el capítulo 4 se analizan en profundidad este tipo de transformaciones. Base radial: Las transformaciones de base radial son métodos de interpolación de datos dispersos donde la transformación geométrica viene dada por la combinación lineal de funciones de base radial simétricas centradas en un punto de control (CP) particular [19, 41, 45]. Las funciones de base radial proporcionan deformaciones suaves con un comportamiento fácilmente controlable. En 2D, una función de base radial consta típicamente de una componente global (habitualmente, una transformación afín) que corrige posibles traslaciones, rotaciones, cambios de escala, etc. (primer término de la expresión (3.7)) y una componente local capaz de modelar distorsiones no lineales. Matemáticamente: xi= m X j=0 j X k=0 ajk (x0 i)j−k(y0 i)k+ n X j=1 Ajg(rj) yi= m X j=0 j X k=0 bjk (x0 i)j−k(y0 i)k+ n X j=1 Bjg(rj) (3.7) donde rj= (x0 i, y0 i)−x0 j, y0 j (3.8) 3.3. COMPARATIVA EXPERIMENTAL 47 yges una función de base radial, habitualmente: Thin−plate−spline g(rj) = r2 jlog r2 j Multi−quadric g(rj) = r2 j+δ±µ, δ > 0, µ 6= 0 Gaussiana g(rj) = e(−r2 j/σ), δ > 0 Shifted−LOG g(rj) = log r2 j+δ3 2, δ ≥0 Cubic−spline g(rj) = krjk3 La función base determina la influencia de cada CP sobre cada píxel de la imagen registrada (segundo término de la expresión (3.7)), esto es, el alcance de cada CP. Algunas funciones base tienen una influencia más global (por ejemplo, thin-plate-spline) mientras otras tiene una influencia mucho más local (por ejemplo, Gaussiana). El lector puede dirigirse a los trabajos de Goshtasby [42, 51] para un estudio completo de este tipo de transformaciones en términos de precisión y complejidad computacional. 3.3 Comparativa experimental En la práctica, debido una amplia variedad de factores como el tamaño y resolución espacial de las imágenes, la geometría de la escena, las posiciones de adquisición, etc., una transformación que ofrece buenos resultados con ciertas imágenes puede no producir resultados aceptables con otras, requiriéndose técnicas capaces de corregir distorsiones más complejas (típicamente, no lineales). Por ejemplo, en el campo de la teledetección, las funciones globales polinomiales muestran habitualmente un buen comportamiento con imágenes de bajay media-resolución (por ejemplo, Landsat, IRS, Spot, etc.) adquiridas en el nadir, pero pueden no ser suficientemente efectivas para registrar imágenes de alta-resolución (por ejemplo, QuickBird, Ikonos o OrbView) ya 48 3. REGISTRO BASADO EN PUNTOS que son capturadas desde posiciones muy dispares1, lo que da lugar a imágenes con diferencias geométricas muy importantes. (ver fig. 3.4). Observación fuera del nadir Figura 3.4: Imágenes QuickBird de la misma escena adquiridas desde distintos ángulos de observación. Obsérvese las diferencias geométricas (y radiométricas) entre ambas imágenes. En la literatura se han propuesto numerosos métodos elásticos que abordan el registro de imágenes afectadas de importantes diferencias geométricas, incluyendo: funciones lineales [39] o cúbicas [40] por trozos, funciones multiquadric [29], funciones thin-plate-spline [19], funciones B-spline [59], etc. La gran mayoría de estos métodos estiman sus coeficientes a partir de un conjunto de puntos de control identificados en ambas imágenes, ya que, en términos generales, son menos costosas computacionalmente que los procedimientos basados en intensidad, que se formulan como procesos iterativos de optimización. Esta tesis contribuye con una evaluación experimental de tres de estos métodos y lo hace respecto a las siguientes cuestiones fundamentales: a) cómo modelan las diferencias relativas no-lineales inducidas por la observación de una escena no-plana desde distintos ángulos y b) cómo influye en la precisión del registro el número y distribución de los CPs utilizados en la estimación. Está claro que, si un método de registro simple (por ejemplo, funciones globales polinomiales) logra la exactitud requerida para una aplicación particular, 1Para ofrecer periodos de re-visita más cortos, estos satélites puede observar la escena desde órbitas y ángulos muy dispares (observación fuera del nadir). 3.3. COMPARATIVA EXPERIMENTAL 49 no es necesario perder tiempo seleccionando los CPs extras requeridos por una técnica mucho más potente. El objetivo de esta evaluación es, por tanto, discernir en que circunstancias se aconseja una u otra de las transformaciones analizadas y el número de correspondencias idóneo para obtener resultados aceptables. Las funciones de transformación evaluadas en esta sección son las siguientes: •una técnica global basada en funciones polinomiales (de diversos órdenes) (expresión (3.5)), •un método local basado en funciones lineales por trozos (expresión (3.6)), y •una técnica híbrida basada en la combinación de una afinidad global y funciones thin-plate-spline locales (expresión (3.7)). En esta evaluación experimental se emplean imágenes de alta-resolución del satélite QuickBird (0.6 meters/píxel) de una región con una orografía variada y adquiridas desde distintos ángulos de observación, lo que da lugar a importantes distorsiones, especialmente, en las regiones montañosas. Para cuantificar la precisión de cada método se utilizan dos medidas de consistencia: el error cuadrático medio (RMSE) (expresión (2.1)) y el error circular con 90% de confianza (CE90) (expresión (2.2)). 3.3.1 Conjuntos de datos En esta sección se describe brevemente las imágenes de prueba, así como los conjuntos de correspondencias utilizados en esta comparativa. Imágenes de prueba Se han considerado tres imágenes pancromáticas (etiquetadas como i1, i2 e i3) de la ciudad del Rincón de la Victoria (Málaga) adquiridas en diferentes fechas y desde diferentes puntos de vista. La figura 3.5(a) muestra dos de estas imágenes de prueba, en las que aparece marcadas las regiones de interés: 50 3. REGISTRO BASADO EN PUNTOS terreno casi plano (área urbana) y accidentado (área montañosa). Para definir dichas regiones se han utilizado las curvas de nivel mostradas en la figura 3.5(b). La gráfica mostrada en la parte superior de la figura 3.5(a) muestra el perfil del terreno2a lo largo de la flecha de la imagen de la derecha (desde la montaña hasta la costa) donde los niveles de altitud van desde los 300 m. hasta los 0 m. Los pares de imágenes considerados para el registro son <i1-i2>, <i1i3> y <i2-i3>. Estos pares presentan ángulos de observación relativos muy dispares que, en combinación con el relieve del terreno (el mismo para todas ellas), dan lugar a distorsiones geométricas como las observadas en la figura 3.4. Conjuntos de puntos de control Para cada par de imágenes se extraen, de manera automática mediante el procedimiento descrito en la sección 4.5.1, dos conjuntos diferentes de correspondencias: puntos de control (CP) para estimar los coeficientes de la función de transformación, y puntos independientes de control (ICP) para evaluar la precisión del registro. Para garantizar su distribución uniforme sobre las imágenes, se selecciona un punto por cada celda de una retícula rectangular dispuesta sobre las imágenes (ver fig. 3.6). Actuando sobre el tamaño de la celda se obtienen diferentes densidades, en particular, se han definido los siguientes: 50 píxeles de ancho para los conjuntos de ICPs y 100, 200 y 400 píxeles para los conjuntos de CPs, que para el caso concreto de las imágenes utilizadas, de 0.6 m./píxel de resolución, se corresponden con 30, 60, 120 y 240 m., respectivamente. De esta manera se dispone de 900, 225 y 64 CPs y un conjunto de 3600 ICPs por cada par de imágenes. Aparte de estos conjuntos de datos, que nos permiten estudiar la precisión de los métodos en regiones en la que se combinan zonas planas con otras con una geometría más compleja (geometría mixta), se consideran las regiones de interés marcadas en la figura 3.5(a), 2La información de elevación ha sido obtenida de un DEM con una resolución espacial de 20 ×20 m. proporcionado por la “Consejería de Medio Ambiente” de la “Junta de Andalucía”. Este es el DEM más preciso de la zona disponible actualmente. 3.3. COMPARATIVA EXPERIMENTAL 51 Terreno mixto Costa Tereno montañoso Terreno casi plano Costa MontañasMontañas 2 5 25 25 25 25 25 25 50 50 50 50 50 50 50 50 5 0 50 50 75 75 75 75 75 75 75 75 75 75 75 75 75 75 100 100 100 100 100 100 1 00 100 100 100 100 100 125 125 125 125 125 125 125 150 150 150 175 200 225 distance (m) distance (m) Elevation contours 3.85 3.855 3.86 3.865 3.87 x 10 5 4.064 4.0645 4.065 4.0655 x 10 6 60 80 100 120 140 160 180 200 220 240 a) b) Terrain profile distance (m) elevation (m) Figura 3.5: Imágenes de prueba utilizadas en el estudio comparativo. (a) Dos imágenes de la ciudad del Rincón de la Victoria (Málaga). Estas imágenes cubren aproximadamente un área de 4 km2e incluyen zonas con distintos perfiles de terreno (ver el perfil del terreno a lo largo de la flecha). (b) Información de elevación utilizada para definir apropiadamente las regiones de interés: áreas casi planas y de relieve accidentado. 58 3. REGISTRO BASADO EN PUNTOS Capítulo 4 Registro lineal por trozos La observación de una escena no plana desde puntos de vista distintos produce típicamente imágenes con diferencias geométricas difíciles de modelar. En el capítulo 3 se describen diferentes técnicas elásticas que tratan de abordar su registro y se realiza un análisis experimental de sus comportamientos. Este análisis concluye que, de las técnicas analizadas, sólo las transformaciones lineales por trozos (PWL) y thin-plate-splines (TPS) producen unos resultados aceptables, en términos de precisión, cuando estas diferencias locales existen. Cuando la escena observada es poliédrica (típica en escenarios de interiores y urbanos), el rendimiento de las TPS cae significativamente en favor de las PWL, que continúan ofreciendo un comportamiento adecuado. Obsérvese en la figura 4.1 cómo la transformación lineal por trozos registra con precisión la imagen de entrada, mientras que la transformación thin-plate-spline introduce distorsiones adicionales en la imagen registrada. Recuérdese que las transformaciones lineales por trozos dividen las imágenes a registrar en triángulos conjugados que son transformados individualmente mediante una transformación afín. Sin embargo, para que este método trabaje adecuadamente, cada par de triángulos correspondientes deben caer sobre las proyecciones de una superficie 3D plana, en otro caso, el registro puede generar artefactos indeseados, tales como “líneas quebradas” que reducen la calidad del registro (ver fig. 4.2). 59 60 4. REGISTRO LINEAL POR TROZOS Imagen de entrada registrada a) b) Thin-plate-spline Lineal por trozos Imagen de referencia Imagen de entrada 1 2 3 4 5 6 7 1 2 3 4 5 6 7 Pares de correspondencia Figura 4.1: Registro de un par de imágenes de una escena poliédrica mediante una transformación: (a) thin-plate-spline y (b) lineal por trozos. Obsérvese la precisión de la transformación lineal por trozos, mientras que la thin-plate-spline produce una imagen registrada con importantes distorsiones. Las implementaciones actuales del registro lineal por trozos incluidas en los paquetes científicos de procesamiento de imágenes tales como Matlab [99], Insight Segmentation and Registration Toolkit (ITK) [79], el software de Image Fusion Systems Research [51] o paquetes de teledetección tales como: ENVI/IDL [54], ERDAS [61], etc. generan las redes triangulares conjugadas a partir de un conjunto de correspondencias, utilizando para ello alguna técnica de triangulación, típicamente el método de Delaunay [93]. La triangulación de Delaunay de un conjunto de vértices maximiza el ángulo entre todas las posibles de triangulaciones de ese conjunto de vértices, pero no asegura que éstos se distribuyan sobre las imágenes de forma óptima, esto es, cayendo sobre proyecciones de trozos planos de la escena. En este capítulo y el siguiente se aborda la generación de redes triangu- 4.1. OPTIMIZACIÓN DE REDES 61 1 2 3 4 5 6 7 1 2 3 4 5 6 7 1 2 3 4 5 6 7 1 2 3 4 5 6 7 b) Red incompatible con la geometría de la escena Imagen de referencia Imagen de entrada Imagen registrada f a) Imagen de referencia Imagen de entrada Imagen registrada Red compatible con a geometría de la escena f Figura 4.2: Para que un registro lineal por trozos trabaje adecuadamente, los triángulos correspondientes deben ser proyecciones de una única superficie plana de la escena, como el triángulo {1,2,5}en (b); en otro caso, se produce el “quebramiento” de líneas y la consistencia del registro decrece, como se aprecia en (a). lares óptimas para mejorar la consistencia del registro lineal por trozos. En este primero, se realiza una revisión de las distintas técnicas de optimización de redes triangulares propuestas en la literatura, poniendo especial énfasis en aquellas propuestas en el campo del registro de imágenes; también formaliza el registro lineal por trozos. En el capítulo 5, se describe el procedimiento de registro propuesto, se analizan diversas estrategias de optimización y se presentan algunos resultados experimentales. 4.1 Optimización de redes La generación de redes triangulares óptimas es crucial en una amplia variedad de campos, como el modelado de objetos, aproximación de superficies, compresión y transmisión de imágenes, etc. Aunque el objetivo final es aproximar, tanto como sea posible, algún conjunto de datos (ya sea 2D o 3D) 62 4. REGISTRO LINEAL POR TROZOS mediante superficies triangulares por trozos, el objetivo concreto de las técnicas de optimización de redes varia con el tipo de problema. Así, en el modelado de objetos la optimización está dirigida a generar redes triangulares que representen apropiadamente la forma 3D del cuerpo o escena, utilizando para ello el menor número de triángulos posible. En procesamiento de imágenes, por el contrario, se pueden encontrar enfoques interesantes tales como la triangulación basada en datos (DDT) [28] que se utiliza, entre otros, para aproximar el contenido de una imagen (por ejemplo, reconstrucción de imágenes [112]), para reducir la cantidad de información redundante (por ejemplo, compresión de imágenes [60]) o para aproximar un conjunto discreto de muestras de una imagen mediante una superficie lineal por trozos (problema conocido como interpolación de imágenes [103]). Las técnicas de optimización de redes se pueden clasificar de acuerdo con los siguientes aspectos: •el mecanismo utilizado para modificar la red (esto es, el tipo y alcance de las acciones), •la métrica empleada para evaluar la bondad de una configuración de red dada (esto es, la función de energía o coste) y •el procedimiento para realizar el refinamiento de la red (es decir, el método de optimización). Atendiendo al tipo y alcance las acciones aplicadas para modificar la red (ver fig. 4.3), en la literatura se pueden encontrar técnicas donde: 1. se modifica la realización topológica mediante el intercambio de aristas (edge swapping) [77, 78, 82, 112], 2. se modifica la realización geométrica mediante el refinamiento de las coordenadas de los vértices [69, 89, 92, 111] y 3. se modifican simultáneamente la realización topológica y geométrica mediante la división/eliminación/intercambio de aristas y el refinamiento de las coordenadas de los vértices [27, 49, 50, 105]. 4.1. OPTIMIZACIÓN DE REDES 63 Rotura de aristas Supresión de aristas Intercambio de aristas i i illl l kkk k h h jj j Red inicial Arista interna acciones Figura 4.3: Acciones basadas en aristas empleadas para modificar la topología/geometría de una red dada. Muchos de estos métodos se han desarrollado en el campo del modelado geométrico para simplificar y/o refinar una red 3D inicialmente obtenida a partir de un conjunto muy denso de vértices proporcionado por un sensor 3D, típicamente, un escáner láser [49, 50]. Métodos similares han sido utilizados más recientemente para reconstruir escenas 3D [77, 78, 105] y en diferentes aplicaciones del enfoque DDT para procesamiento de imágenes [15, 60, 112]. A diferencia de estos métodos, que trabajan sobre redes 3D, en el registro lineal por trozos de imágenes se parte de dos redes 2D conjugadas que se modifican en un intento de maximizar la consistencia del registro. Un ejemplo de esto es el trabajo de Servais [92], que refina la localización de los vértices (la topología de la red permanece fija) para compensar el movimiento afín en secuencias de video. Aunque permite modelar distorsiones suaves, el refinamiento de las coordenadas de los vértices no proporciona la flexibilidad suficiente para modelar las importantes diferencias geométricas que se producen cuando las imágenes son adquiridas desde ángulos muy dispares. Típicamente, las técnicas de optimización se formulan como procesos de minimización (o maximización) que van desde búsquedas aleatorias [78], a procedimientos más complejos basados en el “recocido simulado” [91], mo- 64 4. REGISTRO LINEAL POR TROZOS delos estocásticos bayesianos [105], enfoques variacionales1[79], etc. Independientemente de la técnica de optimización utilizada, uno de los puntos clave en este proceso es la definición de la función de energía o coste idónea para evaluar la mejora obtenida al aplicar cierta acción sobre una configuración de red dada. Mientras que en otros campos (por ejemplo, modelado de objetos, ajuste de superficies, interpolación de imágenes, etc.) medir la calidad de una red puede realizarse sobre los puntos 3D o las muestras de la imagen disponibles, en registro de imágenes se ha de confiar en la similitud radiométrica de la imagen de referencia y la registrada. Hasta ahora, diversas métricas de similitud han sido utilizadas para este propósito: SSD [77], NCC [82], plantillas basados en diferencias de imágenes [78], etc. Como se detalló en la sección 2.3.1, ninguna de estas medidas son invariantes a las diferencias radiométricas no-lineales en las imágenes, como puede ocurrir en el caso de imágenes capturadas con diferentes cámaras, o en imágenes adquiridas con el mismo sensor, pero en condiciones de iluminación muy diferentes que pueden acarrear sombras, saturaciones, reflexiones, etc. En este capítulo se propone un proceso dirigido por una función de coste basada en la información mutua (MI). Aunque la MI ha sido utilizada cómo medida de consistencia del registro en numerosos trabajos [23, 24, 52, 57, 109], esta es la primera vez que se integra en un procedimiento de optimización de redes para el registro lineal por trozos. 4.2 Planteamiento del problema En numerosas aplicaciones, la proyección de la escena sobre el sensor puede aproximarse mediante una transformación para-perspectiva [47, 70, 110], también denominada afín o paralela. Este modelo de transformación (que se ilustra en la fig. 4.4) asume que los puntos de un objeto se proyectan primero sobre un plano de profundidad promedio paralelo al plano imagen, para seguidamente, ser transferidos al plano imagen. Los puntos se proyectan sobre el plano promedio usando rayos paralelos al eje imaginario OG, que une el 1En www.itk.org, el lector puede encontrar una amplia variedad de código abierto en el que se implementa muchas de estas técnicas: campos potenciales, cuerpos elásticos, etc. 4.2. PLANTEAMIENTO DEL PROBLEMA 65 centro óptico de la cámara Oy el centro del objeto G, preservándose de este modo el paralelismo. Esta simplificación es válida en aquellas configuraciones donde los efectos de la distorsión perspectiva son pequeños, esto es, las líneas paralelas en el espacio continúan siéndolo en las imágenes. 1 O2 O Pose 1 Pose 2 Plano promedio Plano promedio Escena 3D Plano imagen Plano imagen 1 S 2 S f > @ |Rt Haces paralelos Figura 4.4: Proyección para-perspectiva de una cámara. Cuando el campo de visión de la cámara es pequeño y el tamaño del objeto respecto a la distancia que lo separa del sensor también lo es, la proyección perspectiva (una transformación no lineal) puede ser aproximada mediante una proyección para-perspectiva, esto es, una transformación lineal. Las cámaras afines dan lugar a una importante reducción en la complejidad de numerosos problemas de visión. En particular, para el problema de registro, esta simplificación significa que tres correspondencias (no-colineales), en lugar de las cuatro necesarias en su forma general, son suficientes para estimar la transformación afín (homografía) que transfiere puntos de un trozo de imagen a otro [47]. Así, cuando se realiza un mapeo afín entre dos regiones triangulares conjugadas de ambas imágenes, estas deben alinearse perfectamente (esto es, la consistencia del registro es máxima); en caso contrario, las regiones triangulares son proyecciones de una superficie no plana. En la figura 4.5 se muestra el resultado de registrar, mediante una transformación lineal por trozos (PWL), imágenes de una escena plana afectadas de deformación perspectiva y afín, respectivamente. Obsérvese como, mien- 66 4. REGISTRO LINEAL POR TROZOS tras la primera presenta diversos artefactos (indicados con flechas en la imagen) que deterioran la calidad del registro, la segunda se registra con absoluta precisión. Como se indica arriba, suponer que la transformación que alinea con precisión dos regiones triangulares conjugadas es una afinidad sólo es válido en aquellas configuraciones cámara-escena que dan lugar a imágenes en las que la distorsión perspectiva es inapreciable. 1 2 3 4 1 2 3 41 2 3 4 Imagen de entrada (dist. perspectiva) Imagen de entrada (dist. afín) a) b) Imagen de referencia Imagen registrada Imagen registrada PWL PWL Figura 4.5: Registro lineal por trozos (PWL) de imágenes afectadas de distorsión (a) proyectiva y (b) afín. Los triángulos empleados para la transformación PWL son el {1,2,3}y el {1,3,4}. Obsérvese en (a) los artefactos (indicados con flechas) que se producen al aplicar este tipo de transformación a imágenes con distorsiones perspectivas. Resumiendo, si se logra modificar la topología y/o geometría de la red triangular utilizada en un proceso genérico de registro lineal por trozos mediante acciones dirigidas a maximizar el número de elementos triangulares que se registran con precisión, esto es, que caen sobre proyecciones de superficies 3D planas: a) se eliminan configuraciones indeseadas como las ilustradas en la figura 4.2 y b) se mejora la precisión global del registro. En esta tesis se explota la relación entre las redes triangulares 2D que se utilizan en el registro PWL y su correspondiente malla triangular 3D que modela la escena observada. En base a esto se formaliza, no sólo la generación 4.3. FORMALIZACIÓN DEL PROBLEMA 67 de redes triangulares 2D óptimas (desde el punto de vista de la consistencia del registro), sino también la reconstrucción a partir éstas de la superficie 3D observada, yendo de este modo, más allá de un puro enfoque 2D del registro. A continuación se introducen los conceptos y definiciones utilizados en la formulación del registro lineal por trozos y del método de optimización propuesto. 4.3 Formalización del problema Supóngase una escena cualquiera cuya geometría es aproximada mediante una superficie lineal por trozos, esto es, un conjunto de triángulos conectados unos con otros conformando una red triangular en R3(ver fig. 4.6). Aunque existen diversas maneras de formalizar mallas triangulares [60, 77, 78, 105], en esta tesis se emplea el complejo simplicial, una estructura matemática empleada en topología algebraica para definir espacios topológicos, sus invariantes, las relaciones entre símplices, etc. El complejo simplicial aporta ciertas ventajas respecto a otras formalizaciones, entre otras: •Modela con una misma estructura puntos, segmentos, triángulos y conjuntos de estos, así como sus relaciones topológicas (conexión, adyacencia, etc.) y sus invariantes topológicos (huecos, túneles, etc.). •Separa las realizaciones topológica y geométrica de la estructura, de modo que se puede actuar sobre la realización topológica con independencia su realización geométrica, y viceversa. •Facilita el almacenamiento, la actualización y el acceso eficiente a los diferentes elementos de la estructura del complejo mediante matrices de adyacencia y conexión. Gracias a estas características, los complejos simpliciales pueden modelar estructuras topológicas con independencia del espacio en el que se encuentren definidas: el espacio o el plano. Aunque, a primera vista, esto puede parecer de un interés marginal en el contexto del registro de imágenes, mantener la 74 4. REGISTRO LINEAL POR TROZOS referencia; finalmente, 3. se toma el complejo simplicial Kde M2 iy se genera la red M2 j= (K, V 2 j), mediante la sustitución de los puntos en V2 ipor sus correspondientes en V2 j. Aunque la realización geométrica φV2 igenerada mediante el procedimiento de triangulación siempre es un encaje3, su correspondiente φV2 jpuede no serlo, lo que se traduce en la aparición de ciertas inconsistencias denominadas patch reversals. Los patch reversals se producen típicamente cuando la escena es observada desde puntos muy diferentes (ver fig. 4.10(a)), y tienen graves implicaciones: 1. φV2 i(|K|)yφV2 j(|K|)no se pueden registrar, ya que, la transformación fno es un homeomorfismo, y 2. φV3(|K|)no se puede reconstruir, ya que, aunque Pies un homeomorfismo, Pjno lo es, luego no se verifica la expresión (4.2). Para eliminar estas inconsistencias, se revisa la topología inicial aplicando las siguientes modificaciones: 1. si la arista que produce el patch reversal es frontera, se elimina el triángulo que la contiene, como se muestra en la figura 4.10(b). 2. si la arista que lo produce es compartida, se intercambia la arista, como se ilustra en la figura 4.10(c). Este proceso se repite iterativamente hasta que se eliminan la totalidad de las inconsistencias. El procedimiento de generación de redes triangulares descrito tiene dos inconvenientes (ver fig. 4.11): 1. no garantiza que los triángulos conjugados caigan sobre proyecciones de superficies planas. Por ejemplo, los triángulos {1,7,3}y{10,12,15} de K. 3Dados dos espacios topológicos XeY,f:X7→ Yes un encaje (embedding), sí y sólo sí, fes un homeomorfismo entre Xy su imagen f(X). 4.4. GENERACIÓN DE LAS REDES 75  2 1 VK I  2 2 VK I a) b) Eliminando patch reversals en aristas compartidas Eliminando patch reversals en aristas frontera Patch reversals c) 1,1 x 1,2 x1,9 x 1,7 x1,8 x 1,3 x 1,10 x 1,16 x 1,15 x 1,4 x1,11 x1,12 x 1,1 x 1,2 x1,9 x 1,7 x1,8 x 1,3 x 1,10 x 1,16 x 1,15 x 1,4 x1,11 x1,12 x 1,1 x 1,2 x1,9 x 1,7 x1,8 x 1,3 x 1,10 x 1,16 x 1,15 x 1,4 x1,11 x 1,12 x 2,7 x 2,8 x 2,3 x 2,1 x 2,11 x2,12 x 2,10 x 2,9 x 2,2 x 2,15 x2,16 x 2,4 x 2,15 x 2,7 x 2,8 x 2,3 x 2,1 x 2,11 x2,12 x 2,10 x 2,9 x 2,2 x 2,15 x2,16 x 2,4 x 2,7 x 2,8 x 2,3 x 2,1 x 2,11 x2,12 x 2,10 x 2,9 x 2,2 x 2,15 x2,16 x 2,4 x Figura 4.10: (a) Patch reversals (PR) producidos por la existencia de caras ocultas cuando la escena (descrita con una red triangular) es observada desde posiciones diferentes. Para eliminar estas inconsistencias se proponen los siguientes cambios topológicos en la redes iniciales: (b) Si la arista es frontera: se elimina el triángulo que la contiene, esto es, {10,12,16},{2,9,10}. (c) Si la arista es compartida: se intercambia, esto es, {9,15}por {10,11}. 2. generan estructuras topológicas regulares 1-conectadas, esto es, todos los símplices máximos de Kson triángulos y todos los triángulos tienen al menos un adyacente, dando lugar a nuevos triángulos sin correspondencia real con la escena observada. Por ejemplo, los triángulos {2,9,4} y{8,11,15}de K. En relación a la optimización de la red, en esta tesis se propone maximizar el número de triángulos que caen sobre proyecciones de superficies planas de la imagen, para lo cual se modifica la topología y/o geometría de las redes M2 i= (K, V 2 i)yM2 j= (K, V 2 j). Las acciones implementadas no 76 4. REGISTRO LINEAL POR TROZOS 1,1 x 1,2 x1,9 x 1,7 x1,8 x 1,3 x 1,10 x 1,16 x 1,15 x 1,4 x1,11 x1,12 x 1 ITriángulos que no caen sobre proyecciones de superficies planas Triángulos sin correspondencia real con la escena Figura 4.11: Triángulos no consistentes con las proyecciones de la escena observada: (a) triángulos que no caen sobre proyecciones de superficies planas y (b) triángulos sin correspondencia real con la escena observada. contemplan la posibilidad de “desconectar” triángulos para generar símplices de dimensión menor (aristas y puntos), así como tampoco, eliminar triángulos sin correspondencia real con la escena. Este tipo de acciones se dejan como trabajo futuro. 4.5 Selección del conjunto de correspondencias Son muchas las aplicaciones de visión por computador que requieren de técnicas robustas de búsqueda de correspondencias para extraer información relevante de dos o más imágenes, ya sea para registrarlas con precisión, para recuperar la disposición relativa del sensor o sensores en el momento de la captura, para reconstruir tridimensionalmente una escena, etc. Una característica es un elemento distintivo o representativo detectado manual, o preferiblemente, de forma automática en la imagen, por ejemplo: regiones, bordes, contornos, segmentos, intersecciones de rectas, esquinas, etc. Las características poseen propiedades que las describen (esto es, descriptores), en forma vectorial d = (d1, . . . , dm)>∈Rm, y vectores de coordenadas, x = (x, y)>∈R2, que las localizan en la imagen. Una correspondencia (o características conjugadas) es el conjunto formado por los vectores de coordenadas que localizan una misma característica 4.5. SELECCIÓN DEL CONJUNTO DE CORRESPONDENCIAS 77 en las imágenes que intervienen en el análisis. Así, dadas dos imágenes Iie Ij, la correspondencia k-ésima seleccionada en ambas imágenes vendría dada por la tupla (xi,k,xj,k). Esta sección se centra exclusivamente en las técnicas de detección y emparejamiento robusto de puntos de interés (esquinas). Para una revisión exhaustiva de los diferentes métodos propuestos en la literatura, el lector puede dirigirse a los siguientes trabajos [20, 42, 113]. 4.5.1 Detección y emparejamiento de esquinas Una esquina es un punto o localización de la imagen que tiene, para una escala determinada, grandes gradientes en todas las direcciones. La búsqueda de esquinas ha recibido, y continúa recibiendo, gran atención en la literatura. A continuación se revisan brevemente algunos de los trabajos más representativos. En 1982, Kitchen y Rosenfeld [58] propusieron explotar las derivadas parciales de segundo orden para detectar esquinas. A pesar de los prometedores resultados, este detector de esquinas demostró ser excesivamente sensible al ruido. Con objeto de minimizar esta limitación, Förstner y Gulch proponen en [35] un detector sensiblemente más robusto basado en la primera derivada, aunque computacionalmente muy costoso. En 1980, Moravec presenta un método basado en las derivadas parciales de primer orden e introduce el concepto de la matriz de auto-correlación. En 1988, Harris y Stephens extienden el trabajo de Moravec en [46], proponiendo una función detección de esquinas muy eficiente. En [94], Shi y Tomasi introducen una nueva función de detección que explota igualmente la matriz de auto-correlación para extraer esquinas fáciles de seguir en secuencias de imágenes. Merece especial atención el trabajo Lowe [63], donde se propone una novedosa técnica para detección de esquinas basada en la búsqueda de extremos (máximos y mínimos locales) en un espacio de escalas construido a partir de diferencias de Gaussinas (DoG). Este enfoque produce, a diferencia de otros detectores como los de Harris o Shi y Tomasi, esquinas en diferentes escalas. 78 4. REGISTRO LINEAL POR TROZOS Descriptores de esquinas Los vectores de propiedades (o descriptores) describen las esquinas o, más formalmente, el entorno de éstas con un único objetivo: lograr que dicho vector sea lo suficientemente distintivo como poder diferenciarlo inequívocamente del resto de esquinas. El procedimiento más sencillo consiste en utilizar los niveles de gris de un entorno de la esquina como descriptor. Esta representación, aunque muy utilizada tiene sus limitaciones: es sensible a los cambios en el brillo y/o contraste, la presencia de ruido y las distorsiones geométricas. A esto hay que sumar problemas asociados específicamente al tamaño del descriptor, un entorno grande es mucho más distintivo pero más sensible a los problemas descritos con anterioridad, y viceversa. Aunque en la literatura se pueden encontrar técnicas que tratan de compensar estas limitaciones mediante procedimientos de normalización radiométrica, o que consideran que los vecindarios pueden verse afectados por un conjunto limitado de posibles diferencias geométricas, la utilización de este tipo de técnicas suele restringirse a aplicaciones en los que las diferencias entre imágenes no son excesivamente importantes, como por ejemplo, para el seguimiento de características en secuencias de imágenes [65]. Numerosos autores han tratado de suplir estas limitaciones dotando a las esquinas extraídas con detectores tradicionales como los de Harris, Shi y Tomasi o Lowe, con descriptores adicionales. Este es el caso de Schmid y Mohr [90], que emplea una versión extendida del detector de Harris para la búsqueda de contenidos en imágenes. Otros autores como Baumberg [16], Mikolajczyk y Schmid [71], Schaffalitzky y Zisserman [88] o Brown y Lowe [21], con su famoso SIFT, mejoran la estabilidad de los descriptores al dotarlos de invarianza a un mayor número de diferencias geométricas, incluidas las transformaciones afines. En la literatura también se pueden encontrar variantes de estos procedimientos, como el denominado PCA-SIFT, propuesto por Ke y Sukthankar [56], que realiza una análisis PCA para generar descriptores equivalentes a los SIFT pero de menor dimensión, o técnicas relativamente recientes como 4.5. SELECCIÓN DEL CONJUNTO DE CORRESPONDENCIAS 79 el SURF [17], un descriptor invariante a la rotación y la escala, basado en el uso de imágenes integrales [100] y en la matriz Hessiana de la imagen. En [74], el lector puede encontrar una comparativa de distintos detectores de esquinas y descriptores evaluada sobre imágenes de objetos tridimensionales. Emparejamiento de esquinas Una vez localizadas las esquinas en ambas imágenes se procede a la generación de correspondencias mediante su emparejamiento. En general, dada una esquina localizada en una imagen, su emparejamiento consiste en localizar, en el conjunto de esquinas detectadas en la otra imagen, aquella con el descriptor más “parecido”, esto es, el par de esquinas conjugadas cuyos descriptores asociados minimizan, por ejemplo, la distancia Euclídea. Ésta búsqueda se puede abordar de forma exhaustiva (esto es, todos con todos) o restringirla al conjunto de esquinas detectadas en un determinado área de la segunda imagen, con el consiguiente ahorro computacional. Por ejemplo, en [110], Xu and Zhang restringen la búsqueda de correspondencias dentro de un determinado área denominada ventana de búsqueda, otros explotan la restricción epipolar, para restringir la comparación a aquellas esquinas que caen sobre sus correspondientes líneas epipolares [75, 76]. Naturalmente, este enfoque requiere un conocimiento a priori de las posibles diferencias geométricas que presentan las imágenes o del desplazamiento relativo de los sensores en el momento de la adquisición (por ejemplo, en estéreo). 4.5.2 Selección robusta de correspondencias En la práctica, el empleo de una u otra técnica de detección de esquinas dependerá fundamentalmente de las diferencias geométricas que presenten las imágenes que intervienen en el análisis y del tipo de aplicación. Así, el detector de Harris es sensible a los cambios de escala, pero tiene un coste computacional pequeño, lo que lo hace especialmente útil en aplicaciones que requieran de rapidez de computo y las imágenes no presenten importantes cambios de escala. Por otro lado, el detector de Lowe tiene un elevado 80 4. REGISTRO LINEAL POR TROZOS coste computacional (especialmente, la etapa de búsqueda de extremos) pero proporciona características altamente estables a los cambios de escala y otras deformaciones geométricas afines, lo que lo hace especialmente apropiado en aplicaciones en las que las imágenes han sido adquiridas desde puntos de vista muy dispares, y que presentan obviamente, importantes diferencias geométricas. En esta tesis, al igual que otros trabajos encontrados en la literatura [32, 76], se explotan los beneficios de ambas técnicas en pos de obtener un mejor rendimiento computacional y aumentar el porcentaje de correspondencias válidas en la etapa de emparejamiento. Más concretamente, el proceso de búsqueda y emparejamiento se aborda del siguiente modo: 1. se localizan las esquinas con el detector Harris4, 2. se calculan sus correspondientes descriptores SIFT, y finalmente, 3. se establecen los emparejamientos utilizando la distancia Euclídea. Además, para dotar de robustez al procedimiento, se recupera la geometría epipolar afín intrínseca a ambas imágenes utilizando para su estimación el conjunto de correspondencias establecidas anteriormente. Esto proporciona tres beneficios: 1) detectar los posibles pares espurios (outliers), es decir, pares no consistentes con la restricción epipolar; 2) refinar la localización de los pares válidos (inliers) y 3) obtener un estimación robusta de la matriz fundamental afín, FA, que será utilizada en el siguiente capítulo para dirigir el proceso de división de aristas en la red triangular. Para estimar la matriz fundamental afín se emplea un procedimiento basado en el algoritmo RANdom SAmple Consensus (RANSAC) [31]. Esta técnica explota la redundancia de información de un conjunto de muestras para proporcionar una estimación robusta del modelo que se ajusta con precisión a la mayoría de ellas. En el apéndice B se discute en detalle la estimación robusta de parámetros y se describen éste y otros métodos propuestos en la literatura. 4El lector puede encontrar en el apéndice D una descripción detallada de los detectores de Harris y Lowe. 4.5. SELECCIÓN DEL CONJUNTO DE CORRESPONDENCIAS 81 En este caso particular, el modelo es la matriz fundamental afín, FA, y la función de error que mide la consistencia de una correspondencia dada, (xi,xj), es el error geométrico de primer order (o distancia de Sampson), esto es, la distancia de una esquina a su correspondiente línea epipolar, que se expresa como: d(xj,FAxi) = x> jFAxi2 (FAxi)2 1+ (FAxi)2 2+F> Axj2 1+F> Axj2 2 (4.3) donde los términos de la forma (v)2 krepresentan el cuadrado de la coordenada k-ésima del vector v. Para ganar en precisión, en su lugar, el procedimiento utilizado en esta tesis emplea el error epipolar simétrico, que considera la distancia del par (xi,xj)a sus respectivas líneas epipolares. El error epipolar simétrico se deriva de (4.3) y se define como: d(xj,FAxi)2+dxi,F> Axj2= =  x> jFAxi2 (FAxi)2 1+ (FAxi)2 2 +x> jFAxi2 F> Axj2 1+F> Axj2 2   (4.4) El paso final de algoritmo RANSAC es la re-estimación del modelo considerando sólo las correspondencias válidas (inliers). Esta re-estimación se realiza mediante un proceso de minimización del que se deriva la estimación de máxima verosimilitud (MLE) de la matriz fundamental afín y una estimación de las coordenadas óptimas de las correspondencias. La MLE asume que los errores cometidos en la identificación de las esquinas siguen una distribución Gaussiana. En tal caso, la MLE de la FAes aquella que minimiza 82 4. REGISTRO LINEAL POR TROZOS el error de re-proyección. Formalmente, se formula como: min {FA,ˆxi,k,ˆxj,k} n X k=1 d(xi,k,ˆxi,k)2+d(xj,k,ˆxj,k)2(4.5) donde (xi,k,xj,k)son las correspondencias iniciales y (ˆxi,k,ˆxj,k)las refinadas, siendo nel número de correspondencias válidas. El lector puede encontrar en el texto de Hartley y Zisserman [47] (pag. 349) el desarrollo de teórico de este problema de minimización. La figura 4.12 ilustra el proceso de selección robusta de correspondencias descrito en esta sección. 4.6 Conclusiones En este capítulo se ilustran las ventajas de emplear funciones lineales por trozos para el registro de imágenes conjugadas de una escena poliédrica. Este procedimiento divide las imágenes en regiones triangulares conjugadas que se registran individualmente mediante afinidades. Sin embargo, para que este método trabaje adecuadamente, las regiones triangulares deben caer sobre proyecciones de superficies planas de la escena, aspecto que no abordan las diferentes implementaciones de esta técnica que se pueden encontrar en el software disponible en el mercado. En este capítulo se ha analizado esta problemática y se propone una formalización del registro lineal por trozos basado en complejos simpliciales que abarca, no sólo la generación de redes triangulares 2D óptimas (desde el punto de vista de la consistencia del registro), sino también la reconstrucción a partir éstas de la superficie 3D observada. También se ha propuesto un procedimiento automatizado para la generación de redes triangulares conjugadas libres de patch reversals. La descripción del procedimiento de optimización propuesto se aborda en el siguiente capítulo. 4.6. CONCLUSIONES 83 12 3 45 6 78 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 12 3 45 6 78 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 2 3 45 6 78 9 10 11 12 13 14 15 17 18 19 21 22 23 24 1 16 20 25 2 3 45 6 78 9 10 11 12 13 14 15 17 18 19 21 22 23 24 1 16 20 25 12 3 45 6 8 9 10 11 12 13 14 15 17 18 19 20 22 23 24 25 12 3 45 6 78 9 10 11 12 13 14 15 17 18 19 20 21 22 23 24 25 17 16 20 25 Imagen de referencia Imagen de entrada a) b) c) 21 7 Figura 4.12: Selección robusta de correspondencias: a) detección y emparejamiento, b) detección de los outliers y c) refinamiento de las coordenadas de los inliers (los círculos denotan las coordenadas refinadas). Los fragmentos ampliados de las imágenes muestran, en detalle, las actuaciones (b) y (c) sobre las correspondencias 16 y 25, respectivamente. 90 5. OPTIMIZACIÓN DE REDES TRIANGULARES Definición 5.5 (collapse) Dado un vértice v∈e, se define la acción eliminar (collapse) arista ˆ K= collapse (e, v, K)(5.5) como el conjunto de operaciones sobre K(ver fig. 5.4): a) eliminar el conjunto de símplices star({i}, K)∪star({j}, K), b) insertar el vértice v, c) insertar una arista {p, v}por cada vértice {p} ∈ bound (e, K)e d) insertar un triángulo {p, q, v}por cada arista {p, q} ∈ bound (e, K). a) Acción “eliminar arista”: b) d)c) ^ ` i ^ ` i ^ ` k ^ ` l ^ ` k ^ ` l ^ ` k ^ ` l ^ ` i ^ ` k ^ ` l ^ ` ^ `   collapse , , ,ij j K Figura 5.4: Pasos de la acción collapse(e, {j}, K): (a) eliminar star({i}, K)∪ star({j}, K); (b) insertar el vértice {i}; (c) insertar una arista {p, i}por cada vértice {p} ∈ bound (e, K)y (d) insertar un triángulo {p, q, i}por cada arista {p, q} ∈ bound (e, K). Nótese que para la arista ese pueden generar dos acciones, una por cada vértice de e. Se podría haber optado, sin embargo, por una única acción que 5.2. APLICACIÓN DE LAS ACCIONES 91 eliminase simultáneamente los dos vértices, lo que, aunque con menor coste computacional, puede obviar interesantes configuraciones de la red. 5.2 Aplicación de las acciones Es importante subrayar que, al igual que ocurre en la generación automática de las redes conjugadas (ver sec. 4.4), la aplicación de una acción puede producir un patch reversal. En este caso, la función de transformación fya no es un homeomorfismo y, por tanto, las realizaciones geométricas de Kni se pueden registrar ni reconstruir. Para evitar este tipo de circunstancias, antes de aplicar cualquier acción, conviene evaluar si ésta produce este tipo de inconsistencias. Dadas dos redes M2 1= (K, V 2 1)yM2 2= (K, V 2 2), una arista compartida e={i, j} ∈ Ky el complejo simplicial s= star(e, K), se dice que: 1. el swap(e, K)es aplicable, si y sólo si, la arista opuesta ¯eno produce un patch reversal; 2. el split(e, v, K)es aplicable, si y sólo si, φV2 1(|v|)∈φV2 1(|s|)yφV2 2(|v|)∈ φV2 2(|s|), es decir, la realización geométrica del nuevo vértice φV2 1(|v|) está dentro de la región cuadrangular definida por la realización geométrica φV2 1(|s|)(análogamente para φV2 2(|v|)); 3. el collapse(e, v, K)es aplicable, si y sólo si, vno es frontera en K. A diferencia de la acción intercambiar arista, donde la generación de un patch reversal impide su aplicación, las inconsistencias que se producen en φV2 1(|ˆ K|)oφV2 2(|ˆ K|)al eliminar aristas se pueden resolver de forma análoga a la descrita en la sección 4.4. La figura 5.5 ilustra este proceso. Especial interés recibe la división de aristas. A diferencia de las acciones intercambiar yeliminar arista que se aplican directamente, la acción dividir arista necesita determinar el punto de división idóneo, de modo que, dependiendo de la localización de éste, la acción puede conllevar una mejora de la consistencia o nó. A esto hay que sumar su mayor coste computacional, 92 5. OPTIMIZACIÓN DE REDES TRIANGULARES Patch reversals  2 1 VK I  2 3 VK I a) b) c) Eliminando patch reversals ^ ` ^ `  collapse 4,8 , 8 ,K 1,1 x 1,2 x1,9 x 1,7 x1,8 x 1,3 x 1,10 x 1,16 x 1,15 x 1,4 x1,11 x1,12 x 1,1 x 1,2 x1,9 x 1,7 x 1,3 x 1,10 x 1,16 x 1,15 x 1,4 x1,11 x1,12 x 1,1 x 1,2 x1,9 x 1,7 x 1,3 x 1,10 x 1,16 x 1,15 x 1,4 x1,11 x1,12 x 3,7 x3,8 x 3,1 x 3,3 x3,4 x3,11 x3,12 x 3,10 x 3,9 x 3,2 x 3,15 x3,16 x 3,7 x 3,1 x 3,3 x3,4 x3,11 x3,12 x 3,10 x 3,9 x 3,2 x 3,15 x3,16 x 3,7 x 3,1 x 3,3 x3,4 x3,11 x3,12 x 3,10 x 3,9 x 3,2 x 3,15 x3,16 x Figura 5.5: (a) Patch reversal (PR) producido por la eliminación de la arista {4,8} ∈ K. Para eliminar estas inconsistencias se procede forma análoga a lo descrito en la sección 4.4, esto es, se intercambia la arista compartida por los triángulos afectados, esto es, {4,12}por su opuesta {11,15}. ya que al no disponerse de información sobre la geometría de la escena, la división debe abordarse en el dominio de la imagen. Por ejemplo, la figura 5.6 muestra una realización geométrica típica: diversas aristas que conectan puntos situados en proyecciones de superficies planas diferentes, por ejemplo la arista e={4,10}. Observando la configuración de la red se aprecia claramente que si se intercambiara o eliminara ela consistencia no mejoraría mucho. Por el contrario, si se lograra dividir een el borde del poliedro (ver flecha en la figura), las realizaciones geométricas de los nuevos triángulos, {9,10, p}y{10,12, p}, caerían sobre proyecciones de superficies planas, obteniéndose una clara mejora de la consistencia en la región afectada y, por ende, en la consistencia global. El procedimiento de división propuesto en este trabajo explota la infor- 5.2. APLICACIÓN DE LAS ACCIONES 93 Acción “dividir arista” 1 I3 I 1 I3 I 3 p ^ ` ^ `   split 4,10 , , p K 1 p 1,1 x 1,2 x1,9 x 1,7 x1,8 x 1,3 x 1,10 x 1,16 x 1,15 x 1,4 x1,12 x 1,1 x 1,2 x1,9 x 1,7 x1,8 x 1,3 x 1,10 x 1,16 x 1,15 x 1,4 x1,12 x 3,7 x3,8 x 3,1 x 3,3 x3,4 x 3,12 x 3,10 x 3,9 x 3,2 x 3,15 x3,16 x 3,7 x3,8 x 3,1 x 3,3 x3,4 x 3,12 x 3,10 x 3,9 x 3,2 x 3,15 x3,16 x Figura 5.6: Inserción óptima de nuevos vértices en la red triangular. Supóngase, por ejemplo, que se desea dividir la arista e={4,10}. Una posición idónea para dividirla sería la intersección de ésta con el borde del poliedro, ya que las realizaciones geométricas de dos de los nuevos triángulos, en particular {9,10, p} y{10,12, p}, caerían sobre una de las caras del poliedro, obteniéndose una clara mejora de la consistencia en la región de influencia de e. mación de bordes de las imágenes (a partir de la cual se extraen segmentos) y la geometría epipolar afín intrínseca a ambas vistas para determinar la localización idónea del punto de división. En aras de una mayor claridad, la descripción del procedimiento se apoya en el caso ilustrado en las figuras 5.6 y 5.7. Considérense las imágenes I1yI3, las redes M2 1= (K, V 2 1)yM2 3= (K, V 2 3) y la arista compartida e={4,10} ∈ Kde la figura 5.6. Para la división de ese procede del siguiente modo: 1. Extracción de segmentos en ambas imágenes: se extraen bordes en ambas imágenes y se aproximan segmentos. Los segmentos de longitud inferior a un umbral dado se descartan. En la implementación de esta etapa se pueden utilizar el detector de Canny [22] y el método de aproximación propuesto por Kovesi en [83]. 94 5. OPTIMIZACIÓN DE REDES TRIANGULARES 2. Extracción de los puntos conjugados de división (ver fig. 5.7(b)): se determinan los puntos de intersección de φV2 1(|e|)con los segmentos de I1: {p1,i = (x1,i, y1,i)>, i = 1, . . . , n} A continuación, se determinan los puntos de intersección de las líneas epipolares l3,i = FA×p1,i con los segmentos de I3: {p3,j = (x3,j, y3,j)>, j = 1, . . . , m} descartándose los p3,j /∈φV2 3(|s|), con s= star(e, K), es decir, los puntos que no están en el interior de la región cuadrangular definida por la realización geométrica de s. Para el resto se calculan sus descriptores SIFT y se emparejan utilizando la distancia Euclídea. Finalmente, los pares cuya distancia entre descriptores esté por encima de un umbral dado también se descartan. El procedimiento descrito puede generar varias correspondencias por arista, en cuyo caso, se selecciona la que presenta una mayor disparidad, ya que se ha comprobado experimentalmente que es la correspondencia con mayor probabilidad de estar situada en el borde que “conecta” proyecciones de diferentes superficies (ver fig. 5.7(c)). Se pueden considerar otras estrategias de selección como, por ejemplo, generar tantas acciones de división como pares hayan sido localizados (con el consiguiente incremento en el coste computacional). Si no se obtiene ninguna correspondencia válida, la acción se marca como no aplicable. La estimación de la matriz fundamental afín, la extracción de segmentos, así como la identificación de los puntos de división se realiza antes de iniciar el proceso de optimización. Una vez iniciado éste, cada vez que se aplica una acción se extraen puntos de división en las nuevas aristas y se eliminan los pares correspondientes a las aristas que ya no existen. El coste computacional de esta acción es superior a las acciones eliminar eintercambiar arista debido, fundamentalmente, al cálculo de los SIFTs. Para ganar en eficiencia, 5.2. APLICACIÓN DE LAS ACCIONES 95 Acción “dividir arista” 1,10 x 1,4 x Segmentos 1,10 x 1,4 x Intersecciones b) c) d) 3,4 x 3,10 x Intersecciones 1,9 x 1,10 x 1,12 x 1,4 x 1 p 3,9 x 3,12 x 3,4 x 3,10 x 3 p 3,4 x 3,10 x a) 3 I 1 I Bordes de Bordes de Segmentos 1 I3 I Figura 5.7: Proceso de inserción de nuevos vértices en las redes triangulares: (a) detección de bordes y extracción de segmentos, (b) extracción de puntos candidatos, (c) emparejamiento y selección de la “mejor” correspondencia (aquella con la disparidad máxima) y (d) configuración de la red una vez finalizada la división. 96 5. OPTIMIZACIÓN DE REDES TRIANGULARES los descriptores SIFT se pueden obtener con la librería SiftGPU [95], una implementación que explota las prestaciones de las modernas unidades de procesamiento gráfico (GPU). En ocasiones, la mejora obtenida al dividir una arista es idéntica1a la obtenida con una operación de intercambio (ver fig. 5.8), en cuyo caso, siempre se aplica la operación de intercambio, evitándose la inserción de vértices superfluos y que ralentizarían el proceso de optimización. 5.3 Descripción del método El objetivo del método presentado en este capítulo es mejorar la precisión del registro lineal por trozos. Para este propósito se modifica iterativamente la topología/geometría de las redes triangulares (conjugadas) iniciales mediante la aplicación de distintas acciones, que se han descrito en la sección 5.1.2. Este proceso se formula como una búsqueda greedy [25] que, en cada iteración, toma una arista concreta, verifica el conjunto de acciones aplicables y mide la mejora en la consistencia de cada una ellas, seleccionando finalmente, la acción que produce el mayor incremento. Para cuantificar la mejora introducida se emplea una función de coste local que mide cuanto mejora el registro en la zona afectada, antes y después de aplicar cada acción. 5.3.1 La función de coste local En este trabajo se explota la robustez de la MI, más concretamente del ECC, para cuantificar la mejora del registro introducida por una acción. Dadas dos imágenes I1yI2; dos redes M2 1= (K, V 2 1)yM2 2= (K, V 2 2); una arista compartida e={i, j} ∈ Ky una acción act sobre e, se definen los siguientes elementos: Definición 5.6 La región de influencia rde la acción act(e, K)en M2 1 (análogamente para M2 2) es el conjunto de coordenadas r⊂N2que se ven 1Puede haber pequeñas diferencias debidas a la interpolación de los niveles de intensidad. 5.3. DESCRIPCIÓN DEL MÉTODO 97 Acción “dividir arista” / “intercambiar arista” ^ ` ^ `   split 10,15 , , p K ^ `   swap 10,15 ,K 1 p2 p 1,1 x 1,2 x1,9 x 1,7 x1,8 x 1,3 x 1,10 x 1,16 x 1,15 x 1,4 x1,11 x1,12 x 1 I2 I 1 I2 I 1 I2 I 1,1 x 1,2 x1,9 x 1,7 x1,8 x 1,3 x 1,10 x 1,16 x 1,15 x 1,4 x1,11 x1,12 x 1,1 x 1,2 x1,9 x 1,7 x1,8 x 1,3 x 1,10 x 1,16 x 1,15 x 1,4 x1,11 x1,12 x 2,7 x 2,8 x 2,3 x 2,1 x 2,11 x2,12 x 2,10 x 2,9 x 2,2 x 2,15 x2,16 x 2,4 x 2,7 x 2,8 x 2,3 x 2,1 x 2,11 x2,12 x 2,10 x 2,9 x 2,2 x 2,15 x2,16 x 2,4 x 2,7 x 2,8 x 2,3 x 2,1 x 2,11 x2,12 x 2,10 x 2,9 x 2,2 x 2,15 x2,16 x 2,4 x Figura 5.8: Comparación entre “intercambiar arista” y “dividir arista”. En ocasiones, una operación de división e intercambio proporcionan la misma mejora en la consistencia, obsérvese por ejemplo la arista e={10,15}, en cuyo caso la acción de división se marca como no aplicable. 98 5. OPTIMIZACIÓN DE REDES TRIANGULARES afectadas por su aplicación, tal que r⊂φV2 1(|s|). La región de influencia viene definida por el complejo simplicial sy depende del tipo de acción, en particular: •s= star(e, K), para las acciones swap y split, y •s= (star({i}, K)∪star({j}, K)), para la acción collapse. Definición 5.7 La mejora ∆ωintroducida por la acción act(e, K)es la variación en la consistencia del registro de su región de influencia antes (ω1) y después (ω2) de aplicarla, esto es: ∆ω(e, act) = ω2(ˆs)−ω1(s)(5.6) donde ω1(s) = ECC (I1(r),I2(fs(r))) , ω2(ˆs) = ECC (I1(r),I2(fˆs(r))) ,(5.7) I1(r)representa los píxeles de I1dados por la región de influencia r,I2(fs(r)) yI2(fˆs(r)) los píxeles de I2dados por la región de influencia transformada de acuerdo con los dos posibles complejos syˆs=act(e, s), respectivamente. La acción se acepta, si y sólo si, ∆ω(e, act)≥δ > 0, esto es, cuando la acción da lugar a una mejora por encima de un umbral δdado. Este umbral debe ser pequeño, pero suficiente, para impedir que se apliquen acciones cuyas regiones de influencia caen sobre proyecciones de superficies planas y que, debido a la interpolación de los niveles de intensidad, dan lugar a pequeñas variaciones positivas de ω. Se ha determinado experimentalmente que un δ≈0,01 resuelve esta problemática. Como se discutió en la sección 2.4, la estimación de la MI y de sus variantes normalizadas (como el ECC), depende de la estimación de la distribución de probabilidad conjunta de las intensidades de las imágenes. Puesto que la fiabilidad de esta estimación se ve comprometida cuando el número de píxeles es pequeño, aquellas acciones cuyas zonas de influencia tienen un tamaño (en 5.3. DESCRIPCIÓN DEL MÉTODO 99 píxeles) inferior a un umbral dado son descartadas (en nuestros experimentos, 1000 píxeles). Por cuestiones de eficiencia, se ha utilizado el método de estimación paramétrico con 32 acumuladores. 5.3.2 Método de optimización El método de optimización se formula como una búsqueda greedy [96], que comienza con dos redes triangulares M2 1= (K, V 2 1)yM2 2= (K, V 2 2)generadas como se detalla en la sección 4.4, y finaliza con sus correspondientes ˆ M2 1yˆ M2 2 optimizadas que maximizan la consistencia global del registro, formalmente: nˆ K, ˆ V2 1,ˆ V2 2o= arg max {K,V 2 1,V 2 2}ω(K)(5.8) donde ω(K) = ECC (I1(m),I2(fK(m))) (5.9) siendo m⊂φV2 1(|K|)las coordenadas de los píxeles localizados en el interior de la red triangular dada por K. La búsqueda greedy, detallada en el pseudocódigo de la figura 5.9, empieza creando una lista ordenada (en orden descendente) de ∆ωindexada por (e, accion). Generar esta lista conlleva un importante costo computacional, pero sólo se realiza una vez, al inicio del procedimiento. En cada iteración del proceso de búsqueda se toma el primer elemento de la tabla, esto es, la acción con el mayor ∆ω, y se aplica si ∆ω > δ. Seguidamente se actualiza la lista ordenada con el conjunto de acciones aplicables sobre las aristas de la región de influencia de la acción y se evalúan sus correspondientes precondiciones, así como la posibilidad de que su aplicación de lugar a la aparición una configuración de red previa (ciclo). El algoritmo finaliza cuando la lista de acciones aplicables está vacía o ninguna de las acciones aplicables produce una mejora por encima del δdado. Por último, para ilustrar el método de optimización propuesto se presenta un experimento consistente en registrar dos imágenes de una escena de inte- 106 5. OPTIMIZACIÓN DE REDES TRIANGULARES •procedimientos de optimización similares que utilizan como funciones de costo medidas de similitud basadas en intensidad: SSD [77] y diferencias de imágenes [78]. •métodos de registro no-rígidos tradicionales: polinomial [37], lineal por trozos (sin optimización) [39], cúbico por trozos [41] y thin-plate-splines [19]. 5.4.1 Datos y metodología En la realización de estos experimentos se han empleado imágenes del repositorio ALOI [36] (que incluye una amplia variedad de objetos), escenas urbanas (por ejemplo, fachadas de edificios) e imágenes de satélite de alta resolución. El propósito de elegir esta amplia variedad de imágenes es validar el comportamiento del método con distintos tipos de iluminación, contenidos de la imagen y ángulos de observación. El conjunto de puntos de control (CP) utilizado para generar las redes triangulares iniciales se identifica automáticamente mediante el procedimiento descrito en la sección 4.4. A pesar de que este método proporciona un importante número de puntos correctamente emparejados, muchos de ellos son identificados en diferentes planos de la escena, los cuales, al ser triangulados mediante Delaunay [93], dan lugar a gran cantidad de triángulos no consistentes. El objetivo es, por tanto, detectar y corregir estas configuraciones indeseadas y generar otra red que maximice el número de triángulos que caen sobre proyecciones de planos 3D de la escena, y por ende, mejore la consistencia del registro. Para cuantificar esta mejora se emplean los siguientes procedimientos: 1. Midiendo la consistencia global del registro. Para ello se utilizan las siguientes métricas: (a) el ECC de las imágenes de referencia y registrada, aunque para ser preciso, sólo se evalúa la región de la imagen delimitada por la envolvente convexa de la red (expresión (5.9)), y 5.4. PRUEBAS EXPERIMENTALES Y COMPARATIVAS 107 Imágenes y redes iniciales Red final a) b) e) f) d) h) c) g) Figura 5.12: (a-d) Imágenes reales de escenas poliédricas y sus correspondientes redes triangulares de Delaunay. (e-h) Redes triangulares optimizadas por el método propuesto (sólo acciones de intercambio). 108 5. OPTIMIZACIÓN DE REDES TRIANGULARES Tabla 5.1: Datos de las redes triangulares, número de acciones de intercambio aplicadas, consistencia global y tiempo empleado en el proceso de optimización. Test Morris & Kanade Nakatsuji et al. Nuestro método Fig. 5.12 # aristas # acciones Tiempo2ECC # acciones Tiempo ECC # acciones Tiempo ECC a 56 13 11.42 0.239 14 9.21 0.235 13 1.42 0.239 b 47 13 7.39 0.226 10 8.11 0.321 8 1.65 0.332 c 275 102 133.78 0.423 65 72.82 0.398 27 3.93 0.425 d 65 18 38.43 0.263 11 9.75 0.211 16 3.15 0.264 2En los experimentos se ha empleado Matlab sobre un Pentium Core 2 Duo 6400. Para ganar en velocidad, algunas de las rutinas más costosas computacionalmente, se han desarrollado en C como extensiones de Matlab, así como el software de búsqueda y emparejamiento de correspondencias. El tiempo se expresa en segundos. 5.4. PRUEBAS EXPERIMENTALES Y COMPARATIVAS 109 (b) el RMSE de un conjunto de ICPs seleccionados en ambas imágenes (sólo en imágenes de satélite). 2. Evaluando las inconsistencias de la reconstrucción 3D (sin escala) generada a partir de las redes triangulares, antes y después de ser optimizadas. Puesto que no se dispone de un DTM de la zona en el caso de las imágenes de satélite ni se conoce la geometría exacta de la escena observada en el resto, la evaluación se realiza mediante inspección visual. 5.4.2 Comparativa con métodos similares En esta sección se compara el comportamiento del método propuesto con los métodos de Morris y Kanade [77] y Nakatuji et al. [78]. Estos dos métodos detectan inconsistencias topológicas en redes triangulares 3D que aproximan la geometría de una escena y tratan de eliminarlas mediante el intercambio de aristas. Más concretamente, el primero emplea una búsqueda greedy y una función de coste global basada en SSD, mientras que el segundo propone una búsqueda aleatoria y una función de detección de inconsistencias basada en la diferencia de imágenes. Esta función emplea una plantilla de tamaño fijo que se aplica a la diferencia de los cuadriláteros afectados por la acción, los cuales han sido mapeados mediante transformaciones afines. Como ambos métodos solamente emplean acciones de intercambio de aristas, con objeto de que éstos se puedan comparar con el método propuesto, el proceso de optimización sólo considera esta acción. La figura 5.12(a-d) muestra algunas de las imágenes empleadas en este trabajo, así como las redes iniciales generadas a partir de los conjuntos de CPs identificados en ellas. La tabla 5.1 reúne los datos de las redes triangulares iniciales, así como el número de acciones aplicadas, el tiempo de procesamiento y la consistencia del registro correspondientes a cada uno de los métodos comparados. De estos resultados se pueden destacar diversos aspectos: el bajo coste computacional del método propuesto, su efectividad (consistencia vs. número de acciones) y su robustez a los cambios en la iluminación, como se aprecia en los resultados obtenidos para el par de imágenes de la figura 5.12(b) (dos imágenes 110 5. OPTIMIZACIÓN DE REDES TRIANGULARES adquiridas en condiciones de iluminación muy diferentes). Obsérvese la consistencia del registro del método de Morris y Kanade para el mismo par de imágenes, a lo que hay que sumar su coste computacional, prohibitivo cuando el número de aristas crece significativamente (ver los resultados para el par de imágenes de la figura 5.12(c)). El método de Nakatsuji et al. tiene un comportamiento similar a éste último a pesar de que emplea una función de detección local, que es costosa ya que requiere la estimación de 4 afinidades por arista. 0 5 10 15 20 25 30 0.395 0.4 0.405 0.41 0.415 0.42 0.425 0.43 ECC Acción nº 0 2 4 6 8 10 12 14 0.195 0.2 0.205 0.21 0.215 0.22 0.225 0.23 0.235 0.24 ECC Acción nº a) c) b) d) Consistencia del registro Consistencia del registro 1 2 3 4 5 6 7 8 0.24 0.25 0.26 0.27 0.28 0.29 0.3 0.31 0.32 0.33 0.34 ECC Acción nº 0 2 4 6 8 10 12 14 16 0.19 0.2 0.21 0.22 0.23 0.24 0.25 0.26 0.27 0.28 ECC Acción nº Figura 5.13: Evolución de la consistencia global (medida con el ECC) durante el proceso de optimización de los pares de imágenes mostradas en las figuras 5.12(ad)). En general, se observa una mejora significativa de la consistencia con todos los pares analizados. Las figuras 5.12(e-h) muestran las redes optimizadas por el método propuesto correspondientes a los pares de imágenes de las figuras 5.12(a-d). Como se aprecia en todos los ejemplos, el proceso de optimización finaliza con redes triangulares consistentes con la geometría de la escena y en todos 5.4. PRUEBAS EXPERIMENTALES Y COMPARATIVAS 111 los casos, en menos de 25 iteraciones, como recogen las gráficas de la figura 5.13(a-d). Nótese que durante el proceso de optimización hay acciones que, aparentemente, no mejoran la consistencia global del registro (como revelan los pequeños trozos planos en las curvas de la figura 5.13(b,c). Obsérvese que se dice aparentemente, ya que, como se menciona arriba, las gráficas muestran la consistencia global del registro, no la consistencia de la región afectada por la acción, que siempre mejora3. Estas acciones, a pesar de que apenas mejoran la consistencia global, dan lugar a configuraciones topológicas que son explotadas en iteraciones posteriores para generar redes consistentes con las proyecciones de la escena. Los métodos analizados proponen el intercambio de aristas como única acción para la búsqueda de nuevas configuraciones topológicas. A continuación se muestra la red optimizada para las imágenes de la figura 5.12(d) (la más densa de las utilizadas) obtenida con el método propuesto, pero incluyendo el resto de las acciones descritas en este capítulo. Como se detalla en la sección previa, en primer lugar se simplifica la red mediante la aplicación iterativa de acciones de eliminación (ver fig. 5.14(b)), una vez simplificada, se lanza el proceso optimización con la totalidad de las acciones (ver fig. 5.14(c)). Comparando este resultado con el obtenido con la sola aplicación de acciones de intercambio (ver fig. 5.12(g)), la red continúa siendo consistente pero mucho menos densa. Se ha comprobado experimentalmente que esta estrategia ofrece mejores resultados (en términos de precisión y eficiencia) que un único proceso de optimización, sin simplificación previa. 5.4.3 Comparativa con métodos tradicionales En esta sección se compara el comportamiento del método propuesto con diversas técnicas de registro no-rígidas utilizadas en la comunidad científica e incluidas en una amplia variedad de software: Matlab [99], Insight Segmentation and Registration Toolkit (ITK) [79], el software de Image Fusion Systems Research [51] o paquetes de procesamiento de imágenes para tele3Como se detalla en la sección 5.3.2, una acción se aplica si y sólo si, ∆ω≥δ > 0lo que asegura un incremento monótono de la consistencia global, aunque cuanto menor sea la variación menor relevancia tendrá sobre ésta. 112 5. OPTIMIZACIÓN DE REDES TRIANGULARES Red inicial Red simplificada Red optimizada b) c)a) Figura 5.14: Red triangular optimizada con el método propuesto: (a) red inicial, (b) red simplificada (sólo acciones de eliminación) y (c) red optimizada (todas las acciones). Obsérvese como, dada la disposiciones de las correspondencias, no ha sido necesaria la inclusión de nuevos vértices. detección como: ENVI/IDL [54], ERDAS [61], etc. Más concretamente, las técnicas analizadas aquí son funciones polinomiales (de 2ohasta 4ogrado) (POL2 - POL4), lineales por trozos (PWL), cúbicas por trozos (PWC) y thin-plate-spline (TPS). Para realizar los experimentos se han utilizado dos imágenes de altaresolución del satélite QuickBird de aprox. 3000 ×1500 píxeles ≡1800 × 900 m2, de la ciudad costera del Rincón de la Victoria (Málaga) (ver fig. 5.15). Las imágenes, adquiridas en invierno y primavera, presentan diferencias radiométricas muy significativas (brillo y contraste, sombras, cambios en la cobertura terrestre, etc.) e importantes distorsiones inducidas por la adquisición fuera del nadir (con diferencias de 32.1 grados). A parte de las imágenes completas, en esta comparativa se analizan tres áreas de estudio extraídas de éstas (de aprox. of 900 ×750 píxeles): •área urbana, que contiene edificios de diferentes alturas sobre un terreno casi plano, •área rural, una región prácticamente plana pero, sin elementos de importante altura (excepto algunos árboles) y •área residencial, una región accidentada que contiene edificios de baja altura (viviendas unifamiliares de una o dos plantas). Al seleccionar estas regiones de estudio se pretende validar el rendimiento 5.4. PRUEBAS EXPERIMENTALES Y COMPARATIVAS 113 Área rural Área residencial Línea de costa Montañas 2 5 25 25 25 25 25 25 50 50 50 50 50 50 50 50 5 0 50 50 75 75 75 75 75 75 75 75 75 75 75 75 75 75 100 100 100 100 100 100 1 00 100 100 100 100 100 125 125 125 125 125 125 125 150 150 150 175 200 225 distance (m) distance (m) Elevation contours 3.85 3.855 3.86 3.865 3.87 x 10 5 4.064 4.0645 4.065 4.0655 x 10 6 60 80 100 120 140 160 180 200 220 240 a) b) Área urbana Terrain profile distance (m) elevation (m) Figura 5.15: Imágenes de prueba utilizadas en el estudio comparativo. (a) Dos imágenes de la ciudad del Rincón de la Victoria (Málaga). Las regiones de interés utilizadas en los tests aparecen marcadas en la imagen de referencia. (b) Información de elevación utilizada para definir apropiadamente las regiones de estudio. del método propuesto con escenas de distinta geometría y contenidos. Los resultados de esta comparativa aparecen resumidos en la tabla 5.2. A la vista de estos resultados se pueden subrayar diversos aspectos: a) las funciones globales polinomiales, como era esperado, ofrecen peores resultados que las técnicas locales e híbridas, ya que las imágenes presentan distorsiones relativas de importancia, especialmente, en las zonas con relieve accidentado; 114 5. OPTIMIZACIÓN DE REDES TRIANGULARES Tabla 5.2: Comparación del método propuesto con otros métodos de registro bien conocidos, concretamente: un método global polinomial (de 2oa 4ogrado); dos métodos locales (lineal y cúbico por trozos); y un método híbrido (thin-plate-spline). RMSE (m.) ECC Función Urbana Residencial Rural Img. completa Urbana Residencial Rural Img. completa POL 2oorden 3.20 2.71 2.18 8.12 0.330 0.362 0.356 0.217 POL 3er orden 3.27 2.05 1.82 7.67 0.310 0.404 0.383 0.219 POL 4oorden 2.86 2.07 1.75 6.00 0.338 0.406 0.393 0.243 PWL 2.00 1.26 1.06 1.82 0.447 0.527 0.472 0.449 PWC 2.24 1.33 1.27 2.15 0.425 0.494 0.444 0.435 TPS 2.02 0.92 1.03 1.66 0.437 0.556 0.487 0.464 Nuestro método 1.80 0.93 0.96 1.45 0.465 0.549 0.494 0.467 5.4. PRUEBAS EXPERIMENTALES Y COMPARATIVAS 115 b) entre los tres métodos tradicionales con alcance local, TPS es el más preciso de los analizados, lo que viene a confirmar lo observado en el capítulo 3, y c) el método propuesto en este trabajo iguala e incluso mejora los resultados del TPS. La figura 5.16 muestra la evolución de la consistencia global al aplicarlo sobre las diferentes áreas de estudio, incluida la imagen completa. El lector puede observar que la consistencia mejora en todos los casos analizados, independientemente de la medida utilizada (RMSE o ECC). Estas mejoras se acentúan en las regiones urbana y residencial ya que, aunque la zona rural es casi plana, las redes de Delaunay generadas inicialmente son razonablemente precisas y no requieren de la aplicación de muchas acciones. Para las tres primeras áreas (a, b y c) el algoritmo requiere entre 20 y 25 iteraciones, mientras que para la imágenes completas require en torno a 71. El número de iteraciones, y por tanto, el coste computacional, depende fuertemente del número de aristas compartidas, como se observa en los tiempos de ejecución mostrados en la tabla 5.3. Tabla 5.3: Datos de la red, número de acciones aplicadas y tiempo empleado en el proceso de optimización. Región # triángulos # aristas # acciones Tiempo Urbana 168 258 23 7.30 Residencial 155 239 17 6.80 Rural 162 248 24 7.73 Toda la imagen 502 759 71 39.04 Nótese que durante el proceso de optimización hay acciones que mejoran la consistencia del registro, según se observa en la curva del ECC (incrementando su valor), mientras que el RMSE permanece fijo (como revelan los pequeños trozos planos en las curvas). Esta discrepancia entre ambas métricas se debe al hecho que los triángulos implicados en las acciones no contienen ICPs en su interior y, consecuentemente, el RMSE no se ve alterado por la aplicación de la acción. Es importante recordar que, aunque los ICPs fuesen 122 6. CONCLUSIONES Y TRABAJOS FUTUROS bles defectos, etc. Estos son sólo algunos ejemplos pero que hacen entrever la importancia de este proceso y que justifica el interés en desarrollar nuevos métodos y procedimientos en este campo. Fruto de este interés son este trabajo y los diferentes artículos científicos publicados durante su realización. Esta tesis aborda la problemática asociada al proceso de registro, así como los distintos elementos que intervienen en él y que, de uno u otro modo, afectan a la precisión de los resultados obtenidos. Por su importancia, una parte signitificativa de este trabajo se centra en el análisis de dos aspectos claves de este proceso: la medida de la consistencia, crucial para determinar cómo de bueno ha sido el “encaje” de una imagen sobre la otra, y el tipo y alcance de la función de transformación utilizada y que determina la complejidad de las diferencias geométricas puede corregir. De esta evaluación se puede subrayar lo siguiente: •La robustez de la información mutua (MI) (y sus diferentes variantes normalizadas) para medir la similitud radiométrica de dos imágenes, incluso cuando las imágenes presentan cambios radiométricos no funcionales. Su dependencia de la correcta estimación de las distribuciones de probabilidad de los niveles de intensidad de las imágenes. Estimaciones que depende, a su vez, de aspectos tales como el tamaño de las imágenes, la función de interpolación utilizada y la correcta estimación del histograma conjunto. Entre los procedimientos de estimación analizados destaca la ventana de Parzen, un método capaz de producir estimaciones poco ruidosas a partir de conjuntos pequeños de muestras. •Los funciones de transformación locales (lineales por trozos) e híbridas (thin-plate-spline) exhiben un comportamiento superior a las funciones globales (polinomial) ya que pueden explotar la información, relativa a las diferencias geométricas de ambas imágenes, proporcionada por conjuntos de correspondencias muy densos, decayendo su rendimiento con conjuntos dispersos. De todas las transformaciones evaluadas, sólo las funciones lineales por trozos proporcionan un rendimiento aceptable cuando la geometría de la escena es poliédrica, típica en entornos de interiores y urbanos. Para que este 6.2. TRABAJOS FUTUROS 123 método trabaje adecuadamente las regiones triangulares conjugadas deben caer sobre proyecciones de superficies planas de la imagen (triángulos consistentes con la escena). Sin embargo, los procedimiento utilizado para generar automáticamente las regiones conjugadas en el software disponible actualmente no garantiza este hecho, lo que introduce importantes errores en el registro. Este trabajo analiza en detalle esta problemática proponiendo un enfoque novedoso para generar redes triangulares conjugadas que se disponen sobre las imágenes maximizando el número de triángulos consistentes. Para ello, plantea, mediante complejos simpliciales, una formalización del registro lineal por trozos que va más allá de un mero enfoque 2D, puesto que abarca no sólo la generación de redes triangulares 2D óptimas (desde el punto de vista de la consistencia del registro), sino también formaliza la reconstrucción 3D de la superficie de la escena observada. También se propone un procedimiento de optimización sustentado en esta formalización que modifica topológica y geométricamente la red triangular (complejo simplicial) utilizada en el registro lineal por trozos para mejorar su precisión. Este método emplea el coeficiente de correlación de entropía como medida de consistencia. El método propuesto se ha comparado experimentalmente con diversas funciones de transformación no-rígidas y procedimientos similares propuestos en los campos del registro de imágenes y de la reconstrucción 3D. Los resultados de estas comparativas reflejan la precisión del método propuesto, con errores en todos los casos inferiores a los otros métodos evaluados. Respecto a su coste computacional, su tiempo de ejecución es sensiblemente menor que los otros métodos iterativos analizados, aunque, lógicamente, continúa siendo excesivamente costoso comparado con los métodos de registro no-rígido basado en funciones polinomiales, lineales por trozos (sin optimización) y thin-plate-splines, lo que limita su uso en cierto tipo de aplicaciones. 6.2 Trabajos futuros Una de las contribuciones de esta tesis es la formalización del proceso de registro lineal por trozos mediante complejos simpliciales. Esta formalización 124 6. CONCLUSIONES Y TRABAJOS FUTUROS asume que los complejos simpliciales utilizados en el registro son proyecciones de una superficie lineal por trozos que aproxima la escena observada. De este modo, si se determina la topología y/o geometría de red óptima, el proceso de registro no sólo mejora la precisión del registro, sino que, de generarse una reconstrucción 3D a partir de estas proyecciones, ésta estaría libre de inconsistencias y se ajustaría con precisión a la escena observada. Pese a las posibilidades este novedoso enfoque, los complejos simpliciales manejados en nuestra implementación son redes 1-conectadas (esto es, todos los símplices máximos son triángulos y todos los triángulos tienen al menos un adyacente). Esta circunstancia da lugar a que, una vez finalizado el proceso de optimización, continúe habiendo triángulos que no se corresponden con proyecciones reales de la escena. Sería conveniente, por tanto, incluir en el proceso de optimización nuevas acciones capaces de conectar/desconectar coherentemente componentes de los complejos simpliciales ganando así en flexibilidad y potencia. La búsqueda greedy asegura un coste computacional bajo, pero puede caer en mínimos locales. Existen mecanismos que pueden mejorar significativamente este aspecto aunque a costa de incrementar su coste computacional. Por ejemplo, una interesante opción consistiría en implementar un procedimiento de lookahead que permitiría anticipar como evolucionaría la consistencia del registro antes de aplicar una determinada acción, de forma que sería posible descartar una acción que, aunque produce la mayor mejora, da lugar a una configuración tal que no es posible seguir mejorando la consistencia. Uno de los principales inconvenientes del método de registro propuesto, en comparación con otras técnicas de registro, es su elevado coste computacional. Quizás sea éste el mayor lastre para su incorporación en sistemas que requieran un tiempo de respuesta mucho más crítico como por ejemplo las aplicaciones robóticas. Pese a esto, la utilización de redes triangulares ofrece la posibilidad de paralelizar operaciones: sería relativamente fácil, por ejemplo, paralelizar el proceso de estimación de la función de coste local. Este hecho que adquiere mayor relevancia dada la reciente aparición de plataformas especialmente preparadas para realizar estas tareas. Este es el caso, por ejemplo, de librerías como RapidMind R [84] u OpenVidia [81] que pro- 6.2. TRABAJOS FUTUROS 125 porciona soporte para procesadores de múltiples núcleos (Intel R y AMD R ), GPUs de última generación (nVidia R y ATI R ), y procesadores Cell R . 126 6. CONCLUSIONES Y TRABAJOS FUTUROS Apéndice A Entropía, entropía relativa e información mutua El concepto de información es muy amplio para ser recogido en una única definición. Sin embargo, para cualquier distribución de probabilidad es posible definir una cantidad que posee muchas de las propiedades que una medida de información debería tener, la entropía. La noción de entropía se extiende para definir la información mutua, que mide la cantidad de información que una variable aleatoria contiene sobre otra. La información mutua es un caso especial de una cantidad más general denominada entropía relativa, que mide la distancia entre dos distribuciones de probabilidad. La entropía, entropía relativa e información mutua están íntimamente relacionadas y comparten numerosas propiedades. Este apéndice define todas y cada una de estas cantidades, deriva algunas de sus propiedades e introduce la mayoría de las ideas básicas necesarias para entender algunos de las conceptos desarrollados en esta tesis. El lector puede referirse al libro de Cover y Tomas [26] para un estudio detallado sobre los elementos de la teoría de la información. A.1 Entropía La entropía es una medida de información, o más concretamente, de la incertidumbre de una fuente de información. Para poder entender el concepto 127 128A. ENTROPÍA, ENTROPÍA RELATIVA E INFORMACIÓN MUTUA que encierra obsérvese el siguiente ejemplo: Ejemplo A.1 Dada una fuente binaria de información, F, que proporciona el ganador de un partido de fútbol, esto es: el equipo A, con una probabilidad p(A)de 3/4, o el B, con probabilidad p(B)de 1/4, de tal manera que la situación que se tiene es la siguiente: F? ABBABBAB A B La transmisión de esta información a través del canal se puede abordar de distintas formas. Supóngase, por ejemplo, que se envían codificados en binario los resultados de tres partidos dando lugar a una codificación como la que sigue: Cadena Prob. Código Longitud AAA 27/64 0 1 ABA 9/64 100 3 BAA 9/64 101 3 AAB 9/64 110 3 BBA 3/64 11100 5 BAB 3/64 11101 5 ABB 3/64 11110 5 BBB 1/64 11111 5 Observando el código se comprueba que se transfieren una mayor cantidad de bits en aquellas cadenas que tienen menor probabilidad de ocurrir. La longitud media del código es: L=27 64 ×1 + 9 64 ×3 + 3 64 ×3 + 1 64 ×1 = 3,47 bits De acuerdo con la nueva codificación, se tiene que p(0) = 0,36 yp(1) = A.1. ENTROPÍA 129 0,63. La entropía es una medida de la cantidad de información que se recibe al final del canal y se mide típicamente en bits. Definición A.1 Dada una variable aleatoria discreta Xque tiene una determinada distribución de probabilidades, p(x), se define la entropía H(X) como: Hb(X) = −X x∈X p(x) logbp(x) =−X x∈X p(x) logbp(x) =−Elogbp(X)(A.1) La unidad en que se mide depende de la base del logaritmo utilizada, estas son: BITS si es base 2, DITS si es base 10 o NATS si son logaritmos neperianos. Ejemplo A.2 De acuerdo con el ejemplo anterior, la entropía de la fuente F sería: H2(X) = 3 4log2 3 4+1 4log24 = 0,81 bits El resultado esperado más probable es el que, obviamente, menos contribuye al valor de la entropía; luego ésta puede entenderse como una ponderación de las contribuciones de cada suceso. A continuación se enumeran algunas propiedades de la entropía: Lema A.0.1 Hb(X)≥0 Prueba: 0≤p(x)≤1implica log (1/p (x)) ≥0 Lema A.0.2 Hb(X) = Ha(X) logba 130A. ENTROPÍA, ENTROPÍA RELATIVA E INFORMACIÓN MUTUA Prueba: logbp= logbalogap Como se desprende del lema A.0.2, la entropía puede cambiarse de una base a otra (esto es, de una unidad a otra) multiplicando por el factor apropiado. Ejemplo A.3 Si F es una fuente de información binaria tal que p(0) = py p(1) = 1 −p=qla entropía (en base 2) vendría dada por la expresión: H(X) = H(p, q) = −plog2p−qlog2q 00.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1  H X p Figura A.1: Comparativa H(X)versus p. La figura A.1 muestra la representación gráfica de H(X)en función de p. De la observación de esta curva se derivan algunas de las propiedades de la entropía: si p(x) = 0 óp(x)=1, esto es, no hay incertidumbre sobre el valor que tomará X, entonces H(X) = 0. Por el contrario, si la incertidumbre es máxima p(x) = 1/2, entonces el valor de H(X)también lo es. A.2 Entropía conjunta y entropía condicional En esta sección se extiende el concepto de entropía a un par de variables aleatorias XeY. A.2. ENTROPÍA CONJUNTA Y ENTROPÍA CONDICIONAL 131 Definición A.2 Dado el par de variables aleatorias discretas (X, Y )con distribución conjunta de probabilidades, p(x, y), se define la entropía conjunta H(X, Y )como: H(X, Y ) = −X x∈XX y∈Y p(x, y) log2p(x, y) =−Elog2p(X, Y )(A.2) igualmente, Definición A.3 Se define la entropía condicional H(Y|X)de la la variable aleatoria Ydada Xcomo: H(Y|X) = X x∈X p(x)H(Y|X=x) =−X x∈X p(x)X y∈Y p(y|x) log2p(y|x) =−X x∈XX y∈Y p(x, y) log2p(y|x) =−Elog2p(Y|X)(A.3) La relación natural entre entropía conjunta y condicional queda patente en la denominada regla de la cadena, esta es: la entropía conjunta de un par de variables aleatorias es la entropía de una de ellas más la entropía condicional de la otra. Formalmente: Teorema A.1 (Regla de la cadena) H(X, Y ) = H(X) + H(Y|X)(A.4)