Visión Computacional

1.	Introducción

El Puntillismo aparece por vez primera en 1884, encabezado por el pintor Neo - Impresionista Georges Seurat, y contando con entre sus más fieles seguidores tales como Henri - Edmond Cross y Vlaho Bokovac.Su procedimiento empleado por estos artistas, consistente en poner puntos de colores puros en vez de pinceladas sobre tela; logrará una técnica que fue el resultado de los estudios cromáticos llevados a cabo por Georges Seurat (1859 - 1891), pintor francés, quien en 1884 llegó a la división de tonos por la posición de toques de colores que, mirados a cierta distancia, crean en la retina las combinaciones deseadas [Wikipedia].

a.	Manejo de colores RGB

Definimos un pixel como la unidad básica de una imagen. Es un punto en cuyo valor se almacena la información conserniente a los colores. Los colores básicos que manejaremos son Rojo, Verde y Azul. Así pixel lo podemos ver como

Manipulate[{{r,g,b},Image[{{{r,g,b}}}, "Byte"]}, {r, 0, 255}, {g, 0, 255}, {b, 0,255}]

En la siguiente figura podemos ver el como se almacena la información conserniente a los 32 bit o 4 Bytes

Una imagen discreta, la podemos definir como un arreglo bidimensional en el cual se almacena la información de cada punto de la imagen (256^3) = 16,777, 216 diferentes colores por pixel. Dado que la información se maneja en Bytes para almacenar un byte necesitaremos de 4 bytes, tres correspondientes a los colores Rojo, verde y Azul y un cuarto que tiene información del nivel de transparencia. En la cámaras es común hablar de su resolución en Megabytes, ¿Que relación existe con esto?. Bueno una camara de 3 Mega pixeles, tiene aproximadamente 314 pixeles de lado. Esto es

N[Sqrt[3*1000*1000]]

En Forma matricial podemos hacer una imagen con colores verdes si definimos, un arreglo bidimensional donde ponemos un vector de colores {R,G,B}

b = {{{0,255, 0}, {0, 0, 0}}, {{0,0,0}, {0, 255, 0}} } ;
 {MatrixForm[b],
Image[b, "Byte"]}

Para este ejemplo si nuestro arreglo de imagen se llama A los colores almacenados en el renglón i columna j, la notación para referirnos a este punto es (Subscript[A, i,j])

b.	La escala de gris

En el caso de una imágen en tono de gris solamente se maneja una sola matriz de datos, así el color rojo será igual que el verde y el azul. La unión de los colores dara como resultado blanco y la ausencia dará negro. Nuestra imagen solamente tendrá colores que van desde el negro al blanco en 256 tonos diferentes.

Manipulate[Image[{{r}}, "Byte"], {r, 0, 255}]

Utilizando imáganes en tono de gris, solamente será necesario almacenar la información conserniente a un canal de color, o un promedio de los tres canales de color. Esto reduce en 1/3 la información que debe ser almacenada.

c.	Interfase para desplegar imágenes en Java.

Ver Algoritmos y ejemplos en la página. ver inicial.java

2.	Operadores Sobre imágenes

a.	Operadores puntuales sobre Imágenes

Una imagen discreta, la podemos definir como un arreglo bidimensional en el cual se almacena la información de cada punto de la imagen. Como definimos en la sección anterior cada unidad es un pixel y en el caso de operadores puntuales consideraremos que cada uno de estos se aplica sobre cada pixel de la imagen.

Linea Recta

Consideremos el caso de la ecuación de una línea recta y = mx + b donde m es la pendiente y b el cruze por cero

Manipulate[Plot[m*x+b,{x,0,255}],{m, -1, 1},{b,-5,5}]

Así para una imágen dadas estas condiciones podemos cambiar la apariencia haciendo una proyección del espacio original de la imagen a un nuevo espacio dado por una línea recta

Manipulate [
A={{1, 0},{0,1}};
B= A;
Amax = Max[A];
Amin = Min[A];
m = (fmax -fmin)/(Amax - Amin);
b = fmin- m*Amin;
B = m*A+b;
{Image[B, "Byte"]}, {{fmax,255}, -255, 255}, {{fmin, 0}, -255, 255}
]

En el caso de una imagen con mayor información esta lucirá como:

