Estimaci´on de homograf´ıas Visi´on en Rob´otica 1er cuatrimestre de 2013

1.

Introducci´ on del problema

Una homograf´ıa es una transformaci´on proyectiva que determina una correspondencia entre puntos. El problema que se plantea resolver es el siguiente: dado un conjunto de puntos xi ∈ Pn y su conjunto de puntos correspondientes x0i ∈ Pn , calcular la transformaci´on proyectiva que lleva xi a x0i . Es decir, calcular h : Pn → Pn tal que h (x) = x´ = Hx. El objetivo es estimar la matriz H.

2.

Algoritmo de transformaci´ on lineal directa (DLT) Este m´etodo sirve para estimar la matriz H antes mencionada.

Algoritmo 1 Dados un conjunto de n puntos de correspondencia (n ≥ 4), {xi ↔ x´i }, determinar H (homograf´ıa) tal que x´i = Hxi 1. Para cada correspondencia xi ↔ x´i , calcular la matriz Ai de dimensi´on 2×9  Ai =

−wi0 xTi 0T

0T wi0 xTi

yi0 xTi −x0i xTi



donde xi = (xi , yi , wi ) y 0T = (0, 0, 0) 2. Generar una matriz A, de dimensi´on 2n × 9, con las matrices Ai 3. Descomponer A seg´ un el SVD: A = U DV t ,

1

4. Luego, la matriz H que estamos buscando es la siguiente:   h1 h2 h3 H =  h4 h5 h6  h7 h8 h9 donde h es la u ´ltima columna de V . Notar que h es un vector de 9 elementos. Ver la demostraci´on correspondiente en el ap´endice A.

3.

DLT con normalizaci´ on de datos

Se hace una traslaci´ on y un escalamiento de tal manera que el algoritmo queda invariante con respecto a la elecci´on arbitraria de la escala y del origen de coordenadas. Se elige una sistema de coordenadas can´onico donde el algoritmo DLT es invariante a las transformaciones por similaridad. Algoritmo 2 DLT con normalizaci´ on de los datos Dado un conjunto {xi → x0i } de n puntos de correspondencia, calcular las e yH e 0 que consisten en una traslaci´ transformaciones de similaridad H on y un escalamiento. 1. Traslaci´ on de las coordenadas en cada imagen: para los datos xi y x0i , se ¯ yx calculan el centroide x ¯0 de cada conjunto de puntos correspondiente, seg´ un: n

¯= x

1X xi n i=1

¯0 = x

1X 0 x n i=1 i

n

2. El centroide ser´ a el nuevo origen de coordenadas en cadacaso. Es decir, se genera un nuevo conjunto de puntos de correspondencia xei → xei 0 donde ei = s (xi − x ¯) y x e0i = s0 (x0i − x ¯0) x con s y s0 escalares que se obtienen como se describe a continuaci´ on. ei y de los x e0i . Hallamos un escalar 3. Calculamos la norma media de los x 0 s y s tal que al multiplicar √ cada uno por la norma media de los puntos trasladados, el resultado sea 2, es decir: Si

n

