Métodos multimalla en problemas lineales de flujo óptico

XX Congreso de Ecuaciones Diferenciales y Aplicaciones X Congreso de Matem´ atica Aplicada Sevilla, 24-28 septiembre 2007 (pp. 1—8) M´etodos multimalla en problemas lineales de flujo ´optico G.Asensio1,2, P.M.Gonz´ alez1,2,C.Platero 1,3, J.M.Poncela 1,2,J.Sanguino 1,2,M.C.Tobar1,2 1Applied BioEngineering Group (ABE-UPM). 2Dpto. Matem´atica Aplicada EUITI, Universidad Polit´ecnica de Madrid, C/Ronda de Valencia 3, 28012, Madrid. E-mail: [email protected]. 3Dpto. Electr´onica, Autom´atica e Inform´atica Industrial,Universidad Polit´ecnica de Madrid, C/Ronda de Valencia 3, 28012, Madrid. E-mail: [email protected]. Palabras clave: Soluci´on num´erica de Ecuaciones en Derivadas Parciales , Flujo ´ Optico Resumen Los m´etodos variacionales son especialmente ´utiles para resolver la ecuaci´on del flujo ´optico OFC (optical flow constraint) u(x)·∇f(t, x)+∂tf(t, x)=0, porque preservan las discontinuidades y funcionan incluso cuando hay variaciones de iluminaci´on. Sin embargo, las distintas regularizaciones implican funcionales cuya minimizaci´on tiene un alto costo computacional. En esta comunicaci´on se formula el problema lineal, se discretiza y se resuelve ´este usando distintos m´etodos multimalla. Se comparan los resultados con los m´etodos iterativos cl´asicos y con los m´etodos directos con almacenamiento disperso. Se comprueba la mayor eficiencia de los m´etodos multimalla. 1. Introducci´on Uno de los problemas m´as interesantes de la visi´on por computador es obtener informaci´on del movimiento de objetos a partir de im´agenes de v´ıdeo. El movimiento es representado por una familia uniparam´etrica de funciones vectoriales fdefinida sobre Ω×[0,T]⊂R2×R, apartirdelacualsedefine una familia de funciones escalares I(−→ x,t)=kf(−→ x,t)k,que recoge la variaci´on de la intensidad (brillo). Sabemos que el movimiento de patrones de intensidad puede tener una relaci´on con el movimiento de objetos en la escena que recoge la secuencia de im´agenes, pero esta relaci´on 1 G. Asensio, C.Platero no est´a del todo clara ya que la percepci´on del movimiento no se comprende todav´ıatan bien como la percepci´on de patrones est´aticos. Basta, por ejemplo pensar en los antiguos anuncios de las barber´ıas, en lo que ocurre al enfocar un objeto fijo cuando viajamos en coche y en la proyecci´on de una h´elice rotando, para ver que la relaci´on entre movimiento aparente y movimiento real es muy compleja [7].Estostresejemplossoncasosdeflujo ´optico, paralaje del movimiento y profundidad cin´etica respectivamente, en los que movimiento aparente y real no coinciden. Al estudiar el movimiento de patrones en el brillo, adem´as de los problemas t´ıpicos de baja resoluci´on o contraste, gradaci´on de sombras debido a la iluminaci´on variable, falta de n´ıtidez en los bordes, agujeros, etc., nos encontramos con problemas de diferentes rangos de movimientos y movimientos esp´ureos de sombras y ruidos [6]. En cualquier caso, en una secuencia t´ıpica se aprecia un fondo est´atico y unos patrones de intensidad que se mueven a distintas velocidades con movimientos lineales, err´aticos, etc. Por eso, el primer paso es segmentar el fondo est´atico, usando por ejemplo la t´ecnica b´asica de diferencias acumulativas, y despu´es, mediante operadores morfol´ogicos, se detectan los patrones din´amicos a estudiar1.Por´ultimo se siguen estos patrones mediante t´ecnicas estad´ısticas (correlaci´on,..etc.)2. M´as modernas son las t´ecnicas basadas en ecuaciones en derivadas parciales [1], como la ecuaci´on OFC, formulada por primera vez por Horn y Schunck en 1981 [4]. La resoluci´on num´ericahaevolucionadodesdelosm´etodos indirectos basados en el gradiente [5] hasta los m´etodos multimalla [2]. En esta comunicaci´on se ponen de manifiesto las ventajas de estos ´ultimos. 2. Ecuacion OFC. Restricci´on y regularizaci´on Una simplificaci´on de los problemas citados anteriormente es el llamado problema del flujo ´optico OFC, consistente en encontrar una distribuci´on de velocidades aparentes del movimiento de patrones (flujo ´optico) −→ v=(u, v)definido sobre Ω×[0,T],apartirdedos im´agenes consecutivas. Para garantizar que existe una forma de relacionar flujo ´optico con el movimiento de los objetos de la imagen hay que hacer una serie de hip´otesis: 1) Los objetos de la escena son planos y est´an iluminados uniformemente. De esta forma nos aseguramos que el brillo es proporcional a la reflectancia del objeto y no existen movimientos de sombras. 2) La reflectancia de los objetos var´ıa suavemente y no hay solapamiento de objetos. Esto garantiza la diferenciabilidad del brillo. 3) El flujo ´optico se debe al movimiento relativo de objetos o al movimiento del observador. Esto hace que los patrones de brillo se muevan como s´olidos r´ıgidos o deformables. 4) El flujo ´optico se define suponiendo que en cada imagen de la secuencia se conservan algunas propiedades. Asumimos que el que el nivel de intensidad se conserva en cada imagen, es decir, suponemos que I(−→ x,t)−I(−→ x+∆−→ x,t+∆t)=0.Esto significa que seguimos la intensidad la intensidad en un p´ıxelcomosideunpatr´on trivial se tratase. Basta entonces aplicar un desarrollo de Taylor, dividir por ∆t, y considerar el l´ımite 1www.dai.ed.ac.uk/homes/s0129556/softwareGLOMO/doc.pdf 2www.cs.ubc.ca/˜deaton/remarks ncc.html 2 M´etodos multimalla en problemas lineales de flujo ´optico ∆t→0, para obtener la ecuaci´on OFC y plantear el siguiente problema: Conocido I(−→ x,t), hallar −→ v(−→ x,t)tal que −→ ∇I(−→ x,t)·−→ v(−→ x,t)+∂I(−→ x,t) ∂t =0.(1) Como es bien sabido, las curvas de nivel del brillo son perpendiculares al gradiente, y por cada pixel s´olo pasa una curva de nivel. Sin embargo, las curvas de nivel de flujo ´optico forman con el gradiente un ´angulo que depende de la variaci´on temporal del brillo ylaecuaci´on OFC permite asociar un conjunto infinito de curvas soluci´on a cada p´ıxel de la imagen inicial. En efecto, cualquier vector que se apoye sobre la recta rperpendicular al gradiente −→ ∇I(−→ x,t) y que diste del origen de coordenadas el cociente entre la derivada temporal y el m´odulo del gradiente de la intensidad, es soluci´on de la ecuaci´on OFC. Estos vectores son vectores autoconjugados respecto de la forma bilineal asociada a la forma cuadr´atica semidefinida positiva definida en el espacio proyectivo (u, v, 1), con matriz asociada J= ⎛ ⎜ ⎜ ⎝¡∂I ∂x¢2∂I ∂x ∂I ∂y ∂I ∂x ∂I ∂t ∂I ∂x ∂I ∂y ³∂I ∂y´2∂I ∂y ∂I ∂t ∂I ∂x ∂I ∂t ∂I ∂y ∂I ∂t ¡∂I ∂t ¢2 ⎞ ⎟ ⎟ ⎠.(2) Existen muchas formas de restringir las posibles curvas soluci´on del problema (apertura). Una primera es a˜nadir una nueva condici´on, por ejemplo que el gradiente tambi´en se conserve y resolverla junto con la ecuaci´on OFC. Pero este m´etodo utiliza derivadas de segundo orden y es muy sensible al ruido. Otra posibilidad es resolver la ecuaci´on en el sentido de m´ınimos cuadrado. Este m´etodo da buenos resultados pero no garantiza que la soluci´on tenga regularidad global. Una forma adecuada de dar apertura y al mismo tiempo regularidad (u, v ∈H2(Ω)) al problema OFC, es reformularlo como un problema de minimizaci´on de un funcional Ψ(u, v) que incluya un t´ermino de fidelizaci´on F(−→ v)yunt´ermino de penalizaci´on, regularizaci´on o suavizado S(−→ v). Por ejemplo, m´ınΨ(u, v) donde Ψ(u, v)=ZΩµk(u, v, 1)k2 J+αµ° ° °−→ ∇u° ° ° 2 +° ° °−→ ∇v° ° ° 2¶¶dx =F(−→ v)+S(−→ v),(3) que garantiza una soluci´on −→ vdel problema OFC con m´odulo del gradiente vectorial m´ınimo. No se puede probar el caracter convexo del funcional Ψ(u, v) por el car´acter semidefinido de J, y por esta misma raz´on es un exceso de notaci´on escribir k·kJ.Elt´ermino S(−→ v) se escoge de acuerdo al tipo de secuencia; por ejemplo, si proviene de una traslaci´on del observador o de una rotaci´on de objetos alrededor de un eje con campo de velocidades de laplaciano nulo, lo adecuado es penalizar el laplaciano. En lo que sigue nos centraremos en el funcional de Horn (3) , a pesar de que tiene un efecto suavizador isotr´opico sobre el gradiente de las componentes del flujo que no es recomendable. Otras propuestas son regularizar relajando el efecto suavizador en p´ıxeles donde el m´odulo del gradiente espacial de la imagen sea mayor, hacer que la regularizaci´on sea anisotr´opicaparaevitarsuavizarelflujo a lo largo de las discontinuidades de la imagen, o incluso hacer que la regularizaci´on est´e conducida por flujo en vez de por el gradiente para conseguir suavizar los bordes del flujo. 3 G. Asensio, C.Platero 3. Ecuaciones de Euler-Lagrange Las ecuaciones de Euler-Lagrange del funcional de Horn son J11u+J12v+J13 −α4u=J12u+J22v+J23 −α4v=0.(4) Estas ecuaciones lineales se discretizan mediante un esquema de diferencias finitas, a partir de una partici´on del dominio Ωen m×nsubrect´angulos (p´ıxeles) y del intervalo temporal en kpasos de tiempo. Si se sustituye el laplaciano por una de sus aproximaciones, por ejemplo 4uijk =4(uijk −uijk), uijk =1 4(ui+1,jk +ui−1,jk +ui,j+1,k +ui,j−1,k),(5) y las componentes del tensor de movimiento se aproximan promediando los valores en un cubo del espacio-tiempo, por ejemplo µ∂I ∂y¶ijk =1 4(Ii,j+1,k−Ii,j,k+Ii+1,j+1,k−Ii+1,j,k+Ii,j+1,k+1−Ii,j,k+1+Ii+1,j+1,k+1−Ii+1,j,k+1), (6) se obtienen las siguientes ecuaciones en diferencias (dos por cada p´ıxel): (J11)ijuij +(J12)ijvij =α(ui+1,j +ui−1,j +ui,j+1 +ui,j−1−4uij)−(J13)ij (J12)ijuij +(J22)ijvij =α(vi+1,j +vi−1,j +vi,j+1 +vi,j−1−4vij)−(J23)ij .(7) Estas ecuaciones junto con las condiciones de contorno nulo uij =vij =0sii/∈{1,...,n} ´o j/∈{1,...,m}definenunsistemaAα h −→ uh=−→ Fhde h=2mn ecuaciones y otras tantas inc´ognitas, siendo la matriz de coeficientes Aα h=µJ11 −α£J12 J12 J22 −᣶∈Mh,£= ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ e £nInOn··· Ine £nIn ... OnIn ...... . . ....... ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ∈Mm×n, e £n= ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ −41 0··· 1−41... 01 ...... . . ....... ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ∈Mn,Jrs =⎛ ⎜ ⎝ (Jrs)11 00 0...0 00(Jrs)mn ⎞ ⎟ ⎠∈Mm×n, (8) y los vectores de Rh −→ uh=(u11,..,u 1n,u 21,.., u2n,..,u m1,...,u mn,v 11,..,v 1n,v 21,.., v2n,..,v m1,...,v mn), −→ Fh=¡(J13)11 ··· (J13)mn (J23)11 ··· (J23)mn ¢(9) El sistema a resolver tiene como matriz de coeficientes una matriz sim´etrica definida positiva de gran talla, pero de una gran dispersi´on (para una resoluci´on de 128×128 s´olo el 0.0182 % de elementos son no nulos), lo que hace que los m´etodos directos e indirectos funcionen muy r´apido, incluso sin reducci´on del ancho de banda [3]. Sin embargo, tienen el gran problema de que una resoluci´on 256 ×256 implica una matriz de 131072 ×131072, muy dif´ıcil de generar y almacenar. 4 M´etodos multimalla en problemas lineales de flujo ´optico 4. M´etodos cl´asicos de resoluci´on En el art´ıculo original de Horn y Schur se propone un m´etodo indirecto consistente en manipular las ecuaciones de Euler-Lagrange discretizadas (7), obtener un sistema de ecuaciones lineales cuya matriz de coeficientes se invierte facilmente por las propiedades de la matriz Jy llegar a µuij −uij vij −vij ¶=−α ¯¯¯e Jα¯¯¯µ(J11)ijuij +(J12)ijvij +(J13)ij (J22)ijvij +(J12)ijuij +(J23)ij ¶,¯¯¯e Jα¯¯¯=α(α+(J11)ij+(J22)ij). (10) El valor del campo en un pixel depende del promedio del valor en los p´ıxeles del entorno. Es posible obtener dicho valor iterativamente: Ãuk+1 ij vk+1 ij !=µuk ij vk ij ¶−α ¯¯¯e Jα¯¯¯µ(J11)ijuk ij +(J12)ijvk ij +(J13)ij (J22)ijvk ij +(J12)ijuk ij +(J23)ij ¶.(11) En las zonas donde el gradiente de la intensidad es nulo, el campo de velocidades es la media de los valores del entorno de vecindad. En estas zonas el gradiente de velocidades no se puede utilizar para restringir el conjunto de soluciones de la ecuaci´on OFC y, de igual modo que ocurre con la distribuci´on de temperaturas en una placa uniforme, la velocidad aparente en el interior depender´adelosvaloresdelamismaenlosp´ıxeles que rodean la zona. En las zonas donde el gradiente de la intensidad es lineal, ocurre algo similar. La restricci´on de soluciones de la ecuaci´on OFC s´olo ser´a efectiva en la direcci´on perpendicular al gradiente de la intensidad. En definitiva, la soluci´on restringida de la ecuaci´on OFC s´olo ser´asignificativa donde el gradiente del brillo var´ıesuficientemente, sin que esta fluctuaci´onseaexcesivaloquedar´ıa origen a pixelaciones indeseables debido a los errores en las estimaciones de las derivadas. Los m´etodos directos, variantes de la eliminaci´on de Gauss, tienen el inconveniente adicional de que acceden a cada uno de los coeficientes de la matriz Aα h.Sin embargo son m´etodos m´as r´apidos, especialmente si se hace previamente una reordenaci´on de filas y columnas para reducir el ancho de banda (algoritmo de Cuthill-McKee) o para reducir el grado de conexi´on del grafo asociado a la matriz (algoritmo de m´ınimo grado). Los m´etodos indirectos, como el de Jacobi, obtienen la soluci´on del sistema sin tener que acceder a los coeficientes para invertir la matriz Aα h. Para ello, parten de la descomposici´on de la matriz de coeficientes en una suma de una matriz diagonal Dhy dos triangulares Uh yLh, obteniendo la siguiente expresi´on: −→ uh=Kh−→ uh+−→ chsiendo Kh=I−D−1 hAα hy−→ ch=D−1 h −→ Fh,(12) lo que justifica la resoluci´on mediante el esquema iterativo −→ uk+1 h=Kh−→ uk h+−→ ch,k= 0,1,...,r . Sin embargo, los algoritmos m´as utilizados son versiones amortiguadas (damped) que incluyen un par´ametro de control 0 <ω≤1: Kω h=I−ωD−1 hAh,−→ cω h=ωD−1 h −→ Fh.(13) 5 G. Asensio, C.Platero Las propiedades del espectro de Kω hgarantizan que estos algoritmos sean consistentes y convergentes, pero sobre todo garantizan la propiedad de ”suavizado”del error en cada iteraci´on. Esto se comprueba f´acimente al discretizar una versi´on muy simplificada de la ecuaci´on OFC con J=0,α=1yn=1.Enestecasoelespectrodelamatrizdeiteraci´on es σ(Kω h)=©µω i=1−ω¡1−cos iπ h¢ªh i=1 <1. La convergencia est´a garantizada, aunque se puede deteriorar si disminuimos el n´umero de p´ıxeles. Esto se debe a que el radio espectral de la matriz de iteraci´on tiende a uno cuando htiende a cero. El error se puede expresar como combinaci´on lineal de los autovectores (−→ vi)j= sin(ijπ/h)queseclasifican en, autovectores de baja frecuencia tales que 1 ≤i≤h/2 y autovectores de alta frecuencia con h/2≤i≤h. El par´ametro de control ωse escoge tal que el factor de suavizado, mayor autovalor de asociado a las altas frecuencias, sea lo m´as peque˜no posible. Con esta elecci´on, se consigue que en unas pocas iteraciones las altas frecuencias desaparezcan y el error se convierta en un error ”suave”, combinaci´on lineal de autovectores de bajas frecuencias. En el caso que nos ocupa, el factor de suavizado alcanza el valor m´ınimo µ∗=1/3si ω∗'2/3,independientemente del valor de h. Si escogemos el m´etodo de Jacobi amortiguado con ω=2/3,garantizamos que el error disminuya m´as lentamente en cada iteraci´on, pero se suaviza (se eliminan las altas frecuencias) m´as r´apidamente, pudiendo adem´as disminuir o aumentar la resoluci´on sin que esto afecte a la convergencia. 5. M´etodos multimalla La principal ventaja de los m´etodos indirectos del tipo del m´etodo de Jacobi, es que son algoritmos vectorizables. Adem´as, si una vez obtenida una soluci´on para el flujo asociado alasim´agenes tyt+1,esta soluci´on se usa como soluci´on inicial para obtener el flujo asociado a las im´agenes t+1yt+2,la la eficiencia del m´etodo aumenta. El mayor inconveniente es la imposibilidad de eliminar del error las bajas frecuencias asociadas a autovectores −→ vicon 1 ≤i≤h/2. La idea fundamental de los m´etodos multimalla es hacer hpeque˜no (dowsamplig, disminuir el n´umero de p´ıxeles) para eliminar mayor cantidad de bajas frecuencias, a´un a costa de que aparezcan frecuencias altas indeseadas (aliasing, pixelado), que luego pueden ser filtradas con un nuevo suavizado. Aunque se sabe que la fundamentaci´on te´orica est´a poco desarrollada, los m´etodos multimalla funcionan bien para problemas de PDE’s el´ıpticas (lineales o no) con contorno suave, necesitan menos memoria porque no almacenan la matriz Aα h,son f´acilmente vectorizables (aunque requieren m´as ´ordenes de comunicaci´on que otros solvers), precisan una cuidadosa “sintonizaci´on”, elecci´on del tipo de ciclos, de algoritmo suavizador, etc. El algoritmo de resoluci´on multimalla (2-malla) m´as sencillo, consiste en los siguientes tres pasos. Primero se aplican riteraciones del m´etodo Jacobi amortiguado para obtener una soluci´on aproximada −→ ur h,unerror−→ er h=−→ uh−−→ ur hy un residual−→ rr h=−→ Fh−Aα h −→ ur h, tales que Aα h −→ er h=Aα h¡−→ uh−−→ ur h¢=−→ rr h. En segundo lugar, el residual se restringe Ih→2h(−→ rr h)=−→ rr 2hyseresuelveelsistema lineal Aα 2h −→ er 2h=−→ rr 2h. Por ´ultimo y en tercer lugar, se interpola el residual I2h→h(−→ er 2h)=−→ er hy se corrige 6 M´etodos multimalla en problemas lineales de flujo ´optico la soluci´on sustituyendo −→ ur hpor −→ ur h+−→ er h. El algoritmo se representa mediante el grafo a) de la figura siguiente, donde el c´ırculo azul representa iteraciones de Jacobi (suavizado), el c´ırculorojorepresentalaobtenci´on del error por el m´etodo directo y el c´ırculo verde la correcci´on de la soluci´on. La flecha descendente representa la restricci´on del residual y la flecha ascendente la interpolaci´on del error. Apartirdeladefinici´on anterior podemos definir diferentes algoritmos n-malla de tipo V y W, representados respectivamente por los grafos b) y c) de la figura 1. En el gr´afico del algoritmo tipo W aparece un c´ırculo de dos colores verde y azul que representa la correcci´on de la soluci´on y el posterior suavizado de la misma. Pero los m´as utilizados son los llamados m´etodos multigrid total, como el que se representa por el grafo d) de la figura 1 y que es el utilizado en los experimentosqueseproponeenlasiguientesecci´on.La ventaja es que partiendo de una soluci´on a baja resoluci´on, con un juego de interpolaci´onrestricci´on, se llega a una soluci´onaaltaresoluci´on, evitando resolver el sistema salvo a baja resoluci´on y por un m´etodo directo. Fig. 1 Representaci´on de los m´etodos multimalla 6. Experimentos num´ericos Los m´etodos de resoluci´on anteriormente descritos se han implementado en c´odigo Matlab que se ha ejecutado en un Pc Pentium IV 3.2 GHz. En el experimento 1 se estudia el conocido test de la secuencia del parque Yosemite3, tanto a alta 256×256 como abajaresoluci´on 32×32. En esta secuencia s´olo se mueve el observador. Se puede ver que, a pesar del suavizado de las l´ıneas de flujo,abajaresoluci´on tambi´en se aprecia el flujo ´optico. En el experimento 2, la secuencia recoge el movimiento de una esfera. En este caso, a baja resoluci´on el flujo ´optico pierde informaci´on acerca del eje de giro de la esfera, s´ıapreciable en el caso de alta resoluci´on. Fig.2Experimento1aaltayabajaresoluci´on 3http://www.informatik.uni-ulm.de/ni/staff/PBayerl/homepage/ECVP2003/exampleYos.html: 7 G. Asensio, C.Platero Fig.3Experimento2aaltayabajaresoluci´on Tabla 1 Tiempos de c´alculoparaelExperimento2aaltaresoluci´on 7. Conclusiones Los tiempos de c´alculo recogidos en la Tabla 1 ponen de manifiesto la mayor eficiencia del m´etodo multimalla. Aunque los c´odigos programados no est´an optimizados se obtienen tiempos de c´alculo que permitir´ıan analizar en tiempo real secuencias de 21 im´agenes por segundo, con una resoluci´on muy aceptable de 256×256 p´ıxeles. Esta t´ecnica multimalla de c´alculo del flujo ´optico ser´a utilizada para analizar secuencias de im´agenes obtenidas por tomograf´ıa. Referencias [1] G. Aubert, P. Kornprobst, Mathematical Problems in Image Processing, Partial Diffential Equations and the Calculus of Variations, sec. 5.3, Springer 2006. [2] A.Bruhm,J.Weikert,T.Kohlbeger,C.Sch¨orr, A Multigrid Platform for Real-Time Motion Computation with Discontinuity-Preserving Variationaal Methods, preprint 136, Universtit¨at Saarlandes, Fachrichtung 6.1-Mathematik, Saarbr¨ucken 2005. [3] I.S. Duff,A.M.Erisman,J.K.Reid,Direct Methods for Sparse Matrices, Monographs on Numerical Analysis, Oxford Science Publications 1989. [4] B. Horn, B. Schunck, Determining Optical Flow, Artificial Inteligence 17 (1981) 185-203, North- Holland. [5] D. Robinson, P. Milanfar, Fast Local and Global Projection-Based Methods for Affine Motion Estimation, Kluwer Academic Publishers, Journal of Mathematical Imaging and Vision 18:35-54, 2003. [6] M.K.Titsias,Unsupervised Learning of Multiple Objects in Images ,PhDthesis,SchoolofInformatics, University of Edinburgh, 2005. [7] C. Ware, Information Visualization. Perception for design, Morgan Kaufmann series in interactive technologies, 1999. 8