Sobre un problema de optimización no-lineal en visión estereoscópica

Sobre un problema de optimizacion no-lineal en vision estereoscopica. Asp ectos Computacionales Luis Alvarez Leon 1 Javier Sanchez Perez 1 Resumen Presentamos un meto do no lineal para la estimacion de la geometra 3-D de una escena a partir de 2imagenes esteroscopicas. El problema principal consiste en calcular la p osicion relativadelas2camaras apartirdeun n umero de puntos que se corresp onden en ambas camaras. La p osicion relativadelas2 camaras viene dada p or un vector de 7parametros: { X =( s l  m n t x t y t z ){. Para calcular estos parametros hay que minimizar una energa no-lineal del tip o E ( X )= k Aq ( X ) j 2 donde A es una matriz 9 x 9y q ( X )esunvector funcion de X . En este traba jo presentamos un algoritmo para la busqueda de mnimos lo cales de E ( X ) basado en una mo dicacion del meto do de gradiente de paso optimo. Presentamos algunas exp eriencias comparativas con otros meto dos clasicos. Intro duccion En los ultimos a ~nos se han investigado diferentes tecnicas que p ermiten determinar la estructura 3 ; D a partir de dos imagenes proyectivas. Para sentar las bases del problema, presentamos en la gura 1 un mo delo proyectivo utilizado normalmente en traba jos sobre imagenes en p ersp ectiva (Ver  ? ], ? ]y  ? ]). Consideramos un conjunto de puntos 3 ; D que se proyectan en cada camara. Para cada camara se utiliza un sistema de co ordenadas distinto. Denotamos por ( x y  z )las co ordenadas 3 ; D de un punto ypor ( u v )las co ordenadas de la imagen del punto proyectado en una camara, y denotamos p or ( x 0 y 0 z 0 ) las co ordenadas 3 ; D de un punto y p or ( u v ) las co ordenadas de la imagen del punto en la otra camara. Consideramos que se cono cen los parametros intrnsecos de las dos camaras. Esto signica, en particular, que el sistema de co ordenadas de la imagen puede estar normalizado de tal manera que el origen de cada camara esta lo calizado en el punto de la imagen corresp ondiente a la intersecccion del punto fo cal con el plano de la imagen el punto fo cal esta en direccion del eje z (eje z 0 en la otra camara), y los vectores unitarios ^ u ^ v (^ u 0  ^ v 0 resp.) estan alineados y con la misma magnitud que los vectores unitarios ^ x ^ y (^ x 0  ^ y 0 resp.). En el sistema de referencia 3 ; D , p o demos asumir tambien, que se cono cen las distancias entre los fo cos y los planos de las imagenes, (denotamos por D y D 0 estas distancias). Con esta normalizacion, obtenemos que el sistema de co ordenadas de la imagen de un punto 3 ; D ( x y  z )proyectado en la imagen, viene dado por ( u v )=( D x=z  D y =z ) ( ( u 0 v 0 )=( D 0 x 0 =z 0 D 0 y 0 =z 0 ) resp ect. ). Figura 1: Mo delo Proyectivo. Consideramos un conjunto de puntos 3 ; D que denotamos por ( x i y i z i ) en un sistema de referencia 3 ; D ypor ( x 0 i y 0 i z 0 i ) en el otro sistema de referencia 3 ; D . La transformacion entre los dos sistemas de referencia 3 ; D esta denido por una traslacion y una rotacion rgida, y se puede expresar como 0 @ x 0 i y 0 i z 0 i 1 A = 0 @ r 11 r 12 r 13 r 21 r 22 r 23 r 31 r 32 r 33 1 A 0 @ x i y i z i ; t x t y t z 1 A (1) donde los r ik son elementos de una matriz de rotacion R y el vector t = ( t x t y t z ) representa la traslacion del primer fo co al segundo. Asumimos que para cualquier punto 3 ; D z i z 0 i > 0 : Asignamos (  i  i ) = ( u i =D  v i =D ) y (  0 i  0 i ) = ( u 0 i =D 0 v 0 i =D 0 ) : Utilizando la ecuacion anterior, obtenemos 4 ecuaciones que envuelven las co ordenadas escaladas de la imagen (  i  i ) y (  0 i  0 i ) : Resolviendo estas ecuaciones se obtiene una sola ecuacion que se puede expresar como (  0 i  0 i  1) Q 0 @  i  i 1 1 A =0 (2) donde Q se puede expresar como 0 @ r 13 t y ; r 12 t z r 11 t z ; r 13 t x r 12 t x ; r 11 t y r 23 t y ; r 22 t z r 21 t z ; r 23 t x r 22 t x ; r 21 t y r 33 t y ; r 32 t z r 31 t z ; r 33 t x r 32 t x ; r 31 t y 1 A (3) Por lo tanto, para cada par de puntos en corresp ondencia (  i  i )y(  0 i  0 i ) obtenemos una ecuacion dada por (2). Con N puntos lo calizados en ambas imagenes, el sistema de ecuaciones resultante se puede escribir como A 0 B B B B @ q 11 q 12 : : q 33 1 C C C C A = 0 B B B B @ 0 0 0 0 0 1 C C C C A (4) donde A es una matriz Nx 9. La ecuacion anterior es homogenea y, por lo tanto, su solucion solo se puede calcular en funcion de una constante multiplicativa indeterminada. Cuando N  8 y los puntos 3 ; D no estan en alguna conguracion geom etrica esp ecial, la solucion de la ecuacion (4) se puede encontrar minimizando la siguiente energa E ( q )= k Aq k 2 = q T A T Aq (5) con la condicion k q k = 1  donde q = ( q 11 q 12  ::::::: q 31 ) es un vector 9 x 1 con los elementos de la matriz Q: Ya se sab e que la solucion del problema de minimizacion anterior viene dado por el autovector aso ciado al autovalor mas peque~no de la matriz B = A T A: Una vez que se ha obtenido Q , se utiliza alguna tecnica estandar para calcular, a partir de Q la matriz de rotacion R y el vector de traslacion t (en funcion de un parametro de escala). Por ejemplo, se puede obtener t como el autovector aso ciado al autovalor mas p eque ~no de la matriz Q T Q y la matriz de rotacion R se puede calcular utilizando la siguiente expresion r 1 = q 1  t + q 2  q 3 k t k 2 r 2 = q 2  t + q 3  q 1 k t k 2 (6) r 3 = q 3  t + q 1  q 2 k t k 2 donde r i =( r i 1 r i 2 r i 3 )y q i =( q i 1 q i 2 q i 2 ) : El principal problema con esta aproximacion es que cuando existe ruido en las co ordenadas de la imagen, la expresion anterior no genera, normalmente, una matriz de rotacion. Esto signica que R no es una matriz ortonormal. En tal caso, se pueden realizar algunas op eraciones adicionales para transformar R en una matriz de rotacion real. Estas op eraciones son de tip o algebraico y no tienen en cuenta la precision en la ecuacion que p one en corresp ondencia los puntos (2). En  ? ]los autores prop onen un meto do que utiliza algunas tecnicas de geometra algebraica que proveen una caracterizacion de las matrices fundamentales, forzando la restriccion de rigidez. En este artculo, presentamos una nuevaaproximacion al problema de la recup eracion de la matriz de rotacion R y del vector de traslacion t: Utilizaremos cuaterniones para representar una matriz de rotacion R se cono ce (ver, por ejemplo,  ? ]) que, al utilizar cuaterniones, la matriz de rotacion R se puede escribir como 0 @ s 2 + l 2 ; m 2 ; n 2 2( lm ; sn ) 2( nl + sm ) 2( lm + sn ) s 2 ; l 2 + m 2 ; n 2 2( mn ; sl ) 2( nl ; sm )2( mn + sl ) s 2 ; l 2 ; m 2 + n 2 1 A donde s 2 + l 2 + n 2 + m 2 = 1 : , el vector ( l m n ) representa el ej e de rotacion y s =cos ( = 2),  representa el angulo de rotacion. Usando la ecuacion anterior y (3), deducimos que el vector q se puede expresar como q ( X )= 0 B B B B B B B B B B B B @ t y (2 nl +2 sm ) ; t z (2 lm ; 2 sn ) t z ( s 2 + l 2 ; m 2 ; n 2 ) ; t x (2 nl +2 sm ) t x (2 lm ; 2 sn ) ; t y ( s 2 + l 2 ; m 2 ; n 2 ) t y (2 mn ; 2 sl ) ; t z ( s 2 ; l 2 + m 2 ; n 2 ) t z (2 lm +2 sn ) ; t x (2 mn ; 2 sl ) t x ( s 2 ; l 2 + m 2 ; n 2 ) ; t y (2 lm +2 sn ) t y ( s 2 ; l 2 ; m 2 + n 2 ) ; t z (2 mn +2 sl ) t z (2 nl ; 2 sm ) ; t x ( s 2 ; l 2 ; m 2 + n 2 ) t x (2 mn +2 sl ) ; t y (2 nl ; 2 sm ) 1 C C C C C C C C C C C C A donde X =( s l  m n t x t y t z ) : Utilizando esta formulacion, reescribimos el problema de optimizacion de energa (5) en funcion de X E ( X )= k Aq ( X ) k 2 = T q ( X ) T Bq ( X ) (7) donde B = A T A con las restricciones s 2 + l 2 + n 2 + m 2 =1 y t 2 x + t 2 y + t 2 z = C 2 : (donde C representa la distancia entre los fo cos de ambas camaras, en el caso en que no se conozca C , jamos C = 1 y obtenemos sus co ordenadas 3 ; D en funcion de un factor de escala). Notese que con esta formulacion las variable son s l  m n t x t y t z , as que tenemos 7 variables en vez de 9 (en el caso de tomar q ) : y 2 condiciones en vez de 1 : Mas a un, con esta formulacion la matriz de rotacion obtenida es p erfecta utilizando la ecuacion que p one en corresp ondencia los puntos (2) Los cuaterniones ya han sido utilizados en  ? ] con el n de calcular la matriz de rotacion R pero en este caso, los cuaterniones se utilizan para obtener una matriz de rotacion real a partir de la matriz Q ,sin tener relacion con la energa (7). En el caso en que los parametros intrnsecos de las camaras sean descono cidos, la situacion es mas compleja, sin embargo, p o demos llegar a una formulacion similar a la presentada en (2), pero en este caso, la matriz Q dep ende tambien de los parametros intrnsecos. En este caso la matriz Q se denomina matriz fundamental, que prop orciona la geometra epip olar de las camaras. Busqueda de un mnimo lo cal de la energa no-lineal E(X). Consideremos el problema de optimizacion no-lineal (7) E ( X )= T q ( X ) T Bq ( X ) con las condiciones s 2 + l 2 + m 2 + n 2 =1 y t 2 x + t 2 y + t 2 z = C 2 : Primero, calculamos las derivadas de la funcion q ( X ) que se pueden calcular facilmente gracias a que los comp onentes de q ( X ) son p olinomios. Denotamos p or r E ( X )elvector gradiente del funcional E ( X ). Al ser B simetrica, r E ( X ) se puede escribir como r E ( X )=2 T q ( X )  B  Jq ( X ) (8) donde Jq ( X )es el Jacobiano 9 x 7de la funcion q ( X ) : Para encontrar un mnimo lo cal de la energa E ( X ) aplicamos un meto do de descenso p or gradiente. Esto signica que utilizando una primera estimacion X 0 del mnimo (en muc has o casiones se puede suministrar una estimacion a priori" sobre la p osicion de las dos camaras a partir de alg un tip o de calculo aproximado), aproximamos el mnimo lo cal mas cercano de E ( X )como el estado asintotico del esquema iterativo: X n +1 = X n ;  n d n (9) donde X n representa la aproximacion del mnimo lo cal de E ( X ) en el paso n . d n es la direccion de descenso en el paso n y  n es un parametro. Por lo tanto, para calcular X n +1 a partir de X n  necesitamos determinar en cada paso el valor de d n y  n : En un problema de minimizacion sin restricciones, la eleccion natural de d n es r E ( X n ) : Para poder incluir la informacion de las restricciones en la direccion de descenso d n  se proyecta el vector r E ( X n ) en la interseccion del espacio tangente a las sup ercies s 2 + l 2 + m 2 + n 2 = 1 y t 2 x + t 2 y + t 2 z = C 2 : De esta manera se minimiza la distorsion con resp ecto a las restricciones de la nueva estimacion X n +1 : En el siguiente lema, se demuestra como se puede calcular esta proyeccion. Lema 1 La proyeccion del vector r E ( X n ) en la interseccion del espacio tangente a las supercies s 2 + l 2 + m 2 + n 2 =1 y t 2 x + t 2 y + t 2 z = C 2 viene dado por d n =( Id ; p  p ; q  q ) r E ( X n ) (10) donde p y q son vectores unitarios p = T ( s n l n m n n n  0  0  0) q = T (0  0  0  0 t n x t n y t n z ) C y p  p , q  q representan la matriz 7 x 7 p T p y q T q: Demostracion: Por un lado, la direccion normal en el punto X n de la supercie s 2 + l 2 + m 2 + n 2 =1 es el vector p y la direccion normal en el punto X n de la supercie t 2 x + t 2 y + t 2 z = C 2 viene dada por el vector q: Observemos que si el vector Y es ortogonal a p y q , es decir, Y es la interseccion de los espacios tangentes, entonces ( Id ; p  p ; q  q ) Y = Y por otro lado, p y q son ortogonales y k p k = k q k =1 ( Id ; p  p ; q  q ) p = p ; p =0 ( Id ; p  p ; q  q ) q = q ; q =0 y, por lo tanto, ( Id ; p  p ; q  q ) representa la matriz de proyeccion en la interseccion de los espacios tangentes a s 2 + l 2 + m 2 + n 2 =1 y t 2 x + t 2 y + t 2 z = C 2 : Una vez que se calcula la direccion de descenso d n , se elige  n minimizando la funcion (  )= E ( X n ; d n ) (11) la condicion de extremo de la funcion (  )= E ( X n ; d n ) es  0 (  )= hr E ( X n ; d n ) d n i =0 (12) Para calcular el valor optimo de  se utiliza una aproximacion de primer orden de r E ( X ), esto es r E ( X n ; d n )  = r E ( X n ) ; HE ( X n ) d n donde HE ( X ) es la matriz Hessiana del funcional E ( X ) : Por lo tanto, si sustituimos la funcion anterior en la ecuacion (12), obtenemos:  n = hr E ( X n ) d n i h d n HE ( X n ) d n i (13) La matriz Hesiana HE ( X n )se puede calcular facilmente utilizando la expresion HE ( X )=2 T Jq ( X )  B  Jq ( X )+2 T q ( X )  B  Hq ( X ) (14) donde q ( X )  B  Hq ( X )es la matriz 7 x 7 dada por ( q ( X )  B  Hq ( X )) ij = q ( X )  B  @ 2 q ( X ) @X i @X j Nota: Alser d n la proyeccion de r E ( X n ) en el espacio ortogonal de p y q entonces existen   tal que r E ( X n )+ p + q = d n entonces hr E ( X n ) d n i = h d n d n i , de tal forma que el numerador en el calculo de  n es igual a cero s y solo s d n = 0 y, por lo tanto en este caso,  y  representan los multiplicadores de Lagrange de un mnimo local del funcional E ( X ) con las restricciones s 2 + l 2 + n 2 + m 2 =1 y t 2 x + t 2 y + t 2 z = C 2 obtenidos por la tecnica de los multiplicadores de Lagrange. En otras palabras, si d n =0 , la solucion asociada X n satisface la condicion del mnimo local suministrado por la tecnica de Lagrange. Resumiendo, el algoritmo completo para encontrar el mnimo lo cal del funcional E ( X )se puede expresar en los siguientes pasos: 1. Se elige un valor inicial para X 0 (p or ejemplo, la rotacion y traslacion obtenidas por el meto do lineal) 2. Hasta la convergencia de X n (a) Se calcula r E ( X n ) utilizando ec. (8) (b) Se calcula d n utilizando ec. (10) (c) Se calcula HE ( X n ) utilizando ec. (14) (d) Se calcula  n utilizando ec. (13) (e) Se calcula X n +1 = X n ;  n d n : (f ) Se normaliza X n +1 para cumplir con las restricciones. Exp eriencias numericas La simulacion que realizamos consistio en lo siguiente: Elegimos 5 puntos (3 ; D )de un cub o inscrito en la esfera de centro (0  0  2) y radio 1 : Situamos el fo co de la primera camara en el origen y su plano proyectivo, tangente a la esfera de centro (0  0  2) y radio 2, en el punto (0  0  4) : El fo co de la segunda camara se sit  ua en el punto (2  0  2) sobre la misma esfera y su plano proyectivo, tangente ala misma, en el punto ( ; 2  0  2) : Proyectamos los 5 puntos (3 ; D ) en ambas camaras. En este caso los valores de los parametros, que determinan la p osicion de la segunda camara con resp ecto al de la primera, vienen dados p or ( s l  m n t x t y t z )=( 1 p 2  0 : 0  1 p 2  0 : 0  2 : 0  0 : 0  2 : 0), que sera el vector X que minimiza la energa E ( X ). Para realizar las exp eriencias numericas tomamos como aproximacion inicial X 0 una p erturbacion de ( X )a~nadiendole un ruido uniformemente distribuido. Para esta aproximacion inicial ( X 0 ), aplicamos los siguientes meto dos de optimizacion no-lineal: 1. El meto do propuesto. 2. El meto do de gradiente paso optimo (normalizando, en cada iteracion, los vectores ( s l  m n ) y ( t x t y t z )). 3. El meto do de gradiente paso alterno, en el que en cada iteracion se toma de forma alternada las direcciones d i =(0 :::  1 |{z} i :::  0). Realizamos 10.000 pruebas distintas para la conguracion anterior teniendo en cuenta dos casos distintos: En el primero a~nadimos un ruido uniformemente distribuido entre  ; 0 : 1  0 : 1] sobre el vector X , y en el segundo un ruido entre  ; 0 : 25  0 : 25]. En las siguientes tablas mostramos los resultados obtenidos para estos dos casos. Calculamos la media del n umero de iteraciones necesarias para converger y su desviacion estandar para cada meto do. Meto do Media Desviacion Meto do propuesto 4814 2 : 178 Gradiente paso optimo 4989 2 : 233 Gradiente paso alterno 12 : 770 2 : 078 Ruido 0.1 sobre X Meto do Media Desviacion Meto do propuesto 6 : 409 2 : 709 Gradiente paso optimo 6 : 695 2 : 693 Gradiente paso alterno 13 : 842 2 : 543 Ruido 0.25 sobre X De los resultados obtenidos se deduce que el meto do propuesto converge, de forma general, mas rapidamente hacia el resultado nal que los otros dos meto dos. Agradecimientos Este traba jo ha sido parcialmente nanciado por la accion integrada HispanoFrancesa HF98-0098 y el proyecto espa ~nol PB95-1225 de la D.G.I.C.Y.T. Referencias 1] O.Faugeras, \3-D computer vision. A geometric viewp oint," MIT Press , 1993. 2] O.Faugeras y S.Maybank \Motion from point matches: multiplicityof solutions," International Journal of Computer Vision ,Vol. 4(3) pp 225-246, 1990. 3] K.Homann,C.Metz y Y.Chen \Determination of 3-D imaging geometry and object congurations from two biplane views: An enhancement of the Metz-Fencil technique.," Med. Phys. ,Vol. 22(8) pp 1219-1227, 1995. 4] H.C.Longuest-Higgins, \A computer algorithm for reconstructing a scene from two pro jections," Nature , Vol. 293, 133, 1981. 5] C.Metz y L.Fencil, \Determination of three-dimensional structure in biplane radiography without prior knowledge of the relationship between the two views: Theory," Med. Phys. , Vol. 16(1), pp. 45-51, 1989. 6] J.Weng,T.S.Huang y N.Ahuja, \Motion and structure from two p ersp ective views: algorithms, error analysis, and error estimation, " IEEE Trans. Pattern Anal. Machine Intel. PAMI, , Vol.11 451(1989) 1. Departamento de Informatica y Sistemas. Universidad de Las Palmas de Gran Canaria. Campus de Tara. 35017 Las Palmas. e-mail: f lalvarez/jsanchez g @dis.ulpgc.es, http://serdis.dis.ulpgc.es/ lalvarez