1X ¯ )k d¯ = k(xi − x n i=1 2

queremos que sd¯ =

√

2

entonces √

2 s= ¯ d En forma an´ aloga se define: √ 2 s0 = ¯0 d 4. Finalmente, los dos conjuntos de puntos hom´ ologos son: √

√

2 ¯) , ei = ¯ (xi − x x d

e0i x

2 ¯ 0 ) con i = 1, ..., n = ¯0 (x0i − x d

e yH e 0 de simiPara aplicar dichas transformaciones se pueden plantear las H laridad correspondientes. De esta forma, los nuevos puntos quedar´ıan definidos seg´ un

e i x ei = Hx e 0 x0i x e0i = H

(1) (2)

Luego de aplicar estas transformaciones, se utiliza el algoritmo DLT sin e0i }. alteraciones para los nuevos conjuntos de n puntos hom´ ologos: {e xi → x Sin embaro, el resultado del algoritmo DLT ser´ a en realidad una transformacion H˙ tal que: x e0i = H˙ x ei

(3)

mientras que lo que se buscaba en realidad es una transformaci´ on H tal que x0i = Hxi

(4)

Asi que, para obtener H aplicaremos las definiciones (1) y (2) de x ei y x e0i en la ecuaci´ on (3) obtieniendo: e 0 x0i = H˙ Hx e i H e 0−1 obtenemos que y multiplicando a ambos lados por H e 0−1 e 0 0 e 0−1 ˙ e |H {z H} xi = H H Hxi I n×n

3

por lo que, utilizando (4), queda que: e 0−1 H˙ H e H=H que era lo que buscabamos.

4.

RANSAC

Si nos restringimos a P1 , el problema de correspondencia ser´ıa el siguiente. Supongamos que queremos hallar una transformaci´on af´ın unidimensional, Ha , entre un conjunto de puntos correspondientes {xi → x0i } que est´an sobre dos l´ıneas.      a b xi axi + b 0 xi = Ha xi = = , 0 1 1 1 entonces xi = axi + b con xi ∈ R.

(5)

El objetivo entonces es hallar una recta que minimice la suma de las distancias ortogonales al cuadrado, d2⊥i , sujeto a la condici´on de que no haya ning´ un punto que se desv´ıe de la recta buscada m´as all´a de una cierta distancia umbral t, que depender´ a del ruido de medici´on. Entonces, dado un conjunto de puntos hom´ ologos debemos 1. estimar una recta de ajuste, y 2. clasificar los puntos como inliers y outliers (seg´ un t). Se propone un estimador robusto para elegir los puntos hom´ologos, denominado RANSAC (por RANdom SAmple Consensus), que puede usarse cuando la proporci´ on de outliers es muy grande.

4.1.

Idea general del algoritmo (en P1 )

1. Se seleccionan aleatoriamente 2 puntos, con los que se calcula una recta r que pasa por ellos. Esta recta es el modelo de ajuste de los puntos (ecuaci´on (5)). 2. Se calcula el soporte de r, definido este como el conjunto de puntos que se encuentran a una cierta distancia del modelo. En esta caso, esa distancia ser´ a la distancia ortogonal de los puntos a la recta. 3. Se repite esta selecci´ on aleatoria (pasos 1 y 2) un n´ umero N de veces. 4. La recta de ajuste con mayor soporte es la recta (o modelo) de ajuste elegida. Los puntos que caen dentro del soporte ser´an los inliers del modelo y definen el conjunto de consenso S (fig. 1). 4

c

b

a d

Figura 1: Dos conjuntos de consenso (demarcadas mediante l´ıneas punteadas) con diferentes soportes. El conjunto de consenso asociado a la recta que pasa por los puntos a,b tiene soporte igual a 10.

4.2.

Algoritmo general (en P2 )

Retornando ahora a P2 , volvemos ahora a la formulaci´on original del problema de correspondencia, planteado al principio de este apunte. En este caso, el modelo ya no ser´ a una recta sino una homograf´ıa planar. Como se mencion´ o previamente, el algoritmo tiene una forma iterativa, con un umbral N en la cantidad de iteraciones. Sin embargo, como se ver´a a continuaci´ on, N ser´ a estimado en forma adaptativa, en base a la proporci´on de outliers obtenidos hasta el momento (inicialmente, en base a un  dado). Algoritmo 3 RANSAC 1. Definir un  ∈ [0, 1] que representa la proporci´ on estimada de outliers asociados al modelo 2. Elegir s pares de correspondencias de puntos no-colineales y estimar la homograf´ıa entre dichos pares (por ejemplo, utilizando DLT). En P 2 , alcanza con s = 4 para fijar los grados de libertad necesarios de H. 3. Para cada par de puntos, calcular la distancia entre estos y los puntos transformados con H. Si bien en P1 esto se hac´ıa directamente con la distancia ortogonal, en P2 la medida de distancia debe ser redefinida. Por ello, se define d2i de la siguiente forma: d2i = d(xi , H −1 xi 0 )2 + d(xi 0 , Hxi )2 donde d(·, ·) corresponde a la distancia eucl´ıdea.

5

4. Clasificar todos los pares de puntos xi → x0i como outliers o inliers en base a un umbral t, definiendo as´ı el conjunto de consenso S para la iteraci´ on actual, que contendr´ a a los inliers. Es decir:  S = xi → x0i / d2i < t2 El umbral t est´ a definido de antemano y se encuentra tabulado. Para el caso de P2 , t2 = 5,99σ 2 (donde se puede tomar σ = 1). 5. Repetir el proceso N veces, qued´ andonos siempre con la H que maximiza la cantidad de inliers (es decir, que maximiza el cardinal del conjunto de consenso S asociado a la H actual). Cada vez que se obtenga un H mejor a la anterior, se debe reestimar  y N , seg´ un: #S n log(1 − p) N← log(1 − (1 − )s ) ←1−

donde p se define de antemano y corresponde a la probabilidad de que al menos una selecci´ on de puntos est´e libre de outliers. En general se toma p = 0,99. 6. Una vez terminado el ciclo, solo se debe re-estimar la H encontrada (nuevamente, puede utilizarse DLT), esta vez usando los puntos pertenecientes al mejor conjunto de consenso S encontrado.

6

A.

Descomposici´ on en valores singulares (SVD)

Sean A una matriz de dimensi´on m × n, con m ≥ n. Entonces se puede descomponer de la forma: A = U DV t donde U es matriz ortogonal de dimensi´on m × n, D es una matriz diagonal de n×n y V es una matriz ortogonal de n×n. Como U tiene columnas ortogonales, entonces U t U = In×n . Adem´as U preserva la norma, es decir, se cumple que kU xk = kxk ∀x pues, kU xk =

q

t

√

(U x) (U x) =

xt U t U x =

√

xt Ix = kxk

Los valores singulares de la matriz A son los valores de la diagonal D (son no negativos), que son las ra´ıces cuadradas de los valores propios (tambi´en llamados autovalores) de la matriz At A. Para ver esto observemos que si: A = U DV t entonces: At A = U DV t


t


t U DV t = V DU U DV t = V D2 V t = V D2 V −1 |{z} In×n

pues, como V t V = In×n entonces V −1 = V t . Luego se tiene que: At A = V D2 V −1 Ahora, si llamamos H = At A, y multiplicamos a izquierda de ambos miembros por V , tenemos que: HV = V D2 . Llamando

   vi =   

v1i v2i . . vni

     

a cada columna de la matriz V , podemos escribir: Hvi = d2i vi ,

7

es decir,  (Hv1

  Hv2 ...Hvn ) =   

v11 . . . vn1

v12 v22 . . vn2

.

.

v1n . . . . vnn

.

= d21 v1 , ..., d2n vn

     

d21 0 . . 0

0 d22

. 0 .

.

.

. . . 0

0 0 . . d2n

   =  




con d21 ≥ d22 ≥ ... ≥ d2n . Estas ecuaciones definen los valores propios de At A, que son los elementos de la diagonal D2 , y vienen dados por d2i con i = 1, .., n. Las columnas de V son los vectores propios (o autovectores) de At A, luego los valores singulares de A son las ra´ıces cuadradas de los valores propios de At A. Como At A es una matriz sim´etrica y definida positiva entonces los valores propios son reales y positivos, y entonces los valores singulares tambi´en son reales y positivos. As´ı escrito tenemos que la u ´ltima columna de la matriz V corresponde al valor singular menor de la matriz A. Proposici´ on 1 Dada la matriz A de m × n, tal que m ≥ n, hallar x, de dimensi´ on n, tal que minimice kAxk sujeto a la condici´ on kxk = 1. La soluci´ on es la u ´ltima columna de la matriz V , donde A = U DV t es la descomposici´ on SVD de la matriz A. Desmostraci´ on: Sea A de m × n, podemos realizamos la descomposici´on SVD, entonces A = U DV t . Luego queremos minimizar: kAxk = kU DV t xk Se tiene que:



U DV t x = DV t x pues U es ortogonal y preserva la norma, luego queremos minimizar kDV t xk bajo la condici´ on kV t xk = 1, esto es cierto pues se cumple que: kV t xk = t (V t x) (V t x) = kxk = 1. Sea y = V t x, hay que minimizarkDyk bajo la condici´on kyk = 1, donde la matriz diagonal   d1 0 . . 0  0 d2 0 . 0    . .  D=  .   . . .  0 . . 0 dn es tal que d1 ≥ d2 ≥ ... ≥ dn , o sea,  d1 0 . . 0  0 d2 0 . 0   . . .   . . . 0 . . 0 dn

     

8

y1 y2 . . yn





    =    

d1 y1 d2 y2 . . dn yn

     

tal que q

y12 + ... + yn2 = 1.   0  0  n X   2  Como kDyk = d2i yi2 ⇒ y =   .  es el vector que minimiza la norma  .  i=1 1 kDyk, entonces      v1n 0 v11 v12 . . v1n      . v22 .   0   .        . . .  .  =  .  x = V y = .   . . . .  .   .  vnn 1 vn1 vn2 vnn que es la u ´ltima columna de la matriz V.

B.

RANSAC

El algoritmo de RANSAC incluye los siguientes t´erminos: : proporci´ on de outliers que se estima existien en el modelo (este par´ametro luego se ir´ a re-estimando en forma adaptativa). t: umbral de clasificaci´on de puntos como inliers o outliers T : tama˜ no del conjunto de consenso, es decir T = #S N : n´ umero de iteraciones o muestras aleatorias para ensayar el modelo.

B.1.

Par´ ametro t

Si definimos α como la probabilidad de que un punto dado sea inlier, es decir: P (x es inlier) = α lo que se busca es encontrar un t tal que satisfaga esta ecuaci´ on. Supongamos que los errores de medici´on, e, siguen una distribuci´on normal
con media 0 y varianza σ 2 , es decir e ∼ N 0, σ 2 , entonces se puede calcular t, a partir de e. Con estas consideracioes, el cuadrado de la distancia ortogonal, d2⊥ , es la suma de variables aleatorias gaussianas, y por lo tanto, d2⊥ ∼ χ2n donde χ2n es la distribuci´ on Chi-cuadrado con n grados de libertad, y n es la codimensi´ on del modelo (en una homograf´ıa planar es n = 2). Para una l´ınea la codimensi´ on es 1 (s´ olo es considerada la distancia ortogonal a la l´ınea).

9

Entonces, como: P χ2n ≤ k


2

= Fn k


2

Zk

2

=

χ2n (ζ) dζ = α

0

tomamos t2 = Fn−1 (α) σ 2 En resumen, consideramos que un par de puntos correspondientes se clasifica como:  inlier si d2⊥ < t2 outlier si d2⊥ ≥ t2 En general se pide α = 0,95, o lo que es lo mismo P (x es inlier) = 0,95; es decir que la probabilidad de rechazar un punto que es inlier sea 0,05. En el caso de la estimaci´ on de una recta de ajuste la codimensi´on ser´a n = 1, mientras que en el caso de la homograf´ıa la codimensi´on ser´a n = 2. En el siguiente cuadro se muestran los valores del umbral t2 para α = 0,95, para cada valor de n, dependiendo del error de medici´on σ. n 1 2 3

B.2.

modelo recta, F homograf´ıa, c´amara tensor

umbral t2 3,84σ 2 5,99σ 2 7,81σ 2

Par´ ametro T

Para determinar un tama˜ no v´alido para el conjunto de consenso S, se puede pensar que deber´ıa ser similar al n´ umero de inliers que se cree hay inicialmente en los datos. Entonces, si  es la proporci´on de outliers y n es la cantidad total de datos, definimos T de la siguiente manera: T = (1 − )N

B.3.

Par´ ametro N

Para asegurar, con probabilidad p, que al menos una muestra aleatoria de s puntos est´ a libre de outliers, se debe obtener un N que sea la cantidad de muestras aleatorias a tomar, por lo ranto, la cantidad de iteraciones del algoritmo. Vamos a obtener N de la siguiente forma: Definimos α = P (el punto seleccionado es inlier) , entonces  = 1 − α = P (el punto seleccionado es outlier) αs = P (los s puntos sean inliers) 10

Recordemos que s es 2 en el caso de que el modelo sea una recta y 4 en el caso de una homograf´ıa planar. Luego 1 − αs = P (haya uno o m´as outliers entre los s puntos) . Con en N selecciones se tiene: N

(1 − αs )

= P (en N selecciones de s puntos haya outliers en todas ellas) .

Si p = P (en N selecciones de s puntos, al menos una no tenga outliers) entonces 1 − p = P (en N selecciones de s puntos, todas tengan outliers) entonces se tiene que: N

1 − p = (1 − αs ) . Tomando logaritmos: log (1 − p) = N log (1 − αs ) , luego N=

log (1 − p) log (1 − p) = s s log (1 − α ) log (1 − (1 − ) )

(6)

En resumen, N es el n´ umero de selecciones que se require para asegurar, con probabilidad p, que al menos una muestra no tiene outliers para un dado tama˜ no de muestra s y una dada proporci´on de outliers . B.3.1.

Determinaci´ on adaptativa del par´ ametro N

Si asumimos ahora  desconocido, para determinar N se puede proceder de la siguiente forma. Algoritmo 4 Algoritmo adaptativo para determinar N 1. Iniciar  = 0,5 y N = ∞. 2. Iniciar un contador de iteraciones (o muestras) i = 0. 3. Repetir, mientras i < N : a) Tomar una muestra aleatoria y contar la cantidad outliers k. b) Asignar  =

k N.

c) Calcular N usando p y , en base a la ecuaci´ on (6). d) Incrementar i. 11

B.3.2.

Funci´ on de costo robusta

P En vez de minimizar la funci´on de costo C = i d2⊥i sobre los inliers otra posibilidad ser´ıa minimizar una funci´on de costo robusta que incluya a todos los datos. Esta funci´ on de costo viene dada por: X D= γ (d⊥i ) i

donde

 γ (d) =

d2 t2

12

si si

d2 < t2 d2 ≥ t2