f= ImageData[ColorConvert[Import["ExampleData/lena.tif"], "Grayscale"]];
f = 255*f;
fmin =Min[f];
fmax = Max[f];
gmax = 150;
gmin = 50;
g=(f-fmin)*(gmax-gmin)/(fmax- fmin)+gmin;
{Image[f, "Byte"], Image[g, "Byte"],
ImageHistogram[Image[f, "Byte"]],ImageHistogram[Image[g, "Byte"]]}

Otro ejemplo de operadores puntuales sobre imagenes es la conversión de colores RGB a escala de gris. Así para convertir la imagen a colores a tono de gris, simplemente sacamos el promedio de cada uno de los pixeles en la 3 bandas aplicando la operación

(Subscript[A, i,j])=((Subscript[R, i,j])+(Subscript[G, i,j])+(Subscript[B, i,j]))/3

Con esta operación tendremos una sola matriz con el promedio de los colores y cuando se crea la imagen a color las matrice R, G, B serán iguales a la matriz A.

Función Sigmoide

Otra manera de hacer esto es considerar una función diferente como la sigmoide, definida como y = (1/(1+(e^(-βx)))). Como luce una imagen en estas condiciones.

A1= ImageData[ColorConvert[Import["ExampleData/lena.tif"], "Grayscale"]];
b= 100;
A2=2A1/Max[A1]-1;
A2 = 1/(1+Exp[-b*A2]);
A2 = (A2+1)*125;
{Image[A1*255, "Byte"], Image[A2, "Byte"],
Plot[1/(1+Exp[-b*x]),{x, -1,1}],
ImageHistogram[Image[A1*255, "Byte"] // ImageAdjust],ImageHistogram[Image[A2, "Byte"]//ImageAdjust]}

En esta figura podemos ver la image original, la imagen procesada, la sigmoide aplicada y los histogramas (en la parte baja) de la imagen antes y después de aplicar la transformación.

b.	Operadores de Ventana sobre imágenes

El objetivo de un operador de ventana es reemplazar el valor de un pixel, por el contenido pesado en una ventana al rededor de esta. Así por ejemplo (Subscript[I, i,j])=(∑_(<k,l>∈W))(Subscript[I, k,l.]) En la imagen tenemos un operador correspondiente a la media, donde cada uno de los pixeles es reemplazado por el promedio en una ventana de tamaña (Subscript[n, k])

nr = 116;
nc = 150;
Img = ColorConvert[Import["ExampleData/lena.tif"], "Grayscale"];
A = ImageData[Img];
ImageDimensions[Img]
B = A;
nk = 1;
For[i = 10, i <= nr - 10, i++,
 For[j = 10, j <= nc - 10, j++, {
    suma = 0;
    For[k = i - nk, k ≤ i + nk, k++,
      For[l = j - nk, l ≤ j + nk, l++,
        suma += A[[k, l]] 
       ]
      ]
     B[[i, j]] = suma;
    }
   ];
 ]

{Image[A], Image[B] // ImageAdjust}

c.	Transformación Afín y Proyectiva.

Una transformación afin es una transformación geométrica que esta constitida por translación, escalamiento, rotación y cizallamiento.

Traslación

Una translación la podemos hacer simplemente asuminedo que las nuevas coordenadas (OverHat[x])=x+tx (OverHat[y])=y+ty  les sumanmos un valor tx o ty según corresponda. En coordenadas homogeneas queda como

(({{(OverHat[x])}, {(OverHat[y])}, {1}}))=(({{1, 0, tx}, {0, 1, ty}, {0, 0, 1}}))(({{x}, {y}, {1}}))

Manipulate [
{TranslationTransform[{tx, ty}], Graphics[Rectangle[], Frame->True],Graphics[GeometricTransformation[Rectangle[],TranslationTransform[{tx, ty}]], Frame->True]}, {tx, 0, 10}, {ty, 0, 10}]

Escalamiento

El escalamiento pude entenderse como hacer que una figura geométrica cambie su tamaño. Un escalamiento en x lo podemos representar por como x=x*sx y en y como y=y*sy.  En coordenadas homogeneas se pude expresar como

ScalingTransform[{sx, sy}]

Manipulate [
{ScalingTransform[{sx, sy}], Graphics[Rectangle[], Frame->True],Graphics[GeometricTransformation[Rectangle[],ScalingTransform[{sx, sy}]], Frame->True]}, {{sx,1}, 0, 10}, {{sy,1}, 0, 10}]

algo

Rotación

Consideremos el caso de un punto que rota respecto a un punto fijo. Las coordenadas x, y en forma polar las podemos obtener como x=r cos(θ) y y=r seno(θ). Si consideramos que esta gira un ángulo ψ entonces podemos representar esta rotación en forma polar como

(({{(OverHat[x])}, {(OverHat[y])}}))=(({{r cos(θ+ψ)}, {r  seno(θ+ψ)}}))=(({{r cos(θ) cos(ψ)-r seno(θ) seno(ψ)}, {r  cos(θ)seno(ψ)+ r seno(θ)cos(ψ)}}))=(({{x cos(ψ)-y seno(ψ)}, {x seno(ψ)+ y cos(ψ)}}))

Una matriz de rotación puede ser definida entonces como

(({{(OverHat[x])}, {(OverHat[y])}}))=(({{cos(ψ) , -seno(ψ)}, {seno(ψ), cos(ψ)}}))(({{x  }, {y}}))

En coordenadas homogeneas puede representarse como

r = RotationTransform[θ]

Manipulate [
{ang, RotationTransform[ang Degree], Graphics[Rectangle[], Frame->True],Graphics[GeometricTransformation[Rectangle[],RotationTransform[ang Degree]], Frame->True]}, {ang, 0, 90}]

Cizallamiento

El cizallamiento es una transformación dada por la matriz, donde cx es el ángulo de cizallamiento respecto al eje x.

ShearingTransform[cx Degree,{1,0},{0,1}]

Manipulate[
{ShearingTransform[cx Degree,{1,0},{0,1}], Graphics[GeometricTransformation[Rectangle[],ShearingTransform[cx Degree,{1,0},{0,1}]], Frame->True]},
{cx, 0, 90}]

En caso de querer aplicar el cizallamiento en la dirección y, la matriz de transformación será.

(({{1, 0, 0}, {tg(cy), 1, 0}, {0, 0, 1}}))

Una vez que unimos todas las transformaciones tenemos una composición dada por una  Rotacion R, un escalamiento S, un Cizallamiento C y una Translación T. Al conjunto de estas operaciones se le conoce como transformación afin.

Composition[TranslationTransform[{tx,ty}],RotationTransform[ang], ScalingTransform[{sx,sy}], ShearingTransform[cx, {1,0},{0,1}]]

Manipulate[
{Composition[
TranslationTransform[{tx,ty}],
RotationTransform[ang],
ScalingTransform[{sx,sy}],
ShearingTransform[cx Degree,{1,0},{0,1}]], Graphics[GeometricTransformation[Rectangle[],
Composition[
TranslationTransform[{tx,ty}],
RotationTransform[ang],
ScalingTransform[{sx,sy}],
ShearingTransform[cx Degree,{1,0},{0,1}]]
], Frame->True]},
{tx, 0, 10}, {ty, 0, 10},{ang, 0, 90},{cx, 0, 90}, {{sx,1}, 0, 10}, {{sy ,1}, 0, 10}]

d.	Interpolación.

Todas las transformaciónes aplicadas anteriormente son para un punto en el espacio de coordenadas [x,y], pero que sucede cuando hablamos de imágenes. Una imagen esta compuesta de colores en una posición con coordenadas enteras [i,j]. Si aplicamos una transformación afín las coordenas enteras pasaran a ser coordenadas en números reales, así por ejemplo en el caso de una rotación tendremos

Manipulate[
{
A={{Cos[0.1], -Sin[0.1], 0},{Sin[0.1], Cos[0.1], 0}, {0, 0,1}};
c = Transpose[{{i,2, 1}}];
x=A.c;
MatrixForm[x], " = ", MatrixForm[A] , MatrixForm[c] },
{i, 0, 10, 1}
]

En estas condiciones en probable descidir hacer un truncamiento, sin embargo cuando se trata de un escalamiento con sx=2 y sy = 1, tendremos que ciertos pixeles no serán visitados.

Manipulate[
A={{2, 0, 0},{0, 1, 0}, {0, 0,1}};
c = Transpose[{{i,2, 1}}];
x=A.c;
{MatrixForm[x], " = ", MatrixForm[A] , MatrixForm[c] },
{i,0,10, 1}]

Esta condición hará que no se tengán cubiertas todos los pixeles de una imagen ya que en el caso del escalamiento (sx = 2) solamente las posiciones pares de la nueva imagen tendrán valor.¿Como corregir este efecto?.En lugar de barrer todos los puntos de la imagen origen, barreremos todos los puntos de la imagen destino y cualquier pixel que no este en posición entera lo calcularemos utilizando interpolación bilineal. Consideremos el caso de una matriz de rotación

A = {{0.9510565162951535, -0.3090169943749474, 0.0},
 {0.3090169943749474, 0.9510565162951535, 0.0}, 
{0.0, 0.0, 1.0}} ;
MatrixForm[A]

al aplicar la transformación para cada uno de los puntos de la imagen de acuerdo con el siguiente código tenemos 

       for(i=0; i<nr; i++)
            for(j=0; j<nc; j++) {
                afin(i, j, x, T);
                
                if(x[0] >=0 && x[0] < nr && x[1] >=0 && x[1] < nc)
                    I2[(int) x[1]][(int) x[0]] = I1[i][j];
            }



Note los puntos negros generados por el redondeo. Esta condición es mas notoria cuando se escala una imagen con un factor superior a 1.

La función de interpolación bilineal en Java es

  static public double bilineal(double f[][], double x0, double y0)
  {
    int ii, jj, nr, nc;
    double a, b, y;

    ii = (int) y0;
    jj = (int) x0;

    a = y0 - ii;
    b = x0 - jj;

    nr = f.length;
    nc = f[0].length;

    if (ii < 0 || jj < 0 || ii > nr - 2 || jj > nc - 2)
      return 0;
    if ( (a == 0) && (b == 0))
      return f[ii][jj];

    if (ii == nr - 1 && jj == nc - 1)
      return f[ii][jj];

    if (b == 0 || jj == nc - 1) {
      y = (1 - a) * f[ii][jj] + a * f[ii + 1][jj];
      return y;
    }

    if (a == 0 || ii == nr - 1) {
      y = (1 - b) * f[ii][jj] + b * f[ii][jj + 1];
      return y;
    }

    y = (1 - a) * (1 - b) * f[ii][jj] + (1 - a) * b * f[ii][jj + 1] +
        a * (1 - b) * f[ii + 1][jj] + a * b * f[ii + 1][jj + 1];
    return y;
  }

Una manera de resolver este problema es aplicar la transformación sobre los puntos de la imagen destino utilizando la inversa de la función y la función de interpolación de acuerdo con el código

  static public void afin(double f[][], double g[][], double T[])
  {
    int i, j, nr, nc;

    double in, jn;

    nr = f.length;
    nc = f[0].length;


    for (i = 0; i < nr; i++)
      for (j = 0; j < nc; j++) { //System.out.println(T[0]);

        jn = T[0] * j  + T[1] * i  + T[2];
        in = T[3] * j  + T[4] * i  + T[5];
        g[i][j] = bilineal(f, jn, in);

      }
  }
  
  La diferencia la podemos ver en la imagen
  
   

ver (afin.java)

Transformación Log Polar

La geometría log polar en imagens fue motivada por la estrcutura de la retina de algunos sistema de vision y esta tiene algunas propiedades de compresion. Cuando comparamos el sistema usual de coordenadas cartesianas en imagenes, las imagenes log - polar permiten muestreo mas rápido en los sistema de visión artificial sin reducción en el tamaño de la vista y la parte centra de la retina (fovea).
  
  The log-polar transformation is a conformal mapping from the points on the cartesian plane (x,y) to points in the log-polar plane (ξ, η):

ξ=Log[(Sqrt[(x^2)+(y^2)])]
η=ArcTan[(y/x)]

Aplicando esto tenemos que una imagen luce como

 ver (logpolar.java)

static public void LogPolar (double datos[][], double res[][])
{
 int nr, nc, i, j, ren, col;
 double alfa, r, x, y;
 nr = datos.length;
 nc = datos[0].length;
 ren = res.length;
 col = res[0].length;
 System.out.println (ren + " " + col);
 double rmax = Math.sqrt (ren*ren + col*col)/2.0;
 for (i = 0; i < ren; i++) 
   for (j = 0; j < col; j++) {
     alfa = (2.0*Math.PI*i)/(ren - 2);
     r = Math.exp (j*Math.log (rmax)/col);
     x = r*Math.cos (alfa) + nc/2.0;
     y = r*Math.sin (alfa) + nr/2.0;
     res[i][j] = bilineal (datos, x, y);
   }
 }

 
static public void InvLogPolar (double datos[][], double[][] res)
{
 int nr, nc, i, j, ii, jj, ren, col;
 double alfa, r;
 double fac = 65.0;
 nr = datos.length;
 nc = datos[0].length;
 ren = res.length;
 col = res[0].length;
 double phimax = Math.log (Math.sqrt (nr*nr + nc*nc)/2.0);
 for (i = 0; i < ren; i++) 
   for (j = 0; j < col; j++) {
     ii = ren/2 - i;
     jj = col/2 - j;
     alfa = (Math.atan2 (ii, jj)/Math.PI + 1.0)*(double) (nr - 2)/2.0;
     r = Math.log (Math.sqrt (ii*ii + jj*jj))/phimax*nc;
     res[i][j] = bilineal (datos, r, alfa);
   }
}

Transformación Proyectiva

Una transformación proyectiva es la generalización de las transformaciones (R^2)->(R^2), en la que las líneas paralelas no son transformadas necesariamene en tales. Se puede expresar de la forma
(({{(OverHat[x])}, {(OverHat[y])}, {(OverHat[z])}}))=(({{A, t}, {(v^T), 1}}))(({{x}, {y}, {1}}))
donde (v^T) es la diferencia entre la transformación proyectiva y la afín. Note para el caso de la transformación proyectiva la coordenada (OverHat[z]) no necesariamente es unitaria lo cual implica que la imagen transfromada se fue a un plano con coordenadas diferentes a esta. Una manera de regresar al sistema de coordenadas es dividir entre (OverHat[z]). Así la matriz de transformación queda

(({{(OverTilde[x])}, {(OverTilde[y])}, {1}}))=(1/(OverHat[z]))(({{A, t}, {(v^T), 1}}))(({{x}, {y}, {1}}))
Esta ecuación es conocida tambien como homógrafía. La aplicación a una imagen de una transformación 

T=(({{1, 0, 0}, {0, 1, 0}, {1e-3, 0, 1}}))  es



El código Java para hacer esta transformación es:

static public void proyecccion (double f[][], double g[][], double T[])
{
 int i, j, nr, nc;
 double in, jn, kn;
 nr = f.length;
 nc = f[0].length;
 for (i = 0; i < nr; i++) 
  for (j = 0; j < nc; j++) { 
   // System.out.println (T[0]);
   jn = T[0]*j + T[1]*i + T[2];
   in = T[3]*j + T[4]*i + T[5];
   kn = T[6]*j + T[7]*i + 1.0;
   jn /= kn;
   in /= kn;
   g[i][j] = bilineal (f, jn, in);
  }
 }

Parametros a partir de puntos correspondientes

En esta sección vamos a calcular los parametros de una homografía, la cual representa un modelo proyectivo,  dada por la ecuación 

(({{(Subscript[x, 2])}, {(Subscript[y, 2])}, {(Subscript[z, 2])}}))=(({{(Subscript[h, 1]), (Subscript[h, 2]), (Subscript[h, 3])}, {(Subscript[h, 4]), (Subscript[h, 5]), (Subscript[h, 6])}, {(Subscript[h, 7]), (Subscript[h, 8]), (Subscript[h, 9])}}))(({{(Subscript[x, 1])}, {(Subscript[y, 1])}, {(Subscript[z, 1])}}))=(({{A, t}, {(v^T), k}}))(({{(Subscript[x, 1])}, {(Subscript[y, 1])}, {(Subscript[z, 1])}}))

donde el vector (Subscript[P, 1])=(([(Subscript[x, 1]),(Subscript[y, 1]), (Subscript[z, 1])])^T) es un punto en el espacio tridimensional, el vector (Subscript[P, 2])=(([(Subscript[x, 2]),(Subscript[y, 2]), (Subscript[z, 2])])^T) es la proyección en un plazo con (Subscript[z, 2])=1 y ambos puntos son correspondientes (Subscript[P, 1])⟷(Subscript[P, 2])

El problema consiste en determinar el vector H=(([(Subscript[h, 1]),(Subscript[h, 2]),(Subscript[h, 3]),(Subscript[h, 4]),(Subscript[h, 5]),(Subscript[h, 6]),(Subscript[h, 7]),(Subscript[h, 8])])^T) a partir de un conjunto de puntos correspondientes. Una manera de resolver este problema es utilizando mínimos cuadrados, sin embargo tiene problema de estabilidad numérica por lo que es preferible utilizar el método siguiente.

Dado los puntos y que estos son correspondientes, entonces podemos decir que (Subscript[P, 2])×H (Subscript[P, 1])=0, dado que ambos vectores on paralelos. Esto no lleva a la siguiente representación

(({{(0^T), -(Subscript[z, 2,i]) (Subscript[P, 1,i]^T), (Subscript[y, 2,i]) (Subscript[P, 1,i]^T)}, {(Subscript[z, 2,i]) (Subscript[P, 1,i]^T), (0^T), -(Subscript[x, 2,i]) (Subscript[P, 1,i]^T)}}))(({{(Subscript[H, 1])}, {(Subscript[H, 2])}, {(Subscript[H, 3])}}))=0

donde (Subscript[H, 1])=(([(Subscript[h, 1]),(Subscript[h, 2]), (Subscript[h, 3])])^T), (Subscript[H, 2])=(([(Subscript[h, 4]),(Subscript[h, 5]), (Subscript[h, 6])])^T) y (Subscript[H, 3])=(([(Subscript[h, 7]),(Subscript[h, 8]), (Subscript[h, 9])])^T) .

La solución del sistema de ecuaciones puede llevarse a cabo utilizando la descompisición en valores singulares SVD. La descomposición en valores singulares calcula tres matrices USV = SVD(A) donde U y V son matrices y S una matriz diagonal. La solución del sistema es el último vector de la matriz V (ver Harley-Zisserman pp 90).

Ejemplo

Dados los puntos P1 = { [10.0 10.0 1.0], [10.0 50.0 1.0],  [50.0 50.0 1.0 ], [50.0 10.0 1.0 ], [25.0 25.0 1.0]} y P2 = {[8.9108,  10.89108,  1.0], [4.95049, 46.53465, 1.0], [39.04761, 48.57142, 1.0], [42.85714, 14.28571, 1.0 ], [20.48780, 25.36585, 1.0 ], determinar la Homografía correspondiente. La siguiente imagen muestra los puntos P1 a la izquierda y los puntos P2 a la derecha.

(ver Proyectiva.java)

De acuerdo con la ecuación 

(({{(0^T), -(Subscript[z, 2,i]) (Subscript[P, 1,i]^T), (Subscript[y, 2,i]) (Subscript[P, 1,i]^T)}, {(Subscript[z, 2,i]) (Subscript[P, 1,i]^T), (0^T), -(Subscript[x, 2,i]) (Subscript[P, 1,i]^T)}}))(({{(Subscript[H, 1])}, {(Subscript[H, 2])}, {(Subscript[H, 3])}}))=0

Tenemos

A ={
{0.000,0.000,0.000,-10.000,-10.000,-1.000, 108.911,108.911,10.891},
{10.000, 10.000, 1.000, 0.000, 0.000, 0.000,-89.109,-89.109,-8.911},
{0.000, 0.000, 0.000,-10.000,-50.000,-1.000, 465.347, 2326.733, 46.535},
{10.000, 50.000, 1.000, 0.000, 0.000, 0.000,-49.505,-247.525,-4.950},
{0.000, 0.000, 0.000,-50.000,-50.000,-1.000, 2428.571, 2428.571, 48.571},
{50.000, 50.000, 1.000, 0.000, 0.000, 0.000,-1952.381,-1952.381,-39.048},
{0.000, 0.000, 0.000,-50.000,-10.000,-1.000, 714.286, 142.857, 14.286},
{50.000, 10.000, 1.000, 0.000, 0.000, 0.000,-2142.857,-428.571,-42.857},
{0.000, 0.000, 0.000,-25.000,-25.000,-1.000, 634.146, 634.146, 25.366},
{25.000, 25.000, 1.000, 0.000, 0.000, 0.000,-512.195,-512.195,-20.488}};
{u, s, v} = SingularValueDecomposition[A]; 

Print["A = ", MatrixForm[A]],
Print["v = ", MatrixForm[v]],
Print["H = ", MatrixForm[v[[All,9]]/v[[9,9]]]]