Métodos variacionales para la estimación del flujo óptico y mapas de disparidad

Departamento de Inform´atica y Sistemas M´etodos Variacionales para la Estimaci´on del Flujo ´ Optico y Mapas de Disparidad Optical Flow and Disparity Maps Estimation using Variational Methods Tesis Doctoral Agust´ın Javier Salgado de la Nuez Las Palmas de Gran Canaria Enero 2010 A mis padres Agradecimientos El desarrollo de una tesis doctoral supone un gran esfuerzo en el que han intervenido un gran n´umero de personas. Aprovecho este momento para agradecer la ayuda y el apoyo recibido durante este periodo. En primer lugar, quisiera agradecer a mis padres el esfuerzo econ´omico realizado durante estos a˜nos para que pudiera estudiar una carrera universitaria y posteriormente continuar con los estudios de doctorado. Hacer una menci´on especial a Luis Alvarez por ofrecerme la posibilidad de realizar los estudios de doctorado en el seno del grupo AMI y continuar la relaci´on iniciada en el proyecto fin de carrera. Adem´as agraceder al tutor de mi tesis, Javier S´anchez su tiempo, apoyo y dedicaci´on durante el desarrollo de la tesis. Durante todos estos a˜nos he tenido la oportunidad de compartir momentos de trabajo y distracci´on con compa˜neros en el laboratorio de investigaci´on AMI de la Universidad de Las Palmas de Gran Canaria. Quisiera recordarlos: Karl, Carlos, David, Jes´us, Pedro, Antonio, Miguel y Laura. Agraceder a los miembros del grupo AMI la ayuda y comentarios recibidos: Agust´ın Trujillo, Carmelo Cuenca, Julio Esclar´ın, Luis Mazorra y Miguel Alem´an. Quisiera agraceder a Joachim Weickert por permitirme realizar una estancia de investigaci´on en el seno de su grupo de investigaci´on. Durante este per´ıodo tuve la posibilidad de asistir a sus clases, ampliar mis conocimientos sobre m´etodos variacionales para la estimaci´on del flujo ´optico e intercambiar ideas con ´el durante much´ısimas horas. Quisiera recordar a Andres Bruhn, Luis Pizarro e Irena Galic por su colaboraci´on, consejos y las ayudas recibidas durante dicha estancia. Por ´ultimo, quisiera agradecer a las instituciones que han financiado parte de los trabajos realizados en el contexto de esta tesis. En primer lugar, al Departamento de Inform´atica y Sistemas y al servicio de investigaci´on de la Universidad de Las Palmas v vi de Gran Canaria por poner los medios suficientes para llevar a cabo este trabajo. A otras instituciones como el DAAD (Deutscher Akademischer Austausch Dienst) que me concedieron una beca para el desarrollo de un proyecto de investigaci´on en el grupo del Profesor Joachim Weickert. A la fundaci´on de la Universidad de Las Palmas de Gran Canaria que gracias a la beca Innova, financiada por Unelco, me di´o la posibilidad de adquirir material utilizado en el desarrollo de esta tesis. Algunos trabajos de esta tesis forman parte de proyectos financiados por Consejer´ıa de Educaci´on Cultura y Deportes del gobierno de Canarias (PI2002/193), Ministerio de Ciencia y Tecnolog´ıa y FEDER (TIC2003-08957). Agredecer tambi´en a la empresa MediaPro la cesi´on y el permiso para la utilizaci´on de una secuencia de f´utbol utilizada en un m´etodo de estimaci´on del mapa de disparidad. ´ Indice general Introducci´on 1 Abstract 11 1. Estado del Arte 15 1.1. Flujo ´ Optico ................................... 15 1.1.1. Clasificaci´on de los M´etodos . . . . . . . . . . . . . . . . . . . . . . . 17 1.1.2. M´etodos Variacionales . . . . . . . . . . . . . . . . . . . . . . . . . . 20 1.2. Visi´on Estereosc´opica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 1.2.1. Geometr´ıa Epipolar . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 1.2.2. Clasificaci´on de los M´etodos . . . . . . . . . . . . . . . . . . . . . . . 28 1.3. Caracter´ısticas de las Im´agenes . . . . . . . . . . . . . . . . . . . . . . . . . 31 1.3.1. Secuencias de Im´agenes . . . . . . . . . . . . . . . . . . . . . . . . . 31 1.3.2. ParesEst´ereo............................... 35 1.4. SecuenciasdePrueba............................... 35 1.4.1. Flujo ´ Optico ............................... 35 1.4.2. Visi´on Estereosc´opica . . . . . . . . . . . . . . . . . . . . . . . . . . 42 1.5. MedidasdeError................................. 50 2. Estimaci´on del Flujo ´ Optico en Secuencias de Im´agenes 53 2.1. Introducci´on.................................... 53 2.1.1. Contribuciones de este Cap´ıtulo . . . . . . . . . . . . . . . . . . . . . 54 2.2. Generalizaci´on de los Modelos de Energ´ıa . . . . . . . . . . . . . . . . . . . 56 2.2.1. Modelos de Energ´ıa Continuos . . . . . . . . . . . . . . . . . . . . . 56 2.2.2. Modelos de Energ´ıa No Continuos Secuencial . . . . . . . . . . . . . 58 2.2.3. Modelos de Energ´ıa No Continuos Aleatorio . . . . . . . . . . . . . . 59 2.3. M´etodo Variacional Multicanal . . . . . . . . . . . . . . . . . . . . . . . . . 61 vii viii ´ INDICE GENERAL 2.3.1. M´etodo Variacional Monocanal . . . . . . . . . . . . . . . . . . . . . 63 2.3.2. Modelo de Energ´ıa . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 2.3.3. Minimizaci´on de la Energ´ıa . . . . . . . . . . . . . . . . . . . . . . . 65 2.3.4. Esquema Num´erico . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 2.3.5. Resultados Experimentales . . . . . . . . . . . . . . . . . . . . . . . 67 2.4. M´etodo Variacional con Regularizaci´on Temporal no Continua . . . . . . . 76 2.4.1. Modelo de Energ´ıa . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 2.4.2. Minimizaci´on de la Energ´ıa . . . . . . . . . . . . . . . . . . . . . . . 80 2.4.3. Esquema Num´erico . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 2.4.4. C´alculo del Flujo Inverso, h∗...................... 82 2.4.5. Resultados Experimentales . . . . . . . . . . . . . . . . . . . . . . . 84 2.5. M´etodos Variacionales basados en el An´alisis Espectral . . . . . . . . . . . . 94 2.5.1. Tensor de Movimiento . . . . . . . . . . . . . . . . . . . . . . . . . . 95 2.5.2. Modelo de Energ´ıa . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96 2.5.3. Minimizaci´on de la Energ´ıa . . . . . . . . . . . . . . . . . . . . . . . 102 2.5.4. Esquema Num´erico . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 2.5.5. Resultados Experimentales . . . . . . . . . . . . . . . . . . . . . . . 111 2.6. Conclusiones ...................................120 2.6.1. M´etodo Variacional Multicanal . . . . . . . . . . . . . . . . . . . . . 120 2.6.2. M´etodo Variacional con Regularizaci´on Temporal no Continua . . . 120 2.6.3. M´etodo Variacional basado en el An´alisis Espectral . . . . . . . . . . 121 3. Estimaci´on del Mapa de Disparidad en Pares Est´ereo 123 3.1. Introducci´on....................................123 3.1.1. Contribuciones de este Cap´ıtulo . . . . . . . . . . . . . . . . . . . . . 123 3.2. Estimaci´on de la Disparidad usando una Secuencia Estereosc´opica . . . . . 125 3.2.1. Notaci´on del Flujo ´ Optico y Mapas de Disparidad . . . . . . . . . . 125 3.2.2. Modelo de Energ´ıa . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 3.2.3. Minimizaci´on de la Energ´ıa . . . . . . . . . . . . . . . . . . . . . . . 129 3.2.4. Resultados Experimentales . . . . . . . . . . . . . . . . . . . . . . . 131 3.3. Mapa de Disparidad combinando M´etodos de Graph–cuts y Variacional . . 143 3.3.1. M´etodo de Correlaci´on . . . . . . . . . . . . . . . . . . . . . . . . . . 144 3.3.2. M´etodo de Graph-cuts . . . . . . . . . . . . . . . . . . . . . . . . . . 146 3.3.3. M´etodo Variacional . . . . . . . . . . . . . . . . . . . . . . . . . . . 148 ´ INDICE GENERAL ix 3.3.4. Combinaci´on del M´etodo de Graph-cuts y Variacional . . . . . . . . 149 3.3.5. Resultados Experimentales . . . . . . . . . . . . . . . . . . . . . . . 150 3.4. Conclusiones ...................................168 3.4.1. Estimaci´on de la Disparidad a partir de una Secuencia Estereosc´opica168 3.4.2. Estimaci´on del Mapa de Disparidad mediante la Combinaci´on de un M´etodo de Graph–cuts y uno Variacional . . . . . . . . . . . . . . . 168 4. Conclusiones 171 4.1. Resumen......................................171 4.2. Trabajofuturo ..................................175 Conclusions 177 Notaci´on 179 Listado de Figuras 183 Listado de Tablas 191 Bibliograf´ıa 195 6´ INDICE GENERAL (2000) [Alvarez00] y Alvarez/Deriche/S´anchez/Weickert (2002) [Alvarez02b]. Ambos m´etodos se basan en una t´ecnica de minimizaci´on de energ´ıa que ofrecen resultados precisos y densos. Se comenzar´a repasando algunos conceptos sobre la geometr´ıa epipolar para ir profundizando en ellos de forma que se vayan fusionando las ideas del m´etodo del flujo ´optico y del est´ereo. La notaci´on definida en esta introducci´on se utilizar´a en la presentaci´on del modelo de energ´ıa. Se continua con la minimizaci´on y el esquema num´erico incluyendo la descripci´on del enfoque multipiramidal que permite detectar los desplazamientos largos. El esquema num´erico utiliza la t´ecnica del descenso del gradiente. Finalizaremos presentando los resultados experimentales obtenidos utilizando secuencias sint´eticas y reales. En el segundo trabajo se propone un m´etodo para la estimaci´on del mapa de disparidad que combina otros dos: uno variacional descrito en Alvarez/Deriche/S´anchez/Weickert (2002) [Alvarez02b] y uno de graph–cuts, [Boykov04]. En primer lugar, haremos menci´on al proceso de rectificaci´on. Cuando la disposici´on de las c´amaras es frontoparalela la formulaci´on de los modelos de energ´ıa se simplifica ya que el desplazamiento de los p´ıxeles queda expresado como un escalar (en una ´unica direcci´on). La convergencia del m´etodo variacional se acelera si se dispone de una buena aproximaci´on inicial. En este trabajo se describen dos t´ecnicas utilizadas para el c´alculo de la aproximaci´on inicial: la correlaci´on ygraph– cuts. Seguidamente se presenta las caracter´ısticas del m´etodo variacional y c´omo se ha combinado con el de graph–cuts para calcular el mapa de disparidad. Se introduce el m´etodo variacional en un enfoque multipiramidal para la detecci´on de largos desplazamientos y evitar su convergencia m´ınimos locales irrelevantes. Por ´ultimo, se muestran los resultados experimentales obtenidos estableciendo una comparaci´on entre el m´etodo variacional con las dos aproximaciones iniciales. Se quiere poner de manifiesto la mejora que supone la inclusi´on del m´etodo de graph–cuts. En una segunda bater´ıa de pruebas que quiere poner a prueba la estabilidad de las distintas t´ecnicas ante la presencia de ruido en los pares est´ereo. Al finalizar este cap´ıtulo se presentan las conclusiones valorando los hechos m´as significativos de los trabajos presentados en este cap´ıtulo. Cap´ıtulo 4, Conclusiones: En la secci´on 4.1 de este cap´ıtulo se exponen las conclusiones finales de la tesis comentando los problemas encontrados durante su desarrollo y destacando las aportaciones realizadas a la literatura en todos los trabajos recogidos en este documento. En la secci´on 4.2 se enumera algunas tareas que por su complejidad y falta de tiempo se ha dejado pendiente para el futuro. Anexo 4.2, Notaci´on: En este anexo se unifica toda la notaci´on utilizada en la tesis. Se trata de una gu´ıa r´apida en la que el lector puede familiarizarse con la nomenclatura empleada en el documento. ´ INDICE GENERAL 7 Las principales aportaciones En este documento se presentan una serie de trabajos que representan las aportaciones novedosas a la literatura en el campo de la estimaci´on del flujo ´optico y est´ereo realizadas en esta tesis. A continuaci´on, se describen brevemente estas contribuciones. M´etodo Variacional Espacial Multicanal: Extensi´on del Modelo de Alvarez/Weickert/S´anchez (2000) [Alvarez00]: Bas´andonos en el funcional de energ´ıa propuesto por [Alvarez00] se le ha aplicado una serie de modificaciones para la inclusi´on de la informaci´on multicanal. Los cambios introducidos en el modelo de energ´ıa son los siguientes: (1) la definici´on de un t´ermino de ligadura como una suma ponderada de la suposici´on lambertiana de cada uno de los canales y (2) la modificaci´on del operador de Nagel-Enkelmann para que su gradiente aglutine informaci´on de varias im´agenes. Se han definido diversas estrategias para la construcci´on de dicho gradiente. Este m´etodo se ha utilizado para la detecci´on de masas nubes en secuencias sat´elites multiespectrales. M´etodo Variacional que incorpora un T´ermino de Regularizaci´on exclusivamente Temporal no Continuo: La principal contribuci´on de este trabajo ha sido el desarrollo de un m´etodo variacional espaciotemporal que incorpora un t´ermino de regularizaci´on exclusivamente temporal. El desacoplamiento del t´ermino de suavizado en dos, uno espacial y otro temporal, pretende evitar algunas de las limitaciones presentes en trabajos con una regularizaci´on temporal continua basada en derivadas parciales temporales, como en [Papenberg06]. Nuestro m´etodo realiza convenientemente la regularizaci´on temporal teniendo en cuenta la presencia de desplazamientos largos sin afectar a la regularizaci´on espacial. M´etodo Variacional basado en el An´alisis Espectral del Tensor de Movimiento: La principal aportaci´on de este trabajo ha sido el desarrollo de un m´etodo variacional espaciotemporal que intercambia informaci´on complementaria entre el t´ermino de ligadura y el de suavizado. El proceso de difusi´on s´olo se aplica en las direcciones ortogonales al flujo dominante. La combinaci´on de varias t´ecnicas existentes posibilita la creaci´on de un framework capaz de definir cualquier modelo de una forma compacta, sencilla, f´acilmente extensible y adaptable. El tensor de movimiento permite representar cualquier invarianza en una notaci´on compacta. La descomposici´on espectral de este tensor asegura la f´acil identificaci´on de las direcciones dominantes del flujo al mismo tiempo que define una notaci´on independiente de las invarianzas definidas en el modelo de energ´ıa. La inclusi´on de t´erminos no-cuadr´aticos no s´olo aporta robustez a las estimaciones frente al ruido sino que la combinaci´on de las funciones de robustificaci´on junto a la descomposici´on espectral del tensor de movimiento permite la construcci´on de 8´ INDICE GENERAL cuatro prototipos capaces de mejorar las estimaciones de m´etodos de caracter´ısticas similares. M´etodo para la estimaci´on del mapa de disparidad utilizando una secuencia de pares est´ereo: La principal aportaci´on de este trabajo es la creaci´on de un m´etodo que combina la informaci´on del flujo est´ereo (Alvarez/Deriche/S´anchez/Weickert (2002), [Alvarez02b]) con la del flujo ´optico ( Alvarez/Weickert/S´anchez (2000), [Alvarez00]) aumentando la robustez de las estimaciones gracias a la fusi´on, dentro del mismo modelo de energ´ıa, de las ideas plasmadas en dos m´etodos de gran precisi´on. La estimaci´on del mapa de disparidad se apoya en la geometr´ıa epipolar para poder establecer las correspondencias de cada par est´ereo. Esas mismas correspondencias se pueden obtener mediante el c´alculo del flujo ´optico a lo largo de la secuencia de ambas c´amaras. Por lo tanto, en una secuencia de pares est´ereo disponemos de hasta cuatro vistas de un punto 3D en dos instantes de tiempo distintos. Este m´etodo aprovecha la informaci´on suministrada por el flujo est´ereo y ´optico, estimados de forma conjunta, para as´ı poder lograr un mapa de disparidad m´as preciso incluso en situaciones donde existen oclusiones. Para asegurar la coherencia del resultado final se ha incluido en el modelo de energ´ıa una restricci´on que fuerza la congruencia de ambos flujos. Combinaci´on de un m´etodo variacional espacial y uno de graph–cuts para la estimaci´on del mapa de disparidad: La segunda aportaci´on en el campo de la estimaci´on del mapa de disparidad consiste en la combinaci´on de un m´etodo variacional (Alvarez/Deriche/S´anchez/Weickert (2002), [Alvarez02b]) y una t´ecnica de graph–cuts ([Kolmogorov01, Boykov04]). El m´etodo de graph–cuts se ha utilizado para obtener una aproximaci´on inicial del mapa de disparidad que el m´etodo variacional se encargar´ıa de refinar. Normalmente, la t´ecnica m´as empleada para estimar la inicializaci´on es una basada en la correlaci´on a ventanas. Dado los buenos resultados que ofrece algunos m´etodos de graph–cuts en el campo de la estimaci´on del mapa de disparidad parece una buena alternativa a los tradicionales m´etodos basados en correlaci´on a ventanas. ´ INDICE GENERAL 9 Publicaciones realizadas En esta secci´on se har´a una breve descripci´on de las publicaciones realizadas en el contexto de esta tesis. Multi-Channel Satellite Image Analysis Using a Variational Approach Alvarez/Casta˜no/Garc´ıa/Krissian/Mazorra/Salgado/S´anchez (2008) [Alvarez08]: En este trabajo se presenta un m´etodo variacional multicanal para hacer frente a algunos de los problemas m´as habituales en el an´alisis de im´agenes por sat´elite, como es la estimaci´on del movimiento de las estructuras nubosas presentes en la atm´osfera. El modelo de energ´ıa propuesto se basa en el trabajo de Alvarez/Weickert/S´anchez (2000) [Alvarez00] y combina la informaci´on de varios canales del sat´elite. Las ventajas y mejoras en las estimaciones que ofrece este nuevo m´etodo se ponen de manifiesto en los experimentos realizadas con dos secuencias de sat´elite. Optical Flow Estimation with Large Displacements: A Temporal Regularizer Salgado/S´anchez (2006, 2007) [Salgado06a], [Salgado06b], [Salgado07b]: En este trabajo se presenta un modelo variacional para la estimaci´on del flujo ´optico en una secuencia de im´agenes. Bas´andonos en el modelo espacial propuesto por Alvarez/Weickert/S´anchez (2000) [Alvarez00] se ha a˜nadido un nuevo regularizador temporal que desacopla la informaci´on espaciotemporal al mismo tiempo que es capaz de detectar los desplazamientos grandes. La disociaci´on entre la regularizaci´on espacial y temporal es necesaria con el objetivo de evitar o resolver algunas incongruencias detectadas en otros trabajos. 3D Geometry Reconstruction from a Stereoscopic Video Sequence Salgado/S´anchez (2005) [Salgado05a], [Salgado05d]: En este trabajo se propone un m´etodo que estima la geometr´ıa 3D de una escena a partir de una secuencia de v´ıdeo tomada desde un par de c´amaras est´ereo. Las c´amaras est´an r´ıgidamente situadas en una posici´on fija, en posici´on frontoparalela, y hay una serie de objetos movi´endose por la escena. El m´etodo propuesto calcula el desplazamiento de los objetos y la estructura 3D de la escena a trav´es de la estimaci´on del flujo ´optico de la secuencia de cada c´amara y de la estimaci´on del mapa de disparidad de cada par est´ereo de dicha secuencia. Para relacionar esta informaci´on se ha establecido una restricci´on temporal que relaciona e impone cierta coherencia entre el flujo ´optico y est´ereo. Esta restricci´on se justifica haciendo uso de la formulaci´on matem´atica com´un que existe entre estos problemas. Combining two Methods to Accurately Estimate Dense Disparity Maps Salgado/S´anchez (2005, 2007) [Salgado05c], [Salgado07a], [Salgado05b]: En este trabajo se combinan dos m´etodos con el objeto de mejorar las estimaciones de la geometr´ıa 3D de una escena. Para ello, se apoya en par de im´agenes estereosc´opicas donde la disposici´on de las c´amaras puede ser arbitraria. Para simplificar la complejidad de los c´alculos se recurre a un proceso de rectificaci´on. La estimaci´on 10 ´ INDICE GENERAL del mapa de disparidad se realiza empleando el m´etodo variacional propuesto en Alvarez/Deriche/S´anchez/Weickert (2002), [Alvarez02b]. A partir de un modelo de energ´ıa se obtiene un PDE cuya soluci´on determina el desplazamiento registrado en cada par est´ereo. Uno de los problemas de este tipo de t´ecnicas es que depende en gran medida de la primera aproximaci´on, es decir, de su inicializaci´on. En este trabajo se compara la influencia sobre el resultado final de dos t´ecnicas que se utilizan para calcular esta aproximaci´on inicial. Por un lado, tenemos una t´ecnica basada en la correlaci´on y, por otro lado, una de graph–cuts. Abstract The overall objective of this thesis is to contribute originally to the development of the a set of variational methods for estimating optical flow and disparity maps from image sequences. Depending on the number of cameras in a scene, we can have two different settings: (i) a sequence of images taken by one camera at different time step and (ii) a sequence of images taken from several cameras at the same time. These two cases represent two key problems in computer vision and have been investigated for decades: (i) the optical flow and (ii) the disparity map estimation. The optic flow is the apparent motion of pixels between images. Given two images, the aim is the retrieval of the pixel displacements from one image to the other. This problem is the base for a broad number of applications such as 3D reconstruction, driver assistance systems, video compression and surveillance systems. There are still open research issues regarding this topic and new important contributions often appear. In order to solve the optic flow problem many techniques have been proposed. Among them, the variational methods have demonstrated to be one of the best techniques to obtain accurate solutions. In the stereo problem, we have two views of the scene taken at the same time. The disparity map is the computation of the pixels displacement from one view to the other. The stereo problem has a very useful tool such as the epipolar geometry. The epipolar geometry allows us to limit the search area of the points in correspondence along a line. The first part of this thesis is focused on the variational optic flow methods using image sequences. In this topic we have developed: 1. a mathematical model, based on a variational approach, to estimate the cloud structure motion by combining information from various satellite channels. We include information of all the channels in a single variational motion estimation model. The initial optical flow technique [Alvarez00] is the base of our multichannel sequential motion tracking algorithm. We extend this variational optical flow method to deal with multichannel sequential data; 2. a model for computing the optical flow in a sequence of images with a spatial– temporal regularizer explicitly designed for large displacements. We study the 11 12 ´ INDICE GENERAL introduction of a temporal regularizer that expands the information beyond two consecutive frames. We propose to decouple the spatial and temporal regularising terms to avoid an incongruous formulation between the data and smoothness term. Our model is based on an energy functional that yields a partial differential equation (PDE). This PDE is embedded into a multi-pyramidal strategy to recover large displacements; 3. an energy model for the optic flow estimation using a sequence of images. It takes advantage of the recently introduced motion tensor within an energy functional. Through a principal component analysis of the motion tensor we construct novel energy functionals. We define a new framework that nicely combines two complementary tensors. This provides a direct mechanism to switch between smoothing the solution along the prominent directions and attracting objects with similar intensity values. Each one is derived from the eigenvalues and eigenvectors of the motion tensor. The main contribution of this work is the inclusion of a modified tensor to steer the diffusion process and the motion tensor decomposition. We use non-quadratic functionals improve the method’s robustness with respect to outliers. This energy model can handle large displacements through the use of warping techniques. The second part of this thesis is dedicated to the variational stereo flow methods using pairs of images. The main contributions in this topic are: 1. a novel method for the reconstruction of the 3D geometry of a scene from a stereoscopic video sequence. There are two video–cameras pointing to the same scene and recording frames at the same time. For every stereoscopic pair of images in the sequence we may compute a disparity map independently from the other pairs, to obtain a set of independent disparity maps. The problem is that, in general, the continuity of the solution is not preserved and it is very sensitive to the presence of noise. If we want to overcome this problem, we have to relate the estimation of disparity maps through the sequence. One way to do it is to compute the displacement of objects on both video–cameras and use this information to constraint the computation of the disparity maps in time. This work is a continuation of previous works on optical flow [Alvarez00] and disparity map estimation [Alvarez02b]. These two methods were also based on energy minimization techniques and proved to be reliable and accurate; 2. a method that combines two techiques for computing disparity maps using a pair of images. The first one is based on graph–cut energy minimization [Kolmogorov01, ´ INDICE GENERAL 13 Boykov04]. This method has demonstrated good results in integer precision which is enough for some applications. If we need a better accuracy, then it is necessary to use a different technique. In this case, we propose to use the variational method described in [Alvarez02b] as a suitable complement. Both methods are based on energy minimization approach. One of the problems of these methods is that they need a good initial approximation in order to obtain an accurate solution. We normally use a cross-correlation technique to compute this initialization. As a conclusion, we have presented in this thesis different approaches to deal with the optical and stereo flow problems. As a result of our research, we have developed new methods that give solutions for different situations where displacements may be large and sub-pixel accuracy is needed. Cap´ıtulo 1 Estado del Arte En los ´ultimos treinta a˜nos se ha desarrollado una intensa actividad investigadora en el campo de la visi´on por computador. En sus or´ıgenes la visi´on por computador estaba estrechamente relacionada con la rob´otica pero gracias al acercamiento de los ordenadores al gran p´ublico y el incremento de la capacidad de c´omputo de las m´aquinas surgieron nuevos problemas a los que la visi´on por computador pod´ıa dar respuesta. Entre estos problemas destacamos la estimaci´on del flujo ´optico en secuencias de im´agenes y el c´alculo del mapa de disparidad en pares est´ereo. Durante todo este tiempo se han desarrollado m´etodos que ofrecen soluciones de gran precisi´on. Sin embargo, todav´ıa siguen existiendo algunas cuestiones sin resolver a la que intentaremos dar respuesta en esta tesis. En este cap´ıtulo se describe el estado del arte referente al problema de la estimaci´on del flujo ´optico y de la carta de disparidad. Se hace un recorrido por los m´etodos m´as importantes de la literatura comentando algunas de las t´ecnicas m´as exitosas utilizadas en los ´ultimos a˜nos y que han supuesto una fuente de inspiraci´on para los m´etodos desarrollados con posterioridad. Por ´ultimo, se describe las secuencias de im´agenes utilizadas en los distintos experimentos y las m´etricas de error en las que nos apoyamos para evaluar cuantitativamente las soluciones de los m´etodos. 1.1. Flujo ´ Optico El c´alculo del flujo ´optico consiste en la estimaci´on del movimiento aparente de los objetos en una secuencia de im´agenes. Dado un conjunto de im´agenes, el objetivo es calcular el desplazamiento de los p´ıxeles entre las distintas im´agenes. Disponemos de una c´amara (o videoc´amara) que capta im´agenes de una escena. En la escena podemos encontrar una serie de objetos est´aticos o din´amicos que, por lo general, se ver´an influenciados por condiciones variables del entorno, tales como fuentes de iluminaci´on, sombras, reflejos y otros efectos luminosos, as´ı como por otras dificultades asociadas a la aparici´on y desaparici´on de objetos en la escena o la oclusi´on de unos objetos con otros. El problema del flujo ´optico se ha convertido en uno de los m´as importantes a resolver en el campo de la visi´on por ordenador. Su importancia radica fundamentalmente en el 15 22 CAP´ ITULO 1. ESTADO DEL ARTE funciones de robustificaci´on se han convertido en una de las herramientas m´as ´utiles para para minimizar los problemas que genera el ruido ya que los outliers son penalizados en menor grado que una funci´on cuadr´atica. El uso de estad´ısticos robustos en los m´etodos de estimaci´on del flujo ´optico los introdujo [Black91, Black96b]. Bas´andose en la investigaci´on realizada por [Hampel86, Huber81] acerca de los estad´ısticos robustos, Black et al. propuso el uso de M-estimators como funciones de robustificaci´on. En este sentido, cabe destacar el trabajo realizado por [M´emin98a] en el cual el problema de optimizaci´on no-convexa se resolvi´o haciendo un ajuste por m´ınimos cuadrados ponderados e iterativo. Existen otras contribuciones como [Black92, Black96a] en la que tambi´en se proponen otras funciones para atenuar el efecto de los outliers. En el trabajo [Haussecker01] se propone un m´etodo para tratar las variaciones de intensidad de las im´agenes. En [Wells96, Viola97, Hermosillo02] se proponen m´etodos m´as robustos basados en la estimaci´on de correspondencias multimodal entre im´agenes para tratar con secuencias donde existe una transformaci´on compleja en las intensidades de los p´ıxeles. En [Bruhn05c, Papenberg06] se ha demostrado que el uso de funciones de robustificaci´on reduce significativamente el error mejorando la precisi´on de las estimaciones. Muchos de los trabajos que hemos comentado utilizan ´unicamente informaci´on espacial para calcular el flujo ´optico. En una secuencia de im´agenes el movimiento de los objetos se propaga m´as all´a de dos frames. Por ello, la informaci´on temporal procedente de los frames anteriores o posteriores nos permite la correcci´on de las estimaciones y mejorar sustancialmente la precisi´on de las mismas. Uno de los primeros trabajos que incluyeron informaci´on espaciotemporal fue el de Nagel [Nagel90], en donde se propone una extensi´on temporal de su operador, [Nagel86]. Black y Anandan [Black91] propusieron un m´etodo, que asume la existencia de una aceleraci´on en el tiempo, en el que el movimiento se calcula incrementalmente. Para hacer este m´etodo m´as robusto frente al ruido se calcula una media de la aceleraci´on en la componente temporal. Pese a que hace m´as de quince a˜nos se propusieron estos m´etodos no ha sido hasta hace unos pocos a˜nos cuando la informaci´on temporal no ha sido ampliamente incluida en los m´etodos. Algunas contribuciones interesantes son [Elad98, Weickert01b, Farneb¨ack01]. Weickert y Schn¨orr [Weickert01b] propusieron un modelo continuo con t´ermino de suavizado espaciotemporal convexo nolineal. En este trabajo tanto las derivadas espaciales como temporales se formularon de una forma homog´enea. En Brox et al. [Weickert04], la funci´on de robustificaci´on conocida como Total Variation se aplica tanto en el t´ermino de ligadura como en el de suavizado. Este trabajo permite detectar los largos desplazamientos e incluye en el t´ermino de ligadura dos invarianzas: la ya tradicional suposici´on lambertiana y el gradiente constante. En el t´ermino de suavizado se asume que el flujo es suave tanto en la direcci´on espacial como temporal. Recientemente, Papenberg et al. [Papenberg06] ha obtenido unos resultados excelentes gracias a la combinaci´on de m´ultiples elementos, como son el uso de funciones de robustificaci´on, la inclusi´on de varias invarianzas dentro del t´ermino de ligadura y la informaci´on espaciotemporal. Los m´etodos m´as sofisticados hasta ahora combinan las t´ecnicas variacionales con los level sets. Amiaz et al. [Amiaz06, Amiaz07] propone un m´etodo que incluye la t´ecnica de segmentaci´on en el modelo de energ´ıa variacional. Brox et al. [Brox06] han propuesto un trabajo muy similar que combina el m´etodo de Papenberg et al. [Papenberg06] e incluyen informaci´on de level sets. Uno de los ´ultimos trabajos m´as relevantes en el campo de la 1.1. FLUJO ´ OPTICO 23 estimaci´on del flujo ´optico es el de [Zimmer09]. En ´el se propone un modelo de energ´ıa que introduce el concepto de complementariedad entre el t´ermino de ligadura y suavizado al mismo tiempo que utiliza alguna de las t´ecnicas m´as innovadoras de la literatura. Otras de las estrategias que se est´an investigando es la inclusi´on de t´ecnicas para la detecci´on de oclusiones. Los m´etodos actuales han llegado a un alto nivel de sofisticaci´on y complejidad. Cada nuevo elemento que se incluye en el modelo de energ´ıa requiere un nuevo par´ametro. Esto conlleva que la optimizaci´on de los par´ametros sea una tarea que consume mucho tiempo. En los ´ultimos a˜nos han surgido algunos m´etodos conocidos como learning optic flow [Roth07, Sun08, Li08] cuya caracter´ıstica principal es la utilizaci´on de un ´unico vector de par´ametros para cualquier secuencia. Para obtener este ´unico vector de par´ametros el m´etodo se entrena con una base de datos de secuencias. En cierto modo, ese vector de par´ametros aglutina la mejor optimizaci´on para cualquier tipo de movimiento recogido en dicha base de datos. Estos par´ametros no son los ´optimos para una determinada secuencia pero de forma gen´erica son los mejores para cualquier secuencia. A la hora de formular el modelo de energ´ıa se asume que el flujo ´optico es asim´etrico, es decir, que el desplazamiento se produce en una sola direcci´on, de una imagen a otra. En la literatura se han realizado una serie de trabajos [Cachier00, Christensen01], Alvarez/Deriche/S´anchez/Weickert (2002, 2007) [Alvarez02a, Alvarez07a], Alvarez/Casta˜no/Garc´ıa/Krissian/Mazorra/Salgado/S´anchez (2007) [Alvarez07b] en el que se describen modelos sim´etricos del flujo ´optico donde la soluci´on se obtiene como una combinaci´on del flujo en ambas direcciones. Como se ha visto en secciones anteriores la ecuaci´on de restricci´on del flujo (OFC) se basa en la idea de la existencia de las derivadas. Esta suposici´on no es v´alida cuando los desplazamientos son largos ya que la condici´on de derivabilidad no se cumple. Para poder abordar este problema y utilizar las t´ecnicas existentes para desplazamientos cortos, en la literatura han surgido lo que se conoce como estrategias multiescala. La idea consiste en crear distintas escalas (versiones) de la secuencia de im´agenes de entrada, de forma que los desplazamientos en cada una de ellas se vaya reduciendo llegando a un punto en el que la OFC es computable. Este desplazamiento corto en la escala inferior representa el m´aximo desplazamiento en la secuencia inicial. El flujo ´optico calculado en la escala inferior ser´a utilizado como inicializaci´on en la escala superior. Dado que en la escala inferior no existe ninguna estimaci´on previa se aplica el m´etodo variacional directamente. Existen varias estrategias para crear cada una de las escalas, entre las que destacamos el enfoque piramidal y gaussiano. La estrategia multipiramidal consiste en crear por cada imagen una familia de subim´agenes, aplic´andole un factor de escalado en cada nivel. Este enfoque ha sido utilizado en gran n´umero de m´etodos como en [Anandan89, Battiti91, Luettgen94, Bornemann96], [Enkelmann88, M´emin02]. Esta estrategia ofrece dos importantes ventajas. Por un lado, este tipo de m´etodos mejora el rendimiento ya que el flujo ´optico se calcular´a m´as r´apidamente en las escalas inferiores y con la inicializaci´on en las superiores la convergencia del m´etodo ser´a m´as r´apida, [Bruhn05b]. Por otro lado, en los funcionales de energ´ıa noconvexo la convergencia al m´ınimo global no est´a asegurada. En las escalas inferiores del esquema multipiramidal la posibilidad de alcanzar m´ınimos locales indeseados desaparece, creando buenas inicializaciones que disminuye el riesgo de alcanzar un m´ınimo 24 CAP´ ITULO 1. ESTADO DEL ARTE local irrelevante en las escalas superiores y mejorando as´ı, las estimaciones obtenidas ([Black96b, M´emin98b, Bruhn05c, Papenberg06]). Cuando se aplica una gaussiana a una imagen se genera un efecto de suavizado difuminando los bordes. Si aplic´asemos sucesivas gaussianas aumentar´ıamos ese efecto sobre la imagen. En cierta forma, los desplazamientos largos se van reduciendo por el suavizado. En esta estrategia las escalas se crean a partir de sucesivas aplicaciones de gaussianas con desviaci´on σ; el tama˜no de las im´agenes no var´ıa. Un trabajo que utiliza esta estrategia es Alvarez/Weickert/S´anchez (2000) [Alvarez00]. 1.2. VISI ´ ON ESTEREOSC ´ OPICA 25 1.2. Visi´on Estereosc´opica En la visi´on estereosc´opica tenemos dos vistas de la misma escena en el mismo instante de tiempo. Cada punto 3D (M) se proyecta en ambas c´amaras, de forma que, existen dos proyecciones mym’ de cada punto 3D. Se define como disparidad al desplazamiento entre el punto 2D msituado en la imagen izquierda (Il) y su punto en correspondencia m’ en la imagen derecha (Ir). Por lo tanto, la visi´on estereosc´opica se reduce a un problema de correspondencias. En este sentido guarda una gran relaci´on con el c´alculo del flujo ´optico. Ambos problemas tienen muchas similitudes en cuanto a las caracter´ısticas y dificultades, de forma que muchas de las ideas que se aplican en un problema se adaptan al otro. El objetivo de la visi´on estereosc´opica es el de reconstruir la escena a partir de la proyecci´on de un par est´ereo. Para poder recuperar la informaci´on 3D es necesario: (1) calibrar el sistema de c´amaras, (2) estimar las correspondencias entre las im´agenes y, por ´ultimo, (3) calcular los puntos 3D a partir de las correspondencias. Salvo este ´ultimo paso los dos anteriores han sido objeto de amplio estudio y existen much´ısimos trabajos que proponen distintas t´ecnicas para su resoluci´on. 1. Calibraci´on de c´amaras: Para la estimaci´on de los mapas de disparidad es necesario que las c´amaras est´en calibradas. En la literatura se definen dos tipos de calibraciones: fuerte, que consiste en hallar las matrices de proyecci´on asociadas a cada c´amara y d´ebil, cuyo objetivo es la estimaci´on de la matriz fundamental, que relaciona la informaci´on entre las c´amaras mediante una transformaci´on lineal. La matriz fundamental ha sido utilizada en trabajos como [Faugeras93a, Faugeras01, Hartley03]. 2. Estimaci´on de las correspondencias: La b´usqueda de correspondencias, como se ha visto en el caso del c´alculo del flujo ´optico, se trata de un problema bastante complejo. Sin embargo, es posible simplificarlo bas´andonos en la geometr´ıa epipolar. Dado que la geometr´ıa epipolar ya ha sido descrita en varios trabajos ([Faugeras01, Hartley03]) en el apartado 1.2.1 s´olo se mencionar´a los conceptos b´asicos para que el lector no tenga que recurrir a documentaci´on adicional. 3. Reconstrucci´on de los puntos 3D: Dado que las im´agenes del par est´ereo se toman en el mismo instante de tiempo, el desplazamiento de los p´ıxeles nos da informaci´on acerca de la profundidad de los objetos respecto a la c´amara. P´ıxeles situados cerca de la c´amara tendr´an un desplazamiento mayor que los situados a mayor distancia. Los mapas de disparidad son unas im´agenes, normalmente expresadas en niveles de grises, que nos indican el desplazamiento de cada uno de los p´ıxeles de la imagen. Esta distancia se suele expresar como la norma del vector desplazamiento √u2+v2, donde uyvrepresentan el desplazamiento horizontal y vertical respectivamente. Si dispusi´esemos de un sistema de c´amaras calibradas fuertemente, ser´ıa posible reconstruir el punto 3D a partir de las de las correspondencias de los puntos 2D. Sin embargo, en la pr´actica observamos que existen muchos factores que dificultan la correcta reconstrucci´on de los puntos 3D. Estos factores se comentar´an en el apartado 1.3.2. 26 CAP´ ITULO 1. ESTADO DEL ARTE Figura 1.2: Geometr´ıa de una c´amara proyectiva. Ces el centro de la c´amara y pes el punto principal. 1.2.1. Geometr´ıa Epipolar El problema est´ereo relaciona la informaci´on captada por dos c´amaras en el mismo instante de tiempo. Esta tarea no es sencilla dado que se utiliza dispositivos independientes situados en posiciones distintas. Aunque dos c´amaras sean exactamente iguales, mismo modelo, marca y configuraci´on, las im´agenes que toman no son id´enticas ya que hay que tener en cuenta que se tratan de dispositivos f´ısicos y el desgaste de sus piezas no es exactamente el mismo. Cualquier c´amara fotogr´afica se puede representar mediante el modelo de pin-hole (figura 1.2) de forma que los puntos de la escena se proyectan a trav´es del foco de la c´amara (C) en el plano proyectivo o imagen, Π. La l´ınea desde el centro de la c´amara hasta el plano de la imagen recibe el nombre de eje principal orayo principal, y el punto donde el eje principal corta al plano de la imagen recibe el nombre de punto principal. El plano que contiene al centro ´optico y paralelo al plano de la imagen es nominado plano principal. En [Faugeras93a] se introduce de manera formal la geometr´ıa proyectiva y se describe la c´amara proyectiva. Una c´amara se puede modelar como una matriz de proyecci´on ˜ P, de tama˜no 3 ×4, que transforma un punto 3D, M= [X, Y, Z]t, en unas coordenadas de imagen 2D, m= [x, y]t. La matriz de proyecci´on incluye informaci´on interna y externa sobre la c´amara. Los par´ametros extr´ınsecos lo forman el vector de traslaci´on y la matriz de rotaci´on de la c´amara respecto al sistema de referencia global. Los par´ametros intr´ınsecos representan ciertos par´ametros que determina c´omo los objetos de la escena son proyectados en el plano de la imagen y lo conforman: origen de la imagen, distancia focal, tama˜no del p´ıxel y desviaci´on de los ejes. En el problema est´ereo intervienen dos c´amaras, cada una de ellas dispone de una matriz de proyecci´on. Cada punto 3D se proyectar´a generando dos puntos, mym0, en las dos im´agenes; estableciendo la siguiente relaci´on entre las matrices de proyecci´on y los puntos m=˜ PM , m0=˜ P0M. (1.9) 1.2. VISI ´ ON ESTEREOSC ´ OPICA 27 M II’ CC’ e’ e m’ m Figura 1.3: Geometr´ıa epipolar. IeI0son las im´agenes del par est´ereo. CyC0son los focos de las c´amaras. eye0los epipolos. mym0son las proyecciones del punto 3D M. La intersecci´on de la recta que une los focos de las c´amaras CyC0con los planos de imagen de cada c´amara crea dos puntos que se denominan epipolos, eye0. La proyecci´on de la recta que pasa por el punto my el epipolo de su c´amara genera una l´ınea sobre el plano de la imagen de la otra c´amara que se conoce como l´ınea epipolar (figura 1.3). La caracter´ıstica m´as relevante es que el punto m’, proyecci´on del punto 3D en la otra c´amara, pertenece a dicha l´ınea. En principio, el ´area de b´usqueda de correspondencias se reduce a una recta y el desplazamiento de cada p´ıxel se podr´ıa expresar en funci´on de un escalar dentro de la l´ınea epipolar. Esto a priori supone una gran ventaja pero requiere un sistema de c´amaras correctamente calibrado. Si no fuera as´ı, la geometr´ıa epipolar no ser´ıa v´alida y la zona de b´usqueda no se reducir´ıa a dicha l´ınea. No siempre es posible disponer de un sistema de c´amaras fuertemente calibrado. Por ello, se recurre a una configuraci´on m´as sencilla que se conoce como calibrado d´ebil. En [Luong96] se define la restricci´on epipolar a trav´es de una matriz, denominada matriz fundamental, mostrando la relaci´on que existe entre las matrices de proyecci´on de las c´amaras y la matriz fundamental. A partir de la matriz fundamental no se puede calcular las matrices de proyecci´on de las c´amaras. ´ Esta s´olo ofrece informaci´on relativa de una c´amara a la otra, no respecto al sistema de referencia global. En la matriz fundamental se encuentran representados los par´ametros intr´ınsecos – origen de la imagen, distancia focal, tama˜no del p´ıxel y desviaci´on de los ejes – y los par´ametros extr´ınsecos – matriz de rotaci´on y vector de traslaci´on – de las c´amaras. A continuaci´on, se describe formalmente c´omo expresar este desplazamiento bas´andonos en la geometr´ıa epipolar y en la matriz fundamental. Vamos a suponer que el sistema de c´amaras est´a calibrado d´ebilmente y que las c´amaras est´an alineadas horizontalmente. Con esta configuraci´on de c´amaras la matriz fundamental ser´ıa muy simple, ec. (1.10). F=  000 001 0−1 0  .(1.10) Se define un conjunto de pares de p´ıxeles A, ec. (1.11), donde puede encontrar la correspondencia. A=m, m0|my=m0 yy 0 ≤m0 x−mx< k,(1.11) 28 CAP´ ITULO 1. ESTADO DEL ARTE Il(x)g(x)Ir(x) Cuadro 1.2: Il(x) y Ir(x) son las im´agenes izquierda y derecha de un par est´ereo. g(x) representa el mapa de disparidad calculado a partir del par est´ereo. donde (x, y) las componentes horizontal y vertical de cada p´ıxel. kes un umbral que define la zona de b´usqueda de la correspondencia dentro del ´area A. En el caso simple, la l´ınea epipolar. La ecuaci´on de restricci´on epipolar m0tFm = 0 establece que dos puntos est´an en correspondencia, m= (x,1) = (x, y, 1) y m0t= (x0,1) = (x0, y0,1), en dos im´agenes est´an definidas por la matriz fundamental, F[Faugeras01, Hartley03]. Esta definici´on permite estimar el flujo est´ereo solamente sobre ciertas l´ıneas. De alg´un modo reduce la zona de b´usqueda de las correspondencias siguiendo las siguientes ecuaciones: a(x) = f11x+f12y+f13 b(x) = f21x+f22y+f23 c(x) = f31x+f32y+f33 Utilizando esta notaci´on, la l´ınea epipolar ∆ se puede escribir como a(x, y)x0+b(x, y)y0+c(x, y)=0.(1.12) En el cuadro 1.2, llamamos al flujo est´ereo como g(x)=(u(x), v(x))t. El flujo est´ereo depende de una funci´on escalar λ(x) y define una distancia en la l´ınea epipolar de la siguiente forma: u(x) = −λ(x)b(x) √a2(x)+b2(x)−a(x)x+b(x)y+c(x) a2(x)+b2(x)a(x) v(x) = λ(x)a(x) √a2(x)+b2(x)−a(x)x+b(x)y+c(x) a2(x)+b2(x)b(x) En Alvarez/Deriche/S´anchez/Weickert (2002) [Alvarez02b] podemos encontrar descrito con m´as detalle la parametrizaci´on anterior. 1.2.2. Clasificaci´on de los M´etodos Los m´etodos de est´ereo se pueden agrupar en cuatro categor´ıas: (i) basados en caracter´ısticas, (ii) basados en ´areas, (iii) basados en frecuencia y (iv) basados en energ´ıas. 1.2. VISI ´ ON ESTEREOSC ´ OPICA 29 En [Brown03] podemos encontrar un resumen de las distintas t´ecnicas existentes y los ´ultimos avances en visi´on estereosc´opica. Basados en Caracter´ısticas: Los m´etodos que pertenecen a este grupo establecen las correspondencias bas´andose en alg´un tipo de caracter´ıstica extra´ıda de las im´agenes [Grimson85], como curvas [Brint90, Robert91, Nasrabadi92], l´ıneas [Medioni85, Ayache87, McIntosh88] o los bordes de los objetos [Ohta85, Pollard85]. Por un lado, con este tipo de m´etodos es posible estimar mapas de disparidad con muy poca informaci´on (caracter´ısticas extra´ıdas de la imagen). Por otro lado, las soluciones obtenidas no son densas. Basados en ´ Areas: Los mapas de disparidad se calculan a partir de la correlaci´on de ciertas zonas de la imagen asumiendo que existe alg´un tipo de similaridad [Scharstein02, Devernay94, Faugeras93b, Fua93, Nishihara84]. Con este tipo de m´etodos se obtienen muy buenos resultados si se aplican sobre pares est´ereo altamente texturados. Basados en Frecuencias: Los m´etodos basados en frecuencia utilizan la informaci´on de las im´agenes en el dominio de Fourier [Froehlinghaus96, Barron94, Fleet90, Fleet93, Jenkin94, Kuglin75, Wiklund92]. Basados en la Energ´ıa: En este tipo de m´etodos la disparidad se calcula a partir de la minimizaci´on de una energ´ıa que penaliza las desviaciones respecto a las restricciones impuestas en el modelo. El mapa de disparidad se obtiene tras la resoluci´on de un sistema de ecuaciones diferenciales parciales (EDP’s), normalmente del tipo difusi´on-reacci´on. Este tipo de ecuaciones diferenciales se componen de dos t´erminos: uno de reacci´on (lo que hasta hora hemos llamado t´ermino de ligadura) y otro de difusi´on (de regularizaci´on o suavizado). Algunos autores proponen distintas formas de clasificar los m´etodos basados en energ´ıas. Al igual que Barron et al. [Barron94] para el caso del flujo ´optico unos autores agrupan los m´etodos en funci´on de la densidad de los mapas de disparidad: locales y globales. Otros lo hacen en funci´on del tipo de informaci´on que utilizan para realizar las estimaciones: probabil´ısticos yvariacionales. Los mapas de disparidad en los m´etodos locales no son densos mientras en los globales s´ı. Al igual que ocurr´ıa con los m´etodos que estimaban el flujo ´optico, los m´etodos locales calculaban la disparidad de cada p´ıxel a partir de la informaci´on dentro de una ventana centrada en dicho p´ıxel. Para ello, utilizaban alguna caracter´ıstica de la imagen ya sea en color o en escala de grises. Los m´etodos m´as representativos son [Kanade94, Yoon05]. Sin embargo, en los m´etodos globales es necesario la inclusi´on de un t´ermino de regularizaci´on para convertir el problema en bien condicionado. A su vez, 30 CAP´ ITULO 1. ESTADO DEL ARTE estos m´etodos se pueden dividir en varios subtipos: graph–cuts [Kolmogorov01, Boykov04], programaci´on din´amica [Lei06], belief propagation [Klaus06] o m´etodos variacionales [Shah93, Robert96, Mansouri98, Alvarez02b, Kim03, Slesareva05]. Los m´etodos probabil´ısticos calculan las estimaciones en funci´on de la disparidad m´as probable, de forma que ´esta surge de la minimizaci´on de una energ´ıa discreta ([Kolmogorov01, Lei06, Klaus06]). Las soluciones obtenidas por este tipo de m´etodos ofrecen muy buenos resultados, aunque debido a su naturaleza discreta, ´estas son en precisi´on entera. Aunque la precisi´on de las soluciones puede ser suficientemente para muchas aplicaciones como segmentaci´on o la estimaci´on del mapa de disparidad, no lo es tanto para otras como reconstrucci´on 3D donde se requiere una precisi´on mayor. Li y Zucker [Li06] comentan que uno de los inconvenientes que tiene este tipo de t´ecnicas es cuando no se cumple la suposici´on de la disparidad constante por trozos, hecho que ocurre cuando la profundidad en la escena var´ıa suavemente. Un ejemplo de este fen´omeno lo podemos observar en los experimentos hechos en esta tesis con la secuencia del pasillo. En los ´ultimos a˜nos se han propuesto nuevos m´etodos que se basan en la teor´ıa de graph–cuts [Roy98, Ishikawa98, Kolmogorov02, Boykov04, Kolmogorov04]. Todos ellos son una extensi´on de la idea m´aximo–m´ınimo presentada originalmente en [Greig89, Wu93, Roy99, Bobick99]. M´etodos Variacionales Los m´etodos variacionales no tienen algunas de las limitaciones que ofrecen los probabil´ısticos. Por un lado, dada la naturaleza continua de los funcionales de energ´ıa que definen, las soluciones obtenidas son densas y con precisi´on subp´ıxel. Por otro lado, todas las restricciones impuestas en el modelo est´an presentes en la definici´on de la energ´ıa y existe una s´olida teor´ıa matem´atica subyacente. El t´ermino de ligadura asume que una determina propiedad en las im´agenes del par est´ereo no var´ıa. En trabajos como [Robert96, Mansouri98, Alvarez02b] se ha utilizado la suposici´on lambertiana. En otros, como [Slesareva05], se ha incluido varias invarianzas dentro de ese t´ermino. [Slesareva05] es una adaptaci´on del m´etodo de [Papenberg06] para el caso est´ereo. Otra contribuci´on a destacar ha sido la de [Ari07] que ha utilizado en el t´ermino de ligadura la norma L1, a diferencia del tradicional t´ermino cuadr´atico ampliamente usado en la literatura. Las similitudes entre las ideas propuestas para el problema del flujo ´optico y est´ereo se reflejan tambi´en en la clasificaciones hechas en los t´erminos de regularizaci´on. Los t´erminos de regularizaci´on se pueden agrupar en funci´on de la informaci´on utilizada en la difusi´on (image–driven,disparity–driven) o en funci´on de la direcci´on en el que se aplica (isotr´opico,anisotr´opico). Los regularizadores anisotr´opicos siempre han ofrecido mejores resultados que los isotr´opicos ya que son capaces de ajustar la difusi´on en los bordes de los objetos. En la literatura nos encontramos algunos ejemplos de regularizadores anisotr´opicos image–driven como [Alvarez02b, Mansouri98], en los que se preservan los bordes utilizando la informaci´on de la imagen. Se basan en la idea que los p´ıxeles con gradiente grandes reflejan una discontinuidad en la escena (objetos distintos). Sin embargo, en el caso de im´agenes texturadas esta suposici´on no siempre se cumple. El segundo tipo de regularizadores son los disparity–driven. En trabajos como [Slesareva05, Ari07] se utilizan 1.3. CARACTER´ ISTICAS DE LAS IM ´ AGENES 31 este tipo de regularizadores que, con un comportamiento similar a los flow-driven en flujo ´optico, suavizan en funci´on de la informaci´on del mapa de disparidad. [Shah93] propuso el uso de una restricci´on de suavizado que preservara las discontinuidades bas´andose en la disparidad. Posteriormente, [Robert96] describi´o un framework para la definici´on de regularizadores disparity-driven. En [Zimmer08] se ha propuesto un m´etodo que incluye un regularizador anisotr´opico disparity–driven que ofrece dos ventajas. Por un lado, se aprovecha de la precisi´on del anisotr´opico y, por otro lado, el disparity-driven evita la sobresegmentaci´on del mapa disparidad como ocurre en los regularizadores image-driven. Las funciones de robustificaci´on se han utilizado en el problema del flujo ´optico con dos fines: para atenuar los efectos del ruido en las im´agenes ([Black91, Black96b, Bruhn05c]) o para controlar el proceso de difusi´on ([Nesi93, Schn¨orr94b, Weickert01a]). Tradicionalmente, en el caso est´ereo, se ha utilizado como t´ermino de regularizaci´on la funci´on de Tikhonov, [Shah93]. En otros trabajos, [Shah91, Kumar97, Slesareva05], se empleado la funci´on conocida como Total Variation (TV) o la propuesta por Perona-Malik, [Perona90]. La funci´on de Total Variation tiene la ventaja que carece de par´ametros y ha demostrado ser muy efectiva ante la presencia de outliers. Sin embargo, la funci´on de Perona-Malik crea discontinuidades m´as n´ıtidas en el flujo. Las oclusiones son unos fen´omenos que se producen por la existencia de p´ıxeles que s´olo son visibles en una imagen del par est´ereo. Hasta hace unos a˜nos los m´etodos no inclu´ıan ning´un tipo de t´ecnica para la detecci´on de las oclusiones [Robert96, Kumar97, Mansouri98, Alvarez02b, Slesareva05]. En [Ari07] se propone un m´etodo que ofrece muy buenos resultados y que aglutina varias t´ecnicas muy exitosas. Est´a basado en el framework de Mumford-Shah [Mumford89], tiene un t´ermino de regularizaci´on disparity–driven y maneja las oclusiones de forma similar a [Shah93]. Las t´ecnicas multiescala utilizadas para la detecci´on de desplazamientos largos desarrolladas para la estimaci´on del flujo ´optico siguen siendo v´alidas en el problema est´ereo. De hecho, existen trabajos como Alvarez/Deriche/S´anchez/Weickert (2002) [Alvarez02b] en el que se describe un m´etodo para la estimaci´on de la disparidad para desplazamientos largos que utiliza un enfoque similar a Alvarez/Weickert/S´anchez (2000), [Alvarez00]. 1.3. Caracter´ısticas de las Im´agenes En este apartado se describen las caracter´ısticas de las im´agenes y la complejidad que puede tener una escena. Conviene tenerlas en cuenta a la hora de dise˜nar los modelos de energ´ıa ya que los m´etodos variacionales propuestos para la estimaci´on del flujo ´optico y del mapa de disparidad se basan en (i) la idea de invarianza sobre alguna propiedad de las im´agenes y en (ii) una restricci´on sobre el desplazamiento de los p´ıxeles. 1.3.1. Secuencias de Im´agenes El sensor de una c´amara captura la luz procedente de una escena. Esta energ´ıa puede proceder directamente de fuentes de luz o reflejada por los objetos, pudiendo alterar la 38 CAP´ ITULO 1. ESTADO DEL ARTE Figura 1.6: De izquierda a derecha, el desplazamiento horizontal y vertical que se registra entre los frames 9 y 10. Los tonos claros representan desplazamientos positivos y los oscuros, negativos. Los p´ıxeles en negro no se tienen en cuenta y corresponden con el fondo de la escena y las aristas de las torres de m´armol. Figura 1.7: Campo de desplazamiento entre los frames 9 y 10 expresado en campo de vectores. Esta secuencia ha sido ampliamente utilizada en la literatura debido a la diversidad de movimientos presentes en la misma. En la evaluaci´on de los m´etodos se utiliza el flujo ´optico calculado entre los frames 8 y 9. S´olo se dispone del campo de desplazamiento con bastante precisi´on entre estos dos frames. En el resto de la secuencia los campos de desplazamientos disponibles tienen una precisi´on menor. Secuencias Reales En este apartado se describe las secuencias reales utilizadas en los experimentos. A diferencia de las secuencias sint´eticas, ´estas presentan una serie de perturbaciones, como efecto de entrelazado,cambios de iluminaci´on,sombras, etc., que dificultan la correcta estimaci´on del flujo ´optico. Otro de los inconvenientes es la ausencia de ground truth por lo que no es posible obtener resultados cuantitativos. Taxi La secuencia del Taxi de Hamburgo se trata de una de las m´as famosas y m´as utilizadas 1.4. SECUENCIAS DE PRUEBA 39 Figura 1.8: Frames 0, 7 y 14 de la secuencia de Yosemite con nubes. Figura 1.9: Mapas de desplazamiento horizontal y vertical que se registra entre los frames 8 y 9. Los tonos claros representan desplazamientos positivos y los oscuros, negativos. en los trabajos sobre el c´alculo del flujo ´optico. Fue creada por H.H. Nagel (KOGS/IAKS, Universidad de Karlsruhe, Alemania, http://i21www.ira.uka.de/image sequences/). Est´a compuesta por cuarenta frames de tama˜no 256 ×190 p´ıxeles y, en ella, se puede observar el efecto del entrelazado,cambios de iluminaci´on yoclusiones. En la figura 1.11 podemos ver tres frames que componen esta secuencia. En ella, se aprecia como un taxi ejecuta un giro hacia la derecha para entrar en una calle y un coche oscuro situado en la esquina inferior izquierda circula hacia la derecha a gran velocidad. En la esquina inferior derecha, circula un cami´on justo detr´as del taxi y, por ´ultimo, en la esquina superior izquierda un peat´on se desplaza por la acera a una velocidad muy inferior a la que lo hacen los veh´ıculos. El resto de objetos presentes en la secuencia no se mueven. Rheinhafen Otra de las secuencias reales utilizadas en los experimentos es la conocida como Rheinhafen. Fue creada tambi´en por Nagel (http://i21www.ira.uka.de/ image sequences). Se trata de una secuencia en escala de grises compuesta por mil frames de tama˜no 688×565 p´ıxeles. Una c´amara situada a cierta altura captura el tr´afico rodado que circula por una v´ıa. Se observan varios veh´ıculos que circulan en distintas direcciones y velocidades. Pr´oximo a la c´amara se capta una furgoneta que circula a gran velocidad. Al fondo, se ven varios veh´ıculos que se mueven a menor velocidad y que se disponen a cambiar de direcci´on. En esta secuencia se percibe con total claridad el efecto de entrelazado en los contornos de la furgoneta situada cerca de la c´amara. En la figura 1.12 se muestra dos frames no consecutivos en los que se aprecia con mayor claridad el movimiento registrado en la escena. 40 CAP´ ITULO 1. ESTADO DEL ARTE Figura 1.10: Mapa de desplazamiento entre los frames 8 y 9 expresado en campo de vectores. Figura 1.11: Frames 0, 10 y 19 de la secuencia del Taxi de Hamburgo. Secuencias del Meteosat Como parte de esta tesis se ha desarrollado un m´etodo variacional multicanal para la estimaci´on del movimiento de estructuras nubosas. Para evaluar este nuevo m´etodo se han utilizado dos secuencias sat´elites captadas por el Meteosat. El Meteosat es un sat´elite de meteorolog´ıa provisto de una serie de sensores que captan la energ´ıa reflejada por la Tierra en distintos rangos de frecuencia. Estas secuencias han sido tomadas por la segunda generaci´on de sat´elites Meteosat que son capaces de capturar im´agenes de diez bits de cuantificaci´on cada 15 minutos, con un tama˜no de p´ıxel de tres kil´ometros cuadrados para cada uno de los once canales con lo que dispone, desde el infrarrojo hasta el visible. Figura 1.12: Dos frames no consecutivos de la secuencia Rheinhafen. 1.4. SECUENCIAS DE PRUEBA 41 ID del Canal Longitud de Onda Aplicaci´on Principal VIS 0.6 0.63 µm Detecci´on y seguimiento de nubes, identificaci´on de superficie VIS 0.8 0.81 µm Detecci´on y seguimiento de nubes, identificaci´on de superficie NIR 1.6 1.64 µm Discriminaci´on entre nieve, hielo y masas nubosas IR 3.9 3.92 µm Detecci´on de nubes bajas y nieblas durante la noche WV 6.2 6.25 µm Masas de vapor de agua situadas a media altura WV 7.3 7.35 µm Masas de vapor de agua situadas a baja altura IR 8.7 8.70 µm Distinci´on entre agua e hielo IR 9.7 9.66 µm Detecci´on de ozono en la parte inferior de la estratosfera IR 10.8 10.80 µm Estimaci´on de la temperatura de las nubes y de la superficie terrestre IR 12 12.00 µm Estimaci´on de la temperatura de las nubes y de la superficie terrestre IR 13.4 13.40 µm Estimaci´on de nubes a gran altura HRV 0.75 µm Alta resoluci´on espacial Cuadro 1.3: Caracter´ısticas y principales aplicaciones de los canales del Meteosat. En el canal visible cuenta con un sensor de alta resoluci´on que captura im´agenes con un tama˜no de p´ıxel de un kil´ometro cuadrado [Schmetz02]. Cada uno de los once sensores que el sat´elite dispone suministran informaci´on que es utilizada en distintas aplicaciones. En la tabla 1.3, podemos ver resumido la aplicaci´on asociada a cada canal. Sin embargo, la mayor´ıa de las aplicaciones meteorol´ogicas combinan la informaci´on de varios canales, principalmente de estos cuatro, VIS 0.8, WV 6.2, WV 7.3 y IR 10.8 [Schmetz02]. El canal VIS 0.8 capta informaci´on del espectro visible. Este canal permite la identificaci´on y seguimiento de las nubes en la atm´osfera, la monitorizaci´on de la vegetaci´on y de la superficie terrestre. Los canales WV 6.2 y WV 7.3 permiten observar el vapor de agua presente en la atm´osfera as´ı como las corrientes de aire. Tambi´en es posible la localizaci´on de nubes semitransparentes situadas a gran altura [Schmetz93]. Por ´ultimo, el canal IR 10.8 capta informaci´on referente a la temperatura de la superficie terrestre, los oc´eanos y de la parte superior de las nubes [Inoue87]. Secuencia del Hurac´an Vince Esta secuencia capta la presencia del hurac´an Vince en el oc´eano Atl´antico. Este hurac´an supone un efecto ins´olito ya que se form´o en una zona situada demasiada al Este de donde se suele producir habitualmente. Este hurac´an categor´ıa 1 (seg´un la escala SaffirSimpson) se form´o a las 18:00 el 9 Octubre del 2005 al noroeste de Funchal (Islas Madeira). A partir de entonces, empez´o a perder fuerza hasta convertirse en una tormenta tropical casi seis horas despu´es [Franklin06]. En la figura 1.13, se puede observar las im´agenes 42 CAP´ ITULO 1. ESTADO DEL ARTE Figura 1.13: Secuencia de Vince. De izquierda a derecha y de arriba hacia abajo, se muestra el canal visible 0,81µm, los canales de vapor de agua 6,25µm y 7,35µm y a la derecha el canal de infrarrojo, 10,80µm. de cuatro canales tomadas por el Meteosat. Las zonas de inter´es de esta secuencia son dos: el hurac´an Vince y la masa nubosa presente en el Atl´antico Norte (pr´oxima a las Islas Brit´anicas). Los movimientos que principalmente se registran en esta secuencia son rotacionales. Secuencia del Atl´antico Norte(June 5th, 2004) Esta secuencia fue tomada por el Meteosat el 5 de Junio del 2004. En ella se muestra los fen´omenos atmosf´ericos acontecidos aquel d´ıa en el hemisferio norte en la franja europeonorteafricana. En la figura 1.14, se puede observar las im´agenes de la atm´osfera de cuatro canales captadas por el Meteosat. En esta secuencia cabe destacar las dos grandes masas nubosas. Por un lado, la procedente del Atl´antico Norte aproxim´andose hacia las islas Brit´anicas y, por otro lado, la gran masa nube que se aproxima desde el oeste hacia la pen´ınsula Ib´erica. El movimiento predominante en estas dos grandes masas nubosas es rotacional. En el resto de la secuencia el movimiento es menos severo y principalmente traslacional. 1.4.2. Visi´on Estereosc´opica Al igual que hac´ıamos en el apartado anterior vamos a comentar las distintas secuencias de pares est´ereo que se han utilizado en los experimentos. Durante las pruebas se han 1.4. SECUENCIAS DE PRUEBA 43 Figura 1.14: Secuencia del Atl´antico Norte. De izquierda a derecha y de arriba hacia abajo, se muestra el canal visible 0,81µm, los canales de vapor de agua 6,25µm y 7,35µm y a la derecha el canal de infrarrojo, 10,80µm. empleado tanto im´agenes reales como sint´eticas. Algunas de ellas han sido ampliamente utilizadas en la literatura y otras han sido creadas por miembros del grupo AMI (ULPGC). Un par est´ereo se compone de dos im´agenes tomadas en el mismo instante de tiempo. Para establecer las correspondencias entre los puntos es necesario que el sistema de c´amaras est´a calibrado. Existen m´ultiples configuraciones pero, normalmente, para simplificar el problema de la estimaci´on de los mapas de disparidad se suele situar ambas c´amaras de forma frontoparalela. De este modo, los focos est´an situados en el mismo plano de proyecci´on y las l´ıneas epipolares son paralelas. En las c´amara se ha seleccionado una configuraci´on horizontal. El desplazamiento de los objetos s´olo ocurre en una sola direcci´on. La matriz fundamental asociada a las c´amaras en este caso es F=  000 001 0−1 0  (1.13) En todas las secuencias sint´eticas que hemos utilizado en los experimentos los focos de las c´amaras est´an situadas en el mismo plano horizontal. En las secuencias reales se ha recurrido a un proceso de rectificaci´on para los desplazamientos se produjesen en un s´olo eje. 44 CAP´ ITULO 1. ESTADO DEL ARTE Figura 1.15: Distintos pares est´ereos de la secuencia del Cilindro. En la primera columna, las im´agenes de la c´amara izquierda. En la segunda columna, las im´agenes de la c´amara derecha. Secuencias Sint´eticas En nuestros experimentos se han utilizado seis secuencias sint´eticas: (1) Cilindro, (2) Cilindro y Esfera, (3) Venus, (4) Map, (5) Sawtooth y (6) Corridor. Las dos primeras han sido creadas por el grupo AMI, las tres siguientes forman parte de la base de datos 2001 de Middlebury (http://vision.middlebury.edu/stereo/data/scenes2001/) y, por ´ultimo, el par est´ereo del pasillo que ha sido creado por el grupo de visi´on por ordenador de la Universidad de Bonn. Cilindro La secuencia conocida como Cilindro consta de trece pares de im´agenes de tama˜no 800 ×600 p´ıxeles. En ella, se observa c´omo se desplaza de izquierda a derecha un objeto cil´ındrico. Este objeto tiene una textura compuesta por unas hebras de lana entrelazadas. En el fondo de la imagen se ha planchado una textura de tono claro para que cause contraste con el cilindro. La velocidad del cilindro es de unos ocho p´ıxeles por imagen. En la figura 1.15 se muestra distintos pares est´ereo de la secuencia. En el mapa de disparidad (fig. 1.16) se observa que el ´unico movimiento est´a en el cilindro, el fondo es est´atico. Dado que los p´ıxeles del cilindro no est´an a la misma 1.4. SECUENCIAS DE PRUEBA 45 Figura 1.16: De izquierda a derecha los mapa de disparidad de los pares est´ereo mostrados en la figura 1.15. Los p´ıxeles que registren un mayor desplazamiento tendr´an un tono claro y los que se muevan a menor velocidad un tono m´as oscuro (gris´aceo). En el fondo de la escena no se registra movimiento alguno por lo que el tono de los p´ıxeles es negro. profundidad el desplazamiento var´ıa en funci´on de la distancia a la c´amara (menor distancia mayor desplazamiento y viceversa). Cilindro y Esfera La secuencia del Cilindro y Esfera se trata de una versi´on algo m´as compleja que la del Cilindro. En ella, se observan dos objetos en movimiento a distinta velocidad: un cilindro y una esfera (figura 1.17). El desplazamiento se produce de derecha a izquierda. La superficie de los objetos est´a recubierta por una serie de texturas. En primer plano encontramos el cilindro que se desplaza a una velocidad de once p´ıxeles por imagen. Justo detr´as del cilindro se observa una esfera. El movimiento de este objeto es de alrededor de seis p´ıxeles. Cerca del fondo de la escena se aprecia un panel est´atico con una textura clara que ocupa casi la totalidad de la imagen. En la figura 1.18 se muestra los mapas de disparidad de la secuencia est´ereo. En esas im´agenes se aprecia, en funci´on de la distancia a la c´amara, los tres objetos con distinta tonalidad. Venus La base de datos de Middlebury incluye una serie de secuencias de im´agenes y pares est´ereo tanto reales como sint´eticas (http://vision.middlebury.edu/stereo/data/). Venus es un par est´ereo sint´etico incluido en dicha base de datos. En la figura 1.19 podemos ver las im´agenes de este par. En ellas se observa un lienzo al fondo y dos texturas situadas en primer plano: una p´agina de un peri´odico de deportes y otra de un p´oster. En la figura 1.20, se aprecia el mapa de disparidad de Venus. Como es de esperar el desplazamiento de los p´ıxeles de los objetos m´as pr´oximos a las c´amaras es mayor. Map El par est´ereo conocido como Map tambi´en forma parte de la base de datos de Middlebury. En este par se observa un objeto rectangular situado en primer plano sobre un fondo. El objeto tiene una textura oscura mientras que el fondo est´a compuesto por un 46 CAP´ ITULO 1. ESTADO DEL ARTE Figura 1.17: Distintos pares est´ereos de la secuencia del Cilindro y la Esfera. En la primera columna, las im´agenes de la c´amara izquierda. En la segunda columna, las im´agenes de la c´amara derecha. Figura 1.18: De izquierda a derecha los mapa de disparidad de los pares est´ereo mostrados en la figura 1.17. Los p´ıxeles que registren un mayor desplazamiento tendr´an un tono claro y los que se muevan a menor velocidad un tono m´as oscuro (gris´aceo). En los mapas de disparidad se aprecia en distinta tonalidad los tres objetos en movimiento: el cilindro (el m´as pr´oximo a la c´amara), la esfera (en una posici´on intermedia) y el panel (situado detr´as de los dos objetos anteriores). La parte del fondo de la imagen que no queda oculta por el panel no se aprecia ning´un desplazamiento por lo que el tono de los p´ıxeles es negro. 1.4. SECUENCIAS DE PRUEBA 47 Figura 1.19: Par est´ereo de Venus. Figura 1.20: El mapa de disparidad de Venus. patr´on que se repite pero con distintas orientaciones. En las figuras 1.21 y 1.22 tenemos el par est´ereo de Map y su mapa de disparidad, respectivamente. Sawtooth El par est´ereo conocido como Sawtooth est´a integrado en Middlebury. Dos lienzos (uno en tonos claros y otro oscuro) conforman el fondo de la escena. En primer plano, tenemos una textura gris´acea en cuya parte superior se asemeja a dientes de sierra, de ah´ı el nombre del par est´ereo. En las figuras 1.23 y 1.24 podemos observar el par est´ereo de Sawtooth y el mapa de disparidad asociado. Secuencia del Pasillo Figura 1.21: Par est´ereo de Map. 54 CAP´ ITULO 2. ESTIMACI ´ ON DEL FLUJO ´ OPTICO EN SECUENCIAS DE IM ´ AGENES Como parte principal de esta tesis hemos querido dedicar un cap´ıtulo completo al problema de la estimaci´on del flujo ´optico. Siguiendo la tendencia de los ´ultimos a˜nos todos los m´etodos propuestos en este cap´ıtulo incluyen informaci´on procedente de varias im´agenes y manejo de largos desplazamientos. Antes de entrar en detalle en cada uno de ellos se comentar´an las contribuciones hechas a la literatura y se mostrar´a una generalizaci´on de la estructura del modelo de energ´ıa que es compartida por los tres m´etodos descritos en este cap´ıtulo. 2.1.1. Contribuciones de este Cap´ıtulo Las contribuciones de los m´etodos descritos en este cap´ıtulo son las siguientes: M´etodo Variacional Multicanal: Extensi´on del Modelo de Alvarez/Weickert/S´anchez (2000) [Alvarez00]: Bas´andonos en el funcional de energ´ıa propuesto por [Alvarez00] se ha realizado una serie de modificaciones para la inclusi´on de la informaci´on multicanal. Este m´etodo se ha utilizado para la detecci´on de masas nubes en secuencias sat´elites multiespectrales. Los cambios introducidos en el modelo de energ´ıa son los siguientes. El t´ermino de ligadura se ha definido como una suma ponderada de la suposici´on lambertiana de cada uno de los canales y, el t´ermino de suavizado, el operador de Nagel-Enkelmann modificado cuyo gradiente aglutina informaci´on de todos los canales. Se han definido diversas estrategias para la construcci´on de dicho gradiente. Con las modificaciones realizadas se ha constatado una mejora en las estimaciones respecto al modelo original. M´etodo Variacional que incorpora un T´ermino de Regularizaci´on exclusivamente Temporal no Continuo: En este trabajo se ha desarrollado un m´etodo variacional espaciotemporal que desacopla la regularizaci´on temporal de la espacial. El desacoplamiento del t´ermino de suavizado en dos, uno espacial y otro temporal, pretende evitar algunas de las limitaciones presentes en trabajos con una regularizaci´on temporal continua basada en derivadas parciales temporales, como en [Papenberg06]. Nuestro m´etodo realiza convenientemente la regularizaci´on temporal teniendo en cuenta la presencia de desplazamientos largos sin afectar a la regularizaci´on espacial. M´etodo Variacional basado en el An´alisis Espectral del Tensor de Movimiento: En este trabajo se ha desarrollado un m´etodo variacional espaciotemporal que intercambia informaci´on complementaria entre el t´ermino de ligadura y el de suavizado. El proceso de difusi´on s´olo se aplica en las direcciones ortogonales al flujo dominante. La combinaci´on de varias t´ecnicas existentes posibilita la creaci´on de un framework capaz de representar cualquier modelo de una forma compacta, sencilla, f´acilmente extensible y adaptable. El tensor de movimiento permite representar cualquier invarianza en una notaci´on compacta. La descomposici´on espectral de este tensor asegura la f´acil identificaci´on 2.1. INTRODUCCI ´ ON 55 de las direcciones dominantes del flujo al mismo tiempo que define una notaci´on independiente de las invarianzas definidas en el modelo de energ´ıa. La inclusi´on de t´erminos no-cuadr´aticos no s´olo aporta robustez a las estimaciones frente al ruido sino que la combinaci´on de las funciones de robustificaci´on junto a la descomposici´on espectral del tensor de movimiento permite la construcci´on de cuatro prototipos capaces de mejorar las estimaciones de m´etodos de caracter´ısticas similares. Todas las aportaciones hechas en este cap´ıtulo hacen uso intensivo de la informaci´on procedente de m´ultiples im´agenes, ya sean del mismo o de distintos canales. Por ello, antes de entrar en detalle con las distintas aportaciones conviene definir un modelo de energ´ıa general que aglutine todas las ideas multicanal descritas en cada uno de los m´etodos. 56 CAP´ ITULO 2. ESTIMACI ´ ON DEL FLUJO ´ OPTICO EN SECUENCIAS DE IM ´ AGENES 2.2. Generalizaci´on de los Modelos de Energ´ıa Los m´etodos variacionales se basan en la idea de que la soluci´on se obtiene mediante la minimizaci´on de un funcional de energ´ıa. En esta secci´on se describe un modelo de energ´ıa que recoge pr´acticamente todas las variantes presentadas en la literatura. A partir de la ecuaci´on de restricci´on de flujo ´optico vamos generalizando dicho modelo gracias a la incorporaci´on de una serie de t´ecnicas ampliamente utilizadas en la literatura. Para facilitar la comprensi´on al lector comenzaremos con la obtenci´on de un modelo de energ´ıa en el dominio continuo que luego extenderemos al dominio discontinuo. 2.2.1. Modelos de Energ´ıa Continuos Los modelos de energ´ıa empleados en el problema de la estimaci´on del flujo ´optico normalmente se componen de dos t´erminos: (1) el de ligadura y (2) el de regularizaci´on (o suavizado). Como ya se coment´o en la secci´on 1.1 el t´ermino de ligadura asume que cierta propiedad en la imagen no var´ıa a lo largo del tiempo mientras que el de suavizado impone cierta restricci´on de suavidad en el flujo. La suposici´on m´as simple expresada por el t´ermino de ligadura se puede representar tal que f(x)−f(x+h(x)) = 0,(2.1) donde fes una propiedad en la imagen, x:= (x, y, t) las coordenadas en la secuencia de im´agenes y h(x) := (u(x), v(x),1)>es una funci´on que representa el desplazamiento de los p´ıxeles en la secuencia de im´agenes. En la ecuaci´on (2.1) se aprecia una no linealidad que se suele sortear con del desarrollo de Taylor y la eliminaci´on de los t´erminos de orden superior f(x+h(x)) = f(x+δx, y +δy, t +δt) = f(x, y, t) + ∂f ∂xδx +∂f ∂y δy +∂f ∂t δt +O(kh(x)nk),(2.2) Sustituyendo la ecuaci´on (2.2) en la expresi´on (2.1) obtendr´ıamos la ecuaci´on de restricci´on del flujo ´optico (ec. (2.3)). Esta expresi´on s´olo es v´alida cuando el desplazamiento de los objetos es suave, o dicho de otra forma, las derivadas de la expresi´on son computables. f(x)−f(x+h(x)) ≈∂f ∂x δx δt +∂f ∂y δy δt +∂f ∂t δt δt = 0, que tambi´en se suele representar como fx(x)u(x) + fy(x)v(x) + ft(x) = 0,(2.3) donde u(x) := δx δt v(x) := δy δt yfx, fy, ftindican derivadas parciales. La soluci´on de la expresi´on anterior se puede estimar mediante m´ınimos cuadrados (fxu+fyv+ft)2= (h(x)∇f)2= 0, 2.2. GENERALIZACI ´ ON DE LOS MODELOS DE ENERG´ IA 57 Una forma de extender la idea expresada en la ecuaci´on anterior consistir´ıa en la inclusi´on de varias invarianzas en el t´ermino de ligadura con el objetivo de detectar ciertos movimientos que s´olo son perceptibles para algunas invarianzas. De esta forma, aunque se viole la suposici´on de una invarianza el desplazamiento puede ser identificado por otra. En (2.4) expresamos el t´ermino de ligadura como un conglomerado ponderado de varias invarianzas Nc X c=1 γcfc x(x)u(x) + fc y(x)v(x) + fc t(x)2= Nc X c=1 γc(h(x)∇fc(x))2,(2.4) Ncel n´umero total de invarianzas definidas en el t´ermino de ligadura y cel ´ındice que especifica cada una de ellas. γces un peso que indica la importancia de cada invarianza. En principio, hemos supuesto la forma cuadr´atica como la ´unica forma de minimizar el t´ermino de ligadura pero existen otras, como la norma L1. Por ello, para generalizar a´un m´as la ecuaci´on anterior (2.4) la expresaremos como ZΩ Nc X c=1 γcDc(∇fc(x),h(x)) ,(2.5) donde D(.) podr´ıa variar dependiendo de la naturaleza de las im´agenes y Ω el dominio de la imagen. Por ejemplo, las funciones de robustificaci´on (propuestas por algunos autores [Black92, Black96a, M´emin98a]) ser´ıan una alternativa a la forma cuadr´atica para atenuar el efecto de los outliers. Normalmente, los funcionales de energ´ıa est´an compuestos por un segundo t´ermino denominado de regularizaci´on osuavizado. Los primeros t´erminos de regularizaci´on utilizaban ´unicamente informaci´on espacial. Unos de los primeros fue el de Horn-Schunck [Horn81]. Supone que el campo de desplazamiento es suave y este t´ermino suaviza en funci´on del m´odulo del gradiente del flujo al cuadrado, ec. (2.6). k∇2hk2=k∇2uk2+k∇2vk2.(2.6) A partir de este t´ermino de regularizaci´on surgieron otros que se pueden representar mediante la ecuaci´on (2.7) ZΩR(∇I(x),∇h(x)) .(2.7) Todav´ıa es posible extender la ecuaci´on anterior considerando, al igual que se hizo con el t´ermino de ligadura, ec. (2.5), la posibilidad de usar otras propiedades de la imagen distintas a la intensidad de los p´ıxeles, ec.(2.8). ZΩ Nc X c=1 αcRc(∇fc(x),∇h(x)) ,(2.8) donde cse trata de un ´ındice asociado a cada caracter´ıstica de la imagen, αcun peso de cada t´ermino de suavizado y R(.) una funci´on de robustificaci´on. 58 CAP´ ITULO 2. ESTIMACI ´ ON DEL FLUJO ´ OPTICO EN SECUENCIAS DE IM ´ AGENES La ecuaci´on (2.9) resume para el caso continuo el funcional de energ´ıa gen´erico que estima el flujo ´optico para una secuencia de im´agenes teniendo en cuenta la informaci´on de los frames vecinos. E(h) = ZΩ Nc X c=1 γcDc(∇fc,h)dx+ZΩ Nc X c=1 αcRc(∇fc,∇h)dx.(2.9) 2.2.2. Modelos de Energ´ıa No Continuos Secuencial En este apartado vamos a tratar de definir un modelo de energ´ıa gen´erico que englobe todas las situaciones para el caso discontinuo en el que se tenga en cuenta la informaci´on aportada por los frames vecinos. La ecuaci´on de restricci´on del flujo ´optico s´olo es v´alida cuando el desplazamiento de los objetos es peque˜no. Cuando los desplazamientos entre frames son grandes esta ecuaci´on deja de ser v´alida. Por lo tanto, es necesario formular una nueva expresi´on que sea v´alida para cualquier situaci´on, tanto para desplazamientos cortos como largos. En el caso discontinuo podr´ıamos aplicar la siguiente expresi´on f1(˜x)−f2(˜x +h1(˜x)) = 0,(2.10) donde ˜x := (x, y), ht(˜x) := h(˜x, t) := (u(˜x, t), v(˜x, t),1)>yfi:= f(˜x, i). En la ecuaci´on (2.10) se asume que cierta propiedad de la imagen en dos p´ıxeles en correspondencia en dos im´agenes consecutivas es constante. En algunos trabajos como [Papenberg06] se han incluido varias invarianzas en el t´ermino de ligadura. Al igual que se desarroll´o en el caso continuo podemos expresar el t´ermino de ligadura como un conglomerado ponderado de varias invarianzas ahora en el plano discontinuo. Nc X c=1 γc(fc 1(˜x)−fc 2(˜x +h1(˜x))2. Ahora debemos utilizar una secuencia de im´agenes y el flujo ´optico es calculado entre cada par de im´agenes de la secuencia. En vez de utilizar esta invarianza s´olo entre dos frames podemos ampliar esta idea a toda la secuencia N−1 X i=1 Nc X c=1 γcfc i(˜x)−fc i+1(˜x +hi(˜x)2.(2.11) Un t´ermino de ligadura gen´erico deber´ıa quedar expresado por la definici´on de m´ultiples invarianzas y la funci´on de robustificaci´on que las englobe, ec. (2.12). N−1 X i=1 Nc X c=1 γcDc(fi,c(˜x), fi+1,c(˜x),hi(˜x)) .(2.12) En la literatura podemos observar c´omo los distintos regularizadores propuestos tienen en com´un la utilizaci´on del gradiente del flujo para suavizar la estimaci´on hecha por el 2.2. GENERALIZACI ´ ON DE LOS MODELOS DE ENERG´ IA 59 t´ermino de ligadura. La ecuaci´on (2.13) representa un t´ermino de regularizaci´on gen´erico que utiliza el gradiente del flujo y de la imagen de toda la secuencia N−1 X i=1 R(∇fi(˜x),∇hi(˜x)) .(2.13) La ecuaci´on anterior se puede ampliar considerando la posibilidad de usar cualquier propiedad de la imagen, ec. (2.14). N−1 X i=1 Nc X c=1 αcRc(∇fc i(˜x),∇hi(˜x)) .(2.14) La ecuaci´on (2.15) resume el funcional de energ´ıa que estima el flujo ´optico para una secuencia de im´agenes teniendo en cuenta la informaci´on de los frames vecinos. E(h1,...,hN−1) = ZΩ N−1 X i=1 Nc X c=1 γcDc(fi,c, fi+1,c,hi)d˜x +ZΩ N−1 X i=1 Nc X c=1 αcRc(∇fc i,∇hi)d˜x. (2.15) 2.2.3. Modelos de Energ´ıa No Continuos Aleatorio En la definici´on de los modelos de energ´ıa no continuos se asume ´ımplicitamente que los flujos ´opticos se deben estimar s´olo a partir de la informaci´on de los frames vecinos. Esto es una limitaci´on que se impone al modelo s´olo con el fin de simplificar la notaci´on y reducir la complejidad del funcional de energ´ıa. En este apartado vamos a describir la estructura de un modelo de energ´ıa no continuo que tiene en cuenta la informaci´on de cualquier frame de la secuencia de im´agenes para la estimaci´on del flujo ´optico. En un modelo no continuo secuencial el t´ermino de ligadura se defin´ıa como N−1 X i=1 Nc X c=1 γcDc(fi,c(˜x), fi+1,c(˜x),hi(˜x)) .(2.16) La generalizaci´on del t´ermino de ligadura se puede llevar a un paso m´as all´a, de forma, que las invarianzas no se tengan que cumplir ´unicamente entre frames consecutivos sino para cualquier frame de la secuencia. En el caso m´as simple el t´ermino de ligadura quedar´ıa expresado como Nc X c=1 γcfc i(˜x)−fc j(˜x +hij(˜x))2.(2.17) donde hij(˜x) es el flujo ´optico calculado entre dos frames, i, j, cualesquiera de una secuencia de im´agenes. Si ampliamos esta idea para cualquier par de im´agenes de la secuencia obtendr´ıamos la ecuaci´on N X i=1 N X j=1,j6=i Nc X c=1 γcDc(fi,c(˜x), fj,c(˜x),hij(˜x)) .(2.18) 60 CAP´ ITULO 2. ESTIMACI ´ ON DEL FLUJO ´ OPTICO EN SECUENCIAS DE IM ´ AGENES Como es de esperar el incremento de los grados de libertad en el t´ermino de ligadura acarrea un aumento de las inc´ognitas a resolver en el sistema de ecuaciones derivado a partir del modelo de energ´ıa. Con esta generalizaci´on hemos querido poner de manifiesto las m´ultiples combinaciones existentes en la definici´on del t´ermino de ligadura en los funcionales de energ´ıa. Aunque a efectos pr´acticos, por simplicidad, se recurrir´a a versiones m´as simples de (2.18), como puede ser (2.11). Si nos fijamos en los t´erminos de ligadura de los m´etodos propuestos en la literatura son casos particulares de la ecuaci´on (2.18). En un modelo no continuo secuencial el t´ermino de suavizado se defin´ıa como N−1 X i=1 Nc X c=1 αsRc(∇fc i(˜x),∇hi(˜x)) . Si ampliamos la ecuaci´on anterior considerando la posibilidad de regularizar el flujo con informaci´on de frames no consecutivos obtendr´ıamos N X i=1 N X j=1,j6=i Nc X c=1 αcRc(∇fc i(˜x),∇hij(˜x)) .(2.19) A partir de las generalizaciones de los t´erminos de ligadura, ec. (2.18), y suavizado, ec. (2.19), se ha definido un modelo de energ´ıa que a diferencia de (2.15) utiliza la informaci´on procedente de frames cualesquiera. E(h12,...,hN(N−1)) = ZΩ N X i=1 N X j=1,j6=i Nc X c=1 γcDc(fi,c, fj,c,hij)d˜x +ZΩ N X i=1 N X j=1,j6=i Nc X c=1 αcRc(∇fc i,∇hij)d˜x.(2.20) 2.3. M ´ ETODO VARIACIONAL MULTICANAL 61 2.3. An´alisis de Im´agenes Sat´elite Multicanal usando M´etodos Variacionales El an´alisis de las im´agenes sat´elites se ha convertido en un campo de estudio muy activo en los ´ultimos a˜nos. Los avances en la tecnolog´ıa de los sensores ha hecho posible la obtenci´on de im´agenes de mayor resoluci´on y en un mayor rango de frecuencias. Esto ha impulsado el campo de la teledetecci´on y de la meteorolog´ıa. La estimaci´on del movimiento de las nubes a partir de las im´agenes sat´elites tiene mucha aplicaci´on en la climatolog´ıa y meteorolog´ıa ([Hasler90]). En la literatura se han propuesto distintas clases de t´ecnicas para la estimaci´on del movimiento de las nubes, como es local cross-correlation ([Leese71, Phillips72, Schmetz93]) ycross-correlation combinado con relaxation labeling ([Wu95, Evans06]), estimaci´on de movimiento usando im´agenes est´ereo ([Young90, Kambhamettu95]), redes neuronales ([Cˆot´e95]) t´ecnicas block-matching ([Brad02]), ajuste local ([LZ01]) o m´etodos variacionales ([Corpetti02]). A pesar de la variedad de las t´ecnicas utilizadas predomina el uso de la correlaci´on. Tiene la ventaja de ser robusta frente a cambios de intensidad global, pero el inconveniente del alto coste computacional que requieren las estimaciones y la no integraci´on de una restricci´on de regularizaci´on global. A la hora de buscar las correspondencias de los p´ıxeles se busca la similaridad entre ventanas o patrones alrededor de un p´ıxel. Para obtener un resultado denso ser´ıa necesario la aplicaci´on de esta operaci´on en cada p´ıxel de la imagen. Dado lo costoso computacionalmente hablando de este tipo de t´ecnicas se selecciona un conjunto de puntos en la imagen en los que se aplica la correlaci´on. A priori se desconoce cu´ales son los puntos de inter´es en la imagen por lo que se seleccionada un conjunto de puntos equidistantes entre ellos alineados representando una rejilla (ver figura 2.1). Para obtener un mapa de desplazamiento denso es necesario aplicar alguna t´ecnica de interpolaci´on a partir de los valores obtenidos en la rejilla. Por un lado, la densidad de la rejilla debe ser lo suficientemente peque˜na para que todos los desplazamientos sean detectados. Por otro lado, esta densidad debe ser lo suficientemente grande para que la soluci´on sea obtenida en el menor tiempo posible. Ser´ıa deseable que este tipo de t´ecnicas incluyese alg´un tipo de restricci´on que asegurase la coherencia de los resultados, es decir, que el desplazamiento sea una funci´on suave. El local cross-correlation puede resultar ´util cuando los movimientos son r´ıgidos pero no funciona tan bien en el caso de desplazamientos r´apidos no r´ıgidos. Los sat´elites captan la radiaci´on luminosa en distintos rangos de frecuencias creando as´ı un amplio abanico de im´agenes. Tradicionalmente para el seguimiento de las nubes se han utilizado las im´agenes del canal infrarrojo (IR: 10.5-12.5 µm, [Leese71, Schmetz93]). El canal visible permite la detecci´on y seguimiento de nubes de baja altura ([Ottenbacher97, LZ01]). Dada la riqueza de informaci´on multiespectral ofrecida por los sat´elites, ser´ıa l´ogico pensar que la combinaci´on de toda esta informaci´on permitir´ıa mejorar el seguimiento de las masas nubosas. Recientemente, se han propuesto algunos m´etodos que combinan datos multicanal, como es una t´ecnica cross-correlation multicanal que usa el canal visible e infrarrojo ([Evans06]) o una estimaci´on multicapa de las nubes tomando informaci´on de varios canales ([H´eas06]). Los m´etodos variacionales han tenido bastante ´exito en el campo de la estimaci´on 62 CAP´ ITULO 2. ESTIMACI ´ ON DEL FLUJO ´ OPTICO EN SECUENCIAS DE IM ´ AGENES Figura 2.1: Un ejemplo de una imagen sat´elite en la que se ha resaltado en rojo los puntos de la rejilla en los que se aplica la correlaci´on. del flujo ´optico. Estos m´etodos imponen una restricci´on de suavizado global creando resultados coherentes. Las soluciones son densas por lo que no es necesario ning´un tipo de interpolaci´on una vez hecho los c´alculos. La combinaci´on de varios canales para la estimaci´on del flujo ´optico no es nada novedoso en los m´etodos variacionales. En la literatura existe un gran n´umero de m´etodos que utilizan secuencias de color (tres canales, rojo, verde y azul). Sin embargo, el m´etodo que presentamos en esta secci´on combina la informaci´on de los canales de un modo diferente a como se suele hacer en el caso de color. En secuencias de color la informaci´on est´a repartida por los tres canales por lo que es imprescindible la combinaci´on de todos ellos para calcular los desplazamientos. En las secuencias multicanal de los sat´elites la informaci´on se concentra en cada canal, es posible estimar el movimiento de las estructuras nubosas en cada canal independiente. Pero hay determinadas estructuras que se detectan mejor en unos canales que en otros. Por ese motivo, la combinaci´on de los canales sat´elites enriquece las estimaciones, de forma que cuando un canal no ofrece una buena estimaci´on, ´esta ser´ıa aportada por otro donde el movimiento se detectara con mayor nitidez. En esta secci´on presentamos un m´etodo variacional multicanal aplicado al an´alisis de im´agenes sat´elites multiespectrales. Este m´etodo combina la informaci´on de varios canales para mejorar la estimaci´on del desplazamiento de las nubes. Se trata de una variante del modelo de Nagel-Enkelmann, [Nagel86]. Comenzaremos, en el apartado 2.3.1, presentando un m´etodo variacional monocanal [Alvarez00], al que le aplicaremos una serie de peque˜nas modificaciones para que pueda utilizar informaci´on multicanal, apartado 2.3.2. Dado que la informaci´on de relevancia no se reparte equitativamente entre los distintos canales ser´ıa deseable establecer algunas estrategias para priorizar unos canales frente a otros. Por ´ultimo, mostraremos en el apartado 2.3.5 los resultados cuantitativos y cualitativos obtenidos con las dos secuencias de sat´elites comentadas en el apartado 1.4.1. En estos experimentos se refleja la importancia y mejora que ofrece el m´etodo multicanal frente a 2.3. M ´ ETODO VARIACIONAL MULTICANAL 63 monocanal. 2.3.1. M´etodo Variacional Monocanal La estimaci´on del movimiento de las masas nubosas por la atm´osfera se trata de un problema de c´alculo del flujo ´optico. Disponemos de unas im´agenes tomadas por un sat´elite donde las estructuras nubosas son unos p´ıxeles que se desplazan por las im´agenes. Para estimar este movimiento hemos utilizado un m´etodo variacional. Los m´etodos variacionales m´as simples utilizan ´unicamente un s´olo canal, por ejemplo una imagen en escala de grises. Otros m´as sofisticados combinan la informaci´on de varios canales. En este apartado se presenta el modelo de energ´ıa utilizado por el m´etodo monocanal. Este modelo corresponde con el descrito en [Alvarez00]. Se trata de un m´etodo 2D compuesto por dos t´erminos: uno de ligadura y otro de suavizado. El t´ermino de ligadura incluye la suposici´on lambertiana mientras que el de suavizado el operador de NagelEnkelmann [Nagel86] con algunas mejoras. El funcional de energ´ıa a minimizar es E(h) = ZΩ (I1(˜x)−I2(˜x +h(˜x)))2d˜x +αZΩ tr(∇htD∇h)d˜x,(2.21) donde ˜x es un punto que pertenece al dominio Ω, αes una constante que pondera el t´ermino de suavizado, ∇es el operador del gradiente y, Des una matriz de proyecci´on regularizada en la direcci´on ortogonal a ∇I1. Operador de Nagel–Enkelmann La matriz de difusi´on Dse define como: D(∇I1) = 1 k∇I1k2+ 2ζ2ξξt+ζ2Id,(2.22) donde ξ= (∂I1 ∂y ,−∂I1 ∂x )tes un vector ortogonal a ∇I1,k.kindica la norma, Id es la matriz identidad y, ζes un coeficiente que determina el comportamiento isotr´opico del suavizado e inhibe los bordes borrosos cuando la magnitud del gradiente es alta: k∇I1k  ζ. Los par´ametros de entrada de este m´etodo son Cyλ∈(0,1). A partir de estas variables se calculan dos pesos utilizados en el modelo de energ´ıa como son C y ζdonde α=C max (|(∇Gσ∗I1)(x)|2),(2.23) λ=Zζ 0 H|∇Gσ∗I1|(z)dz, (2.24) donde Gσ∗I1representa la convoluci´on de I1con una gaussiana de desviaci´on est´andar σ,H|∇Gσ∗I1|(z) representa el histograma normalizado de |∇Gσ∗I1|.λse conoce como fracci´on isotr´opica. Cuando λ→0, la difusi´on aplicada es anisotr´opica mientras λ→1, lo hace isotr´opico. Esta normalizaci´on de αyζpermite a la energ´ıa ser invariante frente a cambios de intensidad del tipo (I1, I2)→(kI1, kI2). 70 CAP´ ITULO 2. ESTIMACI ´ ON DEL FLUJO ´ OPTICO EN SECUENCIAS DE IM ´ AGENES Figura 2.2: Secuencia de Vince. De izquierda a derecha y de arriba hacia abajo, se muestra el canal visible 0,81µm, los canales de vapor de agua 6,25µm y 7,35µm y a la derecha el canal de infrarrojo, 10,80µm. canal. La mejora ha sido de 27.34 % y 27.58 % en el AEE y AAE, respectivamente. Existe una peque˜na diferente de 1.61 % (AEE) y 1.92 % (AAE) entre las dos estrategias del m´etodo multicanal. En la figura 2.5 se muestra una imagen a color que representa el AEE entre el ground truth y la mejor estimaci´on del m´etodo multicanal. Como se puede observar la soluci´on obtenida es bastante buena excepto en la costa africana y en los bordes de la imagen. M´etodo/Canal AEE AAE SF, VIS 0.8 0.2608 5.2055 SF, WV 6.2 0.5280 12.3891 SF, WV 7.3 0.4992 11.9396 SF, IR 10.8 0.3413 7.8710 MC Avg. 0.1926 3.8432 MC Max. 0.1895 3.7696 Cuadro 2.1: Secuencia de Vince: AEE y AAE obtenidos por las distintas versiones del m´etodo monocanal y multicanal. 2.3. M ´ ETODO VARIACIONAL MULTICANAL 71 Figura 2.3: Secuencia 1. Arriba: el campo desplazamiento de los canales VIS 0.8 y WV 6.2. Abajo: el campo desplazamiento de los canales WV 7.3 y IR 10.8. Figura 2.4: Secuencia 1. Campo de desplazamiento obtenido con el m´etodo multicanal (estrategia gradiente M´aximo). 72 CAP´ ITULO 2. ESTIMACI ´ ON DEL FLUJO ´ OPTICO EN SECUENCIAS DE IM ´ AGENES Figura 2.5: El Error Eucl´ıdeo Medio (AEE) entre el ground truth y la mejor estimaci´on del m´etodo multicanal. Figura 2.6: Secuencia del Atl´antico Norte. De izquierda a derecha y de arriba hacia abajo, se muestra el canal visible 0,81µm, los canales de vapor de agua 6,25µm y 7,35µm y a la derecha el canal de infrarrojo, 10,80µm. 2.3. M ´ ETODO VARIACIONAL MULTICANAL 73 Secuencia 2. Atl´antico Norte (5 Junio 2004) La segunda secuencia utilizada en los experimentos corresponde a unos datos del Meteosat del d´ıa 5 de Junio del 2004. Se trata de una vista de la zona afroeuropea del hemisferio norte. En la figura 2.6 podemos observar las im´agenes de cuatro canales del sat´elite. Dado el tama˜no de las im´agenes se ha seleccionado una regi´on de 559×575 p´ıxeles en donde, a nuestro entender, se registran los fen´omenos atmosf´ericos m´as importantes (figura 2.8). En esa regi´on se localizan dos grandes masas nubosas. Una de ellas es una borrasca que se aproxima del Atl´antico hacia las islas brit´anicas por la parte norte. La segunda masa nubosa se sit´ua cerca de la pen´ınsula Ib´erica. El movimiento predominante en ambas masas nubosas es rotacional. En la figura 2.7 se muestra la estimaci´on del desplazamiento de las nubes en cada uno de los canales del sat´elite. Estos resultados son muy parecidos en magnitud. Sin embargo, en el canal visible se detecta con m´as claridad las vorticidades presentes en la imagen; mientras que en el resto de canales ese fen´omeno se detecta como casi traslacional. Por ello, la contribuci´on de los datos de este canal debe ser mayor. En la figura 2.8 se representa el campo de desplazamiento obtenido con el m´etodo multicanal utilizando los pesos con los valores comentados anteriormente. Estos pesos otorgan un valor predominante al canal visible. En las figuras 2.7 y 2.8 se han mostrado los distintos campos de desplazamiento calculados con el m´etodo variacional monocanal y multicanal. Visualmente no es posible percibir las peque˜nas diferencias entre los resultados. Por este motivo, en la tabla 2.2 se recogen unos datos cuantitativos que permiten comparar con m´as exactitud los distintos resultados. La mejora de la estimaci´on hecha por el m´etodo multicanal frente a las del monocanal es de al menos 20.77 % para el AEE y de 23.47 % para el AAE. Al igual que ocurr´ıa en la secuencia anterior las dos estrategias utilizadas en el m´etodo multicanal ofrecen resultados parecidos. En la figura 2.9 se muestra una imagen a color que representa el AEE entre el ground truth y la mejor estimaci´on del m´etodo multicanal. M´etodo/Canal AEE AAE SF, VIS 0.8 0.1704 4.4748 SF, WV 6.2 0.5813 15.1845 SF, WV 7.3 0.4776 12.6222 SF, IR 10.8 0.3064 7.7821 MC Avg. 0.1593 4.1831 MC Max. 0.1350 3.4246 Cuadro 2.2: Secuencia del Atl´antico Norte: AEE y AAE obtenidos por las distintas versiones del m´etodo monocanal y multicanal. 74 CAP´ ITULO 2. ESTIMACI ´ ON DEL FLUJO ´ OPTICO EN SECUENCIAS DE IM ´ AGENES Figura 2.7: Secuencia 2. Arriba: el campo desplazamiento de los canales VIS 0.8 y WV 6.2. Abajo: el campo desplazamiento de los canales WV 7.3 y IR 10.8. 2.3. M ´ ETODO VARIACIONAL MULTICANAL 75 Figura 2.8: Secuencia 2. Campo de desplazamiento obtenido con el m´etodo multicanal (estrategia gradiente M´aximo). Figura 2.9: El Error Eucl´ıdeo Medio (AEE) entre el ground truth y la mejor estimaci´on del m´etodo multicanal. 76 CAP´ ITULO 2. ESTIMACI ´ ON DEL FLUJO ´ OPTICO EN SECUENCIAS DE IM ´ AGENES 2.4. M´etodo Variacional que incorpora un T´ermino de Regularizaci´on exclusivamente Temporal no Continuo En esta secci´on se describe un m´etodo para la estimaci´on del flujo ´optico en una secuencia de im´agenes. Se trata de una modificaci´on del modelo de energ´ıa propuesto por Nagel-Enkelmann [Nagel86] al que se le ha a˜nadido un t´ermino de regularizaci´on temporal. En la ecuaci´on (2.31) podemos observar el modelo de energ´ıa propuesto en [Papenberg06]. En ese trabajo se analiz´o el impacto de la combinaci´on de varias invarianzas en el t´ermino de ligadura. El t´ermino de ligadura utilizado se compone de dos invarianzas: gradiente constante y la suposici´on lambertiana, y el de suavizado se trata de una versi´on espaciotemporal robustificada del propuesto por Horn-Schunck. E(h) = ZΩ Ψ(I1(x)−I2(x+h))2+γ(∇I1(x)−∇I2(x+h))2dx +αZΩ Ψk∇3uk2+k∇3vk2dx.(2.31) Si tomamos como referencia el trabajo de [Papenberg06], su modelo de energ´ıa incluye un t´ermino de regularizaci´on basado en el gradiente espaciotemporal del flujo. Este t´ermino impone una restricci´on de continuidad en el flujo en todas direcciones, es decir, que el desplazamiento sea suave tanto en la direcci´on espacial como temporal. Este m´etodo ha demostrado ser uno de los m´as precisos que existen. En la evaluaci´on de los m´etodos se utilizan secuencias sint´eticas que omiten parte de la complejidad del mundo real. En situaciones reales la restricci´on de Papenberg et al. no es un modelo satisfactorio para la estimaci´on de ciertos tipos de flujos. Pese a que el tratamiento de forma conjunta de la informaci´on espacial y temporal ofrece una serie de ventajas, obliga a una continuidad en el dominio temporal que no es siempre posible. Por ello, como parte de esta tesis se ha desarrollado un m´etodo que a˜nade un t´ermino de regularizaci´on temporal independiente del espacial. De esta forma, se desacopla la informaci´on del flujo espacial del temporal evitando la limitaci´on anteriormente comentada. Antes de entrar en detalle en la descripci´on de nuestro modelo de energ´ıa conviene repasar alguno de los t´erminos de regularizaci´on m´as conocidos en la literatura. Uno de los primeros m´etodos variacionales que se propusieron para la estimaci´on del flujo ´optico fue el de Horn–Schunck [Horn81]. El funcional de energ´ıa se defin´ıa mediante dos t´erminos: (1) la ecuaci´on de restricci´on del flujo ´optico, ec. (1.2), y (2) un t´ermino que consiste en la norma del gradiente del flujo ´optico. k∇hk2=k∇uk2+k∇vk2,(2.32) donde k.kes el la norma de un vector y ∇es el operador gradiente. Uno de los inconvenientes que tiene este t´ermino de regularizaci´on es que no preserva las discontinuidades en el flujo. El modelo de energ´ıa propuesto por Nagel-Enkelmann [Nagel86] supone otra de las grandes aportaciones en el campo del flujo ´optico. Este modelo, similar al definido por Horn-Schunck, incluye un t´ermino de suavizado anisotr´opico. Este t´ermino var´ıa la cantidad de difusi´on aplicada en funci´on del gradiente de la imagen. En las regiones donde 2.4. M ´ ETODO VARIACIONAL CON REGULARIZACI ´ ON TEMPORAL NO CONTINUA 77 el gradiente es peque˜no act´ua de forma isotr´opica y en aqu´ellas donde el gradiente es alto de forma anisotr´opica, suavizando a lo largo de los contornos y no a trav´es de ellos. Inicialmente, el modelo de Nagel-Enkelmann [Nagel86] se dise˜n´o para dos frames, ec. (2.33); pero luego se cre´o una versi´on espaciotemporal del mismo, [Nagel90]. El m´etodo descrito en esta secci´on se basa en la versi´on espacial. E(h) = ZΩ (I1(˜x)−I2(˜x +h))2d˜x +αZΩ tr(∇h>D∇h)d˜x,(2.33) donde ˜x es un punto en el dominio Ω, (.)>se trata del operador de trasposici´on, I1y I2son las dos im´agenes de entrada, tr(.) es el operador conocido como trace,αes un peso constante que pondera al t´ermino de suavizado y Des una matriz de proyecci´on regularizada en la direcci´on ortogonal al ∇I1. La matriz Dse define como D(∇I1) = 1 k∇I1k2+ 2ζ2ξξ>+ζ2Id,(2.34) donde ξ= (∂I1 ∂y ,−∂I1 ∂x )>es el vector ortogonal a ∇I1,Id es la matriz identidad, y ζes un coeficiente que determina el comportamiento isotr´opico aplicado en el suavizado e inhibe la difusi´on a trav´es de los bordes: k∇I1k  ζ. La estructura de esta secci´on se divide de la siguiente forma. En la secci´on 2.4.1 se introducir´a el modelo de energ´ıa y se justificar´a la divisi´on de la regularizaci´on en dos t´erminos distintos, uno espacial y otro temporal. En la secci´on 2.4.2, se describe con detalle el esquema num´erico utilizado a partir de la minimizaci´on del funcional de energ´ıa. As´ı mismo, en la secci´on 2.4.5, se mostrar´a una serie de experimentos que justificar´an la inclusi´on del t´ermino de regularizaci´on temporal. Para ello, se ha evaluado nuestro m´etodo utilizando secuencias tanto reales como sint´eticas. 2.4.1. Modelo de Energ´ıa Los primeros m´etodos variacionales propuestos para la estimaci´on del flujo ´optico inclu´ıan t´erminos de regularizaci´on espacial. Estos t´erminos suavizan la soluci´on utilizando la informaci´on procedente de dos im´agenes. Cuando se dispone de una secuencia de im´agenes ´estos no hacen uso de la informaci´on de los frames vecinos. Para aumentar la precisi´on de las estimaciones y aprovechar la transferencia de informaci´on de los frames consecutivos se han dise˜nado nuevos t´erminos de regularizaci´on espaciotemporales. Muchos de ellos, s´olo son una extensi´on temporal de su hom´ologo espacial ([Weickert01b, Nagel90]). En el trabajo de Weickert-Schn¨orr [Weickert01b] se propone un m´etodo que no es m´as que una extensi´on espaciotemporal del m´etodo de Horn-Schunck. El t´ermino de regularizaci´on es de la forma R(|∇3u|2+|∇3v|2), donde Rse trata de una funci´on de robustificaci´on que penaliza en menor grado las desviaciones en el flujo. Este t´ermino trata de forma conjunta las derivadas espaciales y temporal. Por un lado, la formulaci´on del modelo queda expresada de una forma homog´enea. Por otro lado, esta formulaci´on continua est´a condicionada a que el desplazamiento de los p´ıxeles sea peque˜no. 78 CAP´ ITULO 2. ESTIMACI ´ ON DEL FLUJO ´ OPTICO EN SECUENCIAS DE IM ´ AGENES A partir del m´etodo [Weickert01b] se cre´o una extensi´on para desplazamientos largos, [Weickert04]. El nuevo modelo de energ´ıa propuesto inclu´ıa como t´ermino de ligadura, la ecuaci´on de restricci´on de flujo ´optico para largos desplazamientos, ec. (2.10) y, para el de suavizado se utiliz´o –R(|∇3u|2+|∇3v|2)–, el mismo que en [Weickert01b]. Este funcional de energ´ıa tiene dos limitaciones. En primer lugar, cuando se producen desplazamientos largos la regularizaci´on temporal puede influir negativamente en la regularizaci´on espacial debido al acoplamiento que existe entre las derivadas espaciotemporales. En segundo lugar, se da una ligera incongruencia entre el t´ermino de ligadura y el de suavizado. Uno permite discontinuidades en el tiempo y, el otro, exige que el flujo sea continuo en el tiempo. Una forma de evitar estas dos limitaciones consiste en reemplazar el t´ermino R(|∇3u|2+ |∇3v|2) por R(|∇u|2+|∇v|2)+T(|ut|2+|vt|2) en el que las derivadas espaciales y temporales est´en separadas. Todav´ıa se sigue teniendo el segundo inconveniente, la continuidad del flujo y la imposibilidad de manejar largos desplazamientos. Esto se soluciona sustituyendo la derivada temporal por estimaci´on que tenga en cuenta los desplazamientos largos. El t´ermino T(.) debe permitir la detecci´on de los desplazamientos largos en la direcci´on temporal por lo que la derivada temporal del flujo no es v´alida. Si se asume que la velocidad del objeto es suave en toda la secuencia, los flujos ´opticos entre frames vecinos deben ser parecidos. Esta semejanza se puede expresar como hi(˜x) y hi+1(˜x +hi(˜x)) para el flujo en el frame iy el consecutivo, respectivamente. Por ello, el t´ermino T(.) estar´ıa expresado en funci´on de T(hi(˜x),hi+1(˜x +hi(˜x))). De acuerdo a la generalizaci´on descrita en la secci´on 2.2, nuestro modelo de energ´ıa se ajusta a un caso particular de la ecuaci´on (2.15). Se trata de un m´etodo espaciotemporal que utiliza la informaci´on de toda la secuencia. Se compone de dos partes: un t´ermino de ligadura muy simple que s´olo incluye una sola invarianza y un t´ermino de regularizaci´on algo m´as complejo que suaviza por separado el flujo espacial y temporal. El modelo de energ´ıa propuesto tiene la estructura siguiente E(h) = N−1 X i=1 Z Ω D(Ii, Ii+1,hi)d˜x +α N−1 X i=1 Z Ω R(∇Ii,∇hi)d˜x +β N−2 X i=1 Z Ω T(hi,hi+1)d˜x.(2.35) Este funcional se compone de tres t´erminos: (i) la ecuaci´on de restricci´on del flujo ´optico para largos desplazamientos (ec. (2.11) con una sola invarianza, la suposici´on lambertiana), (ii) el operador de Nagel-Enkelmann y (iii) un t´ermino temporal T(hi,hi+1) = 2.4. M ´ ETODO VARIACIONAL CON REGULARIZACI ´ ON TEMPORAL NO CONTINUA 79 Φ(khi−hi+1(˜x +hi)k2). E(h) = N−1 X i=1 Z Ω (Ii−Ii+1 (˜x +hi))2d˜x +α N−1 X i=1 Z Ω trace(∇hT iD(∇Ii)∇hT i)d˜x +β N−2 X i=1 Z Ω Φ(khi−hi+1(˜x +hi)k2)d˜x,(2.36) donde Φ(x2) = 1 −γe−x2 γ. Cuando el desplazamiento es peque˜no, hi−hi+1(˜x +hi) es una aproximaci´on a la derivada temporal. Un inconveniente que tiene este t´ermino es que la transferencia de informaci´on se hace desde los ´ultimos frames a los primeros. Existe una fuerte dependencia respecto al ´ultimo flujo de la secuencia; si ´este es malo, esa estimaci´on se propagar´a negativamente al resto de la secuencia. Una forma de compensar este problema consiste en la inclusi´on de otro t´ermino, tambi´en temporal, que favorezca la transferencia inversa desde los primeros frames a los ´ultimos, ec. (2.37). De esta forma, nuestro modelo de energ´ıa incluye informaci´on temporal del flujo sin premiar/penalizar a ninguno en especial. La transferencia inversa del flujo (h∗) desde los primeros frames a los ´ultimos se explica con detalle en el apartado 2.4.4. E(h) = N−1 X i=1 Z Ω (Ii−Ii+1 (˜x +hi))2d˜x +α N−1 X i=1 Z Ω trace(∇hT iD(∇Ii)∇hT i)d˜x +β N−2 X i=1 Z Ω Φ(khi−hi+1(˜x +hi)k2)d˜x +β N−1 X i=2 Z Ω Φ( hi−hi−1(˜x +h∗ i−1) 2)d˜x (2.37) Los dos ´ultimos t´erminos de la ecuaci´on (2.37) asumen un modelo en la velocidad de los objetos. Esto obliga a que las funciones hisean similares en direcci´on y magnitud. Por ello, este modelo de energ´ıa es m´as adecuado cuando los desplazamientos de los objetos por la escena son constantes y en la misma direcci´on. Cuando el desplazamiento no se ajuste a esta suposici´on, Φ(.) decrecer´a por el efecto de T(.) perdiendo importancia la regularizaci´on temporal convirtiendo al modelo en puramente espacial. Definir un modelo gen´erico para T(.) que tenga en cuenta todas las posibles situaciones que pueden darse en secuencias reales no es tarea sencilla. Otras alternativas para el t´ermino temporal que podr´ıan favorecer un movimiento constante independientemente de la direcci´on podr´ıa ser Φ(khik2− khi+1(˜x +hi)k2) + Φ(khik2− hi−1(˜x +h∗ i−1) 2). La 86 CAP´ ITULO 2. ESTIMACI ´ ON DEL FLUJO ´ OPTICO EN SECUENCIAS DE IM ´ AGENES Figura 2.13: AEE de la secuencia del Cuadrado. Las diferencias entre las estimaciones del m´etodo espacial y espaciotemporal son m´ınimas aunque si se observa la estabilidad que aporta la informaci´on temporal. Los cambios en la magnitud del error entre los distintos frames es menos abrupta que en el caso espacial. Cuadro 2.3: Lista de par´ametros utilizados por los m´etodos Spatial,Temporal y Bi–Temporal en la secuencia del cuadrado. Escalas = n´umero de escalas en el enfoque multipiramidal. C= par´ametro de regularizaci´on espacial. β= par´ametro de regularizaci´on temporal. γ= par´ametro de control de la funci´on Φ del t´ermino de regularizaci´on temporal. λ= factor de isotrop´ıa del operador de Nagel–Enkelmann. Cuadrado M´etodo Escalas C β γ λ Spatial 4 0.6 - - 0.1 Temporal 4 0.6 0.01 0.1 0.3 Bi–Temporal 4 0.6 0.01 0.5 0.1 soluciones si que se puede percibir la aportaci´on del t´ermino de regularizaci´on temporal que ha permitido controlar el proceso de difusi´on respecto a la estimaci´on espacial. El desplazamiento del fondo de la escena se aprecia m´as homog´eneo en la soluci´on temporal que en la espacial. En la tabla 2.3 mostramos los par´ametros utilizados en los experimentos. Dado que el desplazamiento del objeto es de 15 p´ıxeles es necesario emplear al menos cuatro escalas en el enfoque multipiramidal. A la hora de estimar correctamente el flujo ´optico en la escena y evitar que el proceso de difusi´on act´ue m´as all´a de las discontinuidades del cuadrado debemos utilizar un valor de λbajo, anisotr´opico, mientras que Ctendr´ıa un valor intermedio. Como se ha visto, el t´ermino de regularizaci´on temporal permite mejorar la estimaci´on hecha por el espacial. Cuando el movimiento en la secuencia es traslacional y constante se manifiesta la eficacia de este nuevo t´ermino de regularizaci´on temporal. En la siguiente secuencia veremos c´omo se comporta el m´etodo ante las aceleraciones de los objetos en 2.4. M ´ ETODO VARIACIONAL CON REGULARIZACI ´ ON TEMPORAL NO CONTINUA 87 Figura 2.14: De arriba a abajo los campos de desplazamiento correspondiente a los frames 0, 4 y 8 de la secuencia del Cuadrado. En la primera columna los desplazamientos reales. En las siguientes columnas las estimaciones obtenidas con los m´etodos Spatial yBi–Temporal. Las estimaciones obtenidas por los distintos m´etodos son muy parecidas. Si comparamos las soluciones temporales con la espacial observamos que la aportaci´on del t´ermino de regularizaci´on temporal ha permitido controlar o acotar el proceso de difusi´on. En la estimaci´on espacial se aprecia una regularizaci´on algo descontrolada. determinados frames. Marble Blocks Marble Blocks se trata de una secuencia donde la c´amara se mueve horizontalmente mientras una serie de torres de m´armol permanecen est´aticas (figura 2.16). La secuencia se compone de treinta frames. En nuestros resultados experimentales s´olo se ha mostrado la estimaci´on de quince de ellos (los situados en la parte central, del frame 5 al 20). El movimiento aparente de la escena indica que todas la torres se mueven a distinta velocidad, en funci´on de la distancia a la c´amara. En principio, el movimiento registrado en esta secuencia se asemeja mucho a la anterior. Sin embargo, la velocidad de los objetos en cada uno de los frames no es constante como se pudiera pensar. En los frames 4, 8, 13, 18, 23 y 28 existen unas aceleraciones de los objetos respecto a los frames contiguos. Estos cambios en la velocidad se pueden percibir sutilmente a simple vista. Para verificar esta hip´otesis hemos calculado la velocidad media entre frames, utilizando el flujo real suministrado por Nagel, y hemos observado que este valor se incrementa notablemente en los frames especificados anteriormente. La velocidad media en los distintos frames de la secuencia es de 1,33 con una desviaci´on t´ıpica de 0,03. En los frames indicados previamente la velocidad media se incrementa hasta los 1,58 con una desviaci´on t´ıpica de 0,02. En la secuencia de Marble Blocks las superficies de los objetos est´an cubiertas por texturas. Los regularizadores del tipo image-driven, como es el caso del operador de Nagel– Enkelmann, no se comportan muy bien ya que podr´ıan dar lugar a la creaci´on de flujos 88 CAP´ ITULO 2. ESTIMACI ´ ON DEL FLUJO ´ OPTICO EN SECUENCIAS DE IM ´ AGENES Secuencia del Cuadrado M´etodo AAEµAAEσAEEµAEEσ Spatial 6,2116o0,3455o1,3393 0,0830 Temporal 5,8100o0,1347o1,3819 0,0902 Bi-Temporal 5,7423o0,1612o1,3428 0,0612 Cuadro 2.4: AAE y AEE en la secuencia del cuadrado. Los sub´ındices µyσdenotan la media y la desviaci´on t´ıpica, respectivamente. Seg´un vemos en los resultados cuantitativos el m´etodo Bi–Temporal mejora las estimaciones respecto a los otros dos m´etodos. La mejora del temporal unidireccional respecto al espacial no es tan evidente debido a la dependencia que existe con la estimaci´on del ´ultimo frame; aunque como podemos ver en la tabla, la desviaci´on t´ıpica es muy baja lo que refleja la estabilidad ofrecida por la informaci´on temporal. Figura 2.15: AAE de la secuencia del Cuadrado. Las estimaciones obtenidas por los m´etodos temporales son ligeramente m´as precisas. Nuevamente, la estimaci´on del m´etodo espacial se vuelve inestable debido a la ausencia de la informaci´on temporal. hipersegmentados al confundir las fuertes variaciones en la magnitud del gradiente en la textura con discontinuidades en la imagen. La presencia de las aceleraciones y las texturas en los objetos no supone el escenario ideal para nuestro m´etodo. A pesar de estos dos inconvenientes se ver´a que los resultados experimentales son relativamente buenos. En la figura 2.17 podemos ver una comparaci´on visual entre el ground truth y los flujos obtenidos por los m´etodos Spatial yBi–Temporal. Si le echamos un vistazo con m´as detenimiento observamos que las estimaciones obtenidas por la versi´on Spatial yBi– Temporal del m´etodo son bastante parecidas aunque existen peque˜nas diferencias. En la estimaci´on temporal el suelo se percibe m´as suavizado y el fondo de la escena est´a mucho m´as definido. Otras diferencias a destacar est´an en la torre situada en la parte derecha de la imagen que parece ser m´as continua, sin grandes artificios en el flujo. Por otro lado, el flujo de las dos torres situadas en la parte izquierda alejadas de la c´amara parecen estar subestimadas en su parte central. 2.4. M ´ ETODO VARIACIONAL CON REGULARIZACI ´ ON TEMPORAL NO CONTINUA 89 Figura 2.16: Frames 5, 10, 15, y 20 de la secuencia Marble Blocks. Marble Blocks M´etodo AAEµAAEσAEEµAEEσ Spatial 6,695o2,698o0,2480 0,0963 Temporal 5,402o1,327o0,2081 0,0638 Bi–Temporal 4,731o1,330o0,1848 0,0661 Cuadro 2.5: AAE y AEE para la secuencia Marble Blocks. En esta secuencia se aprecia la ventaja que ofrece la informaci´on espaciotemporal respecto a los m´etodos espaciales. La mejora del temporal unidireccional es notable y en una proporci´on similar al Bi–Temporal. En las figuras 2.18 y 2.19 disponemos de dos gr´aficas que nos muestran los AEE y AAE de cada frame de la secuencia. Al igual que ocurr´ıa con la secuencia de los cuadrados la estimaci´on hecha por el m´etodo puramente Spatial es m´as inestable y menos precisa. Las versiones temporales del m´etodo mejoran sustancialmente respecto al espacial. El comportamiento del Bi–Temporal es similar al Temporal pero con un error en magnitud relativamente inferior. Si observamos las figuras 2.18 y 2.19 existen unos picos en los valores de AEE y AAE. Estos frames coinciden con aqu´ellos en los que hab´ıa aceleraciones de los objetos. Por este motivo, la estabilidad de las estimaciones se ve ligeramente comprometida, es decir, que el t´ermino de regularizaci´on temporal no identifica estas aceleraciones en la imagen y la estimaci´on no es tan buena como en los frames vecinos donde si existe cierta continuidad. En la tabla 2.5 se muestra los AAE y AEE globales junto con las desviaciones t´ıpicas 90 CAP´ ITULO 2. ESTIMACI ´ ON DEL FLUJO ´ OPTICO EN SECUENCIAS DE IM ´ AGENES Figura 2.17: De arriba a abajo, el campo de desplazamiento correspondiente a los frames 5, 10, 15 y 20 de la secuencia del Marble-Blocks. En la primera columna los desplazamientos reales. En la columna central las estimaciones obtenidas con el m´etodo Spatial y en la columna de la derecha, las estimaciones del m´etodo Bi–Temporal. Existe bastante parecido entre las estimaciones. Sin embargo, si nos fijamos con m´as detenimiento la estimaci´on del m´etodo Bi–Temporal es mucho m´as suave en las zonas homog´eneas y no se aprecian grandes artificios. En el caso de la estimaci´on espacial, la presencia de artificios es manifiesta y los bordes de los objetos est´an borrosos y menos definidos. 2.4. M ´ ETODO VARIACIONAL CON REGULARIZACI ´ ON TEMPORAL NO CONTINUA 91 Figura 2.18: AEE para la secuencia de Marble Blocks. En esta gr´afica s´olo se muestra los errores desde el frame 5 hasta el 20. Se observa el efecto de las distintas aceleraciones detectadas en la secuencia. En los frames anteriormente se˜nalados (8, 13 y 18) hay unos cambios notables de velocidad. El m´etodo temporal, en sus dos variantes, penaliza los cambios bruscos entre los flujos ´opticos por lo que esos frames el error en la estimaci´on ser´a mayor. La ´ultima estimaci´on del m´etodo temporal unidireccional debe coincidir con la del espacial. Sin embargo, los datos mostrados en la gr´afica s´olo recogen los errores del frame 5 al 20. Por ´ultimo, cabe decir que se han seleccionado estos frames porque, por un lado, reflejan la estabilidad de los m´etodos temporales y, por otro lado, muestra el aumento del error en aquellos frames donde se producen las aceleraciones en los objetos. Figura 2.19: AAE para la secuencia de Marble Blocks. En esta gr´afica s´olo se muestra los errores desde el frame 5 hasta el 20. Al igual que se observaba en la gr´afica 2.18 el AAE se incrementa en determinados frames que coinciden con los que se produce ciertas aceleraciones detectadas. El comportamiento de los m´etodos es similar al comentado en la gr´afica 2.18. 92 CAP´ ITULO 2. ESTIMACI ´ ON DEL FLUJO ´ OPTICO EN SECUENCIAS DE IM ´ AGENES Cuadro 2.6: Lista de par´ametros utilizados por los m´etodos Spatial,Temporal yBi– Temporal en la secuencia del Marble Blocks. Escalas = n´umero de escalas en el enfoque multipiramidal. C= par´ametro de regularizaci´on espacial. β= par´ametro de regularizaci´on temporal. γ= par´ametro de control de la funci´on Φ del t´ermino de regularizaci´on temporal. λ= factor de isotrop´ıa del operador de Nagel–Enkelmann. Marble Blocks M´etodo Escalas C β γ λ Spatial 2 0.6 - - 0.5 Temporal 2 0.6 0.3 0.5 0.5 Bi–Temporal 2 0.6 0.3 0.5 0.5 asociadas. Seg´un los datos de esta tabla el m´etodo Temporal obtiene respecto al Spatial en torno al 19,31 % de mejora del AAE y un 16,09 % en el caso de AEE. Si comparamos la mejora del Bi–Temporal respecto al Spatial, ´esta es del 29,34 % y 25,48 % para el AAE y AEE, respectivamente. Tambi´en cabe destacar la reducci´on de la desviaci´on t´ıpica que se produce en las versiones temporales. En la tabla 2.6 mostramos los valores de los par´ametros empleados en los experimentos. En general, la inclusi´on del t´ermino temporal influye positivamente en la precisi´on de las estimaciones hechas por el m´etodo. Hay que tener en cuenta que esta secuencia no era, a priori, la m´as id´onea debido a los inconvenientes comentados anteriormente. Taxi de Hamburgo La ´ultima secuencia utilizada en los experimentos es la conocida como Taxi de Hamburgo (fig. 2.20). En el apartado 1.4.1 del cap´ıtulo anterior podemos encontrar m´as informaci´on acerca de ella. Se trata de una secuencia real en la que se observan varios objetos en movimiento: tres veh´ıculos y una persona; el resto de objetos de la secuencia permanecen est´aticos. Al tratarse de una secuencia real no disponemos de ground truth por lo que solamente se mostrar´an resultados visuales. En la figura 2.21, se muestran las soluciones de los m´etodos Spatial yBi–Temporal. Figura 2.20: Frames 0, 10 y 19 de la secuencia del Taxi de Hamburgo. 2.4. M ´ ETODO VARIACIONAL CON REGULARIZACI ´ ON TEMPORAL NO CONTINUA 93 Figura 2.21: A la izquierda, la estimaci´on obtenida por el m´etodo espacial. A la derecha, la soluci´on ofrecida por el m´etodo temporal bidireccional. Los valores de los par´ametros son: N´umero de escalas = 2, C= 0,6, β= 0,1, γ= 0,5, λ= 0,1. Los dos veh´ıculos situados en la parte inferior de la secuencia se mueven a velocidades constantes en la misma direcci´on pero en sentido opuesto. El tercer veh´ıculo, un taxi, realiza un peque˜no giro a baja velocidad y un transe´unte camina por la acera. Este fen´omeno es casi imperceptible a simple vista. Nuestro modelo de energ´ıa no incorpora ninguna t´ecnica de robustificaci´on lo que le hace altamente sensible al ruido. En el caso espacial esta sensibilidad se hace m´as notable. Para una correcta identificaci´on del movimiento de los tres veh´ıculos se ha optado por un regularizador espacial con un comportamiento anisotr´opico (λ= 0,1) y un peso para ambos t´erminos de suavizado de C= 0,6 y β= 0,1 primando el espacial sobre el temporal. Como vemos en la figura 2.21 el t´ermino de regularizaci´on temporal suaviza la soluci´on espacial eliminando en gran medida los falsos movimientos detectados en esta soluci´on. Teniendo en cuenta todo esto, existen importantes diferencias entre las estimaciones de ambos m´etodos y que comentamos a continuaci´on: 1. Fondo de la escena. Salvo los cuatro objetos mencionados anteriormente en la escena no se registra ning´un tipo de movimiento. El flujo ´optico en esos p´ıxeles deber´ıa ser nulo. La soluci´on temporal es bastante suave aunque se observan algunos artificios producto de los efectos del ruido local. Por el contrario, la soluci´on espacial dista mucho de la realidad, el n´umero de movimientos fantasmas es elevado y el flujo no es homog´eneo. 2. Peat´on. El desplazamiento del peat´on por la acera es un movimiento casi imperceptible a simple vista. El m´etodo temporal es capaz de detectarlo mientras que este movimiento al espacial se le solapa con los falsos positivos registrados en el fondo de la escena. 3. Movimiento de los veh´ıculos. Los desplazamientos del taxi y del cami´on son similares en ambos m´etodos. El movimiento del coche situado a la izquierda parece estar mejor estimado en el caso espacial, pero en general, la estimaci´on obtenida con el m´etodo temporal es m´as estable y precisa. 94 CAP´ ITULO 2. ESTIMACI ´ ON DEL FLUJO ´ OPTICO EN SECUENCIAS DE IM ´ AGENES 2.5. M´etodos Variacionales basados en el An´alisis Espectral del Tensor de Movimiento En los ´ultimos a˜nos los m´etodos variacionales presentados han incrementado a´un m´as la precisi´on y robustez de sus estimaciones. Esta mejora se ha debido principalmente gracias a la incorporaci´on de nuevos elementos al modelo de energ´ıa, a la utilizaci´on de nuevos esquemas num´ericos y al uso de discretizaciones m´as exactas. Los modelos de energ´ıa normalmente se componen de dos t´erminos. El t´ermino de ligadura atrae a las estructuras en movimientos mientras que el de regularizaci´on act´ua en aquellas zonas donde el t´ermino de ligadura no dispone de informaci´on suficiente para realizar una buena estimaci´on. El funcionamiento id´oneo de un m´etodo consistir´ıa en la detecci´on del desplazamiento en los bordes de los objetos (regiones donde se aprecia claramente el movimiento) y, a partir de esa informaci´on, rellenar el flujo para el resto de p´ıxeles. El proceso de rellenado de los objetos s´olo actuar´ıa en el interior de los objetos y se detendr´ıa exactamente en las discontinuidades, no se propagar´ıa m´as all´a de los bordes de los objetos. A pesar de que en la literatura se han propuesto todo tipo de regularizadores que permiten guiar el proceso de difusi´on, sigue presente la sensaci´on que los t´erminos del modelo colaboran pero no se complementan del todo. El mayor problema radica en la identificaci´on y elecci´on de la direcci´on o direcciones a suavizar. En esta secci´on se presenta un m´etodo para la estimaci´on del flujo ´optico usando secuencias de im´agenes. El modelo de energ´ıa se apoya en t´ecnicas recientes, como es el tensor de movimiento. A trav´es del an´alisis espectral del tensor de movimiento es posible definir un modelo de energ´ıa que sea capaz de intercambiar informaci´on complementaria entre el t´ermino de ligadura y el de regularizaci´on. De esta forma, se establece un mecanismo directo que suaviza en las direcciones no dominantes del flujo. Para mejorar la precisi´on de las estimaciones se utiliza t´erminos no-cuadr´aticos que han demostrado ser robustos frente a los outliers. Nuestro modelo es capaz de detectar los desplazamientos largos gracias al uso de t´ecnicas de warping. Esta secci´on se divide en cinco apartados. En primer lugar, se describe la notaci´on y las t´ecnicas utilizadas para la definici´on de nuestro modelo de energ´ıa gen´erico. Veremos como gracias a estas t´ecnicas podemos establecer una relaci´on de complementariedad entre ambos t´erminos de la energ´ıa. En el apartado 2.5.2, comenzamos presentando el modelo de energ´ıa de Nagel–Enkelmann. La capacidad de este modelo para direccionar el proceso de difusi´on encierra la idea base de nuestro m´etodo. La combinaci´on de las t´ecnicas anteriores y del modelo de Nagel–Enkelmann dar´a como resultado la estructura general de nuestro framework. Una vez descrito las distintas variantes del modelo general, en los apartados 2.5.3 y 2.5.4, se minimizar´an y se comentar´a el esquema num´erico utilizado. En el apartado 2.5.5, se muestran los resultados experimentales obtenidos con secuencias reales y sint´eticas. Comparamos nuestros resultados con los mejores m´etodos de la literatura demostrando que las ventajas de este nuevo framework no son s´olo apreciables desde el punto de vista te´orico. 2.5. M ´ ETODOS VARIACIONALES BASADOS EN EL AN ´ ALISIS ESPECTRAL 95 2.5.1. Tensor de Movimiento Para detectar las correspondencias de los p´ıxeles entre dos im´agenes se suele suponer que alguna propiedad de la imagen no var´ıa a lo largo del tiempo. Esta suposici´on se puede representar tal que f(x+u, y +v, t + 1) −f(x, y, t)=0,(2.44) donde fes alg´un tipo de propiedad en la imagen y, t+ 1 y t, representan dos im´agenes de la escena en distintos instantes de tiempo. Esta no linealidad se suele solventar a trav´es del desarrollo de Taylor de la ecuaci´on (2.44) y eliminado los t´erminos de orden superior con lo que obtendr´ıamos fxu+fyv+ft= 0, donde los sub´ındices indican derivadas parciales. Para estimar el flujo ´optico, h(x) := (u(x), v(x),1)>, normalmente se recurre a la minimizaci´on por m´ınimos cuadrados de la expresi´on anterior (fxu+fyv+ft)2=h>∇f∇f>h=h>Jh.(2.45) La expresi´on (2.45) se puede descomponer de tal forma que obtengamos una matriz, que por definici´on es, semidefinida positiva, [Bruhn06a]. Esta matriz se conoce como tensor de movimiento,J:= ∇f∇f>. Entre las numerosas ventajas que ofrece destaca: (i) la posibilidad de representar cualquier invarianza en una estructura compacta y, (ii) la inclusi´on m´ultiples suposiciones dentro de esa misma estructura. En tal situaci´on obtendr´ıamos un tensor de movimiento acumulado Jmc = Nc X i=1 γi∇fi∇f> i, donde Nces el n´umero de invarianzas y γiun peso asociado a cada invarianza. A continuaci´on, se muestra un ejemplo de dos invarianzas, suposici´on lambertiana y gradiente constante, en el que se puede comprobar c´omo quedan expresadas varias invarianzas mediante un tensor de movimiento acumulado. ((f(x+u, y +v, t + 1) −f(x, y, t))2+γ|∇f(x+u, y +v, t + 1) −∇f(x, y, t)|2≈ (fxu+fyv+ft)2+γ(fxxu+fxyv+fxt)2 +(fxyu+fyyv+fyt)2= h>∇f∇f>h+γh>(∇fx∇f> x+∇fy∇f> y)h= h>Jfh+γh>J∇fh= h>Jmc h(2.46) donde Jmc =∇f∇f>+γ(∇fx∇f> x+∇fy∇f> y). Una caracter´ıstica de las matrices semidefinidas positivas es que todos sus autovalores son mayores o iguales a cero. Adem´as, este tipo de matrices se puede expresar como una combinaci´on lineal del producto de sus autovalores y autovectores: J=µ1r1r> 1+µ2r2r> 2+···+µNrNr> N= N X i=1 µirir> i,(2.47) 102 CAP´ ITULO 2. ESTIMACI ´ ON DEL FLUJO ´ OPTICO EN SECUENCIAS DE IM ´ AGENES El modelo de energ´ıa de mayor grado de libertad de todos los presentados en este apartado corresponde al F. En este prototipo la funci´on de robustificaci´on se aplica por separado sobre cada autovector de cada invarianza. La ubicaci´on coincide con el punto (2) de la ecuaci´on (2.54). Est´a colocado a la derecha de ambos sumatorios sin ejercer ninguna influencia sobre ellos, o dicho de otra forma, actuando s´olo en cada direcci´on de cada invarianza. De igual modo, el t´ermino de suavizado act´ua inversamente proporcional sobre las direcciones penalizadas por el t´ermino de ligadura. La estructura de cada t´ermino del modelo de energ´ıa se descompone en tres niveles: (i) en el nivel superior se sit´ua el sumatorio de las invarianzas creando el tratamiento diferenciado de cada invarianza, (ii) en un nivel intermedio, el sumatorio de la descomposici´on del tensor, que determina las direcciones dominantes y la magnitud del tensor asociado y, por ´ultimo, (iii) en el nivel inferior, la funci´on de robustificaci´on que se aplica ´unicamente a cada direcci´on de cada invarianza definida en el modelo de energ´ıa. EF(h) = ZΩ Nc X c=1 γc 3 X i=1 ψµic(h>ric)2!dx dy dz +ZΩ Nc X c=1 αc 3 X i=1 ψg(µic)(∇u>ric)2+ (∇v>ric)2!dx dy dz. Para finalizar la presentaci´on de los distintos modelos de energ´ıa de este trabajo quisiera hacer menci´on a las principales aportaciones de este framework. (i) El uso de tensores en todo el modelo de energ´ıa y (ii) su descomposici´on en autovalores y autovectores. Gracias a la combinaci´on del tensor de movimiento y las funciones de robustificaci´on ha sido posible el desarrollo de un framework que establece una relaci´on de complementariedad entre el t´ermino de ligadura y el de regularizaci´on. Adem´as, los prototipos propuestos son f´acilmente adaptables y escalables para cualquier tipo de invarianza. 2.5.3. Minimizaci´on de la Energ´ıa En este apartado se describe la minimizaci´on de los seis prototipos presentados en el apartado anterior. Comenzaremos con la minimizaci´on del prototipo Adado su sencillez. Iremos derivando cada uno de los restantes modelos incrementando la complejidad de los sistemas de ecuaciones obtenidos. En la derivaci´on de los prototipos Cal Faparecen t´erminos no-lineales que hacen que la minimizaci´on no sea trivial. Vamos a seguir una estructura muy parecida a la del apartado anterior. Se presentan cada una de las ecuaciones de Euler–Lagrange asociadas a los modelos de energ´ıa comentando las peculiaridades de cada uno de ellos. Para facilitar la comprensi´on al lector comenzaremos con el prototipo m´as b´asico (A, lineal con una sola invarianza) y terminaremos con el modelo m´as complejo y flexible (F, m´ultiples invarianzas con t´erminos no cuadr´aticos separados). Prototipo A/B: 2.5. M ´ ETODOS VARIACIONALES BASADOS EN EL AN ´ ALISIS ESPECTRAL 103 El modelo de energ´ıa es el m´as simple que nos podemos encontrar ya que no incluye ninguna t´ecnica de robustificaci´on. El t´ermino de ligadura puede incluir una o varias invarianzas y el t´ermino de regularizaci´on que incorpora el tensor inverso al definido en el de ligadura. Como vimos en el apartado anterior los modelos de energ´ıa de los prototipos A y B son id´enticos. La ´unica diferencia est´a en la interpretaci´on de la descomposici´on espectral. Por este motivo hemos unificado la minimizaci´on de ambos prototipos. Las ecuaciones de Euler-Lagrange asociadas a estos prototipos son 3 X i=1 µi(h>ri)ri1−α div (D(u, v, ri)∇u) = 0, 3 X i=1 µi(h>ri)ri2−α div (D(u, v, ri)∇v) = 0, donde D(u, v, ri) = 3 X i=1 g(µi)rir> i. En las ecuaciones de Euler-Lagrange se aprecia el formalismo y elegancia del framework que definimos en este trabajo. Independientemente de la invarianza que se incluya en el t´ermino de ligadura ´esta se expresar´a en funci´on de los autovalores y autovectores del tensor de movimiento. Esta forma de representaci´on de los tensores ofrece la m´axima flexibilidad a nuestro modelo ya que una vez hecha la descomposici´on espectral de las invarianzas el resto del esquema permanece inalterado. Prototipo C: Las ecuaciones de Euler-Lagrange de los siguientes prototipos crecen en complejidad y su minimizaci´on se convierte en una tarea no trivial. Estos prototipos incluyen t´erminos no cuadr´aticos. Adem´as los tensores incluidos en estos t´erminos permiten varias invarianzas lo que complica a´un m´as la formulaci´on de las ecuaciones de Euler-Lagrange. Las ecuaciones de Euler-Lagrange asociadas a este modelo son ψ0 3 X i=1 µi(h>ri)2!3 X i=1 µi(h>ri)ri1−α div (D(u, v, ri)∇u)=0, ψ0 3 X i=1 µi(h>ri)2!3 X i=1 µi(h>ri)ri2−α div (D(u, v, ri)∇v)=0, donde D(u, v, ri) = ψ0 3 X i=1 g(µi)(∇u>ri)2+ (∇v>ri)2!3 X i=1 g(µi)rir> i. El prototipo Cincorpora una funci´on de robustificaci´on que trata de forma conjunta a las invarianzas y a todas las direcciones de cada tensor. Cada ecuaci´on de Euler-Lagrange 104 CAP´ ITULO 2. ESTIMACI ´ ON DEL FLUJO ´ OPTICO EN SECUENCIAS DE IM ´ AGENES se compone de la derivada de la funci´on de robustificaci´on que act´ua como un peso sobre la estimaci´on del t´ermino de ligadura. Ese peso se agrega a la estimaci´on obtenida de fusionar todas las direcciones a trav´es del sumatorio situado a la derecha de la funci´on de robustificaci´on. En el t´ermino de regularizaci´on ocurre exactamente lo mismo pero lo ´unico que cambia es el valor de los autovalores, que son los inversos a los definidos en el t´ermino de ligadura. Prototipo D: Partiendo del modelo anterior e intercambiando la funci´on de robustificaci´on y el sumario que expresa las tres direcciones en la que se ha descompuesto el tensor de movimiento obtendr´ıamos el prototipo D. Ahora cada autovector es penalizado independientemente. Las ecuaciones de Euler-Lagrange asociadas a este modelo son 3 X i=1 ψ0µi(h>ri)2µi(h>ri)ri1−α div (D(u, v, ri)∇u) = 0, 3 X i=1 ψ0µi(h>ri)2µi(h>ri)ri2−α div (D(u, v, ri)∇v) = 0, donde D(u, v, ri) = 3 X i=1 ψ0g(µi)(∇u>ri)2+ (∇v>ri)2g(µi)rir> i. Al igual que ocurr´ıa en el modelo anterior todas las invarianzas se tratan conjuntamente. La derivada de la funci´on de robustificaci´on se convierte ahora en un peso asociado a cada direcci´on del autovector. Por un lado, nuestro modelo de energ´ıa se compone ´unicamente dos t´erminos, independientemente del n´umero de invarianzas que incorpore. Por otro lado, la complejidad de cada t´ermino no var´ıa cuando se incrementa el n´umero de invarianzas. Lo que aumenta es la complejidad de la estructura del tensor pero una vez hecha su descomposici´on espectral su manipulaci´on es sencilla. Las direcciones expresadas por los autovectores representan una combinaci´on de todas las invarianzas. Prototipo E: Los dos ´ultimos modelos tratan de forma separada cada una de las invarianzas. Estos prototipos crean funcionales de energ´ıa que disponen por cada invarianza de un par de t´erminos de ligadura y regularizaci´on. Ahora el sumatorio que engloba a todas las invarianzas queda situado fuera de la robustificaci´on creando esa diferenciaci´on entre invarianzas. 2.5. M ´ ETODOS VARIACIONALES BASADOS EN EL AN ´ ALISIS ESPECTRAL 105 Las ecuaciones de Euler-Lagrange asociadas a este modelo son Nc X c=1 γcψ0 3 X i=1 µic(h>ric)2! 3 X i=1 µic(h>ric)ric1!− Nc X c=1 αcdiv (Dc(u, v, ri)∇u) = 0, Nc X c=1 γcψ0 3 X i=1 µic(h>ric)2! 3 X i=1 µic(h>ric)ric2!− Nc X c=1 αcdiv (Dc(u, v, ri)∇v) = 0, donde Dc(u, v, ri) = ψ0 3 X i=1 g(µic)(∇u>ric)2+ (∇v>ric)2! 3 X i=1 g(µic)ricr> ic!. Los sumatorios de las invarianzas est´an situados a la izquierda de cada t´ermino generando un nuevo elemento por cada invarianza. Ahora por cada invarianza tenemos un par´ametro de regularizaci´on asociado, αc, pudi´endose controlar de una forma m´as precisa e independiente el proceso de difusi´on. Prototipo F: Como ya se ha comentado, el prototipo F es el m´as flexible de todos los presentados en esta secci´on. La funci´on de robustificaci´on est´a situada de tal forma que ofrece total grado de libertad tanto a nivel de invarianza como de autodirecci´on. Este caso recoge dos de las ideas presentadas en los prototipos D yE. Del Dse toma la idea de aplicar la funci´on de robustificaci´on sobre cada autodirecci´on por separado. El prototipo E desacopla las estimaciones de cada invarianza en distintos t´erminos. Las ecuaciones de Euler-Lagrange asociadas a este prototipo son Nc X c=1 γc 3 X i=1 ψ0µic(h>ric)2µic(h>ric)ric1!− Nc X c=1 αcdiv (Dc(u, v, ri)∇u) = 0, Nc X c=1 γc 3 X i=1 ψ0µic(h>ric)2µic(h>ric)ric2!− Nc X c=1 αcdiv (Dc(u, v, ri)∇v) = 0, donde Dc(u, v, ri) = 3 X i=1 ψ0g(µic)(∇u>ric)2+ (∇v>ric)2g(µic)ricr> ic. Al igual que ocurr´ıa en el prototipo anterior, E, en el t´ermino de regularizaci´on se define un par´ametro de regularizaci´on asociado cada invarianza pudi´endose controlar independientemente la cantidad de difusi´on aplicada por cada invarianza. 106 CAP´ ITULO 2. ESTIMACI ´ ON DEL FLUJO ´ OPTICO EN SECUENCIAS DE IM ´ AGENES 2.5.4. Esquema Num´erico En este apartado se describe el esquema num´erico del prototipo F, que por su flexibilidad consideramos el caso m´as completo de todos los presentados. S´olo habr´a peque˜nas diferencias respecto al resto de prototipos. Para hallar las componentes del flujo ´optico (u, v) vamos a resolver el siguiente sistema el´ıptico 0 = Nc X c=1 αcdiv (Dc(u, v, ri)∇u)− Nc X c=1 γc 3 X i=1 ψ0µic(h>ric)2µic(h>ric)ric1! 0 = Nc X c=1 αcdiv (Dc(u, v, ri)∇v)− Nc X c=1 γc 3 X i=1 ψ0µic(h>ric)2µic(h>ric)ric2!,(2.57) donde Dc(u, v, ri) = 3 X i=1 ψ0g(µic)(∇u>ric)2+ (∇v>ric)2g(µic)ricr> ic. El tensor de movimiento se define en un dominio continuo donde se producen desplazamientos peque˜nos y la existencia de las derivadas est´a asegurada. Dado que nuestro modelo se basa en la descomposici´on espectral de dicho tensor, cuando nos encontramos ante desplazamientos largos de los objetos toda la teor´ıa alrededor de esta t´ecnica no se cumple. A continuaci´on vamos a describir c´omo se puede extender esta estructura para que pueda detectar los desplazamientos largos conservando la notaci´on y propiedades comentadas hasta ahora. La descripci´on es similar al desarrollo hecho en [Bruhn06a]. (fxu+fyv+ft)2=h>∇f∇f>h=h>Jh.(2.58) El tensor de movimiento es una forma de representar cualquier invarianza dentro de una estructura compacta. En la ecuaci´on (2.58) se obtiene su estructura a partir de la ecuaci´on de restricci´on del flujo ´optico. Para manejar los desplazamientos largos ser´a necesario deducir el tensor a partir de la ecuaci´on no lineal, ec. (2.44). Esto es posible hacerlo f´acilmente en el momento de la discretizaci´on. Siguiendo la misma estrategia que [Bruhn06a] se ha combinado un enfoque multipiramidal con una discretizaci´on fixed point iteration. Para eliminar los t´erminos nolineales presentes en el sistema de ecuaciones se ha optado por la linealizaci´on del t´ermino de ligadura y por dividir el flujo ´optico a estimar en la siguiente iteraci´on, hk+1 = (uk+1,vk+1,1)>, como la suma del flujo calculado en la iteraci´on actual, hk= (uk,vk,1)>, m´as un peque˜no incremento, dhk= (duk,dvk,1)>. Fijando como referencia el flujo actual queremos estimar su valor en la siguiente iteraci´on hallando un flujo desconocido que necesitamos para ir desde (uk,vk) a (uk+1,vk+1). Dicho de otra manera uk+1 =uk+duk, vk+1 =vk+dvk, hk+1 =hk+dhk. 2.5. M ´ ETODOS VARIACIONALES BASADOS EN EL AN ´ ALISIS ESPECTRAL 107 Si expresamos la ecuaci´on del flujo ´optico para largos desplazamientos, ec. (2.44), siguiendo la notaci´on utilizada en nuestro esquema num´erico tendr´ıamos f(x+hk)−f(x)=0. Para eliminar la no-linealidad de la expresi´on anterior y establecer una formulaci´on en funci´on de las nuevas inc´ognitas del flujo ´optico (duk, dvk) linealizamos f(x+hk+1) obteniendo f(x+hk+1) = f(x+hk) + fx(x+hk)duk+fy(x+hk)dvk.(2.59) Si sustituimos la ecuaci´on (2.59) en la discretizaci´on de cualquier invarianza podremos obtener una expresi´on cuadr´atica en la que tengamos un tensor de movimiento que sea capaz de reflejar los largos desplazamientos de la escena. (f(x+hk+1)−f(x))2≈ (f(x+hk) + fx(x+hk)duk+fy(x+hk)dvk−f(x))2= (fx(x+hk)duk+fy(x+hk)dvk+ (f(x+hk)−f(x)))2= (dhke ∇f(x+hk))2= (dhk)>e ∇f(x+hk)e ∇f>(x+hk) | {z } J(x+hk) (dhk),(2.60) dhk= (duk, dvk,1)>representa el incremento espaciotemporal del movimiento y e ∇es una variante del operador gradiente donde la ´ultima componente es una aproximaci´on de la derivada temporal mediante una diferencia. De acuerdo a la ecuaci´on (2.60) el tensor de movimiento se puede utilizar para representar cualquier invarianza que tenga en cuenta los desplazamientos largos. La descomposici´on espectral de este tensor seguir´a conservando esta caracter´ıstica. Por lo tanto, nuestro modelo est´a preparado para detectar cualquier movimiento presente en la imagen. Sin embargo, el modelado de los desplazamientos largos no se hace expl´ıcitamente en el funcional de energ´ıa. Esto se pospone al momento de la discretizaci´on. Adem´as, nuestro esquema num´erico ser´a embebido en un enfoque multipiramidal. Si fusionamos las ideas descritas en las ecuaciones (2.55) y (2.60), el esquema num´erico quedar´ıa expresado en funci´on del incremento del flujo dhk, ec. (2.61). dhk>e ∇f(x+hk)e ∇f>(x+hk)dhk=dhk>J(x+hk)dhk= dhk> 3 X i=1 µirir> i!dhk= 3 X i=1 µidhk>rir> idhk= 3 X i=1 µi(dhk>ri)2.(2.61) Esta notaci´on es f´acilmente extensible para m´ultiples invarianzas. M´etodo de Resoluci´on del Sistema de Ecuaciones La soluci´on del sistema de ecuaciones de Euler-Lagrange descrito en el apartado anterior, ec. (2.57), se obtiene mediante un m´etodo num´erico iterativo. Successive 108 CAP´ ITULO 2. ESTIMACI ´ ON DEL FLUJO ´ OPTICO EN SECUENCIAS DE IM ´ AGENES Overrelaxation (SOR, [Young71]) se trata de un m´etodo iterativo basado en una extrapolaci´on de los resultados de Gauss-Seidel. Es un compromiso entre simplicidad y eficiencia. Aunque el coste por iteraci´on sea algo superior que el m´etodo de Gauss-Seidel necesita menos iteraciones para converger lo que lo convierte en una buena opci´on. La estructura general de m´etodo SOR se expresa en la siguiente ecuaci´on xk+1 i= (1 −β)xk i+β(Ai,i)−1 bi− i−1 X j=1 Ai,j xk+1 j− 2N X j=i+1 Ai,j xk j , | {z } Iteraci´on de Gauss-Seidel (2.62) donde xpuede ser una de las componentes del flujo ´optico (u, v) y i= 1, ..., 2N;Aes una matriz que almacena los coeficientes del sistema de ecuaciones a resolver. Este m´etodo num´erico ser´a utilizado en la implementaci´on del nuestro m´etodo. Discretizaci´on del Sistema de Ecuaciones de Euler–Lagrange Una vez definido el modelo de energ´ıa y minimizarlo el siguiente paso consiste en la resoluci´on del sistema de ecuaciones de Euler–Lagrange, ec. (2.57). Para ello debemos resolver un sistema de ecuaciones no-lineal del tipo Am(xm) = bm,(2.63) donde Am(xm) es un operador no-lineal y bmes la parte derecha del sistema. Am(xm) se puede descomponer en Am(xm) = Bm(xm)xm+cm(xm),(2.64) donde Bm(xm) y cm(xm) son operadores no-lineales pero para cada valor de xm,Bm(xm) es una matriz sim´etrica y definida positiva; cm(xm) es un vector. Los m´etodos Lagged-Diffusivity ([Kacur68, Fucik73, Chan99]) se basan en la idea de la resoluci´on de un sistema de ecuaciones nolineal, como es el caso de (2.63), a trav´es de su descomposici´on en un conjunto de problemas lineales. La resoluci´on de estos problemas lineales se puede llevar a cabo con t´ecnicas est´andar, como son el m´etodo de SOR o Gauss–Seidel. La descomposici´on hecha en la ecuaci´on (2.64) permite aprovechar las propiedades de la matriz Bm(xm) y el vector cm(xm) para convertir el sistema inicial, ec. (2.63), en lineal, gracias a la evaluaci´on de los operadores Bm(xm) y cm(xm) en el instante anterior k xm,k+1 = (Bm(xm,k))−1(bm−cm(xm,k)), Bm(xm,k)xm,k+1 = (bm−cm(xm,k)).(2.65) Para la resoluci´on del sistema de ecuaciones descrito en la ecuaci´on (2.57) se ha optado por la combinaci´on de las t´ecnicas de fixed point iterations y el m´etodo Lagged-Diffusivity. A continuaci´on, vamos a desarrollar el esquema num´erico resultante. Para ello, vamos a descomponer el proceso en tres niveles: 2.5. M ´ ETODOS VARIACIONALES BASADOS EN EL AN ´ ALISIS ESPECTRAL 109 1. Enfoque Multipiramidal. Para detectar los desplazamientos largos nuestro m´etodo se apoya en un enfoque multipiramidal. En este enfoque se crean un cierto n´umero de escalas s1, s2, . . . , sn, donde cada una representa im´agenes de distinto tama˜no. Durante el cambio de escala, los valores del flujo deben ser actualizados y adaptados a la nueva escala. Para hallar las componentes del flujo ´optico {us, vs}tenemos que resolver en cada escala el sistema de ecuaciones descrito en la ecuaci´on (2.57). El ´ındice srefleja la escala actual. El tensor de movimiento se recalcula y descompone espectralmente en cada escala antes de volver a resolver el sistema de ecuaciones no lineal. 2. Lagged-Diffusivity. En cada escala se resuelve el sistema de ecuaciones de Euler-Lagrange. Este sistema se convierte en lineal gracias al Lagged-Diffusivity. El sistema se resuelve un cierto n´umero de iteraciones. El´ındice lrepresenta el punto de iteraci´on en el que est´a fijado cada elemento del sistema, p.e. (us,l, vs,l). En este nivel se resuelve un problema convexo nolineal a trav´es de un conjunto de problemas convexos lineales manteniendo la nolinealidad en el t´ermino de ligadura y regularizaci´on. En cada iteraci´on se recalculan las derivadas de las funciones de robustificaci´on y el tensor de difusi´on. En nuestro caso, el sistema resultante es 0 = Nc X c=1 αcdiv Dc(us,l +dus,l, vs,l +dvs,l, rs i)∇(us,l +dus,l) − Nc X c=1 γc 3 X i=1 ψ0µs ic((dhs,l)>rs ic)2µs ic((dhs,l)>rs ic)rs ic1, 0 = Nc X c=1 αcdiv Dc(us,l +dus,l, vs,l +dvs,l, rs i)∇(vs,l +dvs,l) − Nc X c=1 γc 3 X i=1 ψ0µs ic((dhs,l)>rs ic)2µs ic((dhs,l)>rs ic)rs ic2,(2.66) donde Dc(us,l +dus,l, vs,l +dvs,l, rs i) es Dc= 3 X i=1 ψ0g(µs ic)(∇(us,l +dus,l)>rs ic)2+ (∇(vs,l +dvs,l)>rs ic)2g(µs ic)rs icrs> ic. 3. Resoluci´on del sistema de ecuaciones. En el nivel inferior se resuelve el sistema de ecuaciones lineales calculando los nuevos valores del incremento del flujo. Para ello, se utiliza el m´etodo num´erico SOR. En este nivel es necesario la inclusi´on de un nuevo ´ındice kque plasme la iteraci´on del SOR en el que nos encontramos. Para simplificar la notaci´on de la discretizaci´on vamos a introducir una serie de expresiones (dhs,l,k)>rs ic= (rs i1cdus,l,k +rs i2cdvs,l,k +rs i3c) ψ0s,l D=ψ0µs ic((dhs,l)>rs i)2 110 CAP´ ITULO 2. ESTIMACI ´ ON DEL FLUJO ´ OPTICO EN SECUENCIAS DE IM ´ AGENES La matriz D, se define en cada p´ıxel icomo: Di=  aibici bidiei cieifi . La discretizaci´on del t´ermino de la divergencia en cada p´ıxel ise realiza de la siguiente manera: div(Di∇u) =   ai∂xu+bi∂yu+ci∂tu bi∂xu+di∂yu+ei∂tu ci∂xu+ei∂yu+fi∂tu =∂x(ai∂xu) + ∂x(bi∂yu) + ∂x(ci∂tu) +∂y(bi∂xu) + ∂y(di∂yu) + ∂y(ei∂tu) +∂t(ci∂xu) + ∂t(ei∂yu) + ∂t(fi∂tu), N∗ ise define como los vecinos alrededor del p´ıxel i. Utilizando un esquema de diferencias est´andar la divergencia se puede expresar como: div(Di∇ui) = X j∈N∗ i wjuj+wiui, para los coeficientes adecuados wj. Para el caso de div(Di∇v) se procede de la misma manera. En el sistema de ecuaciones (2.57) utilizamos un esquema semi-impl´ıcito en el t´ermino de ligadura e impl´ıcito en el de regularizaci´on. El esquema iterativo a implementar surge de la inclusi´on de las ecuaciones (2.66) en el esquema SOR, ec. (2.62). dus,l,k+1 i= (1 −β)dus,l,k i +β Nc X c=1 αc X j∈N∗ i wj(us,l j+dus,l,k j)  Nc X c=1 αc X j∈N∗ i wj + Nc X c=1 γc 3 X i=1 ψ0s,l Dµs irs2 i1! −β Nc X c=1 αc X j∈N∗ i wjus,l j + Nc X c=1 γc 3 X i=1 ψ0s,l Dµs i(rs i2dvs,l,k i+rs i3)rs i1! Nc X c=1 αc X j∈N∗ i wj + Nc X c=1 γc 3 X i=1 ψ0s,l Dµs irs2 i1! 2.5. M ´ ETODOS VARIACIONALES BASADOS EN EL AN ´ ALISIS ESPECTRAL 111 en la direcci´on vertical dvs,l,k+1 i= (1 −β)dvs,l,k i +β Nc X c=1 αc X j∈N∗ i wj(vs,l j+dvs,l,k j)  Nc X c=1 αc X j∈N∗ i wj + Nc X c=1 γc 3 X i=1 ψ0s,l Dµs irs2 i2! −β Nc X c=1 αc X j∈N∗ i wjus,l j + Nc X c=1 γc 3 X i=1 ψ0s,l Dµs i(rs i1dus,l,k i+rs i3)rs i2! Nc X c=1 αc X j∈N∗ i wj + Nc X c=1 γc 3 X i=1 ψ0s,l Dµs irs2 i2! donde βes un par´ametro de relajaci´on cuyo valor est´a en el intervalo (0,2). El sub´ındice irepresenta un p´ıxel en la imagen. Los valores m´as comunes para βsuelen oscilar entre 1.5 y 1.99. 2.5.5. Resultados Experimentales En este apartado se presentan los resultados obtenidos con los distintos prototipos descritos en este trabajo. Se ha establecido una comparaci´on entre todos ellos y los mejores m´etodos propuestos en la literatura. Nuestros prototipos se basan en un modelo anisotr´opico que utiliza informaci´on espaciotemporal y la descomposici´on espectral de los tensores del t´ermino de ligadura y de regularizaci´on. La inclusi´on de t´erminos no cuadr´aticos da robustez a nuestro m´etodo frente al ruido. En estos experimentos se quiere demostrar que la elegante descomposici´on de nuestro m´etodo no s´olo se aprecia desde punto de vista te´orico sino tambi´en pr´actico. En los tests se ha utilizado una secuencia sint´etica, Yosemite con nubes, y una real, Rheinhafen. Antes de comentar los resultados conviene recordar los distintos par´ametros que intervienen en nuestro m´etodo. El par´ametro de regularizaci´on, α, determina el peso del t´ermino de suavizado. Cuando el valor de este par´ametro es alto genera campos de desplazamiento m´as suaves, el t´ermino de suavizado gana importancia respecto al de ligadura. Los prototipos del Aal Dtienen un s´olo par´ametro de regularizaci´on. Para los prototipos EyFexistir´a un par´ametro por cada invarianza definida en t´ermino de ligadura, αcpara c= 1, ..., N siendo Nel n´umero de invarianzas. El ratio de descenso η∈(0,1) para el c´alculo de las escalas. En [Weickert04] se recomienda que el factor de escala ηest´e dentro del rango [0,5,0,95]. En todos los experimentos se ha fijado su valor a η= 0,9. γces un peso que se utiliza en el t´ermino de ligadura para ponderar la importancia de cada invarianza c. 118 CAP´ ITULO 2. ESTIMACI ´ ON DEL FLUJO ´ OPTICO EN SECUENCIAS DE IM ´ AGENES Figura 2.26: Los distintos campos de desplazamiento obtenidos con los prototipos Aal F. En la parte superior, los campos de desplazamiento estimados por el prototipo AyB. En el medio, la estimaci´on del prototipo CyD. En la parte inferior, la soluci´on ofrecida por EyF. 2.5. M ´ ETODOS VARIACIONALES BASADOS EN EL AN ´ ALISIS ESPECTRAL 119 Cuadro 2.10: Lista de par´ametros utilizados en los prototipos Aal Fpara la secuencia de Rheinhafen.Nc= N´umero de invarianzas definidas en el modelo de energ´ıa. γ1= peso asignado a la primera invarianza en el t´ermino de ligadura. γ2= peso asignado a la segunda invarianza en el t´ermino de ligadura. α1= peso asignado al primer t´ermino de suavizado. α2= peso asignado al segundo t´ermino de suavizado. Rheinhafen εd= 10−3,εs= 10−3,εg= 10−3, η= 0,9, β= 1,99, σs= 1,0, σt= 10−3 M´etodo Ncγ1γ2α1α2 Prototipo A 1 1 - 0.00045 - Prototipo B 2 1 0.001 50 - Prototipo C 2 1 0.001 0.01 - Prototipo D 2 1 0.001 0.01 - Prototipo E 2 1 0.05 0.03 0.001 Prototipo F 2 1 1 0.04 0.01 en el t´ermino de ligadura. La soluci´on obtenida por este prototipo se ve afectada por los cambios de iluminaci´on y el ruido presente en la secuencia. Dado que en Rheinhafen el ruido influye m´as que los cambios de iluminaci´on la inclusi´on de una segunda invarianza (prototipo B) no consigue mejorar la estimaci´on del prototipo A. El resto de los prototipos incluyen funciones de robustificaci´on lo que da al m´etodo una mayor insensibilidad frente al ruido. No existe visualmente grandes diferencias entre las versiones joint CyDo entre separate EyF. Si se percibe la influencia de la robustificaci´on junta o separada sobre las invarianzas definidas en el modelo. Cabe destacar que el movimiento de los coches situados al fondo ha sido detectado con bastante nitidez. 120 CAP´ ITULO 2. ESTIMACI ´ ON DEL FLUJO ´ OPTICO EN SECUENCIAS DE IM ´ AGENES 2.6. Conclusiones En este cap´ıtulo se han presentado tres trabajos relativos a la estimaci´on del flujo ´optico con informaci´on procedente de varias im´agenes. En los m´etodos descritos se ha querido plasmar distintos enfoques a la hora de estimar el flujo ´optico en una secuencia de im´agenes: (1) un m´etodo espacial que combina informaci´on multicanal, (2) un m´etodo espacial al que se le ha a˜nadido un t´ermino de regularizaci´on temporal y, (3) un m´etodo espaciotemporal que comparte informaci´on entre el t´ermino de ligadura y suavizado. A continuaci´on, se describe las conclusiones m´as importantes extra´ıdas de estos tres trabajos. 2.6.1. M´etodo Variacional Multicanal La incorporaci´on en los sat´elites de m´ultiples sensores ha multiplicado la cantidad de informaci´on disponible haciendo necesario la automatizaci´on de las tareas de almacenamiento y procesado. La meteorolog´ıa y la climatolog´ıa son campos de estudio que m´as hacen uso de las im´agenes sat´elites. Estas im´agenes facilitan el an´alisis y seguimiento de los distintos fen´omenos que ocurren en la atm´osfera. Tradicionalmente se han utilizado t´ecnicas de correlaci´on para el an´alisis de los fen´omenos atmosf´ericos debido a su simplicidad y buenos resultados. Sin embargo, son computacionalmente costosos y no ofrecen soluciones densas. Los m´etodos variacionales han demostrado un enorme ´exito a la hora de resolver problemas de visi´on por ordenador, como es la estimaci´on del flujo ´optico. En este trabajo se ha propuesto un m´etodo variacional para la detecci´on del desplazamiento de las estructuras nubosas en las im´agenes sat´elites multicanal. Se trata de una extensi´on multicanal de un conocido m´etodo para la estimaci´on del flujo ´optico ([Alvarez00]). El modelo de energ´ıa se compone de un t´ermino de ligadura incluye la suposici´on lambertiana, aplicada a cada canal por separado, mientras que en el t´ermino de suavizado se utiliza el operador de Nagel-Enkelmann ([Nagel86]) con algunas mejoras. Dado que el tensor de difusi´on de Nagel-Enkelmann usa el gradiente de la imagen para determinar la cantidad de difusi´on se han desarrollado dos estrategias para calcular dicho gradiente a partir de la informaci´on multicanal disponible. La combinaci´on de la informaci´on de todos los canales permite mejorar la estimaci´on del movimiento de las nubes. Dado que las secuencias sat´elites son datos reales y no disponemos de los desplazamientos de los p´ıxeles se han creado dos secuencias sint´eticas usando un modelo de movimiento realista. Con este tipo de secuencias es posible realizar una comparaci´on cuantitativa de las estimaciones ofrecidas por los m´etodos. Los resultados experimentales demuestran que el m´etodo multicanal mejora significativamente la precisi´on de las estimaciones si la comparamos con su versi´on monocanal. Los canales visibles e infrarrojo son los que m´as informaci´on aportan acerca del desplazamiento de las nubes. 2.6.2. M´etodo Variacional con Regularizaci´on Temporal no Continua En este trabajo se ha descrito un m´etodo variacional que incorpora un novedoso t´ermino de regularizaci´on exclusivamente temporal. Algunos autores han propuesto 2.6. CONCLUSIONES 121 regularizadores que tratan de forma conjunta la informaci´on espacial y temporal. Estos regularizadores imponen una continuidad en el modelo obligando a que los desplazamientos de los objetos sean peque˜nos tanto en la direcci´on espacial como en la temporal. En situaciones reales un modelo discontinuo, como el nuestro, se ajusta mejor. El funcional de energ´ıa propuesto se trata de una modificaci´on del modelo de NagelEnkelmann [Nagel86] al que se le ha a˜nadido un t´ermino de regularizaci´on temporal. La divisi´on del t´ermino de regularizaci´on en dos, uno espacial y otro temporal, pretende evitar los inconvenientes derivados del acoplamiento de la informaci´on espaciotemporal. Se han dise˜nado dos variantes del t´ermino temporal: (1) Temporal cuya estimaci´on se basa en la informaci´on del flujo del frame siguiente y, (2) Bi–Temporal que utiliza la informaci´on del flujo en ambos sentidos. Un inconveniente que tiene la versi´on Temporal est´a en su fuerte dependencia respecto al ´ultimo flujo de la secuencia, de forma que si la estimaci´on del ´ultimo flujo es mala ese error se propaga a los flujos anteriores. Por el contrario, en la versi´on Bi–Temporal los flujos se compensan entre sus vecinos sin verse influenciados por uno en concreto. Este hecho se aprecia en los experimentos donde el m´etodo bidireccional mejora sustancialmente las estimaciones respecto a su hom´ologo unidireccional. En los resultados experimentales se observa que la incorporaci´on del nuevo t´ermino de regularizaci´on temporal aporta estabilidad a las soluciones del m´etodo espacial. El m´etodo propuesto en esta tesis se trata de unos de los primeros trabajos que considera regularizadores temporales para largos desplazamientos. 2.6.3. M´etodo Variacional basado en el An´alisis Espectral En el ´ultimo trabajo de este tema se propone un framework para la estimaci´on del flujo ´optico usando secuencias de im´agenes. Este framework permite la definici´on de un modelo de energ´ıa gen´erico que gracias a la combinaci´on de algunas de las t´ecnicas m´as innovadoras propuestas en la literatura es f´acilmente adaptable y escalable para cualquier tipo de invarianzas. El tensor de movimiento es una t´ecnica que permite representar cualquier invarianza en una estructura compacta. La descomposici´on espectral de este tensor permite la f´acil identificaci´on de las direcciones dominantes del flujo. Esta informaci´on resulta muy pr´actica a la hora de guiar de una forma m´as precisa el proceso de difusi´on. Al mismo tiempo, se puede establecer una complementariedad entre el t´ermino de ligadura y el de regularizaci´on. La utilizaci´on de funciones de robustificaci´on aportan insensibilidad al m´etodo frente a la presencia de outliers. La combinaci´on de esta t´ecnica junto con la descomposici´on del tensor de movimiento permite la creaci´on de cuatro modelos de energ´ıa distintos, donde las autodirecciones son penalizadas conjuntamente o por separado en funci´on de la colocaci´on de la robustificaci´on. Los experimentos realizados con las distintas variantes de nuestro modelo nos han permitido comprobar que la elegante descomposici´on de nuestro m´etodo no s´olo es apreciable desde punto de vista te´orico sino tambi´en pr´actico. Para ello, se ha utilizado la secuencia sint´etica de Yosemite con nubes para comparar los modelos con los mejores m´etodos propuestos en la literatura. Los resultados son bastante buenos mejorando las estimaciones de otros m´etodos que incorporan t´ecnicas parecidas. Incluso uno de los 122 CAP´ ITULO 2. ESTIMACI ´ ON DEL FLUJO ´ OPTICO EN SECUENCIAS DE IM ´ AGENES modelos obtiene el segundo mejor registro de la literatura. Recordemos que este m´etodo no incluye ninguna t´ecnica de segmentaci´on. En la secuencia real de Rheinhafen tambi´en se pone de manifiesto la calidad de las soluciones obtenidas con los distintos prototipos. Las principales aportaciones realizadas en este trabajo son: (i) el uso de tensores en toda el modelo de energ´ıa y (ii) la descomposici´on espectral de los mismos. Gracias a la combinaci´on del tensor de movimiento y las funciones de robustificaci´on ha sido posible el desarrollo de un framework que establece una relaci´on de complementariedad entre el t´ermino de ligadura y el de regularizaci´on. Cap´ıtulo 3 Estimaci´on del Mapa de Disparidad en Pares Est´ereo 3.1. Introducci´on En la visi´on estereosc´opica disponemos de dos vistas de la misma escena en el mismo instante de tiempo. La estimaci´on del mapa de disparidad consiste en el c´alculo del desplazamiento de los p´ıxeles de una vista a la otra. Este problema se reduce a una b´usqueda de correspondencias entre dos im´agenes por lo que en cierto modo tiene muchas similitudes con la estimaci´on del flujo ´optico. La componente temporal que tienen las secuencias de im´agenes frente a los pares est´ereos supone un impedimento para la adaptaci´on de muchas de las ideas surgidas en el problema del flujo ´optico. Sin embargo, la visi´on estereosc´opica dispone de una herramienta muy ´util como es la geometr´ıa epipolar. La geometr´ıa epipolar permite acotar el ´area de b´usqueda de las correspondencias. Toda esta teor´ıa se basa en la existencia de un sistema de c´amaras perfectamente calibradas. Como parte fundamental de esta tesis se ha dedicado un cap´ıtulo completo a la descripci´on de los distintos m´etodos desarrollados para la estimaci´on del mapa de disparidad. En este cap´ıtulo se describen algunas de las t´ecnicas m´as exitosas surgidas en los ´ultimos a˜nos y c´omo la combinaci´on de ellas permite aumentar la precisi´on de los estimaciones. Antes de entrar en detalle en cada uno de ellos se comentar´an las contribuciones hechas a la literatura. 3.1.1. Contribuciones de este Cap´ıtulo Las contribuciones en el ´ambito cient´ıfico que se han hecho en este cap´ıtulo son las siguientes: M´etodo para la Estimaci´on del Mapa de Disparidad utilizando una Secuencia de Pares Est´ereo: La primera aportaci´on novedosa supone la creaci´on de un m´etodo que combina la informaci´on est´ereo Alvarez/Deriche/S´anchez/Weickert (2002) [Alvarez02b] con 123 124 CAP´ ITULO 3. ESTIMACI ´ ON DEL MAPA DE DISPARIDAD EN PARES EST ´ EREO la del flujo ´optico Alvarez/Weickert/S´anchez (2000) [Alvarez00] aumentando la robustez de las estimaciones al fusionar dentro del mismo modelo de energ´ıa las ideas de dos m´etodos de gran precisi´on. La estimaci´on del mapa de disparidad se apoya en la geometr´ıa epipolar para establecer las correspondencias. Esas mismas correspondencias se pueden obtener mediante el c´alculo del flujo ´optico a lo largo de la secuencia de ambas c´amaras. Por lo tanto, en una secuencia de pares est´ereo disponemos de hasta cuatro vistas de un punto 3D en dos instantes de tiempo distintos. Este m´etodo pretende aprovechar las ventajas de los m´etodos de flujo ´optico y est´ereo para mejorar la estimaci´on de los mapas de disparidad. Al disponer de m´as informaci´on permite que la b´usqueda de las correspondencias sea m´as robusta. Combinaci´on de un M´etodo Variacional Espacial y uno de Graph–cuts para la Estimaci´on del Mapa de Disparidad: La segunda aportaci´on novedosa consiste en la combinaci´on de un m´etodo variacional Alvarez/Deriche/S´anchez/Weickert (2002) [Alvarez02b] y una t´ecnica de graph– cuts [Kolmogorov01, Boykov04]. El m´etodo de graph–cuts se ha utilizado para obtener una aproximaci´on inicial del mapa de disparidad que el m´etodo variacional se encargar´ıa de refinar. Normalmente, la t´ecnica m´as empleada para estimar la inicializaci´on es una basada en la correlaci´on a ventanas. La soluci´on en los m´etodos variacionales se alcanza mediante la minimizaci´on de un funcional de energ´ıa. Para garantizar y acelerar la obtenci´on de una soluci´on pr´oxima al m´ınimo global se suele ayudar al algoritmo mediante una inicializaci´on. O sea, se le ofrece una soluci´on relativamente pr´oxima al m´ınimo global de forma que el algoritmo converja r´apidamente. Es muy importante que la inicializaci´on est´e cerca del m´ınimo porque en el caso contrario la convergencia no est´a garantizada. En los ´ultimos a˜nos se han propuesto m´etodos graph–cuts para la estimaci´on del mapa de disparidad. Con esta t´ecnica se obtienen muy buenos resultados pero en precisi´on entera. Una buena inicializaci´on no tiene porqu´e tener una gran precisi´on s´olo se requiere que est´e pr´oxima a la soluci´on final. Por este motivo, la combinaci´on de un m´etodo variacional [Alvarez02b] y otro de graph– cuts [Kolmogorov01, Boykov04] parece m´as ventajosa que un variacional con una t´ecnica de correlaci´on. 3.2. ESTIMACI ´ ON DE LA DISPARIDAD USANDO UNA SECUENCIA ESTEREOSC ´ OPICA125 3.2. M´etodo para la Estimaci´on del Mapa de Disparidad utilizando una Secuencia de Pares Est´ereo En este apartado se describe un nuevo m´etodo para la recuperaci´on de la geometr´ıa 3D a partir de una secuencia de pares est´ereo. Disponemos de dos c´amaras captando una escena y que est´an fijadas sobre un soporte r´ıgido orientadas en una cierta direcci´on. En esa escena pueden haber objetos est´aticos o din´amicos. Suponemos que el sistema de c´amaras est´a d´ebilmente calibrado y que las c´amaras est´an situadas en posici´on frontoparalela (ver apartado 1.2.1). Debido a las perturbaciones presentes en las im´agenes captadas por un sistema de visi´on estereosc´opico no es posible hacer un perfecto seguimiento del desplazamiento de los p´ıxeles de una vista a la otra. Una forma de minimizar el efecto que producen algunas perturbaciones, como puede ser las oclusiones o el ruido, consistir´ıa en transferir informaci´on desde los flujos ´opticos de cada c´amara a los mapas de disparidad. De esta forma podr´ıamos incluir informaci´on temporal en las estimaciones del flujo est´ereo y mejorar su precisi´on. El m´etodo presentado en este apartado combina la informaci´on est´ereo con la del flujo ´optico para mejorar la estimaci´on de los mapas de disparidad. La estimaci´on del flujo ´optico en ambas c´amaras se utiliza como restricci´on en el c´alculo de los mapas de disparidad haciendo que la b´usqueda de las correspondencias sea m´as robusta. El modelo de energ´ıa permite detectar los desplazamientos largos de los objetos. Para ello, se utiliza un enfoque multipiramidal en el que la soluci´on de la escala inferior se utilizar´a como aproximaci´on inicial de la superior. Este m´etodo es una continuaci´on de los trabajos presentados en Alvarez/Weickert/S´anchez (2000) [Alvarez00] y Alvarez/Deriche/S´anchez/Weickert (2002) [Alvarez02b]. Los m´etodos de ambos trabajos se basan en una t´ecnica de minimizaci´on de energ´ıa que ofrecen resultados precisos y densos. La organizaci´on de este trabajo es la siguiente: en primer lugar, se presenta la notaci´on utilizada para la descripci´on del modelo de energ´ıa. Una vez descrita la notaci´on, se comenta con todo lujo de detalle el modelo de energ´ıa y el esquema num´erico utilizado. Finalmente, se muestra los experimentos realizados con secuencias sint´eticas y reales que pondr´an de manifiesto la aportaci´on de la informaci´on del flujo ´optico en la estimaci´on del mapa de disparidad. 3.2.1. Notaci´on del Flujo ´ Optico y Mapas de Disparidad En este apartado se describe la notaci´on utilizada en el modelo de energ´ıa para combinar la informaci´on del flujo ´optico y del mapa de disparidad. El flujo ´optico representa el desplazamiento de cada p´ıxel de una imagen a la siguiente. Este desplazamiento se suele representar como la funci´on h(x) = (u(x), v(x))t. En nuestro caso, existen dos flujos ´opticos a estimar: el de la imagen derecha e izquierda. Para poder identificar correctamente a cada flujo es necesario introducir ´ındices que hagan referencia a cada c´amara. Un ´ındice ipara especificar un determinado flujo de la secuencia; y otro l, r para indicar si ese flujo pertenece a la c´amara izquierda o derecha, respectivamente. Dado que la secuencia de entrada esta compuesta por Npares est´ereos, existir´an N−1 126 CAP´ ITULO 3. ESTIMACI ´ ON DEL MAPA DE DISPARIDAD EN PARES EST ´ EREO Figura 3.1: Ii,l(x) representa a las im´agenes tomadas por la c´amara izquierda y Ii,r(x) a la de la derecha. hi,l(x) y hi,r(x) son los flujos ´opticos de ambas c´amaras. gi(x) es el mapa de disparidad desde la c´amara izquierda a la derecha del par est´ereo i. flujos ´opticos y Nflujos est´ereos. En la figura 3.1, se muestra un gr´afico que resume los distintos campos de desplazamiento estimados: hi,l(x) para la c´amara izquierda, hi,r(x) para la derecha y gi(x) es el mapa de disparidad. Siguiendo esta notaci´on el flujo ´optico se puede expresar como Ii,l(x)≃Ii+1,l (x+hi,l(x)) Ii,r(x)≃Ii+1,r (x+hi,r(x)) para la imagen izquierda y derecha respectivamente. Ii,l(x) es la imagen tomada en el instante ipor la c´amara izquierda y x= (x, y) es la coordenada de un p´ıxel en la imagen. Ii,r(x) es la imagen de la c´amara derecha. Como se coment´o en el apartado 1.2.1 la b´usqueda de las correspondencias entre los p´ıxeles de un par est´ereo se puede simplificar si disponemos de un calibrado de las c´amaras correcto. Apoy´andonos en la geometr´ıa epipolar es posible reducir la zona de b´usqueda a la l´ınea epipolar. En la figura 3.1 el flujo est´ereo se define como gi(x)=(ui,s(x), vi,s(x))t. Para diferenciar los campos de desplazamiento del flujo ´optico de los del est´ereo se ha incluido la etiqueta s. De forma que el flujo est´ereo se expresar´a en funci´on de dos ´ındices (i, s), donde sindica que se trata del mapa de disparidad e i, la posici´on dentro de la secuencia. El flujo est´ereo depende de un escalar λ(x) que expresa el desplazamiento sobre la recta epipolar. En funci´on de esto, el campo de desplazamiento se puede expresar como ui,s(x) = −λi(x)b(x) √a2(x)+b2(x)−a(x)x+b(x)y+c(x) a2(x)+b2(x)a(x) vi,s(x) = λi(x)a(x) √a2(x)+b2(x)−a(x)x+b(x)y+c(x) a2(x)+b2(x)b(x) (3.1) 3.2. ESTIMACI ´ ON DE LA DISPARIDAD USANDO UNA SECUENCIA ESTEREOSC ´ OPICA127 Figura 3.2: Las cuatro vistas del mismo punto 3D est´an interconectadas a trav´es del flujo ´optico y est´ereo de ambas c´amaras. En un caso ideal, si partimos desde el punto situado en la c´amara izquierda en el instante iy seguimos cualquiera de los dos caminos: (1) a trav´es del flujo ´optico (hi,l) y luego con el flujo est´ereo (gi+1) o (2) a trav´es del flujo est´ereo (gi) y posteriormente el flujo ´optico (hi,r), debemos llegar a la correspondencia de la c´amara derecha en el instante i+ 1. Esta restricci´on temporal se puede expresar como hi,l(x) + gi+1(x+hi,l(x)) ≡gi(x) + hi,r(x+gi(x)). En el trabajo de [Alvarez02b] se pueden encontrar m´as detalles acerca de la parametrizaci´on descrita en la ecuaci´on (3.1). Para establecer las correspondencias entre los p´ıxeles en un par est´ereo es necesario escoger alguna propiedad invariante de las im´agenes. La ecuaci´on de restricci´on del flujo est´ereo puede expresarse como Ii,l(x)≃Ii,r (x+gi(x)) (3.2) Si nos fijamos en la figura 3.2 se puede ver f´acilmente que para cada dos frames consecutivos de la secuencia es posible establecer una relaci´on entre las proyecciones de los puntos que pertenecen al mismo punto 3D mediante el flujo ´optico y el est´ereo hi,l(x) + gi+1(x+hi,l(x)) ≡gi(x) + hi,r(x+gi(x)) (3.3) Esta relaci´on se puede expresar matem´aticamente mediante la ecuaci´on (3.3). En una situaci´on ideal si partimos de un punto y nos desplazamos de acuerdo al flujo est´ereo y ´optico, se alcanzar´ıa el mismo destino independientemente si se ha tomado primero el flujo ´optico o est´ereo. 134 CAP´ ITULO 3. ESTIMACI ´ ON DEL MAPA DE DISPARIDAD EN PARES EST ´ EREO Cilindro AEE Frame M´etodo espacial Nuestro m´etodo Diferencia % Mejora 00.127043 0.098521 0.028521 22.45 % 10.129855 0.103104 0.026751 20.60 % 20.130148 0.099989 0.030159 23.17 % 30.135433 0.110151 0.025282 18.67 % 40.122247 0.097061 0.025186 20.60 % 50.131151 0.099486 0.031665 24.14 % 60.122922 0.094766 0.028156 22.91 % 70.129470 0.106723 0.022747 17.57 % 80.126313 0.114762 0.011551 9.14 % 90.148327 0.103394 0.044933 30.29 % 10 0.132380 0.099746 0.032634 24.65 % 11 0.158627 0.107786 0.050841 32.05 % 12 0.150153 0.118870 0.031283 20.83 % Cuadro 3.1: Secuencia del Cilindro: AEE obtenido en los distintos frames con el m´etodo est´ereo espacial [Alvarez02b] y el m´etodo propuesto en esta secci´on. Cilindro AAE Frame M´etodo espacial Nuestro m´etodo Diferencia % Mejora 00.387795 0.294634 0.093161 24.02 % 10.383048 0.305028 0.078020 20.37 % 20.382562 0.291168 0.091394 23.89 % 30.396488 0.389646 0.006842 1.73 % 40.356934 0.285340 0.071594 20.06 % 50.376084 0.289093 0.086991 23.13 % 60.346269 0.272438 0.073831 21.32 % 70.372812 0.319556 0.053256 14.28 % 80.372630 0.415654 -0.043024 -11.55 % 90.606374 0.308366 0.298008 49.15 % 10 0.390006 0.296871 0.093135 23.88 % 11 0.480648 0.313414 0.167234 34.79 % 12 0.442590 0.348788 0.093802 21.19 % Cuadro 3.2: Secuencia del Cilindro: AAE obtenido en los distintos frames con el m´etodo est´ereo espacial [Alvarez02b] y el m´etodo propuesto en esta secci´on. 3.2. ESTIMACI ´ ON DE LA DISPARIDAD USANDO UNA SECUENCIA ESTEREOSC ´ OPICA135 Figura 3.7: Gr´afica del AAE obtenido en la secuencia del Cilindro. el proceso de difusi´on no se ha detenido en los contornos de cilindro pese a su contraste con el fondo de la escena. Debido al uso del gradiente de la imagen en el operador de Nagel– Enkelmann, si en ciertas zonas del contorno del objeto el gradiente no es lo suficientemente alto puede provocar que la regularizaci´on act´ue m´as all´a de los l´ımites del objeto creando una especie de ´aurea alrededor del mismo. En este caso, la utilizaci´on de funciones de robustificaci´on permitir´ıa mejorar las estimaciones. Los modelos de energ´ıa de [Alvarez00] y [Alvarez02b] en los que se basa nuestro m´etodo tienen como caracter´ıstica la estabilidad de sus par´ametros, peque˜nas variaciones en sus valores ofrecen resultados similares. Esta caracter´ıstica facilita las pruebas y la optimizaci´on de los par´ametros no se convierte en una tarea tediosa. Los par´ametros que configuran nuestro m´etodo variacional son: α, β, ζ y el n´umero de escalas. αes el peso que pondera la importancia de los t´erminos de los flujos ´optico y est´ereo por separado y, βes el peso del t´ermino que combina la energ´ıa de ambos flujos. ζestablece el comportamiento isotr´opico o anisotr´opico del tensor de difusi´on de Nagel–Enkelmann. Con valores peque˜nos el tensor de difusi´on tiene un comportamiento anisotr´opico. Los par´ametros αyζse calculan a trav´es de las ecuaciones descritas en el apartado 2.3.1. Los datos de entrada son λyCα. Los valores de estos par´ametros son Cα= 1,5 (α= 1,86 ×10−4), β= 1 ×10−3y λ= 0,1. Aunque los valores de Cαyβparecen muy dispares una vez que se ha estimado el peso del t´ermino de regularizaci´on en funci´on de la normalizaci´on de las im´agenes de entrada obtenemos valores de similar magnitud. Se ha tomado un valor de λ= 0,1 para evitar que el proceso de difusi´on actuara m´as all´a de los contornos del cilindro. Nuestro m´etodo utiliza un enfoque multipiramidal para recuperar los desplazamientos largos y, al mismo tiempo, evitar converger a un m´ınimo local irrelevante. El n´umero de escalas ´optimo depende del desplazamiento m´aximo registrado en los datos. En nuestro caso, ese desplazamiento m´aximo es de alrededor de 8 p´ıxeles por frame. Dado que el factor de reducci´on aplicado en cada escala se ha fijado a 0,5 el n´umero de escalas utilizadas ha 136 CAP´ ITULO 3. ESTIMACI ´ ON DEL MAPA DE DISPARIDAD EN PARES EST ´ EREO sido tres. En los experimentos qued´o demostrado que un n´umero superior de escalas no mejoraba los resultados obtenidos. Los valores de los par´ametros en el m´etodo est´ereo espacial han sido C= 0,7(α= 3,91 ×10−4), λ= 0,1 y el n´umero de escalas de tres. Dado que nuestro m´etodo se basa en ´este es l´ogico que los valores de los par´ametros sean similares. La mejora aportada por el m´etodo presentado en esta secci´on no se puede cuantificar sino establecemos una comparaci´on con otros m´etodos. En los experimentos esta comparaci´on se ha hecho con el m´etodo est´ereo [Alvarez02b]. De esta forma es posible determinar el grado de influencia del t´ermino de flujo ´optico en el nuevo m´etodo. En las tablas 3.1 y 3.2 podemos encontrar los resultados cuantitativos obtenidos de la comparaci´on de ambos m´etodos. Para facilitar el an´alisis de la informaci´on de ambas tablas hemos representado dicha informaci´on de forma gr´afica (figuras 3.6 y 3.7). Se observa de manera general la reducci´on en el error que nuestro m´etodo ofrece. Sin embargo, existen ciertos frames, m´as concretamente el tercero y el octavo, en el que la estimaci´on es ligeramente peor. Las diferencias apreciables en la figura 3.5 se reflejan en los errores detectados en todos los frames de la secuencia. Analizando los resultados de las tablas podemos sacar dos conclusiones: (1) la estabilidad de las estimaciones del m´etodo y, (2) la contribuci´on de la informaci´on del flujo ´optico en la mejora de la estimaci´on del mapa de disparidad. Cilindro y Esfera La segunda secuencia sint´etica utilizada se trata de una versi´on algo m´as compleja que la secuencia del Cilindro (figura 3.8). Ahora disponemos de dos objetos, un cilindro y una esfera, que se desplazan a distinta velocidad por la escena. Se produce oclusiones entre ambos objetos. Cerca del fondo de la escena se aprecia un panel est´atico con una textura clara que ocupa casi la totalidad de la imagen. Para m´as detalles acerca de la secuencia ver el apartado 1.4.2 de la tesis. A diferencia de los experimentos realizados con la secuencia anterior hemos tenido en cuenta todos los p´ıxeles de la imagen a la hora de obtener los resultados cuantitativos. En la figura 3.9 se muestran los mapas de disparidad densos asociados a los pares est´ereo de la figura 3.8. Como se puede observar se ha producido una correcta detecci´on de los tres objetos presentes en la escena. El cilindro es el objeto m´as pr´oximo a las c´amaras. Por este motivo aparece en un tono claro. En el interior del cilindro se aprecia el efecto de degradado presente en los mapas de disparidad de referencia. Sin embargo, en ambas soluciones se observan algunos artificios y zonas de disparidad subestimadas. La esfera tambi´en es detectada correctamente por ambos m´etodos. La diferencia m´as significativa entre ellos est´a en el contorno de la esfera y sobre todo en la discontinuidad con el cilindro, zonas donde se producen oclusiones. El ´ultimo objeto presente en la escena es el panel situado al fondo, justo detr´as del cilindro y la esfera. En el interior de dicho objeto la disparidad ha sido identificada correctamente y apenas existe diferencia entre ambos m´etodos. Sin embargo, en su contorno se aprecia los errores m´as relevantes del mapa de disparidad. En nuestro m´etodo el proceso de difusi´on se detiene con mayor precisi´on en el contorno del panel. Los par´ametros que configuran nuestro m´etodo variacional son: α, β, ζ y el n´umero de 3.2. ESTIMACI ´ ON DE LA DISPARIDAD USANDO UNA SECUENCIA ESTEREOSC ´ OPICA137 Cilindro y Esfera AEE Frame M´etodo espacial Nuestro m´etodo Diferencia % Mejora 00.079917 0.064413 0.015503 19.40 % 10.082677 0.069917 0.012761 15.43 % 20.078151 0.067839 0.010312 13.20 % 30.099163 0.066084 0.033079 33.36 % 40.075023 0.065726 0.009297 12.39 % 50.085858 0.071507 0.014351 16.71 % 60.071050 0.067640 0.003409 4.80 % 70.085979 0.068504 0.017475 20.32 % 80.082747 0.065617 0.017130 20.70 % 90.085170 0.063675 0.021496 25.24 % 10 0.067150 0.062035 0.005115 7.62 % 11 0.103858 0.073978 0.029880 28.77 % 12 0.097868 0.055609 0.042259 43.18 % Cuadro 3.3: Secuencia del Cilindro y la Esfera: AEE obtenido en los distintos frames con el m´etodo est´ereo espacial [Alvarez02b] y el m´etodo propuesto en esta secci´on. escalas. Los valores de entrada son: Cα= 0,03 (α= 9,48 ×10−3) y β= 1 ×10−3. Dado los fuertes contrastes que se producen en los contornos de los objetos nos interesa que el t´ermino de regularizaci´on tenga un comportamiento anisotr´opico. En nuestro m´etodo λ= 0,1. En esta secuencia el desplazamiento m´aximo es de alrededor de 12 p´ıxeles por frame. Dado que el factor de reducci´on aplicado en cada escala se ha fijado a 0,5 el n´umero de escalas utilizadas ha sido cuatro. En los experimentos qued´o demostrado que un n´umero superior de escalas no mejoraba los resultados obtenidos. Los valores de los par´ametros en el m´etodo est´ereo espacial han sido tambi´en C= 0,03 (α= 9,48 ×10−3), λ= 0,1 y cuatro escalas. El aumento en la precisi´on de las estimaciones ofrecida por nuestro m´etodo no se puede cuantificar sino establecemos una comparaci´on con [Alvarez02b]. En las tablas 3.3 y 3.4 podemos encontrar los AEE y AAE obtenidos por ambos m´etodos. Para facilitar el an´alisis de la informaci´on de ambas tablas hemos representado dicha informaci´on de forma gr´afica (figuras 3.10 y 3.11). En estas gr´aficas podemos observar la estabilidad del error en cada uno de los frames y c´omo la combinaci´on de informaci´on del flujo ´optico y est´ereo permite aumentar la precisi´on de las estimaciones. En general nuestro m´etodo ha permitido la obtenci´on de soluciones m´as suaves en el interior de los objetos y una mejor identificaci´on de la disparidad en las discontinuidades de los objetos, zonas ´estas ´ultimas donde los m´etodos fallan debido a las oclusiones. Secuencia de F´utbol La secuencia de f´utbol corresponde a un partido celebrado en el MiniEstadi del F.C. Barcelona. En ella se observa un ´area del campo desde dos c´amaras situadas en la grada (fig. 3.12). El juego se desarrolla en la parte del campo que no es visible; en el ´area se 138 CAP´ ITULO 3. ESTIMACI ´ ON DEL MAPA DE DISPARIDAD EN PARES EST ´ EREO Figura 3.8: Distintos pares est´ereos de la secuencia del Cilindro y la Esfera. En la primera columna, las im´agenes de la c´amara izquierda. En la segunda columna, las im´agenes de la c´amara derecha. Cilindro y Esfera AAE Frame M´etodo espacial Nuestro m´etodo Diferencia % Mejora 04.059050 2.946890 1.112160 27.40 % 14.056370 3.016360 1.040010 25.64 % 24.041220 3.051960 0.989260 24.48 % 34.666410 2.998850 1.667560 35.74 % 43.943700 2.991910 0.951790 24.13 % 54.201020 3.071630 1.129390 26.88 % 63.871540 3.182070 0.689470 17.81 % 73.921320 3.108290 0.813030 20.73 % 84.467930 3.071850 1.396080 31.25 % 94.372420 3.034130 1.338290 30.61 % 10 3.934800 3.116080 0.818720 20.81 % 11 4.750460 3.330940 1.419520 29.88 % 12 4.050140 3.281430 0.768710 18.98 % Cuadro 3.4: Secuencia del Cilindro y la Esfera: AAE obtenido en los distintos frames con el m´etodo est´ereo espacial [Alvarez02b] y el m´etodo propuesto en esta secci´on. 3.2. ESTIMACI ´ ON DE LA DISPARIDAD USANDO UNA SECUENCIA ESTEREOSC ´ OPICA139 Figura 3.9: En la columna izquierda, los mapas de disparidad reales correspondiente a los frames 1, 4 y 10. P´ıxeles con tonos claros indican desplazamientos mayores y, los oscuros, los menores. En la columna central, los mapas de disparidad estimados con el m´etodo espacial. En la columna derecha, los mapas de disparidad obtenidos con nuestro m´etodo. Figura 3.10: Gr´afica del AEE obtenido en la secuencia del Cilindro y la Esfera. 140 CAP´ ITULO 3. ESTIMACI ´ ON DEL MAPA DE DISPARIDAD EN PARES EST ´ EREO Figura 3.11: Gr´afica del AAE obtenido en la secuencia del Cilindro y la Esfera. puede ver al portero caminando desde el punto de penalti hacia el borde del ´area. Dada las dimensiones de las im´agenes (1900 ×1080 p´ıxeles) se ha seleccionado una regi´on m´as peque˜na, de tama˜no 430×170 p´ıxeles, donde se percibe el ´unico movimiento de inter´es de la secuencia, el caminar del portero en el ´area. Desgraciadamente la posici´on y orientaci´on de las c´amaras no es frontoparalela. Por ello, se ha recurrido a un proceso de rectificaci´on de las im´agenes de ambas c´amaras. La secuencia original est´a compuesta por miles de frames. En nuestros experimentos se han seleccionado once frames con el objetivo de verificar, en un tiempo razonable, el correcto funcionamiento del m´etodo en situaciones reales. En esta secuencia no disponemos de los mapas de disparidad reales, pero dada su simplicidad ser´a muy f´acil verificar visualmente la validez de las estimaciones obtenidas. En funci´on de la distancia a las c´amaras los objetos aparecer´an en un tono claro los m´as pr´oximos y oscuros los m´as lejanos. En la figura 3.14 se muestran los mapas de disparidad obtenidos con el m´etodo est´ereo [Alvarez02b] y el nuestro en los frames 0, 5 y 9 de la secuencia. Los dos objetos presentes en las im´agenes son el portero y el punto de penalti. La estimaci´on de la disparidad de ambos objetos se ha hecho correctamente gracias al contraste que se produce con el c´esped. La figura del portero se aprecia con distintas intensidades dependiendo de la distancia respecto a la c´amara de cada parte del cuerpo. Las partes m´as cercanas a la c´amara como son el guante izquierdo y el punto de penalti aparecen en tono cercano al blanco. Las partes del cuerpo m´as lejanas pierna y brazo derecho tonalidades oscuras. El c´esped tiene una tonalidad constante debido a la actuaci´on del t´ermino de suavizado sobre el de ligadura ya que ´este no tiene informaci´on suficiente para establecer las correspondencias entre las im´agenes. Si comparamos los mapas de disparidad estimados por cada m´etodo podemos apreciar que la soluci´on del nuestro es mucho m´as suave. Como era de esperar la utilizaci´on de informaci´on del flujo ´optico aporta estabilidad a la soluci´on. 3.2. ESTIMACI ´ ON DE LA DISPARIDAD USANDO UNA SECUENCIA ESTEREOSC ´ OPICA141 Figura 3.12: Secuencia de f´utbol. En la primera columna, los frames 0, 5 y 9 tomados por la c´amara izquierda. En la segunda columna, los mismos frames pero captados por la c´amara derecha. Figura 3.13: Im´agenes rectificadas modificadas de las c´amaras izquierda y derecha correspondientes a los frames mostrados en la figura 3.12. 142 CAP´ ITULO 3. ESTIMACI ´ ON DEL MAPA DE DISPARIDAD EN PARES EST ´ EREO Figura 3.14: Mapas de disparidad obtenidos con las im´agenes rectificadas presentadas en la figura 3.13. En la columna de la izquierda, los mapas de disparidad calculado con un m´etodo variacional [Alvarez02b]. En la columna derecha, el mapa de disparidad obtenido con nuestro m´etodo. El modelo de energ´ıa de nuestro m´etodo no incorpora ninguna t´ecnica que lo haga insensible a las perturbaciones t´ıpicas de cualquier secuencia real (p. e. ruido o cambios de iluminaci´on). Esto provoca que las estimaciones hechas por los t´erminos est´ereo y del flujo ´optico se vean influenciadas por estas perturbaciones. Durante las pruebas hemos observado que cuando la relaci´on entre el peso del par´ametro Cαyβera inferior a 100 la estabilidad del m´etodo se ve´ıa seriamente comprometida y los resultados obtenidos no eran correctos. Por este motivo, el valores de los par´ametros son Cα= 0,3 (α= 1,61 ×10−4), yβ= 3x10−3. Para valores de Cαmayores que 0,3 el mapa de disparidad era demasiado suave. Dado en las im´agenes el desplazamiento de los objetos es peque˜no solamente se ha utilizado una ´unica escala en el enfoque multipiramidal, con un factor de reducci´on fijado a 0,5. Para evitar crear una sobresegmentaci´on en el flujo y favorecer una soluci´on suave el par´ametro del regularizador de Nagel–Enkelmann es de λ= 0,5. 3.3. MAPA DE DISPARIDAD COMBINANDO M ´ ETODOS DE GRAPH–CUTS Y VARIACIONAL143 3.3. Estimaci´on del Mapa de Disparidad mediante la Combinaci´on de un M´etodo de Graph–cuts y uno Variacional Los m´etodos variacionales han tenido bastante ´exito en el problema de la estimaci´on del mapa de disparidad y han demostrado poder obtener muy buenos resultados. Para evitar la convergencia de la soluci´on a m´ınimos locales se recurre a un enfoque multipiramidal aunque otra alternativa ser´ıa la utilizaci´on de una inicializaci´on. Para obtener esta inicializaci´on se suele recurrir a una t´ecnica basada en la correlaci´on. Una alternativa que ha surgido en los ´ultimos a˜nos a los m´etodos variacionales han sido los m´etodos de graph–cuts. En esencia comparten muchas caracter´ısticas: basados en la minimizaci´on de energ´ıas, resultados densos, etc. Aunque proceden de enfoques matem´aticos diferentes. Dado las similitudes y los buenos resultados que ofrecen parece razonable que la uni´on del m´etodo de graph–cuts y el m´etodo variacional puede ser m´as ventajosa que la aportada por la t´ecnica basada en correlaci´on. Por este motivo, en esta secci´on se propone un m´etodo para la estimaci´on de los mapas de disparidad que combina uno basado en graph–cuts y otro variacional. El m´etodo de graph–cuts ha sido propuesto por Kolmogorov–Zabih [Kolmogorov01, Boykov04]. Ofrece muy buenos resultados en muchas aplicaciones como pueden ser la segmentaci´on o la estimaci´on de los mapas de disparidad. Es capaz de detectar y manejar las oclusiones pero las soluciones tienen precisi´on entera. El m´etodo variacional es el descrito en el art´ıculo Alvarez/Deriche/S´anchez/Weickert (2002) [Alvarez02b] y se trata de la versi´on est´ereo del m´etodo de flujo ´optico Alvarez/Weickert/S´anchez (2000) [Alvarez00]. El modelo de energ´ıa descrito en Alvarez/Deriche/S´anchez/Weickert (2002) [Alvarez02b] no tiene en cuenta las oclusiones. A continuaci´on se describe con detalle el trabajo desarrollado en esta secci´on. En primer lugar, se define la estructura de la escena y el proceso de rectificaci´on de im´agenes que se aplica a las secuencias cuando la disposici´on de las c´amaras no es frontoparalela. En segundo lugar, se comenta brevemente las t´ecnicas utilizadas como aproximaci´on inicial del m´etodo variacional: en 3.3.1, una implementaci´on b´asica de la correlaci´on a ventanas utilizada en [Alvarez02b] y, en el apartado 3.3.2, el m´etodo de graph–cuts. Una vez comentadas estas dos t´ecnicas se explica c´omo se ha combinado el m´etodo variacional con el del graph-cuts y las similitudes que existen entre ellos, adem´as del enfoque multiescala para la detecci´on de los desplazamientos. En el apartado 3.3.5 se muestra los resultados experimentales obtenidos con el m´etodo propuesto y se establecer´a una comparaci´on entre el m´etodo variacional con las dos aproximaciones iniciales. Con esta comparaci´on se pretende evaluar la influencia de cada una de las inicializaciones en la soluci´on final. Para ello, se utilizar´an pares est´ereos reales y sint´eticos ofrecidos por la base de datos Middlebury (versi´on 2001). Se concluir´a con las observaciones m´as importantes de la combinaci´on de estas dos t´ecnicas.