Extracción de características de imagen: Patrones Binarios Locales con Cython
Introducción
El objetivo común de la extracción de características es representar los datos crudos como un conjunto reducido de características que describan mejor sus atributos principales [1]. De esta forma, podemos reducir la dimensionalidad de la entrada original y usar las nuevas características como entrada para entrenar técnicas de reconocimiento de patrones y clasificación.
Aunque hay varias características que podemos extraer de una imagen, los Patrones Binarios Locales (LBP) son un enfoque teóricamente sencillo, pero eficiente, para la clasificación de texturas en escala de grises e invariante a la rotación. Funcionan porque los patrones más frecuentes corresponden a microcaracterísticas primitivas como bordes, esquinas, puntos y regiones planas [2].
En [2], Ojala et al. mostraron que el histograma de ocurrencias discretas de los patrones uniformes es una característica de textura muy poderosa. La textura de una imagen se define como un fenómeno bidimensional caracterizado por dos propiedades: (1) estructura espacial (patrón) y (2) contraste.
Metodología
Conjunto de vecinos circularmente simétrico
Un conjunto de vecinos circularmente simétrico para un píxel dado gc se define por los puntos con coordenadas (i, j) que rodean al punto central en un círculo de radio R, con un número de elementos P.
Textura
Definimos una textura T como la colección de píxeles en una imagen en escala de grises
donde gp corresponde al valor de gris del vecino local p.
Interpolación
Cuando un vecino no se encuentra en el centro de un píxel, ese valor de gris del vecino debe calcularse por interpolación. Por lo tanto, necesitamos definir una función que, dada una coordenada, devuelva el valor de gris interpolado.
Logrando invarianza a la escala de grises
Considerando una posible pérdida de información, es posible convertir la textura en la diferencia conjunta. Para calcularla, restamos el valor de gris del píxel central a todo el conjunto de vecinos. La distribución de diferencias conjuntas es un operador de textura altamente discriminativo. Registra las ocurrencias de varios patrones en el vecindario de cada píxel en un histograma P-dimensional.
donde gp es el valor de gris del vecino p. Esta distribución es invariante ante desplazamientos de escala de grises.
Patrón Binario Local
El operador LBP_{P,R} es, por definición, invariante ante cualquier transformación monótona de la escala de grises. Mientras el orden de los valores de gris se mantenga igual, la salida del operador LBP_{P,R} permanece constante.
donde
Patrones Binarios Locales Uniformes
En [2], Ojala menciona que, en su experiencia práctica, LBP no es un buen discriminador. Proponen seleccionar únicamente el conjunto de patrones binarios locales tal que el número de transiciones espaciales (cambios de bit 0/1) no supere 2. Por ejemplo, el patrón ‘1111’ tiene 0 transiciones espaciales, el patrón ‘1100’ tiene 1 transición espacial y el patrón ‘1101’ tiene 2 transiciones espaciales. A cada patrón uniforme se le asocia un índice único. La fórmula para crear el índice se tomó de aquí.
Ahora podemos calcular los patrones binarios locales para un píxel central. El siguiente paso es calcular los patrones binarios locales para todos los píxeles.
Pista: Por simplicidad, no estoy considerando el caso en que un índice seleccionado sea negativo (es decir, img_gray[-1][0] devuelve el último píxel de la primera columna). Si quisiéramos tener un cálculo más preciso, deberíamos considerar este caso y tratarlo.
Código en Cython
El código anterior no es perfecto; sin embargo, lo que realmente lo hace lento es que iteramos sobre todos los píxeles de la imagen. Esperar 1 minuto y 10 segundos para calcular nuestras características es mucho si tomamos en cuenta que además tenemos que entrenar una técnica de reconocimiento de patrones. Por lo tanto, necesitamos una implementación alternativa que sea mucho más rápida para los bucles. En este caso, usaremos Cython. El código se presenta en la siguiente imagen; es un bloque de código grande. Algunas partes podrían mejorarse, pero ya es mucho más rápido. Siéntete libre de dejar comentarios si no entiendes algo del código.
El código está escrito de tal forma que la mayor parte se ejecuta enteramente en la API de C. Esta estrategia acelera considerablemente la ejecución, pero también nos permite aprovechar el módulo paralelo de Cython. Vamos a repartir el trabajo entre varios núcleos de la CPU.
from libc.math cimport sin, cos, pi, ceil, floor, pow
from libc.stdlib cimport abort, malloc, free
import numpy as np
cimport numpy as np
cimport cython
from cython.parallel import prange, parallel
cimport openmp
cdef double get_pixel2d(
double *image,
Py_ssize_t n_rows,
Py_ssize_t n_cols,
long x,
long y) nogil:
if (y < 0) or (y >= n_rows) or (x < 0) or (x >= n_cols):
return 0
else:
return image[y * n_cols + x]
cdef double bilinear_interpolation(
double *image,
Py_ssize_t n_rows,
Py_ssize_t n_cols,
double x,
double y) nogil:
cdef double d_y, d_x, top_left, top_right, bottom_left, bottom_right
cdef long min_y, min_x, max_y, max_x
min_y = <long>floor(y)
min_x = <long>floor(x)
max_y = <long>ceil(y)
max_x = <long>ceil(x)
d_y = y - min_y
d_x = x - min_x
top_left = get_pixel2d(image, n_rows, n_cols, min_x, min_y)
top_right = get_pixel2d(image, n_rows, n_cols, max_x, min_y)
bottom_left = get_pixel2d(image, n_rows, n_cols, min_x, max_y)
bottom_right = get_pixel2d(image, n_rows, n_cols, max_x, max_y)
top = (1 - d_x) * top_left + d_x * top_right
bottom = (1 - d_x) * bottom_left + d_x * bottom_right
return (1 - d_y) * top + d_y * bottom
cdef double *joint_difference_distribution(
double *image,
Py_ssize_t n_rows,
Py_ssize_t n_cols,
int x0,
int y0,
int P,
int R
) nogil:
cdef Py_ssize_t p
cdef double *T = <double *> malloc(sizeof(double) * P)
cdef double x, y, gp, gc
if T is NULL:
abort()
gc = get_pixel2d(image, n_rows, n_cols, x0, y0)
for p in range(P):
x = x0 + R * cos(2 * pi * p / P)
y = y0 - R * sin(2 * pi * p / P)
gp = bilinear_interpolation(image, n_rows, n_cols, x, y)
T[p] = gp - gc
return T
cdef int *binary_joint_distribution(double *T, Py_ssize_t T_size) nogil:
cdef int *s_T = <int *> malloc(sizeof(int) * T_size)
cdef Py_ssize_t i = 0
for t in range(T_size):
if T[t] >= 0.0:
s_T[t] = 1
else:
s_T[t] = 0
return s_T
cdef long LBP(double *T, int *s_T, Py_ssize_t T_size) nogil:
cdef long LBP_pr = 0
cdef Py_ssize_t i = 0
for i in range(0, T_size):
LBP_pr = LBP_pr + 2 ** i * s_T[i]
return LBP_pr
cdef int is_uniform_pattern(int *s_T, Py_ssize_t s_T_size) nogil:
cdef Py_ssize_t i = 0
cdef int counter = 0
for i in range(s_T_size - 1):
if s_T[i] != s_T[i + 1]:
counter += 1
if counter > 2:
return 0
return 1
cdef int create_index(int *s_T, Py_ssize_t s_T_size) nogil:
cdef int n_ones = 0
cdef int rot_index = -1
cdef int first_one = -1
cdef int first_zero = -1
cdef int lbp = -1
cdef Py_ssize_t i
for i in range(s_T_size):
if s_T[i]:
n_ones += 1
if first_one == -1:
first_one = i
else:
if first_zero == -1:
first_zero = i
if n_ones == 0:
lbp = 0
elif n_ones == s_T_size:
lbp = s_T_size * (s_T_size - 1) + 1
else:
if first_one == 0:
rot_index = n_ones - first_zero
else:
rot_index = s_T_size - first_one
lbp = 1 + (n_ones - 1) * s_T_size + rot_index
return lbp
cdef int LBP_uniform(int *s_T, Py_ssize_t s_T_size) nogil:
cdef int LBP_pru = 0
cdef Py_ssize_t i = 0
if is_uniform_pattern(s_T, s_T_size):
LBP_pru = create_index(s_T, s_T_size)
else:
LBP_pru = 2 + s_T_size * (s_T_size - 1)
return LBP_pru
@cython.boundscheck(False)
@cython.wraparound(False)
def local_binary_patterns(
double[:, ::1] image,
int P,
int R,
int num_threads=1
):
cdef Py_ssize_t x = 0
cdef Py_ssize_t y = 0
cdef int n_rows = image.shape[0]
cdef int n_cols = image.shape[1]
cdef int[:, ::1] lbp = np.zeros([n_rows, n_cols], dtype=np.int32)
with nogil, parallel(num_threads=num_threads):
for y in prange(n_rows, schedule='static'):
for x in prange(n_cols, schedule='static'):
T = joint_difference_distribution(&image[0][0], n_rows, n_cols, x, y, P, R)
s_T = binary_joint_distribution(T, P)
lbp[y, x] = LBP_uniform(s_T, P)
return np.asarray(lbp)Usando 4 hilos, pudimos calcular los patrones binarios locales para todos los píxeles en menos de 150 ms. Es tanto más rápido que ni siquiera me voy a molestar en calcular cuántas veces.
Comparación con una imagen similar
Tomemos otra imagen de ladrillos, pero esta tendrá una textura distinta.
Ambos histogramas son muy similares, y deberían serlo, ya que al final ambos son ladrillos. Sin embargo, las características del 20 al 40 son muy distintas en ambas imágenes. Esto significa que, con un buen algoritmo de machine learning, podríamos clasificarlas correctamente.
Conclusión
Los patrones binarios locales son características simples pero eficientes. La teoría detrás de ellos no es difícil de entender y son fáciles de programar. Sin embargo, si los programamos enteramente en Python, tendremos algunos problemas de rendimiento. Abordamos el problema con Cython y obtuvimos resultados muy impresionantes. El siguiente paso es recolectar distintas imágenes de textura y entrenar tu algoritmo de machine learning favorito para clasificarlas.
Notebook de Jupyter
Bibliografía
[1] Marques, O. (2011). Practical image and video processing using MATLAB. John Wiley & Sons.
[2] Ojala, T., Pietikäinen, M., & Mäenpää, T. (2002). Multiresolution gray-scale and rotation invariant texture classification with local binary patterns. IEEE Transactions on Pattern Analysis and Machine Intelligence, 24(7), 971–987.