Repositorio Institucional de Documentos

Aplicación de técnicas de alineamiento de imágenes para la georreferenciación automática de imágenes aéreas y satelitales Proyecto Fin de Carrera - Ingeniería Informática - Febrero de 2012 Autor: Pablo Fuertes Correa Director: Pedro R. Muro Medrano ´ Indice 1. Introducci´on 1 1.1. Resumen................................ 1 1.2. Contexto profesional . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.3. Contexto tecnol´ogico . . . . . . . . . . . . . . . . . . . . . . . . . 2 2. Trabajo realizado 4 2.1. Primera aproximaci´on (basada en caracter´ısticas) . . . . . . . . . 4 2.2. M´etodo seleccionado (basado en la intensidad) . . . . . . . . . . 10 2.3. Implementaci´on............................ 14 3. Conclusiones y futuras l´ıneas de investigaci´on 23 4. Anexos 25 4.1. Modelos no lineales . . . . . . . . . . . . . . . . . . . . . . . . . . 25 4.2. C´alculo del Gradiente y la Hessiana . . . . . . . . . . . . . . . . 26 4.3. M´etodo Levenberg-Marquardt . . . . . . . . . . . . . . . . . . . . 27 4.4. Modelos de deformaci´on y estimaci´on de par´ametros . . . . . . . 29 4.5. Newton-Raphson........................... 32 4.6. Descenso del gradiente . . . . . . . . . . . . . . . . . . . . . . . . 33 4.7. LMA(Newton-Raphson y Descenso del gradiente) . . . . . . . . . 33 4.8. Estimaci´on perspectiva mediante aproximaciones afines . . . . . 34 4.9. Transformaci´on af´ın inversa . . . . . . . . . . . . . . . . . . . . . 39 4.10. Transformaci´on singular . . . . . . . . . . . . . . . . . . . . . . . 41 4.11. Diferencias Centrales . . . . . . . . . . . . . . . . . . . . . . . . . 44 4.12. Transformaci´on Log-Polar . . . . . . . . . . . . . . . . . . . . . . 45 4.13. Transformada de Fourier-Mellin . . . . . . . . . . . . . . . . . . . 46 4.14. Pir´amide de multiples resoluciones . . . . . . . . . . . . . . . . . 47 4.15. Levenberg-Marquardt Optimizado . . . . . . . . . . . . . . . . . 49 4.16. Lenguajes contemplados . . . . . . . . . . . . . . . . . . . . . . . 52 4.16.1.C++.............................. 52 4.16.2.R................................ 52 4.16.3.Python ............................ 52 4.16.4.Java.............................. 53 4.16.5. Matlab/Octave . . . . . . . . . . . . . . . . . . . . . . . . 53 4.16.6.OpenCV............................ 53 4.17. Acr´onimos y abreviaturas . . . . . . . . . . . . . . . . . . . . . . 55 ´ Indice de Figuras 56 Referencias 58 1. Introducci´on 1.1. Resumen En el presente proyecto fin de carrera se investigan las t´ecnicas de alineaci´on de im´agenes con la finalidad de establecer la correspondencia de una imagen capturada con una base de im´agenes que sirvan como referencia. Esta base de im´agenes ser´a la aportada por el Plan Nacional de Ortofoto A´erea (PNOA) que ofrece una fotograf´ıa continua de toda Espa˜na perfectamente georreferenciada. Con tal fin se hace una implementaci´on de una t´ecnica basada en caracter´ısticas llegando a la conclusi´on de que esta no ofrece unos resultados aceptables en el ´ambito de aplicaci´on. Teniendo en cuenta estos resultados se han buscado alternativas que pudieran resolver esta problem´atica. Para ello se ha ampliado el abanico de t´ecnicas posibles puesto que se ha considerado que el ´ambito de aplicaci´on no requiere los resultados en tiempo real o cerca del mismo. Esto ha hecho que tambi´en sean contempladas t´ecnicas basadas en la intensidad. As´ı mismo se han realizado pruebas con diferentes lenguajes para determinar cual resulta m´as adecuado en cuanto a velocidad de prototipado y a eficiencia. De estas pruebas se ha decidido realizar la implementaci´on en C++ junto a OpenCV. No obstante el prototipo final ha sido realizado en Matlab, por motivos que ser´an detallados m´as adelante. Una vez escogido el algoritmo [1] que m´as se ajusta a las restricciones del problema y que mejores resultados ofrece, se ha procedido a su implementaci´on. Se ha observado que era demasiado complejo como primera aproximaci´on y por tanto se ha buscado un algoritmo anterior [2], de los mismos autores. Este art´ıculo es la base en el que se apoya el que se pretende implementar, este hecho hace que compartan la tecnolog´ıa subyacente pero con menor n´umero de optimizaciones y por tanto resulta m´as f´acil su implementaci´on e interpretaci´on. Se ha implementado dicho algoritmo, obteniendo resultados no satisfactorios. Esto ha hecho necesario realizar una revisi´on en profundidad del funcionamien- to del m´etodo y de los aspectos t´ecnicos de su implementaci´on. Siendo tambi´en necesario rehacer el an´alisis matem´atico planteado por los autores para buscar alguna errata en su publicaci´on u omisi´on de alg´un paso que provocase el funcionamiento incorrecto. Con el fin de validar que la implementaci´on realizada fuera correcta se ha realizado una implementaci´on alternativa en Matlab con id´entico resultado. Desde este momento se ha continuado trabajando con el prototipo de Matlab y se han explorado m´ultiples opciones para intentar solucionar el problema, delegando en funciones de matlab, transformando el problema en el inverso, realizando el c´alculo num´erico de las derivadas y otras muchas variantes hasta llegar a una soluci´on que aportase un resultado aceptable. 1.2. Contexto profesional El presente proyecto fin de carrera ha sido desarrollado en el Grupo de Sistemas de Informaci´on Avanzados (IAAA) del Departamento de Inform´atica e Ingenier´ıa de Sistemas (DIIS) perteneciente a la Universidad de Zaragoza. La actividad de investigaci´on del grupo de Sistemas de Informaci´on Avanzados (IAAA) est´a enfocada en las tecnolog´ıas de software para la creaci´on 1 de sistemas de informaci´on con datos georreferenciados (tambi´en denominados generalmente datos geogr´aficos y ahora m´as com´unmente datos espaciales) y, especialmente, en la configuraci´on de un nuevo paradigma que se ha convenido en llamar Infraestructuras de Datos Espaciales (IDEs). Se trata de cubrir desde aspectos b´asicos de ingenier´ıa de servicios hasta aspectos metodol´ogicos de la puesta en funcionamiento de sistemas industriales que explotan la informaci´on geogr´afica. Concretamente, la labor de investigaci´on del grupo aborda problem´aticas vinculadas a las ontolog´ıas y meta-datos; interoperabilidad, composici´on y encadenamiento de servicios; visualizaci´on inteligente de informaci´on geogr´afica, meta-datos, procesamiento sem´antico y recuperaci´on inteligente de informaci´on (indexaci´on, interrogaci´on y recuperaci´on); m´etodos y procesos para la creaci´on de informaci´on y el modelado de contenidos heterog´eneos; y, finalmente, un aspecto que se considera fundamental para la realimentaci´on t´ecnica y que es la definici´on, creaci´on y puesta en funcionamiento real de nuevos servicios y aplicaciones. 1.3. Contexto tecnol´ogico El desarrollo de este proyecto se realiza en el ´ambito del conjunto de problemas englobados dentro del t´ermino anglosaj´on image registration que en adelante denominaremos alineamiento geom´etrico, puesto que en el fondo se trata de encontrar una relaci´on entre las distintas im´agenes. De forma m´as precisa, tal como se describe en el estudio de Brown [3]: El alineamiento geom´etrico es el proceso de transformaci´on de diferentes conjuntos de datos a un sistema de coordenadas com´un. Los datos pueden ser m´ultiples fotograf´ıas, datos de diferentes sensores, de diferentes ´epocas o de diferentes puntos de vista. Se utiliza en visi´on artificial, imagen m´edica, reconocimiento autom´atico de objetivo y en la recopilaci´on y an´alisis de im´agenes y datos de los sat´elites. El registro es necesario para poder comparar o integrar los datos obtenidos de estas diferentes mediciones. Clasificaci´on general de los algoritmos: Basados en la intensidad versus basados en caracter´ısticas Los m´etodos basados en la intensidad comparan los patrones de intensidad en las im´agenes a trav´es de m´etricas de correlaci´on, mientras que los m´etodos basados en las caracter´ısticas encuentran la correspondencia entre las caracter´ısticas de las im´agenes, tales como puntos, l´ıneas y contornos. Los m´etodos basados en la intensidad registran im´agenes completas o subim´agenes. Si las sub-im´agenes est´an alineadas, los centros de las correspondientes sub-im´agenes se consideran como puntos de caracter´ısticas correspondientes. Mientras que los m´etodos basados en caracter´ısticas establecen la correspondencia entre varios puntos en las im´agenes. Luego, conociendo la correspondencia entre varios puntos en las im´agenes, se determina una transformaci´on para mapear la imagen objetivo a las im´agenes de referencia, estableciendo punto por punto, la correspondencia entre las im´agenes de referencia y la objetivo. 2 Modelos de transformaci´on La primera categor´ıa de modelos de transformaci´on incluye las transformaciones lineales, que son la traslaci´on, rotaci´on, escalamiento y otras transformaciones afines. Estas son de naturaleza global, por lo tanto, no pueden modelar diferencias geom´etricas locales entre las im´agenes. La segunda categor´ıa de transformaciones permiten que estas sean el´asticas ono r´ıgidas. Estas transformaciones son capaces de deformar localmente la imagen objetivo para alinearla con la imagen de referencia. Las transformaciones no r´ıgidas incluyen las funciones de base radial, modelos f´ısicos continuos y los modelos de grandes deformaciones. M´etodos de dominio espacial versus frecuencia Los m´etodos espaciales operan en el dominio de la imagen, haciendo coincidir los patrones de intensidad o las caracter´ısticas de las im´agenes. Algunos de los algoritmos que emparejan caracter´ısticas son derivaciones de las t´ecnicas tradicionales para realizar el alineamiento manual de la imagen, en el que un operador decide los correspondientes puntos de control en las im´agenes. Cuando el n´umero de puntos de control supera el m´ınimo requerido para definir el modelo de transformaci´on apropiado, los algoritmos iterativos como RANSAC (RANdom SAmple Consensus) se pueden utilizar para estimar de manera robusta los par´ametros de un tipo particular de transformaci´on (por ejemplo, af´ın) para el alineamiento de las im´agenes. Los m´etodos en el dominio de la frecuencia encuentran los par´ametros de transformaci´on para el registro de las im´agenes mientras trabajan en el dominio de transformaci´on. Estos m´etodos trabajan por simple transformaci´on, tales como la traslaci´on, la rotaci´on y el escalado. La aplicaci´on del m´etodo de correlaci´on de fase a un par de im´agenes produce una tercera imagen que contiene un solo pico. La ubicaci´on de este pico corresponde a la traslaci´on relativa entre las im´agenes. A diferencia de muchos algoritmos en el dominio espacial, el m´etodo de correlaci´on de fase es insensible al ruido, las oclusiones y otros defectos t´ıpicos de las im´agenes m´edicas o por sat´elite. Adem´as, la correlaci´on de fase utiliza la transformada r´apida de Fourier (FFT) para calcular la correlaci´on cruzada entre las dos im´agenes, lo que por lo general resulta en grandes mejoras en el rendimiento. Este m´etodo puede ser extendido para determinar la rotaci´on y las diferencias de escala entre dos im´agenes convirtiendo, en primer lugar, las im´agenes a coordenadas log-polar. Debido a las propiedades de la Transformada de Fourier, los par´ametros de rotaci´on y de escalamiento se pueden determinar de forma invariable a la traslaci´on. M´etodos monomodal versus multimodal Los m´etodos monomodales tienden a registrar las im´agenes en la misma modalidad adquirida por el ´unico tipo de esc´aner/sensor, mientras que los m´etodos de multimodales tienden a registrar las im´agenes obtenidas por diferentes tipos de esc´aner/sensor. Los m´etodos de registro multimodal se utilizan a menudo en imagen m´edica ya que las im´agenes de un paciente se obtienen con frecuencia a partir de esc´aneres diferentes. Los ejemplos incluyen el registro de im´agenes del 3 cerebro por tomograf´ıa axial computarizada/resonancia magn´etica o de todo el cuerpo por PET/TAC para la localizaci´on de tumores, el registro de im´agenes TAC con y sin realce de contraste para la segmentaci´on de im´agenes de partes espec´ıficas de la anatom´ıa y el registro de ultrasonido y de im´agenes TAC para la localizaci´on de la pr´ostata en radioterapia. M´etodos autom´aticos versus interactivos Se han desarrollado m´etodos manuales, interactivos, semi-autom´aticos y autom´aticos. Los m´etodos manuales proporcionan herramientas para alinear las im´agenes de forma manual. Los m´etodos interactivos reducen el sesgo del usuario mediante la realizaci´on de ciertas operaciones claves de forma autom´atica mientras que todav´ıa conf´ıan en que el usuario gu´ıe el registro. Los m´etodos semiautom´aticos realizan m´as pasos de registro de forma autom´atica pero dependen del usuario para verificar la correcci´on de un registro. Los m´etodos autom´aticos no permiten la interacci´on del usuario y realizan todos los pasos del registro de forma autom´atica. Medidas de similitud para el registro de la imagen Una medida de similitud de im´agenes cuantifica el grado de similitud entre los patrones de intensidad de dos im´agenes. La elecci´on de una medida de similitud de im´agenes depende de la modalidad de las im´agenes a ser alineadas. Los ejemplos m´as comunes de las medidas de similitud de im´agenes incluyen la correlaci´on cruzada, la informaci´on mutua, la suma de los cuadrados de las diferencias de intensidad y la raz´on de uniformidad de la imagen. La informaci´on mutua y la informaci´on mutua normalizada son las medidas de similitud de imagen m´as populares para el registro de im´agenes multimodales. La correlaci´on cruzada, la suma de los cuadrados de las diferencias de intensidad y la raz´on de uniformidad de la imagen se utilizan para el registro de im´agenes en la misma modalidad. 2. Trabajo realizado 2.1. Primera aproximaci´on (basada en caracter´ısticas) En esta aproximaci´on al problema se implement´o un m´etodo basado en features SURF[4] para encontrar los descriptores de ambas im´agenes y emparejarlos seg´un su similitud. Luego mediante un procedimiento estad´ıstico de RANSAC se descartan los emparejamientos espurios al mismo tiempo que se calcula la homograf´ıa con el conjunto resultante de inliers. Una vez que tenemos la transformaci´on proyectiva entre ambas im´agenes podemos alinearlas o georreferenciarlas. Para realizar este prototipo se ha utilizado como lenguaje C++ y la librer´ıa de visi´on por computador OpenCV. Se ha hecho uso de su funcionalidad para poder extraer los descriptores, emparejarlos y encontrar la transformaci´on realizando el proceso de RANSAC de manera impl´ıcita en la funci´on que encuentra la homograf´ıa. Como referencia para determinar el algoritmo y que funciones utilizar se ha trabajado con los libros Learning OpenCV: computer vision with the OpenCV library [5], OpenCV 2 Computer Vision Application Programming 4 Figura 1: Diferencia de escala Figura 2: Diferencia de escala y rotaci´on Cookbook [6], Computer Vision: Algorithms and Applications [7], y el art´ıculo Image alignment and stitching: A tutorial [8] Una vez se tiene el prototipo implementado se realizan una serie de pruebas con distintas im´agenes para evaluar su comportamiento y en el caso de resultar satisfactorio refinar utilizando m´etodos m´as sofisticados y eficientes. En la figura 1 se aprecia como var´ıan los resultados conforme se hace un cambio de escala. Si ambas im´agenes tienen un factor de escala similar se consiguen muchos emparejamientos correctos y el RANSAC es capaz de desestimar los incorrectos de manera adecuada. Conforme aumenta la diferencia de escala se van obteniendo menos emparejamientos hasta que llega un momento en el que ya no hay emparejamientos con los que trabajar. En la figura 2 se aprecia el mismo proceso de variaci´on en escala pero con el a˜nadido de una rotaci´on de −30◦, el efecto de esta rotaci´on es que la degradaci´on de los emparejamientos se produce antes. En la figura 3 se aprecia el proceso de alineamiento seg´un aumenta la diferencia en perspectiva de ambas im´agenes a una escala similar. En general se obtienen buenos resultados incluso con variaciones bastante grandes. En este caso se aprecia mejor que se trata de una transformaci´on proyectiva entre ambas im´agenes, provocando deformaciones no lineales que son m´as complicadas de emparejar. En las im´agenes donde ya no se puede encontrar una transformaci´on satisfactoria se aprecia que el RANSAC es un m´etodo de montecarlo y por lo 5 Figura 3: Diferencia de perspectiva 6 Figura 4: Casos reales tanto nos devuelve siempre un resultado, sea este correcto o no, entonces nos corresponde a nosotros discriminar verosimilitud del mismo. En las pruebas anteriores se han usado im´agenes obtenidas de la misma fuente. Por lo tanto, aunque sean diferentes tienen otras muchas cosas en com´un como por ejemplo la resoluci´on con la que han sido captadas, niveles de ruido similares pero sobre todo iluminaci´on y textura que son los par´ametros que m´as afectan a los descriptores empleados y en general a cualquier m´etodo. Por este motivo se han realizado una serie de pruebas con im´agenes reales, mostradas en la figura 4. En esta figura se puede observar como en los casos extremos donde hay mucha variaci´on de perspectiva adem´as de im´agenes de fuentes muy distintas no se obtiene un resultado positivo. El resultado negativo en el caso anterior era esperable puesto que se ha planteado un caso muy complicado al introducir una gran variaci´on en el punto de vista con el que est´an tomadas las im´agenes y que ambas proceden de fuentes e instantes distintos. Por ello, en la figura 5 se ha considerado un caso trivial pero con im´agenes obtenidas de distintas fuentes. Lo sorprendente de este caso es que el m´etodo no ha sido capaz de darnos un resultado positivo como se puede apreciar en la figura 6 Este resultado inesperado ha hecho que analicemos este caso particular con m´as detalle en la figura 7, para poder encontrar el origen del problema y poder presentar distintas l´ıneas de actuaci´on que permitan solventarlo. Lo que se aprecia en esta figura 7. El mismo recorte de la misma imagen obtenida de fuentes distintas a la misma escala, es que la resoluci´on es demasiado 7 las im´agenes en tiles del mismo tama˜no para poder calcular la funci´on χ2y se realiza el alineamiento af´ın de todas ellas. En la segunda etapa, dado que hemos conseguido un conjunto de par´ametros afines, se puede resolver un sistema de ecuaciones sobredeterminado para obtener los ocho par´ametros de la transformaci´on perspectiva buscada, siendo posible a˜nadir m´as restricciones a la soluci´on si tenemos en cuenta que conocemos los centros de dichas tiles. 2.3. Implementaci´on Este algoritmo, para hallar la transformaci´on perspectiva, se basa en particionar las im´agenes en un conjunto regular de tiles y alinear estas estimando ´unicamente la transformaci´on af´ın que las relaciona. De esta manera una vez se obtiene este conjunto de datos relacionados, se puede estimar la transformaci´on perspectiva que las relaciona. El n´umero de particiones es arbitrario y depender´a del tiempo que se est´e dispuesto a invertir en su resoluci´on ya que al aumentar el n´umero de tiles se aumenta el n´umero de emparejamientos necesarios y por lo tanto el coste computacional pero a cambio se obtiene una mayor precisi´on en la estimaci´on de los par´ametros de la transformaci´on. La primera fase de la implementaci´on ha consistido en construir el esquema exterior del algoritmo, de esta manera se pospone la implementaci´on del m´etodo de optimizaci´on hasta haber podido realizar ciertas pruebas de rendimiento y usabilidad. En esta parte se ha llevado a cabo la preparaci´on de los datos y se ha decidido como particionar el conjunto de tiles. En concreto se ha decidido no obtener el conjunto como tal, en su lugar se ha preferido ir trabajando sobre la regi´on correspondiente a cada tile en las propias im´agenes puesto que esto se puede realizar de manera eficiente sin duplicar informaci´on de manera innecesaria. Esta parte se ha realizado de tres formas diferentes (C++ con OpenCV, Pyhton con OpenCV y Python con PIL) con el objetivo de escoger la que mejor se adapte a nuestro prop´osito, v´ease el anexo 4.16 para m´as detalles de los lenguajes contemplados. Se ha partido de la base de utilizar C++ por ser el lenguaje de referencia en la mayor´ıa de los trabajos y la documentaci´on asociada que esto nos proporciona. Junto con OpenCV que nos proporciona un conjunto de funciones asociadas con el manejo de im´agenes, ya que esta librer´ıa se encarga de la carga de las mismas sin preocuparnos del formato en el que se encuentren y la forma de manejarlas es muy eficiente ya que puede trabajar sin duplicarlas al trabajar mediante referencias, junto con el conteo y gesti´on de las mismas. Esta caracter´ıstica se da incluso en el acceso a subregiones, lo que en nuestro caso se traduce en un acceso eficiente y en no necesitar estructuras de datos m´as complejas para manejarlas. Se ha mostrado inter´es en Python porque se trata de un lenguaje que ofrece una elevada productividad y expresividad, aunque esto haga que de manera inevitable sea m´as lento que otros lenguajes compilados. No obstante, si en un futuro fuera necesario aumentar la velocidad de ejecuci´on de este lenguaje, se podr´ıa optar por RPython (Restricted Python) que es un subconjunto de Python que se puede portar a c´odigo C de manera autom´atica y genera c´odigo eficiente. O bien por Cython que al convertirlo en un lenguaje tipado permite compilar Python a C y al mismo tiempo permite llamar funciones y m´etodos 14 de C/C++, lo que ofrece la posibilidad de optimizar los cuellos de botella del m´etodo sin necesidad de modificar toda la implementaci´on. Python cuenta con PIL(Python Imaging Library) que se trata de una extensi´on del lenguaje que ofrece primitivas b´asicas para trabajar con im´agenes. Aunque es bastante limitada en principio es suficiente para el prop´osito del algoritmo que se va a implementar puesto que las ´unicas funciones que vamos a delegar en ella son la carga de im´agenes y la interpolaci´on al reescalar. Al realizar la implementaci´on del esquema exterior de esta manera se ha observado que se trabaja de manera muy c´omoda porque el acceso a las subregiones se realiza de manera muy intuitiva y simple. El inconveniente encontrado al hacerlo de esta manera es que el tiempo de ejecuci´on mayor que el de C++ en exceso y es- to teniendo en cuenta que no se ha implementado la parte computacionalmente costosa del algoritmo. Esto ha hecho desestimar esta aproximaci´on, pero no el lenguaje Python que ha permitido prototipar la implementaci´on de una manera bastante r´apida. Para mantener el lenguaje y al mismo tiempo ganar eficiencia se ha hecho uso de la interfaz para Python que nos ofrece OpenCV. Esta interfaz la ofrece de manera nativa y est´a documentada junto a la documentaci´on de C++. Esto es posible por la forma de trabajar de OpenCV, ya que hace uso de referencias y maneja el conteo de las mismas de manera autom´atica, por lo que es compatible con la gesti´on de memoria din´amica de Python. Esta implementaci´on ha resultado ser similar en tiempo de ejecuci´on a la de C++, lo que hace tenerla en cuenta como muy buena opci´on para implementar todo tipo de m´etodos de manera r´apida. El inconveniente es que algunos detalles de la interfaz y la forma de trabajar de Python han hecho que para acceder a algunas propiedades de las im´agenes, en concreto las subregiones, haya sido necesario hacerlo de manera indirecta y de forma no muy intuitiva. Esto ha supuesto una p´erdida de tiempo considerable, lo que ha hecho replantear esta opci´on, por lo menos de momento, porque este inconveniente hace que resulte m´as r´apido el uso de un lenguaje al que se est´a m´as acostumbrado como C++. Una vez que se ha decidido el lenguaje se ha procedido al estudio minucioso del art´ıculo para su implementaci´on. Para ello ha sido necesario recurrir a bibliograf´ıa sobre m´etodos matem´aticos aplicados en el ´ambito de la visi´on por computador [13] [14] [15] [16] [17]. Esto nos ha dotado de las herramientas necesarias para poder comprender el funcionamiento del m´etodo y para poder implementarlo de manera adecuada. Por lo tanto se ha comenzado por implementar el m´etodo de optimizaci´on no lineal LMA (Levenberg-Marquardt Algorithm), que es una combinaci´on del m´etodo de Newton-Raphson y el de descenso del gradiente. La implementaci´on del m´etodo LMA ha sido relativamente r´apida porque en el art´ıculo se especifican todos los detalles de manera adecuada. S´olo se ha encontrado alg´un inconveniente al tener que derivar la funci´on que mide la similitud entre las im´agenes, puesto que se trata de una funci´on multidimensional. Y al elegir como resolver el sistema de ecuaciones lineales porque sobre esto hay demasiados m´etodos para elegir dependiendo de las condiciones concretas de las ecuaciones. En este caso no son ecuaciones que requieran un tratamiento muy complejo y se ha optado por resolverlo usando la pseudoinversa, porque no hace falta calcular la inversa y as´ı se evitan posibles problemas. Como paso previo a construir el resto del algoritmo e integrarlo con ´el, se 15 han realizado pruebas para evaluar su correcto funcionamiento. De estas pruebas no se han obtenido resultados satisfactorios, lo que ha conllevado una revisi´on de la soluci´on implementada. Al no encontrarse ning´un problema en la parte algor´ıtmica se han buscado los detalles de implementaci´on y en dicha revisi´on se ha observado que hab´ıa un error en el acceso a los datos de las im´agenes debido a la forma en que estas son almacenadas y como eran accedidas por los m´etodos. Una vez solucionados estos errores de implementaci´on se han vuelto a obtener resultados no satisfactorios, esto ha hecho necesaria una revisi´on muy detallada de todos los pasos del algoritmo. Lo primero en revisarse ha sido la resoluci´on del sistema de ecuaciones al ser lo que gobierna el funcionamiento del algoritmo. Para ello se han usado distintos m´etodos de resoluci´on de ecuaciones pero de todos ellos se ha obtenido el mismo resultado. Esto nos ha hecho centrar la atenci´on en los datos con los que se resuelve el sistema de ecuaciones, estos son la Hessiana(matriz de derivadas segundas) y B(vector de derivadas), y en su c´alculo se ha encontrado y arreglado alg´un error a causa de la precisi´on de los c´alculos pero no tan grave como para hacer fallar el algoritmo. Al no encontrar ning´un error en el c´odigo, se ha procedido a revisar las demostraciones matem´aticas en las que se basa el algoritmo. Tras un an´alisis detallado se ha llegado a la conclusi´on de que hab´ıa una errata. Esto hac´ıa que el sistema de ecuaciones no dependiera del par´ametro que disminuye o aumenta la resoluci´on de la soluci´on seg´un estemos m´as cerca o m´as lejos del m´ınimo buscado. Una vez solucionado este problema de dise˜no se ha seguido obteniendo una soluci´on similar. Esto ha conllevado a una revisi´on y verificaci´on de todo el c´odigo escrito y funciones utilizadas. Se ha verificado cualquier detalle que pudiera ser una fuente de errores, por ejemplo se ha comprobado que el paso de par´ametros se realizase de manera correcta y que se estuviera trabajando con la matriz de transformaci´on correcta y no con su traspuesta y se ha revisado cada una de las funciones y par´ametros utilizados en busca de cualquier posible error. Al no encontrar la causa del funcionamiento incorrecto del algoritmo, se ha vuelto a revisar como se calculan los datos m´as relevantes para el algoritmo. Como el c´alculo de la Hessiana y el vector B depende de la transformaci´on de la imagen, se ha revisado este aspecto. Hasta este momento se hab´ıa usado una funci´on de OpenCV de transformaci´on perspectiva porque la de transformaci´on af´ın es un subconjunto de esta y al no variar los par´ametros de perspectiva la transformaci´on obtenida es af´ın igualmente y por lo tanto tambi´en se han hecho pruebas para verificar que esto era correcto. Al ser esta la ´unica funci´on que se ha usado de OpenCV (a excepci´on de cargar las im´agenes y mostrarlas) se ha pretendido prescindir de ella y por lo tanto se ha implementado una funci´on de transformaci´on propia. El problema que hemos encontrado aqu´ı es que la interpolaci´on implementada ha sido menos refinada que la que realiza OpenCV de manera interna y esto ha resultado en un empeoramiento del resultado. Al no disponer de una forma alternativa adecuada para calcular dicha transformaci´on, se ha modificado como se calcula la Hessiana y el vector B para que no fuesen dependientes de esta transformaci´on por si esta no estuviese haciendo su funci´on de manera adecuada. Se han construido casos de prueba triviales para verificar el comportamiento del algoritmo frente a ellos pero tambi´en ha fallado, aunque esto nos ha ofre- 16 Figura 12: Variaci´on del resultado en funci´on del incremento/decremento La imagen central se corresponde a I2 transformada en I1 con la matriz de transformaci´on resultante. 17 cido alguna pista ya que se ha observado que no se dan los saltos del tama˜no suficiente para alcanzar la soluci´on y por lo tanto el algoritmo se queda estancado en un m´ınimo local. Por lo tanto se ha variado la precisi´on con la que se tratan las im´agenes y todos los datos intermedios pero tampoco se ha llegado a una mejora significativa. Siguiendo con el an´alisis de la traza se ha observado que el incremento hacia la soluci´on(la soluci´on del sistema de ecuaciones) iba descendiendo aunque no nos estuvi´eramos acercando a la soluci´on. Esto aunque contra intuitivo es un resultado que se ajusta a la descripci´on matem´atica del algoritmo, ya que de manera simplificada este par´ametro afecta “dividiendo” al sistema de ecuaciones y por lo tanto hace que cada vez la soluci´on sea menor. Como los ´unicos par´ametros con los que se puede ajustar el m´etodo son el tama˜no de aumento o disminuci´on del paso y los umbrales de terminaci´on del paso y del error, a lo que se ha procedido es a preparar una serie de pruebas con distintos valores para estos par´ametros. Lo que se ha observado es que el efecto de los umbrales en el resultado final ha sido m´ınimo pero no as´ı en el tiempo de ejecuci´on. El que mayor influencia ha tenido, como se muestra en la figura 12, es el incremento o decremento del paso. El m´etodo es muy dependiente de este valor pero no muestra una dependencia lineal, esto hace que los resultados sean demasiado variables. El principal inconveniente de esto es que lo hace muy dependiente de los datos concretos con los que se trabaja, lo que implica que aunque se consiguiera ajustar un par´ametro adecuado s´olo ser´ıa v´alido para esos datos. Adem´as, esto no nos proporciona informaci´on ´util para ajustar el par´ametro para otros datos. Al seguir sin encontrar una soluci´on satisfactoria se ha decidido intentar validar el algoritmo implement´andolo en matlab. De esta forma si los resultados obtenidos son distintos y correctos se encontrar´ıa de manera inmediata donde se est´a cometiendo el error. Mientras que si se obtienen los mismos resultados se podr´a descartar un problema de implementaci´on y en todo caso denotar´ıa un problema en el dise˜no o interpretaci´on del algoritmo. Esto podr´ıa ser posible porque a los autores del art´ıculo se les hubiera olvidado incluir alguna etapa necesaria para el correcto funcionamiento del m´etodo. Para volver a implementar este algoritmo en matlab ha hecho falta muy poco tiempo puesto que ya se ten´ıa el conocimiento necesario adquirido al haberlo implementado en otro lenguaje. Aunque los resultados obtenidos por esta nueva implementaci´on han seguido siendo negativos, se ha encontrado que en el art´ıculo hab´ıa una interpretaci´on ambigua referida al comienzo de coordenadas, aunque ha resultado ser irrelevante para los resultados. Al no encontrar errores ni en la implementaci´on ni en dise˜no del m´etodo se han estudiado alternativas a la resoluci´on del mismo pero de tal manera que estas no modifiquen la forma b´asica de funcionamiento para as´ı poder seguir utilizando el mismo esquema general. Se han estudiado m´etodos de premio-castigo y otros muchos m´etodos de optimizaci´on no lineal, pero todos ellos requer´ıan modificaciones demasiado importantes en el esquema del algoritmo y esto habr´ıa supuesto volver a realizar todo el trabajo desde el principio. Por lo tanto como alternativa viable se ha optado por intentar resolverlo mediante las herramientas de optimizaci´on que ofrece matlab. El inconveniente de usar estas funciones es que est´an pensadas para resolver un amplio espectro de problemas, esto hace que ofrezcan mucha m´as funcionalidad de la que se necesita en este caso concreto y a su vez demasiados par´ametros para ajustar y seleccionar el comportamiento 18 Figura 13: Resultados con fminsearch de los m´etodos ofrecidos. Al final se ha optado por utilizar fminsearch que encuentra el m´ınimo de una funci´on de variables m´ultiples sin restricciones y sin necesidad de usar las derivadas de la misma. Los resultados usando fminsearch que se observan en la figura 13, son mejores que los que se hab´ıan obtenido con anterioridad. Esto es cierto en el caso de estas im´agenes sint´eticas con mucho contraste y de poco tama˜no, ya que en el caso de im´agenes de mayor tama˜no los resultados vuelven a ser inservibles. Esta experiencia no ha resultado del todo infructuosa ya que se ha observado una peque˜na mejor´ıa en el resultado y al analizar de manera concienzuda por que ha sucedido esto se ha concluido que la posible causa del mal funcionamiento del m´etodo sea debido al c´alculo incorrecto de las derivadas. Tras el an´alisis de las derivadas se ha llegado a la conclusi´on de que es posible que al realizar una proyecci´on se pierda informaci´on. Por este motivo se ha realizado una implementaci´on alternativa donde se trabaja en el espacio inverso4.9. Para ello ha sido necesario realizar el c´alculo de la inversa de la matriz con los par´ametros de transformaci´on de manera anal´ıtica. Una vez obtenida esta nueva matriz de transformaci´on se han calculado de manera anal´ıtica las derivadas (vector B) y las derivadas segundas (hessiana). Una vez resuelto este proceso, el m´etodo resultante es similar pero cambiando el espacio donde se buscan los valores transformados. Esta alternativa tampoco ha resultado en resultados satisfactorios, por lo que ha sido necesario seguir buscando formas alternativas de resolver el problema. Se sigue suponiendo que el problema son las derivadas, puesto que se ha vuelto a testar toda la implementaci´on en busca de posibles errores y para ello se ha realizado una implementaci´on que permite variar par´ametros ajenos al m´etodo en si pero con una fuerte influencia en los resultados. Cabe destacar el formato de trabajo en coma flotante de las im´agenes ya que son almacenadas como valores entre 0 y 255. En busca de errores en este aspecto se han realizado pruebas con todas las posibles combinaciones pero no se ha llegado a ninguna conclusi´on satisfactoria que la hip´otesis inicial con la que se trabajaba. Si bien en alg´un caso particular se han conseguido resultados m´as prometedores, estos no han sido consistentes al realizar m´as pruebas con otro datos. Las funciones tanto de OpenCV (warpPerspective() ywarpAffine()) como de Matlab (imtransform()) no est´an concebidas para tratar con este tipo de transformaciones proyectivas que son m´as comunes en el ´ambito de los gr´aficos por computador puesto que en este ´ambito si que es necesario realizar estas 19 proyecciones de un espacio 3D a un espacio 2D para obtener la vista de una c´amara. No obstante las funciones con las que se ha trabajado no lo soportan y devuelven un resultado t´ecnicamente correcto pero incorrecto en el ´ambito en el que se est´a trabajando. Esto es a causa de la forma en que trabajan estas funciones ya que para hacer su trabajo de una manera m´as eficiente trabajan con la inversa de la matriz de transformaci´on y en este caso esta matriz al ser singular no tiene inversa. Visto de otra manera el problema que tienen estas funciones es la interpretaci´on de la proyecci´on, que en el caso de estas funciones es “encoger” la imagen hasta un vector o punto con valor nulo y por lo tanto esto no aporta ninguna informaci´on. Se ha tratado de solucionar este inconveniente en la funci´on imtransform mediante el uso de la funci´on maketform con la que se pueden generar transformaciones arbitrarias que luego pueden ser interpretadas por imtransform. Esta funci´on trabaja con la inversa de la transformaci´on y por lo tanto ha sido necesario implementar una funci´on de transformaci´on inversa general y una funci´on de transformaci´on inversa para cada una de las derivadas de primer orden y de segundo orden, en funci´on de los 8 par´ametros de una transformaci´on perspectiva. En el caso que se est´a tratando es necesario disponer de esta informaci´on puesto que se trata de la informaci´on b´asica con la que trabaja el m´etodo LMA y esto ha hecho necesario implementar estas funciones. Por un lado hay que decidir como realizar la transformaci´on de una imagen de un espacio a otro, en este caso se ha decidido por simplicidad realizar la transformaci´on directa. Por otro lado ha sido necesario interpretar la proyecci´on, para ello se han estudiado ocho formas distintas de realizarla.4.10 Mientras tanto se han seguido buscando soluciones alternativas para el c´alculo de las derivadas y se ha buscado la manera de prescindir del c´alculo de la derivada segunda en el c´alculo de la hessiana. Esta manera aproximada de obtener la hessiana ha sido posible puesto que la informaci´on que aporta la derivada segunda no es representativa del problema en el que se est´a trabajando y en consecuencia o bien aporta ruido o bien un valor constante que no afecta a la soluci´on. No obstante este valor suele ser despreciable en comparaci´on con el que aporta la primera derivada, por lo tanto no se esperaba un cambio efectivo en el resultado y esto ha sido lo que se ha podido observar al ponerlo a prueba. Por otro lado, tambi´en como esper´abamos se han obtenido resultados menos variables en funci´on de los par´ametros del algoritmo pero igual de incorrectos. Otra soluci´on implementada ha consistido en realizar el c´alculo de las derivadas de manera num´erica como realiza matlab por defecto en la funci´on fminsearch, comentada anteriormente. Esta forma de proceder es m´as costosa al tener que realizar el c´alculo de la transformaci´on para el incremento de cada par´ametro y ofrece valores aproximados por lo que es mucho m´as dependiente de los datos concretos. En concreto es necesario realizar dos transformaciones por cada derivada de primer orden y de tres a cuatro transformaciones por cada derivada de segundo orden. Otro inconveniente que tiene esta forma de proceder es que hay que escoger de manera arbitraria un tama˜no de incremento para los par´ametros. Adem´as hay que tener en consideraci´on que el c´alculo de las derivadas no es el mismo en todos los casos y por lo tanto este incremento aunque proporcional tampoco es el mismo. Todo esto ha llevado a tener m´as par´ametros para modificar y ajustar, por ello se han realizado una serie de pruebas en fun- 20 ci´on del tama˜no del incremento para comprobar cual resultaba m´as adecuado. El problema encontrado con esta forma de proceder es que ha resultado en exceso dependiente del valor asignado al incremento. Como ejemplo, entre incrementos de tan solo una mil´esima se han obtenido resultados completamente diferentes y sin seguir ning´un patr´on de convergencia hacia una soluci´on com´un. Por otro lado en otros casos el valor del incremento ha sido demasiado peque˜no, como consecuencia el m´etodo no obten´ıa valores adecuados haciendo que el camino hacia la soluci´on fuese incorrecto. Se ha seguido buscando solucionar el problema con este ´ultimo m´etodo y para ello se ha decidido implementar una variante del mismo que tomase el valor del incremento de manera din´amica y para ello se han estudiado varias posibilidades. La primera implementada ha sido hacer que el tama˜no del incremento fuese proporcional al valor que tomase el par´ametro en el momento de la consulta. Seg´un fuentes consultadas a tal efecto se recomendaba un valor proporcional a 10−8que es el m´ınimo valor en el que se puede incrementar el 1 en coma flotante, aunque tambi´en se han probado otros valores mayores para que la variaci´on experimentada fuese mayor. No obstante el problema encontrado en esta implementaci´on ha sido que el valor de los par´ametros de transformaci´on era demasiado pr´oximo a cero, o cero, en la mayor´ıa de los casos y esto ha provocado demasiada poca variaci´on en el valor de las derivadas. Todo esto ha hecho que esta aproximaci´on lejos de conseguir resultados satisfactorios haya resultado ser contraproducente. En la misma l´ınea se ha contemplado la variaci´on independiente del incremento de cada par´ametro en funci´on del valor λque gobierna el funcionamiento del algoritmo LMA. Esta aproximaci´on ha resultado demasiado compleja al requerir almacenar una traza de cada par´ametro y tener en cuenta la variaci´on que representaba cada combinaci´on en una funci´on multivariable. Esto ha hecho que esta aproximaci´on se considerase demasiado esfuerzo y se desestimase antes de seguir la implementaci´on al no tener la demostraci´on matem´atica que garantizase su correcto funcionamiento ni la certeza o convicci´on de que esta forma de proceder fuese a reportar una soluci´on mejor que en caso anterior. Por lo tanto se ha continuado explorando la variaci´on del incremento de manera homog´enea para todos los par´ametros involucrados. A tal efecto y con la influencia del estudio anterior, se ha considerado conveniente la variaci´on del incremento en funci´on del par´ametro λ, que gobierna el m´etodo Levenberg- Marquardt. Esta aproximaci´on hace que conforme el m´etodo se acerca a un m´ınimo el intervalo o incremento en el que se calcula la derivada sea menor y por lo tanto m´as preciso. Mientras que cuando estamos lejos de la soluci´on o se considera que se est´a en un m´ınimo local, el valor del incremento aumenta permitiendo una mayor variaci´on en la derivada que hace dirigir la b´usqueda a mayor distancia ofreciendo la direcci´on a seguir en las sucesivas iteraciones. Esta vez los resultados han sido satisfactorios, como se puede observar en la figura 14 y por lo tanto se ha procedido a realizar m´as pruebas con distintos conjuntos de datos. Se ha comprobado en un caso anterior que un resultado prometedor en im´agenes de peque˜no tama˜no no significa que se vaya a obtener un resultado significativamente distinto en im´agenes de mayor tama˜no. Para ello se han utilizado dos im´agenes de 512x512 pixels, siendo una de ellas generada de manera sint´etica a partir de la otra para obtener un caso de prueba que tenga una posible soluci´on. El resultado de esta prueba se puede observar en la figura 21 Figura 14: Resultados con derivadas num´ericas y variaci´on del incremento en funci´on de λ (a) Datos sint´eticos (b) Datos reales Figura 15: Resultados con derivadas num´ericas y variaci´on del incremento en funci´on de λen im´agenes de mayor tama˜no 15a] que no ofrece la precisi´on que se hubiese deseado. No obstante si se tienen en consideraci´on los resultados obtenidos con los m´etodos anteriores se trata de un buen resultado. Para estar m´as seguros de que el m´etodo se va a comportar de manera correcta en la mayor´ıa de las situaciones a las que puede ser sometido se ha considerado oportuno, sin llegar a realizar un estudio exhaustivo, probarlo con datos reales. Para ello se han tomado las im´agenes que han originado la b´usqueda alternativa a los m´etodos basados en caracter´ısticas y los resultados se pueden observar en la figura 15b, que vuelven a obtener un nivel de precisi´on similar al anterior y por lo tanto se puede considerar como satisfactorio. El inconveniente de esta aproximaci´on es el coste computacional, como se ha anticipado anteriormente. Si bien en los m´etodos en los que las derivadas se calculaban de manera anal´ıtica y luego se evaluaba su valor el tiempo de ejecuci´on era alto pero no excesivo, en este caso ha resultado muy elevado teniendo en cuenta que las im´agenes eran de s´olo 512x512 pixels y ha costado de media unas 22h de c´alculo de media. Esta diferencia tan abrumadora hace esperar que se siga una progresi´on cuadr´atica o c´ubica y en el ´ambito de im´agenes mayores resultar´ıa en un tiempo de c´alculo inaceptable. 22 3. Conclusiones y futuras l´ıneas de investigaci´on Como lenguaje de implementaci´on, Matlab ha demostrado ser el m´as pr´actico por su rapidez de prototipado. En su mayor parte esto ha sido as´ı porque la funcionalidad necesaria ha sido implementada y no se han necesitado funciones implementadas por terceros. Por otro lado C++ junto a OpenCV se han mostrado mucho m´as vers´atiles porque este ´ultimo ofrece un conjunto de funcionalidad muy amplio en este dominio de aplicaci´on, aunque este s´olo haya sido utilizado en el prototipo inicial basado en caracter´ısticas por tener una implementaci´on de SURF. El m´etodo basado en caracter´ısticas SURF ha obtenido unos resultados aceptables en dominios controlados. El problema de este tipo de m´etodos es que como ha ocurrido en el caso estudiado, un caso que se consideraba trivial no ha sido capaz de resolverlo e incluso el m´etodo que ten´ıa disponible el experto lo ha conseguido de manera forzada ajustando los par´ametros. No obstante visto los resultados de las otras t´ecnicas habr´ıa que seguir explorando en esta l´ınea. Para ello ser´ıa necesario construir algoritmos m´as robustos y buscar descriptores que sean m´as adecuados para este ´ambito de aplicaci´on. Tambi´en habr´ıa que buscar un m´etodo de validaci´on del resultado estimado para obtener un intervalo de confianza donde poder estimar la verosimilitud del resultado obtenido de manera autom´atica. Esto es necesario puesto que en los casos probados se depende de la validaci´on visual del resultado a excepci´on de casos muy evidentes. La gran ventaja de estos m´etodos es la posibilidad de clasificar los descriptores mediante un k-d tree oScalable Vocabulary Tree (SVT), seleccionar las m´as representativas con un algoritmo tipo PCA y construir una tabla de b´usqueda invertida, de esta manera se pueden obtener de manera muy r´apida los posibles candidatos y realizar un ajuste posterior. Por otro lado el m´etodo implementado, basado en la intensidad, ha ofrecido unos resultados aceptables pero no con la precisi´on esperada. Esto unido al alto coste de c´alculo necesario para llegar a este resultado y sobre todo teniendo en cuenta la progresi´on observada, hace que no resulte viable aplicarlo en casos reales donde se espera trabajar con tama˜nos de imagen mucho mayores a los empleados en nuestros casos de prueba. Para seguir esta l´ınea habr´ıa que implementar el algoritmo descrito en Image Registration Using Log-Polar Mappings for Recovery of Large-Scale Similarity and Projective Transformations [1], ya que debido a que trabaja en el dominio log-polar y las optimizaciones a˜nadidas sobre el m´etodo que se ha implementado promete unos resultados de mayor calidad a los obtenidos. Si esta opci´on se mostrase viable se podr´ıa tener en cuenta una aproximaci´on h´ıbrida, donde un m´etodo basado en caracter´ısticas ofreciese una transformaci´on aproximada y se refinara haciendo uso de este m´etodo. Pues- to que este m´etodo estimar´ıa mejor la transformaci´on una vez se est´a cerca de la soluci´on, pero su coste tambi´en ser´ıa demasiado elevado como para utilizarlo como primera opci´on. Otra posible l´ınea de investigaci´on, relacionada tanto con los m´etodos basados en descriptores pero quiz´as m´as relevante en el caso de los m´etodos basados en intensidad, se trata de mejorar la medida de similitud entre las im´agenes. Esta se ha mostrado poco informativa a la hora de ofrecer indicios de como llegar a la soluci´on. En la parte que s´olo intervienen transformaciones afines, se puede usar el domino log-polar o la transformaci´on al espacio de frecuencias. 23 Transformaciones r´ıgidas Son necesarios 3 par´ametros para definirlas y son traslaci´on y rotaci´on. Se mantienen los ´angulos entre las l´ıneas rectas. Transformaciones afines Son necesarios 6 par´ametros para definirlas y son traslaci´on, rotaci´on, escala y cizalla. Si dos rectas son paralelas lo seguir´an siendo en el espacio transformado. Transformaciones perspectivas Son necesarios 8 par´ametros para definirlas. Se trata de un superconjunto de la transformaci´on af´ın pero es una transformaci´on no lineal que no garantiza que dos rectas paralelas sigan si´endolo en el espacio transformado. La matriz que representa esta transformaci´on se suele llamar homograf´ıa. Transformaciones locales no r´ıgidas En este caso no existe n´umero de par´ametros fijo puesto que se trata de un modelo de deformaci´on que puede ser diferente para cada punto de la imagen. Estos modelos son necesarios cuando se trata de obtener una correspondencia densa para tareas como la fotogrametr´ıa, reconstrucci´on de informaci´on 3D a partir de 2 o m´as im´agenes. En el caso que se contempla se van a alinear im´agenes entre las que existe una ´unica transformaci´on perspectiva. Se puede suponer que s´olo hace falta una transformaci´on perspectiva para alinearlas puesto que las im´agenes estar´an tomadas desde una distancia (z) mucho mayor que las coordenadas (x,y) y por tanto se asemejar´an a im´agenes planas. Al tratarse de im´agenes planas se garantiza que s´olo sea necesaria una ´unica transformaci´on perspectiva. No obstante habr´ıa que tener en cuenta aberraciones causadas por el sensor o cualquier otro tipo de deformaci´on que complicar´ıa el modelo, pero se puede suponer que estas deformaciones no son relevantes puesto que no se requiere una correspondencia densa. Por tanto la funci´on de mapeo para una transformaci´on perspectiva ser´a la siguiente u=a0x+a1y+a2 a6x+a7y+ 1 v=a3x+a4y+a5 a6x+a7y+ 1 Esto define la transformaci´on, en coordenadas homog´eneas, T{I}como   u0 v0 w0  =  a0a1a2 a3a4a5 a6a71    x y 1     u v 1  =   u0 w0 v0 w0 w0 w0    30 Mientras que para una transformaci´on af´ın a6=a7= 0 dando como resultado u=a0x+a1y+a2 v=a3x+a4y+a5 y su transformaci´on correspondiente TA{I}   u v 1  =  a0a1a2 a3a4a5 001    x y 1   Como criterio para medir la similitud entre las im´agenes se utilizar´a la suma de las diferencias al cuadrado o sum of squared differences (SSD) que en un dominio bidimensional se define de la siguiente manera χ2(a) = ZZ u⊂R2 [I1(u)−I0 2(u)]2du =ZZ u⊂R2 [I1(u)−T{I2(x)}]2du =kI1(u)−T{I2(x)}k2 donde T es una transformaci´on geom´etrica aplicada a la imagen I2para mapearla desde su sistema de coordenadas (x,y) al sistema de coordenadas (u,v) de I1 El caso anterior trata un modelo de transformaci´on continuo pero en el caso en el que se trabaja utiliza datos discretos y por tanto un modelo discretizado para dichos datos es el siguiente: χ2(a) = N X i=1 [I1(ui, vi)−I2(ui, vi)]2 = N X i=1 [I1(ui, vi)−T{I2(xi, yi)}]2 = N X i=1 I1(ui, vi)−I2a0xi+a1yi+a2 a6xi+a7yi+ 1 ,a3xi+a4yi+a5 a6xi+a7yi+ 1 2 que en el caso de tratarse de un modelo de transformaci´on af´ın, TA, la discretizaci´on del modelo resultar´a de la siguiente forma 31 χ2(a) = N X i=1 [I1(ui, vi)−I2(ui, vi)]2 = N X i=1 [I1(ui, vi)−TA{I2(xi, yi)}]2 = N X i=1 [I1(ui, vi)−I2(a0xi+a1yi+a2, a3xi+a4yi+a5)]2 4.5. Newton-Raphson El m´ınimo de una funci´on se da en ellos puntos donde la primera derivada de la funci´on es cero. Este se describe para el caso af´ın, puesto que ha sido usado en dicho caso. Por tanto hay que resolver para: Bk(a) = ∂χ2(a) ∂ak =−2 N X i=1 [I1(ui)−I0 2(ui)]∂I0 2(ui) ∂ak= 0 donde I0 2(ui) = TA{I2(xi)}yk= 0,1,...,5 Esto supone un sistema de seis ecuaciones no lineales Bk= 0 k= 0,1,...,5 La parte lineal de la serie de Taylor es la siguiente Bk(a+ ∆a)≈Bk(a) + 5 X l=0 ∂Bk(a) ∂al ∆al Seis inc´ognitas y seis ecuaciones lineales dan como resultado B(a+ ∆a) = B(a) + H(a)∆a Cuando (a+ ∆a) sea la soluci´on, entonces B(a+ ∆a) = 0 y por tanto H(a)∆a=−B(a) donde H(a) es la hessiana, una matriz 6x6 de derivadas de segundo orden sim´etrica, hkl =hlk hkl =∂Bk(a) ∂al =∂2χ2(a) ∂ak∂al = 2 N X i=1 ∂I0 2(ui) ∂ak ∂I0 2(ui) ∂al −[I1(ui)−I0 2(ui)]∂2I0 2(ui) ∂ak∂al 32 Este sistema de ecuaciones se puede resolver para los seis valores de ∆amediante el m´etodo de Gauss, este m´etodo es adecuado porque converge de manera r´apida. Sus inconvenientes son la dificultad de encontrar buenos par´ametros iniciales y la convergencia lenta u oscilaciones si la funci´on no est´a bien balanceada. Para evitar estas oscilaciones, se puede realizar una aproximaci´on de la misma de tal manera que se desprecien los t´erminos relacionados con la segunda derivada. hkl =∂Bk(a) ∂al =∂2χ2(a) ∂ak∂al = 2 N X i=1 ∂I0 2(ui) ∂ak ∂I0 2(ui) ∂al 4.6. Descenso del gradiente En cualquier minimizaci´on iterativa, la iteraci´on (i+1) est´a relacionada con la iteraci´on i por la ecuaci´on ai+1 =ai+ ∆ai Se denomina {ai+1}una secuencia de descenso si χ2(aai+1 < χ2(a)). Una direcci´on que con seguridad produce un descenso es la direcci´on del gradiente negativo. Por tanto ai+1 =ai−C∇χ2(a) donde C controla el tama˜no del paso a lo largo de cada par´ametro en ∆ai=ai+1 −ai =−C∇χ2(a) =−CB(a) La ventaja de este m´etodo es que siempre converge hasta un m´ınimo local. Por otro lado su desventaja es su convergencia lenta hacia el mismo. 4.7. LMA(Newton-Raphson y Descenso del gradiente) Este m´etodo se basa en la combinaci´on del m´etodo de Newton-Raphson y el m´etodo del descenso del gradiente. De Newton-Raphson se tiene H(a)∆a=−B(a) y del descenso del gradiente se tiene ∆a=−CB(a) Si se define 33 Ck=1 λhkk entonces se tiene λhkk∆a=−B(a) Con esto, Marquardt prob´o que el siguiente m´etodo h´ıbrido minimiza χ2(a) de manera m´as robusta que cualquiera de ambos por separado. [H(a) + λI]∆a=−B(a) donde el par´ametro λ≥0 controla el grado en el que el m´etodo de comporta como el de Newton-Raphson o como el del descenso del gradiente. En el caso de que la actualizaci´on anterior haya conseguido reducir χ2(a), entonces λes reducido para realizar la actualizaci´on del m´etodo de Newton-Raphson. En el caso contrario, aumenta χ2(a), se incrementa λpara realizar la actualizaci´on del descenso del gradiente. Por tanto, el m´etodo resultante , conocido como Algoritmo Levenberg-Marquard (LMA), es el que se describe a continuaci´on. Este algoritmo requiere dos umbrales T1yT2para especificar las condiciones de terminaci´on. El algoritmo terminar´a cuando el cambio de χ2(a) caiga por debajo de T1o bien cuando λsupere el umbral T2 begin I n i c i a l i z a r l o s par´ametros como la matriz identidad ( l a transformaci´on obtenida tambi´en es la ide ntidad ) Inicializar λcon un v a l o r dado , por ejemplo λ= 0,001 mientras que ( |∆χ2(a)|> T1| | λ < T2) Calcular l a Hessiana Calcul ar e l v ect or B Resolver e l sistema l i n e a l (H(a) + λI)∆a=−Bpara ∆a Evaluar χ2(a+ ∆a) s i χ2(a+ ∆a)< χ2(a) entonces λ←λ/10 a←a+ ∆a s i n o λ←λ∗10 fin si fi n m ie n tr a s q u e 4.8. Estimaci´on perspectiva mediante aproximaciones afines En este caso se extiende la alineaci´on geom´etrica af´ın de tal manera que pueda tratar con transformaciones perspectivas. Para ello hace falta estimar los ocho par´ametros que definen dicha transformaci´on perspectiva. Se puede proceder de manera directa como en el caso anterior, el problema es que si la posici´on inicial no es muy similar a la soluci´on es muy f´acil quedarse estancado en un m´ınimo local. Esto es mucho m´as grave con las transformaciones perspectivas porque se trata de transformaciones no lineales. Para tratar de evitar este problema se intenta aproximar la transformaci´on perspectiva por transformaciones 34 afines en distintos trozos de la imagen. De esta manera se utiliza la estimaci´on af´ın, que es m´as robusta, y se estima luego la transformaci´on perspectiva en funci´on de estos resultados. Para poder hacer esto, se hace uso de una aproximaci´on af´ın local entorno a un punto mediante la serie de Taylor de primer orden, y se define como u=U(x, y) =U(x0, y0) + ∂U(x0, y0) ∂x (x−x0) + ∂U(x0, y0) ∂y (y−y0) (4.8.1) =A0x+A1y+A2 v=V(x, y) =V(x0, y0) + ∂V (x0, y0) ∂x (x−x0) + ∂V (x0, y0) ∂y (y−y0) (4.8.2) =A3x+A4y+A5 donde A0=∂U(x0,y0) ∂x A1=∂U(x0,y0) ∂y A2=U(x0, x0)−A0x0−A1y0 A3=∂V (x0,y0) ∂x A4=∂V (x0,y0) ∂y A5=V(x0, x0)−A3x0−A4y0 A continuaci´on se va a mostrar como es el proceso para determinar la transformaci´on af´ın entre dos tiles y su uso para inferir los par´ametros de la transformaci´on perspectiva. Se resume la transformaci´on af´ın con la transformaci´on de las cuatro esquinas, de esta manera se reduce el n´umero de puntos a tener en cuenta y se evitan errores en los datos. Figura 16: Mapeo de 4 esquinas Aunque la transformaci´on perspectiva se puede estimar de manera directa una vez conocida la correspondencia entre las parejas cada una de las cuatro esquinas, en este caso se trata de encontrar la mejor transformaci´on af´ın que aproxima este mapeo. Dadas las cuatro esquinas de una tile en la imagen I2y sus correspondientes en la de referencia I1, se puede resolver el mejor ajuste de la transformaci´on af´ın usando un ajuste por m´ınimos cuadrados. Por tanto se puede definir el mapeo de la siguiente manera u=a0x+a1y+a2 v=a3x+a4y+a5 Se resuelve para los par´ametros afines mediante la minimizaci´on de la expresi´on de χ2mostrada a continuaci´on 35 χ2(a) = 4 X i=1 (ui−a0xi−a1yi−a2)2+ (vi−a3xi−a4yi−a5)2 Se pueden relacionar estas correspondencias de la forma U=W A             u1 u2 u3 u4 v1 v2 v3 v4             =             x1y11 0 0 0 x2y21 0 0 0 x3y31 0 0 0 x4y41 0 0 0 0 0 0 x1y11 0 0 0 x2y21 0 0 0 x3y31 0 0 0 x4y41                     a0 a1 a2 a3 a4 a5         y para terminar de inferir los seis par´ametros afines se calcula la soluci´on con el m´etodo de la pseudoinversa A= (WTW)−1WTU La estimaci´on anterior es para una pareja de tiles, pero lo que se busca es inferir los par´ametros afines para cada una de las N parejas en I1eI2y aplic´arsela a ellas en el centro de cada tile. Esto establece la correspondencia de N puntos entre I1eI2, haciendo que los ocho par´ametros de la transformaci´on perspectiva puedan ser inferidos resolviendo el siguiente sistema de ecuaciones u=a0x+a1y+a2−a6ux −a7uy v=a3x+a4y+a5−a6vx −a7vy que define el siguiente sistema de ecuaciones sobredeterminado       ui 0 . . . vi 0 . . .      2N×1 =       xi 0yi 01 0 0 0 −ui 0xi 0−ui 0yi 0 . . . 0 0 0 xi 0yi 01−vi 0xi 0−vi 0yi 0 . . .      2N×8             a0 a1 a2 a3 a4 a5 a6 a7            8×1 (4.8.3) donde (ui 0, vi 0)y(xi 0, yi 0) representan los centros de la tile ien I1eI2para 1≤i≤N. De manera similar al caso af´ın, se calcula la soluci´on utilizando el m´etodo de la pseudoinversa para calcular los ocho par´ametros de la transformaci´on perspectiva. A= (WTW)−1WTU 36 Estos resultados se pueden mejorar si se a˜naden restricciones adicionales a la estimaci´on de los m´ınimos cuadrados. La restricci´on que se puede imponer es la correspondencia entre los centros de las tiles. De las ecuaciones 4.8.1 y 4.8.2 se ten´ıa que la aproximaci´on a la transformaci´on af´ın es de la forma u=A0x+A1y+A2 v=A3x+A4y+A5 donde A0=∂U(x0, y0) ∂x =a0(a6x0+a7y0+ 1) −a6(a0x0+a1y0+a2) (a6x0+a7y0+ 1)2 =a0−a6u0 a6x0+a7y0+ 1 =a0−a6(A0x0+u0)−a7A0y0 A1=∂U(x0, y0) ∂y =a1(a6x0+a7y0+ 1) −a7(a0x0+a1y0+a2) (a6x0+a7y0+ 1)2 =a1−a7u0 a6x0+a7y0+ 1 =a1−a6A1x0−a7(A1y0+u0) A2=U(x0, y0)−A0x0−A1y0=a0x0+a1y0+a2 a6x0+a7y0+ 1 −A0x0−A1y0 =a0x0+a1y0+a2−a6u0x0−a7u0y0−(u0−A2) A3=∂V (x0, y0) ∂x =a3(a6x0+a7y0+ 1) −a6(a3x0+a4y0+a5) (a6x0+a7y0+ 1)2 =a3−a6v0 a6x0+a7y0+ 1 =a3−a6(A3x0+v0)−a7A3y0 A4=∂V (x0, y0) ∂y =a4(a6x0+a7y0+ 1) −a7(a3x0+a4y0+a5) (a6x0+a7y0+ 1)2 =a4−a7v0 a6x0+a7y0+ 1 =a4−a6A4x0−a7(A4y0+v0) A5=V(x0, y0)−A3x0−A4y0=a3x0+a4y0+a5 a6x0+a7y0+ 1 −A3x0−A4y0 =a3x0+a4y0+a5−a6v0x0−a7v0y0−(v0−A5) Hay que tener en cuenta que los t´erminos de A2yA5generan ecuaciones de la forma u0=a0x0+a1y0+a2−a6u0x0−a7u0y0yv0=a3x0+a4y0+a5−a6v0x0− a7v0y0respectivamente. Estas son las mismas correspondencias definidas en la ecuaci´on 4.8.3. De forma compacta se puede expresar las ecuaciones anteriores para obtener la relaci´on de los par´ametros afines en cada una de las N tiles con 37 los par´ametros desconocidos de la transformaci´on perspectiva buscada.                               Ai 0 . . . Ai 1 . . . ui 0 . . . Ai 3 . . . Ai 4 . . . vi 0 . . .                              6N×1 =                               1 0 0 0 0 0 −(Ai 0xi 0+ui 0)−Ai 0yi 0 . . . 0 1 0 0 0 0 −Ai 1xi 0 −(Ai 1yi 0+ui 0) . . . xi 0yi 01 0 0 0 −ui 0xi 0 −ui 0yi 0 . . . 0 0 0 1 0 0 −(Ai 3xi 0+vi 0)−Ai 3yi 0 . . . 0 0 0 0 1 0 −Ai 4xi 0 −(Ai 4yi 0+vi 0) . . . 0 0 0 xi 0yi 01−vi 0xi 0 −vi 0yi 0 . . .                              6N×8            a0 a1 a2 a3 a4 a5 a6 a7           8×1 (4.8.4) donde (ui 0, vi 0)y(xi 0, yi 0) representan los centros de la tile ien I1eI2para 1≤i≤N. Hay que tener en cuenta que la ecuaci´on 4.8.4 es un superconjunto de la ecuaci´on 4.8.3. Y de la misma manera que en el caso anterior, se resuelve utilizando el m´etodo de la pseudoinversa. Otra cosa a tener en cuenta si se utiliza la ecuaci´on 4.8.4 es que la soluci´on de la misma no es estable si se calcula en las coordenadas de los pixels de manera directa. Por lo tanto hay que realizar una normalizaci´on de las coordenadas de los pixels a un dominio de coordenadas [0,1] Con todo esto, se va a ver como es el algoritmo que resuelve esta aproximaci´on a la transformaci´on perspectiva. Consid´erese una imagen a alinear I2y una imagen de referencia I1, ambas divididas en tiles de forma regular y con las mismas dimensiones para poder calcular χ2. Para cada tile Ti 2se realiza el alineamiento geom´etrico af´ın usando el m´etodo Levenberg-Marquardt para encontrar la tile m´as similar de Ti 1. Esto ofrece una colecci´on de par´ametros afines y centros de tiles, con los que se estimar´an los par´ametros de la transformaci´on perspectiva. Particionar I2en N t i l e s : Ti 2, para 1 ≤i≤N i←1 mientras que ( i < N ) Seleccionar la tile Ti 2 para todas l a s p o s i c i o n e s (x, y)enI1 Seleccionar la tile Ti 1 Alinear Ti 2con Ti 1usando LMA s i χ2(a) es m´ınimo entonces xi←x yi←y ai←a fin si f i n p a r a i←i+ 1 fi n m ie n tr a s q u e Resolver l a ecuaci´on 4.8.3 o 4.8.4 para l o s ocho par´ametros de la transformaci´on p e rsp e c tiv a ( m´ınimos cuadrados l i n e a l e s ) 38 4.9. Transformaci´on af´ın inversa M´etodo eficiente del c´alculo de la inversa de una matriz 3x3 a−1=  a0a1a2 a3a4a5 a6a7a8   −1 =1 det(a)  A0A1A2 A3A4A5 A6A7A8   T =1 det(a)  A0A3A6 A1A4A7 A2A5A8   A0= (a4a8−a5a7)A3= (a2a7−a1a8)A6= (a1a5−a2a4) A1= (a5a6−a3a8)A4= (a0a8−a2a6)A7= (a2a3−a0a5) A2= (a3a7−a4a6)A5= (a6a1−a0a7)A8= (a0a4−a1a3) det(a) = a0A0+a1A1+a2A2 Transformaci´on af´ın ⇔a6=a7= 0 y a8= 1 a−1=  a0a1a2 a3a4a5 001   −1 =1 a0a4−a1a3  a4−a1a1a5−a2a4 −a3a0a2a3−a0a5 0 0 a0a4−a1a3   La siguiente expresi´on es la forma compacta de expresar las transformaciones afines porque cuando hace uso de ellas se a˜nade la ´ultima fila, que es siempre el mismo vector constante [0 0 1] a−1=1 a0a4−a1a3a4−a1a1a5−a2a4 −a3a0a2a3−a0a5 A continuaci´on se muestra el c´alculo de las derivadas de primer orden, que forman el vector B, en funci´on de cada par´ametro involucrado en el proceso de la transformaci´on af´ın. ∂a−1 ∂a0 =1 (a0a4−a1a3)2−a2 4a1a4a2a2 4−a1a4a5 a3a4−a1a3a1a3−a2a3 ∂a−1 ∂a1 =1 (a0a4−a1a3)2a4a3−a0a4a0a4a5−a2a3a4 −a2 3a0a3a2a2 3−a0a3a5 ∂a−1 ∂a2 =1 (a0a4−a1a3)20 0 −a4 0 0 a3 ∂a−1 ∂a3 =1 (a0a4−a1a3)2a1a4−a2 1a2 1a5−a1a2a5 −a0a4a0a1a0a2a4−a0a1a5 ∂a−1 ∂a4 =1 (a0a4−a1a3)2−a1a3a0a1a1a2a3−a0a1a5 a0a3a2 0a2 0a5−a0a2a3 ∂a−1 ∂a5 =1 (a0a4−a1a3)20 0 a1 0 0 −a0 39 distancia al mismo. Este muestreo no uniforme es simulado por la escala logar´ıtmica. Aunque las coordenadas log-polar no son exactamente lo mismo, est´an aceptadas como modelo de representaci´on de la retina de los primates. La transformaci´on al dominio log-polar tiene ventajas as´ı como inconvenientes, pero a continuaci´on se van a resaltar las dos ventajas m´as relevantes. Por un lado, las im´agenes resultantes son invariantes a rotaci´on y escala. Por otro lado, la variaci´on espacial del muestreo no uniforme de la retina es la soluci´on para reducir la cantidad de informaci´on que tiene que atravesar el nervio ´optico, manteniendo una alta resoluci´on en la f´ovea y al mismo tiempo un amplio campo de visi´on. Sean dos im´agenes I1eI2y sus transformaciones correspondientes en el dominio log-polar I1peI2p. Si I2es una versi´on rotada de I1, entonces I2p ser´a una versi´on trasladada en el eje θde I1p. Aplicando una correlaci´on cruzada de I1peI2pse puede encontrar el offset entre ambas im´agenes. Hay que resaltar que la t´ecnica de correlaci´on cruzada es usada de manera habitual para encontrar offsets de traslaci´on entre dos im´agenes y este no funciona bien en coordenadas cartesianas, en presencia de rotaciones o escalados. Sin embargo, se ha visto que en el espacio polar encontrar la componente de traslaci´on entre I1peI2p corresponde a encontrar la rotaci´on entre I1eI2. Para determinar el cambio de escala el comportamiento es similar. Consid´erese una ampliaci´on de I1de tal manera que sea cuatro veces mayor que I2. Entonces, todos los puntos (x, y) de I1se corresponder´an a los an´alogos (4x, 4y) en I2. Para determinar el factor de escala se hace uso de las propiedades de los logaritmos, puesto que en el espacio logar´ıtmico (x, y)→(log x, log y), y (4x, 4y)→(log 4x, log 4y)→(log x+ log 4,log y+ log 4). De este modo, se observa de manera expl´ıcita como en el espacio logar´ıtmico la introducci´on de un cambio de escala se corresponde con una traslaci´on en el eje logar´ıtmico (ro logr, seg´un la notaci´on empleada) de la imagen transformada. 4.13. Transformada de Fourier-Mellin El m´etodo de alineamiento geom´etrico de Fourier-Mellin se basa en la correlaci´on de la fase y las propiedades del an´alisis de Fourier. El m´etodo de la correlaci´on de la fase puede encontrar la traslaci´on entre dos im´agenes. Mientras que el m´etodo de Fourier-Mellin extiende la correlaci´on de la fase para poder alinear im´agenes que tienen tanto traslaci´on como rotaci´on [18] [19] [20] [21] [22] [23] [24]. De acuerdo con las propiedades de rotaci´on y traslaci´on de la transformada de Fourier, las transformaciones se representan de la siguiente forma I1(x, y) = I2(u, v) u=xcos θ0+ysin θ0−x0 v=−xsin θ0+ycos θ0−y0 F1(ωx, ωy) = F2(ωu, ωv)e−j(ωxx0+ωyy0) ωu=ωxcos θ0+ωysin θ0 ωv=−ωxsin θ0+ωycos θ0 46 Se observa que la magnitud o amplitud del espectro de frecuencias |F1|es una r´eplica rotada de |F2|y ambos espectros comparten el mismo centro de rotaci´on. Por lo tanto, se puede recuperar esta rotaci´on representando el espectro |F1|y |F2|en coordenadas polares |F1(r, θ|=|F2(r, θ −θ0| La amplitud de Fourier en coordenadas polares se diferencia s´olo por la traslaci´on. Esto hace que se pueda usar el m´etodo de correlaci´on de fase para encontrar esta traslaci´on y estimar θ0. Este m´etodo ha sido extendido para encontrar la escala por medio de mapear la amplitud de Fourier en coordenadas log-polar. Por lo tanto, se encuentra la escala y rotaci´on mediante una correlaci´on de fase, que obtiene la cantidad de desplazamiento en el espacio (log r, θ). La ventaja de este m´etodo es que es tolerante al ruido aditivo. Pero por otro lado s´olo puede obtener buenos resultados con cambios de escala y rotaci´on moderados. Esto es debido a que cuando se realiza un cambio de escala o rotaci´on importante, se generan efectos en los bordes que pueden alterar de manera dram´atica los coeficientes de Fourier [21] [25] [26]. 4.14. Pir´amide de multiples resoluciones Una pir´amide de resoluciones m´ultiples, o espacio de escalas, consiste en un conjunto de im´agenes que representan a una imagen en m´ultiples resoluciones. La imagen original se sit´ua en la base de la pir´amide y es reducida por un factor de escala constante en cada dimensi´on para generar el siguiente nivel. Esto es repetido nivel a nivel hasta que se alcanza la c´uspide de la pir´amide. Lo habitual es usar un factor de escala 2, y por tanto el tama˜no de la imagen en el nivel i es reducido respecto a la original en un factor de 2ien cada dimensi´on. Para que los niveles sean m´as significativos se suele llamar nivel m´as fino al nivel 0 y nivel m´as grueso al ´ultimo nivel. El uso de esta t´ecnica tiene dos ventajas principales. La primera es que cuando se aplica el m´etodo LMA al nivel m´as grueso de la pir´amide, el n´umero de pixels disminuye en un factor de 22(n−1). Esto conlleva a una gran mejora computacional puesto que la mayor´ıa de las iteraciones son ejecutadas en los niveles m´as gruesos. La segunda es que al disminuir el tama˜no de la imagen, esta queda suavizada puesto que es equivalente a realizar un desenfoque gaussiano. Este suavizado de la imagen provoca que χ2sea calculada en im´agenes m´as suaves y esta suavidad evita que la minimizaci´on quede atrapada en m´ınimos locales. Como en el nivel m´as grueso permanecen s´olo las caracter´ısticas de mayor tama˜no, el proceso de alineamiento va desde el nivel m´as grueso hasta el m´as fino de manera progresiva. De manera adicional, esta aproximaci´on en el nivel anterior es pasada al siguiente nivel como estimaci´on inicial de los par´ametros, aunque para ello es necesario escalar los par´ametros de manera correspondiente a lo largo de los niveles. Si el factor de escala entre los niveles es s:xi+1 =sxi yi+1 =syi 47 vi+1 =svi donde ui=ai 0xi+ai 1yi+ai 2 ai 6xi+ai 7yi+ 1 vi=ai 3xi+ai 4yi+ai 5 ai 6xi+ai 7yi+ 1 Sustituyendo las coordenadas del siguiente nivel m´as fino en las ecuaciones anteriores se obtiene ui+1 s=ai 0xi+1 s+ai 1 yi+1 s+ai 2 ai 6xi+1 s+ai 7 yi+1 s+ 1 vi+1 s=ai 3xi+1 s+ai 4 yi+1 s+ai 5 ai 6xi+1 s+ai 7 yi+1 s+ 1 Multiplicando ambos lados por sse obtiene ui+1 =ai+1 0xi+1 +ai+1 1yi+1 +ai+1 2 ai+1 6xi+1 +ai+1 7yi+1 + 1 =ai 0xi+1 +ai 1yi+1 +sai 2 ai 6xi+1 s+ai 7 yi+1 s+ 1 vi+1 =ai+1 3xi+1 +ai+1 4yi+1 +ai+1 5 ai+1 6xi+1 +ai+1 7yi+1 + 1 =ai 3xi+1 +ai 4yi+1 +sai 5 ai 6xi+1 s+ai 7 yi+1 s+ 1 y por lo tanto la relaci´on entre los par´ametros es ai+1 0=ai 0ai+1 1=ai 1ai+1 2=sai 2 ai+1 3=ai 3ai+1 4=ai 4ai+1 5=sai 5 ai+1 6=ai 6 sai+1 7=ai 7 s En el caso habitual, s= 2, por lo tanto los par´ametros de traslaci´on a2ya5 son multiplicados por 2 y los par´ametros a6ya7son divididos por 2. 48 4.15. Levenberg-Marquardt Optimizado Esta modificaci´on est´a basada en el art´ıculo de [27], donde el alineamiento es realizado en im´agenes m´edicas sujetas a transformaciones de rotaci´on, escala y traslaci´on. Por suerte, esta forma de proceder sigue siendo v´alida si se extiende al caso de transformaciones perspectivas. En el m´etodo LMA est´andar, se calcula en vector B8×1y la matriz Hessiana H8×8en cada iteraci´on. En realidad el c´alculo de la Hessiana es innecesario realizarlo en cada iteraci´on, pudiendo ser realizado una ´unica vez al comienzo. Consid´erese la siguiente funci´on objetivo que establece una medida de similitud entre I1eI2 χ2(a) = kI1(u)−I2(TA{x})k2 donde TArepresenta la transformaci´on que provoca la matriz de transformaci´on perspectiva 3 ×3 Se puede suponer que I2se transforma en I1tras una serie de transformaciones perspectivas A. Durante el proceso iterativo, las nuevas estimaciones de Ason calculadas de la siguiente manera Ai+1 =Ai+ ∆Ai donde Ai=I+ ∆A1+ ∆A2+. . . + ∆i−1. Como I2es transformada en cada iteraci´on, la matriz Hessiana necesita ser recalculada porque es una representaci´on del gradiente de I2y por lo tanto la matriz Hessiana es responsable del c´alculo de los t´erminos anteriores ∆A. El prop´osito de esta modificaci´on es eliminar el c´alculo de la matriz Hessiana. Esto se puede conseguir transformando el problema en otro donde I1sea transformada en I2sin tener que modificar I2entre iteraciones sucesivas. Esta aproximaci´on permite calcular la matriz Hessiana una ´unica vez, en la primera iteraci´on. Para poder determinar las nuevas estimaciones de los par´ametros de la transformaci´on perspectiva de Aen el m´etodo modificado, se debe expresar χ2en t´erminos de una transformaci´on tal que mapee I1en I2. La diferencia fundamental entre el m´etodo LMA est´andar y el modificado es que en el m´etodo est´andar se actualiza la estimaci´on actual mediante el cambio de los valores iniciales hacia el m´ınimo global, mientras que en el modificado se trae el m´ınimo global hacia el valor inicial. La consecuencia de esta definici´on se puede resumir con la siguiente regla del LMA modificado A−1 i+1 = (I+ ∆Ai)−1(I+ ∆Ai−1)−1. . . (I+ ∆A2)−1(I+ ∆A1)−1 = (I+ ∆Ai)−1A−1 i Una diferencia importante entre el m´etodo est´andar y el modificado es la forma en la que los par´ametros son actualizados en cada iteraci´on. En el m´etodo est´andar los par´ametros de inicio son elegidos como la matriz identidad Icomo 49 suposici´on inicial para el punto P0. A continuaci´on se calculan las derivadas direccionales de I2,gx=∂I2/∂x ygy=∂I2/∂y. Estos procesos son del orden de O(N×3×3), donde Nes el n´umero de pixels y 3 ×3 es el tama˜no del kernel o matriz de convoluci´on. El m´etodo est´andar proporciona un ∆Aque se usa para a˜nadirlo a la suposici´on inicial Ipara moverse desde el punto P0hasta el P1 en la curva de χ2(a). En la siguiente iteraci´on, a causa de que la imagen I2es transformada por A+ ∆A, es necesario volver a calcular gxygypara encontrar el nuevo ∆A. Por lo tanto en el m´etodo LMA est´andar, la soluci´on ´optima Pise desliza por la curva de χ2(a). Sin embargo, en el m´etodo modificado se cambia la curva de χ2(a) hacia la suposici´on inicial P0. Esto es conseguido remuestreando I1mediante la transformaci´on inversa (I+ ∆Ai)−1A−1 i. Por consecuencia, la imagen I0 1es “acercada” a I2. Ahora, la nueva imagen I0 1y la imagen I2son utilizadas para minimizar χ2(a). El resultado genera un nuevo ∆Aque es a˜nadido siempre a I, la suposici´on inicial en el punto P0. Puesto que I2no cambia, no es necesario calcular gxygy. Es ´util rescribir la funci´on χ2en t´erminos del nuevo mapeo directo, as´ı como el mapeo inverso. Esta descomposici´on va a permitir aplicar una parte importante de la transformaci´on a I1(u). Como resultado, la transformaci´on inversa restante es peque˜na y al aplicarla a I2(x) permite prescindir del c´alculo de la Hessiana. Sup´ongase que se descompone la transformaci´on T{} en dos transformaciones m´as peque˜nas TATI+∆A{}. De esta manera TA{} es la transformaci´on de la iteraci´on previa y TI+∆Aes la transformaci´on, peque˜na, que minimiza χ2(a) en el m´etodo Levenberg-Marquardt. χ2(a) = kI1(u)−I2(TA{TI+∆A{x}})k2(4.15.1) =1 |A|kI1(TA−1{u})−I2(TA−1{TA{TI+∆A{x}}})k2(4.15.2) =1 |A|kI1(TA−1{u})−I2(TI+∆A{x})k2(4.15.3) =1 |A(I+ ∆A)|kI1(T(I+∆A)−1A−1{u})−I2(x)k2(4.15.4) =1 |A(I+ ∆A)|kI1(T−1{u})−I2(x)k2(4.15.5) Estas ecuaciones muestran los pasos necesarios para transformar TATI+∆A{} desde el sistemas de coordenadas de I2en el sistema de coordenadas de I1con la normalizaci´on correspondiente. En vez de minimizar la funci´on 4.15.1 se va a minimizar la 4.15.3 respecto a los par´ametros δA. En el m´etodo LMA modificado es necesario derivar ∂I2(TI+∆A{x})/∂∆ay la regla de actualizaci´on para cada par´ametro de la transformaci´on. Se puede descomponer ∂I2(TI+∆A{x})/∂∆ade la siguiente manera ∂I2(TI+∆A{x}) ∂∆a=∂I2(TI+∆A{x}) ∂x ∂x ∂∆a y puesto que la transformaci´on ∆aes lo suficientemente peque˜na, se pueden realizar las siguientes aproximaciones 50 I2(TI+∆A{x})≈I2 ∂x ∂∆a≈∂u ∂∆a de esto se deduce que ∂I2(TI+∆A{x}) ∂∆a≈∂I2 ∂x ∂x ∂∆a≈∂I2 ∂x ∂u ∂∆a donde ∂I2 ∂x es el gradiente de I2. Entonces, ∂I2(x) ∂∆akpara los ocho par´ametros de la transformaci´on perspectiva resultan ∂I2 ∂∆a0 =∂I2 ∂x ∂u ∂∆a0 =gx x w ∂I2 ∂∆a1 =∂I2 ∂x ∂u ∂∆a1 =gx y w ∂I2 ∂∆a2 =∂I2 ∂x ∂u ∂∆a2 =gx ∂I2 ∂∆a3 =∂I2 ∂y ∂v ∂∆a3 =gy x w ∂I2 ∂∆a4 =∂I2 ∂y ∂v ∂∆a4 =gy y w ∂I2 ∂∆a5 =∂I2 ∂y ∂v ∂∆a5 =gy ∂I2 ∂∆a6 =∂I2 ∂x ∂u ∂∆a6 +∂I2 ∂y ∂v ∂∆a6 =−xgx a0x+a1y+a2 (a6x+a7y+ 1)2+gy a3x+a4y+a5 (a6x+a7y+ 1)2 ∂I2 ∂∆a7 =∂I2 ∂x ∂u ∂∆a7 +∂I2 ∂y ∂v ∂∆a7 =−ygx a0x+a1y+a2 (a6x+a7y+ 1)2+gy a3x+a4y+a5 (a6x+a7y+ 1)2 En el m´etodo LMA est´andar, la regla de actualizaci´on es la siguiente Anew =Aold + ∆A mientras que en el LMA modificado es Anew =Aold(I+ ∆A) =  a0a1a2 a3a4a5 a6a71  ×  ∆a0+ 1 ∆a1∆a2 ∆a3∆a4+ 1 ∆a5 ∆a6∆a71   A continuaci´on se muestra un esquema de como deber´a ser el Algoritmo Levenberg-Marquardt modificado 51 Con strui r l a s pir´amides de r e s o lu c i o n e s para I1eI2 I n i c i a l i z a r l o s par´ametros como la matriz identidad Inicializar λcon un v a l o r dado , por ejemplo λ= 0,001 para i=Lhasta 0 //L e l n i v e l m´as grueso de l a pir´amide Calcular l o s g ra di ent es d i r e c c i o n a l e s : ∂Ii 2 ∂x y∂Ii 2 ∂y Calcular l a matriz Hessiana 8 ×8 mientras que ( |∆χ2(a)|> T1| | λ < T2) Aplicar la transformaci´on TA−1sobre I1en e l n i v e l i Calcul ar e l v ect or B Resolver e l sistema l i n e a l (H(a) + λI)∆a=−Bpara ∆a Evaluar χ2(a+ ∆a) s i χ2(a+ ∆a)< χ2(a) entonces λ←λ/10 a←a+ ∆a s i n o λ←λ∗10 fin si fi n m ie n tr a s q u e f i n p a r a 4.16. Lenguajes contemplados 4.16.1. C++ Este lenguaje se posiciona como el mejor candidato por motivos de eficiencia y la disponibilidad de m´ultiples frameworks para la visi´on por computador, en especial OpenCV. La gran ventaja de este lenguaje es que la gran mayor´ıa de los desarrollos que existen en este ´ambito est´an realizados en C++, lo que proporciona una buena base para trabajar y documentaci´on disponible de la misma. Esto ´ultimo es lo que nos interesa de manera principal de este lenguaje, ya que salvo por el hecho de cargar las im´agenes y su manejo eficiente se pretende implementar el resto de funcionalidad, salvo funciones concretas utilizadas de manera habitual. 4.16.2. R Se trata de un lenguaje cuyo prop´osito es el an´alisis estad´ıstico y dispone de gran parte de la instrumentaci´on matem´atica necesaria en el ´ambito de la visi´on por computador, aunque no de manera espec´ıfica. Los inconvenientes encontrados son que tiene una sintaxis y modo de trabajo diferentes al que estamos acostumbrados y hay poca documentaci´on o trabajos realizados en este ´ambito. 4.16.3. Python Se trata de un lenguaje que permite un prototipado r´apido porque es de muy alto nivel e incorpora manejo de memoria din´amico entre otras cosas. Adem´as existe un wrapper para OpenCV que hace que trabaje de manera muy eficiente manteniendo el manejo din´amico de la memoria. 52 El inconveniente encontrado es que el wrapper tiene algunas particularidades que lo hacen diferir de como trabaja python y esto provoca que para poder utilizar algunas funcionalidades haya que hacerlo de manera especial y se pierda mucho tiempo en estos detalles, y este tiempo es el que se esperaba ganar al utilizarlo como lenguaje de prototipado. Pyhton tambi´en cuenta con una librer´ıa propia para tratamiento de im´agenes, PIL (Python Image Library), esta tiene un modo de trabajo muy consistente y c´omodo pero las pruebas realizadas con esta librer´ıa han sido un orden o dos de magnitud m´as lentos que con OpenCV y no en las tareas m´as costosas computacionalmente. Se ha contemplado tambi´en el uso de Cython, que es una variante de Python compatible pero que a˜nadiendo algunas directivas, sobre todo de tipado, permite compilarlo a c´odigo C haci´endolo comparable en eficiencia a este. No se ha profundizado en esta aproximaci´on puesto que no es un proyecto maduro y se ha preferido evitar tener problemas por este motivo. 4.16.4. Java Java dispone de un wrapper para OpenCV a trav´es de JNI (Java Native Interface) lo que en principio lo har´ıa bastante eficiente. Se ha contemplado el uso de este lenguaje porque es el que se emplea de manera mayoritaria en los desarrollos del grupo y es interesante desde el punto de vista de la integraci´on en otras aplicaciones. 4.16.5. Matlab/Octave Ambos son lenguajes de muy alto nivel que permiten un protototipado r´apido de los algoritmos. As´ı mismo son lenguajes orientados al tratamiento matem´atico, lo que nos es de bastante utilidad. Adem´as en teor´ıa tambi´en deber´ıan ser capaces de acceder a OpenCV mediante un wrapper. Tambi´en hay que tener en cuenta que la sintaxis de ambos es muy similar y salvo algunos detalles se puede decir que son compatibles. La ventaja principal de Octave es que su licencia es GPL y su mayor inconveniente encontrado es que tiene algunos errores en su implementaci´on que hacen desperdiciar tiempo en solucionarlos. Por ejemplo, no cargar una imagen por un fallo en el algoritmo que compone las rutas (se suministr´o a los desarrolladores de Octave la versi´on corregida). Por contra, el entorno y el flujo de trabajo en Matlab est´a mucho m´as pulido, es m´as c´omodo de trabajar y no presenta errores de este tipo. 4.16.6. OpenCV OpenCV (Open Source Computer Vision) es una librer´ıa cuyo prop´osito es proporcionar funciones para la visi´on por computador en tiempo real, bajo licencia BSD. De esta librer´ıa lo que realmente nos interesa es usarla para delegar en ella el manejo en memoria de las im´agenes, puesto que lo hace de manera muy eficiente y ofrece interfaces para C++, C, C#, Ruby, Python y Java manteniendo la eficiencia y permitiendo el manejo de las mismas como memoria din´amica. 53 Tambi´en vamos a hacer uso de algunas de sus funciones como por ejemplo las transformaciones de las im´agenes mediante la matriz de homograf´ıa ya que es una funcionalidad muy probada y no produce artefactos en la imagen. 54 4.17. Acr´onimos y abreviaturas FFT Fast Fourier Transform GMAPS Google MAPS GPL General Public License GPU Graphics processing unit HOG Histograms of Oriented Gradients IDE Infraestructuras de Datos Espaciales JNI Java Native Interface LMA Levenberg-Marquardt Algorithm OpenCV Open source Computer Vision PCA Principal Component Analysis PIL Python Image Library PNOA Plan Nacional de Ortofotograf´ıa A´erea RANSAC RANdom SAmple Consensus SSD Sum of Squared Differences SIFT Scale-invariant feature transform SURF Speeded Up Robust Feature SVT Scalable Vocabulary Tree 55