Siguiendo con la segunda parte de nuestro artículo Transformaciones afines en Python con Numpy, Pillow y OpenCV (Parte 1), vamos en esta ocasión a aplicar las transformaciones afines que aprendimos en la parte anterior al mundo de la visión por computador a través de imágenes.  

En esta segunda parte, pasaremos a demostrar un uso más práctico a través de aplicar transformaciones afines con Python a imágenes vectoriales.

Trabajando las transformaciones afines en Python con imágenes

A estas alturas espero haber sido capaz de construir una cierta intuición acerca de cómo las transformaciones afines se utilizan simplemente para moverse alrededor de los puntos en el espacio 2D, así que con eso fuera del camino me gustaría empezar a trabajar con algunos datos de imágenes reales para dar una demostración más concreta de cómo funciona todo esto.

Esto también me permite cubrir otro tema importante de las transformaciones afines que tiene que ver con la tercera dimensión. La tercera dimensión de los datos de una imagen representa el valor real del píxel, a veces denominado dominio de intensidad, mientras que la ubicación física 2D de los píxeles en las otras dos dimensiones se denomina dominio espacial.

Para empezar voy a leer y mostrar una imagen utilizando matplotlib, que es simplemente una gran letra mayúscula R.

img = plt.imread('letterR.jpg')
img.shape #  (1000, 1000, 4)

Utilizando el método imread(...) puedo leer la imagen JPG, que representa la letra mayúscula R, en un ndarray NumPy. A continuación, muestro las dimensiones de la matriz que son 1000 filas por 1000 columnas, que en conjunto constituyen 1.000.000 de ubicaciones de píxeles en el dominio espacial. Los datos de los píxeles individuales están en forma de una matriz de 4 enteros sin signo que representan un canal rojo, verde, azul y alfa (o muestra) que juntos proporcionan los datos de intensidad de cada píxel.

plt.figure(figsize=(5, 5))
plt.imshow(img)

A continuación, me gustaría aplicar la escala y la rotación anteriores al dominio espacial de los datos de la imagen, transformando así las ubicaciones de los píxeles de forma similar a lo que demostré anteriormente con los datos de puntos. Sin embargo, necesito adoptar un enfoque ligeramente diferente porque los datos de la imagen están organizados de una forma distinta a la de las filas de puntos de datos con las que trabajé anteriormente. Con los datos de la imagen tengo que asignar los índices de cada píxel de los datos de entrada a los índices de salida transformados utilizando la matriz de transformación T, definida anteriormente.

# 2x scaling requires a transformation image array 2x the original image
img_transformed = np.empty((2000, 2000, 4), dtype=np.uint8)
for i, row in enumerate(img):
    for j, col in enumerate(row):
        pixel_data = img[i, j, :]
        input_coords = np.array([i, j, 1])
        i_out, j_out, _ = T @ input_coords
        img_transformed[i_out, j_out, :] = pixel_data

plt.figure(figsize=(5, 5))
plt.imshow(img_transformed)

La representación de la imagen tras aplicar la transformación afín muestra claramente que la imagen original se ha girado 90 grados en el sentido de las agujas del reloj y se ha escalado 2X. Sin embargo, el resultado es ahora obviamente menor, ya que se puede ver fácilmente una discontinuidad en las intensidades de los píxeles.

Para entender la razón de esto, volveré a utilizar un simple gráfico de cuadrícula como demostración. Consideremos un gráfico de 4 cuadrados en una cuadrícula de 2x2 similar al dominio espacial de una imagen de 2x2.

def plot_box(plt, x0, y0, txt, w=1, h=1):
    plt.scatter(x0, y0)
    plt.scatter(x0, y0 + h)
    plt.scatter(x0 + w, y0 + h)
    plt.scatter(x0 + w, y0)
    plt.plot([x0, x0, x0 + w, x0 + w, x0], [y0, y0 + h, y0 + h, y0, y0], color="gray", linestyle='dotted')
    plt.text(x0 + (.33 * w), y0 + (.5 * h), txt)

#             x0, y0, letter
a = np.array((0,  1,  0))
b = np.array((1,  1,  1))
c = np.array((0,  0,  2))
d = np.array((1,  0,  3))

A = np.array([a, b, c, d])
fig = plt.figure()
ax = plt.gca()
for pt in A:
    x0, y0, i = I @ pt
    x0, y0, i = int(x0), int(y0), int(i)
    plot_box(plt, x0, y0, f"{string.ascii_letters[int(i)]} ({x0}, {y0})")

ax.set_xticks(np.arange(-1, 5, 1))
ax.set_yticks(np.arange(-1, 5, 1))
plt.grid()
plt.show()

Ahora observa lo que ocurre cuando aplico una transformación de escala 2X como se muestra a continuación. Recordemos que

Te darás cuenta de que una transformación espacial de este tipo produce... bueno, "huecos" para decirlo en términos sencillos, que he hecho evidentes trazando signos de interrogación junto a las coordenadas. La cuadrícula de 2x2 se transforma en una de 3x3, y los cuadrados originales se reposicionan en función de la transformación lineal aplicada. Esto significa que (0,0) * Ts sigue siendo (0,0) debido a sus propiedades como vector 0, pero todos los demás se escalan por dos, como por ejemplo (1,1) * Ts -> (2,2).

