12.05.2020

Vectores:
 
     [a]
va = [b] , va = [a, b, c, d]^T
     [c]
     [d]

Las matrices en letras mayúsculas 

    [a11 a12 a13]
A = [a21 a22 a23]
    [a31 a32 a33]


Homografía

La homografía es una transformación general
de los puntos de un plano a otro plano

 ----------                  --
|          |                -   --
|          |               -       --  
|  {vp1}   | ---> H --->  -   {vp2}  -
|          |              --        -
 ----------                 --    -
                               --

Una homografía transforma los puntos vp1 a vp2
Un punto en 2D homogéneo vp1 = [x1, y1, 1]^T

vp2 = H  vp1         (1)
3x1  3x3 3x1 
       3x1

H es una matríz de tamaño 3x3

  [x2]   [h11 h12 h13][x1]
l [y2] = [h21 h22 h23][y1]    (2)
  [ 1]   [h31 h32 h33][ 1]

l x2 = h11 x1 + h12 y1 + h13
l y2 = h21 x1 + h22 y1 + h23    (3)
   l = h31 x1 + h32 y1 + h33

l es un factor de escala (por usar puntos 2D
homogéneos). Una homografía no tiene 9 grados de libertad.
No tiene 9 incógnitas. Por el factor de escala, el número
de incógnitas es 8.
Una homografía tiene 8 grados de libertad. 

Vamos a ver un procedimiento para calcular una homografía.
En (3), el factor de escala l lo sustituimos en la primera
equationd de (3):

(h31 x1 + h32 y1 + h33) x2 = h11 x1 + h12 y1 + h13
(h31 x1 + h32 y1 + h33) y2 = h21 x1 + h22 y1 + h23  (4)

Dados un par de puntos, deben ser un par de correspondencias
de puntos, esto es, de la misma características en ambas 
imágenes de entrada, se obtienen dos ecuaciones como en (4).

De (4):
h31 x1 x2 + h32 y1 x2 + h33 x2 = h11 x1 + h12 y1 + h13
h11 x1 + h12 y1 + h13 - h31 x1 x2 - h32 y1 x2 - h33 x2 = 0  (5)

(h31 x1 + h32 y1 + h33) y2 = h21 x1 + h22 y1 + h23 
h31 x1 y2 + h32 y1 y2 + h33 y2 = h21 x1 + h22 y1 + h23 
h21 x1 + h22 y1 + h23 - h31 x1 y2 - h32 y1 y2 - h33 y2 = 0  (6)

Para encontrar los elementos de H, pongo todos sus elementos en
un solo vector:
  vh = [ h11, h12, h13, h21, h22, h23, h31, h32, h33 ]^T

Y contruimos el sistema de ecuaciones:
  A vh = vb

que se resulve como:
  vh = A^{-1} vb

Si vb = 0, es un sistema homogeneo de ecuaciones y no se puede
resolver invirtiendo la matriz A.

Fijamos la escala de la homografía haciendo h33 = 1, entonces
las ecuaciones (5) y (6) quedan como,
de (5):

h11 x1 + h12 y1 + h13 - h31 x1 x2 - h32 y1 x2 - h33 x2 = 0  
Aqui hacemos h33 = 1: 
h11 x1 + h12 y1 + h13 - h31 x1 x2 - h32 y1 x2 - x2 = 0  
h11 x1 + h12 y1 + h13 - h31 x1 x2 - h32 y1 x2 = x2    (7)
De (6):
h21 x1 + h22 y1 + h23 - h31 x1 y2 - h32 y1 y2 - h33 y2 = 0 
Aqui hacemos h33 = 1: 
h21 x1 + h22 y1 + h23 - h31 x1 y2 - h32 y1 y2 - y2 = 0
h21 x1 + h22 y1 + h23 - h31 x1 y2 - h32 y1 y2 = y2    (8)

[x1, y1, 1,  0,  0, 0, - x1 x2, - y1 x2 ][h11]   [x2]
[ 0,  0, 0, x1, y1, 1, - x1 y2, - y1 y2 ][h12]   [y2]
[                                       ][h13]   [  ]
[                                       ][h21] = [  ]
[                                       ][h22]   [  ]
[                                       ][h23]   [  ]
[                                       ][h31]   [  ]
[                                       ][h32]   [  ]

Necesitamos cuatro correspondencias de puntos para poder
calcular una homografía.

Los puntos deben de normalizarse a media cero y varianza uno,
porque si no es es inestable numéricamente.

vp1' = T1 vp1
vp2' = T2 vp2

l vp2' = H' vp1'
l T2 vp2 = H' T1 vp1 
l T2^{1} T2 vp2 = T2^{-1} H' T1 vp1 
l vp2 = T2^{-1} H' T1 vp1 

H = T2^{-1} H' T1

x - pmx
-------- así se transforman las coordenadas x a media 0 y desviación estándar 1
  sx

sx = sigma (la desviación estándar) de todos los valores de x
pmx = promedio de todos los valores de x 

vp1' = T1 vp1
[x1']   [1/sx,  0,   -pmx/sx ][x1]
[y1'] = [  0 , 1/sy, -pmy/sy ][y1]
[ 1 ]   [  0 ,  0r,     1    ][ 1]


------------------------------------------------------- 
Para poder programar este procedimiento nos falta poder
resolver :

  A vh = vb

Se puede resolver usando la descomposición QR
A = Q R
Q es matriz ortogonal, det(Q) = 1, inv(A) = Q^T,
cualquier vector renglón o columna su norma es igual a 1.
R es una matriz triangular superior.

  A vh = vb,
  QR vh = vb,
  Q^{-1} QR vh = Q^{-1} vb,
  Q^T Q R vh = Q^T vb,    Q^T Q = I
  R vh = Q^T vb

Aquí no tenemos que invertir R.
Si tenemos:

  R vh = vc

  [r11  r12  r13][h1]   [c1]
  [  0  r22  r23][h2] = [c2]
  [  0    0  r33][h3]   [c3]
  
  La última incógnita se calcula como:

  r33 h3 = c3
  h3 = c3/r33   (*1)

  La segunda ecuación es:
  r22 h2 + r23 h3 = c2 
  h2 = ( c2 - r23 h3 )/r22  (*2)

  La primera ecuación es:
  r11 h1 + r12 h2 + r13 h3 = c1 
  r11 h1 = c1 - r12 h2 - r13 h3
  h1 = (c1 - r12 h2 - r13 h3)/r11  (*3)

  El procedimiendo para realizar la retrosustitución
  en resolver A vh = vc:

  Entrada: A de tamaño n x n
           El vector vc
  Salida: El vector hv

  Q, R = calcula_QR( A )
  h_n = c_n/r_nn

  for( i=n-1, i>0, i-- ) :
     suma = 0
     for( j=i+1, j<=n, j++ ) :
        suma += r_ij * h_j

     h_i = (c_i - suma)/r_{ii} 
     
 
     # h1 = (c1 - r12 h2 - r13 h3)/r11  (*3)
     #  i      
     # Debemos calcular la suma r12 h2 + r13 h3
     #                           ij  j    ij  j
     #                           j=i+1 hasta n 

Una matriz de tranlación

[ 1 0 dx ]
[ 0 1 dy ]
[ 0 0  1 ]

Una matrix de rotación
[cos t  -sen t  0 ]
[sen t   cos t  0 ]
[  0       0    1 ]