fig = plt.figure()
ax = plt.gca()
for pt in A:
    xt, yt, i = T_s @ pt
    xt, yt, i = int(xt), int(yt), int(i)
    plot_box(plt, xt, yt, f"{string.ascii_letters[i]}' ({xt}, {yt})")

delta_w, delta_h = 0.33, 0.5
plt.text(0 + delta_w, 1 + delta_h, "? (0, 1)")
plt.text(1 + delta_w, 0 + delta_h, "? (1, 0)")
plt.text(1 + delta_w, 1 + delta_h, "? (1, 1)")
plt.text(1 + delta_w, 2 + delta_h, "? (1, 2)")
plt.text(2 + delta_w, 1 + delta_h, "? (2, 1)")

ax.set_xticks(np.arange(-1, 5, 1))
ax.set_yticks(np.arange(-1, 5, 1))
plt.grid()
plt.show()

La cuestión sigue siendo qué hacer con esos huecos que se han introducido. Un pensamiento intuitivo sería simplemente buscar la respuesta en la imagen original. Sucede que si aplicamos la inversa de la transformación a una coordenada en la salida obtendré el lugar correspondiente de la entrada original.

En operaciones matriciales como el mapeo inverso se ve así:

donde x', y' son las coordenadas en la cuadrícula 3x3 transformada anterior, concretamente el lugar que falta, como (2, 1), T-1s (los valores reales se muestran a continuación) es la inversa de la matriz de escala 2x Ts y x, y son las coordenadas que se encuentran en la cuadrícula original 2x2.

Sin embargo, pronto te darás cuenta de que hay un pequeño problema que todavía necesita ser resuelto debido al hecho de que cada una de las coordenadas de la brecha se mapean de nuevo a valores fraccionarios del sistema de coordenadas 2x2. En el caso de los datos de imagen no se puede tener una fracción de píxel. Esto será más claro con un ejemplo de mapeo de la brecha (2, 1) de nuevo a la original 2x2 espacio, así:

En este caso redondearé y' = 1/2 a 0 y diré que corresponde a (1, 0). En general, este método de seleccionar un valor en la cuadrícula original de 2x2 para colocarlo en los huecos de la cuadrícula transformada de 3x3 se conoce como interpolación, y en este ejemplo concreto estoy utilizando una versión simplificada del método de interpolación del vecino más próximo.

Bien, volvamos a los datos de la imagen. Debería estar bastante claro lo que hay que hacer ahora para arreglar esos huecos en la versión escalada y rotada de la letra R. Debo desarrollar una implementación de la interpolación del vecino más cercano basada en el mapeado hacia atrás, usando la inversa de la matriz de transformación T, de las coordenadas de los píxeles en la imagen transformada para encontrar la coincidencia exacta o el vecino más cercano en la imagen original.

T_inv = np.linalg.inv(T)

# nearest neighbors interpolation
def nearest_neighbors(i, j, M, T_inv):
    x_max, y_max = M.shape[0] - 1, M.shape[1] - 1
    x, y, _ = T_inv @ np.array([i, j, 1])
    if np.floor(x) == x and np.floor(y) == y:
        x, y = int(x), int(y)
        return M[x, y]
    if np.abs(np.floor(x) - x) < np.abs(np.ceil(x) - x):
        x = int(np.floor(x))
    else:
        x = int(np.ceil(x))
    if np.abs(np.floor(y) - y) < np.abs(np.ceil(y) - y):
        y = int(np.floor(y))
    else:
        y = int(np.ceil(y))
    if x > x_max:
        x = x_max
    if y > y_max:
        y = y_max
    return M[x, y,]

img_nn = np.empty((2000, 2000, 4), dtype=np.uint8)
for i, row in enumerate(img_transformed):
    for j, col in enumerate(row):
        img_nn[i, j, :] = nearest_neighbors(i, j, img, T_inv)

plt.figure(figsize=(5, 5))
plt.imshow(img_nn)

No está mal, ¿verdad?

Debo señalar que, en la mayoría de los casos, el método del vecino más próximo no será suficiente. Hay otros dos métodos de interpolación más comunes conocidos como bilineal y bicúbica que generalmente proporcionan resultados mucho mejores. Hablaré más sobre estos otros algoritmos de interpolación cuando presente las librerías Pillow y OpenCV en secciones posteriores. El propósito de esta sección es sólo para construir una comprensión intuitiva de cómo funcionan las cosas.

Conclusión

En esta segunda parte hemos seguido con más ejemplos de cómo se usan las transformaciones afines en Python en el caso del tratamiento de imágenes con un ejemplo sencillo que seguía lo que habíamos aprendido en nuestra parte 1 del post.

En la parte 3 del post podrás ver ejemplos prácticos de lo que hemos aprendido en esta parte acerca de cómo usar transformaciones afines en python con las librerías OpenCV y Pillow, así que ¡no te lo pierdas!

¡Nos vemos en clase!