MUI_IIMAI_1415_Analisis_Imagenes

1016 days ago by trinidad1415

Análisis de Imágenes

Introducción

Comencemos recordando que una imagen no es más que una matriz tal que en cada posición tiene almacenado un color (bien como un número o como vector de tres números). Podemos obtener la imagen de una matriz con matrix_plot.

M=Matrix([[1,0,1],[1,1,0],[0,1,0]]);show(M) 
       
matrix_plot(M) 
       

Importando herramientas e imágenes

De modo nativo, Sage no tiene muchos métodos para analizar imágenes, por lo que tendremos que importar algunos paquetes de python.

from numpy import * import scipy.ndimage import matplotlib.image import matplotlib.pyplot 
       

Ahora ya podemos leer una imagen usando matplotlib.image.imread. Esto nos devuelve una imagen, es decir una matriz de nxm con 3 valores en cada posición. La función no nos devuelve una matriz del tipo Matrix de Sage, sino una matriz del tipo array de numpy, por lo que tendremos que usar algunos métodos especiales. Siempre lo podríamos pasar a una lista con el método to_list(), pero entonces es más lento trabajar con ellas.

img=matplotlib.image.imread(DATA+'london.png') matrix_plot(img) 
       

Un array se utiliza de modo análogo a una matriz de Sage. Por ejemplo, las siguientes líneas seleccionan una parte de la imagen y la guardan en otra matriz.

img2=img[150:300,600:800,:] matrix_plot(img) 
       

Análisis pixel a pixel

En cada posición (i,j), la imagen tiene tres valores, correspondientes a los tonos de RGB.

img[0,0] 
       

Podemos dibujar cada uno de los canales de la imagen por separado.

pR=matrix_plot(img[:,:,0]) pG=matrix_plot(img[:,:,1]) pB=matrix_plot(img[:,:,2]) (pR,pG,pB) 
       

Histograma

Vamos a generar el histograma de cada uno de los canales de la imagen. Recordemos que se trata de, para cada valor posible del color (son 256 números entre 0 y 1), contar su frecuencia de aparición.

def histograma(imagenMonocroma): Rl=flatten(imagenMonocroma) escalas = uniq(Rl) n=len(Rl) return [RDF(Rl.count(x))/RDF(n) for x in escalas] 
       

Ahora podemos dibujar el histograma del canal rojo. Para ello convertimos el canal rojo en una lista (realmente en una lista de filas, donde cada fila es una lista de valores) y podemos aplicar la función anterior.

histRojo=histograma(img[:,:,0].tolist()) list_plot(zip(range(256),histRojo),plotjoined=True, color='red') 
       

El método anterior es relativamente lento, por lo que es mejor usar la función matplotlib.pyplot.hist, que devuelve en particular los valores del histograma.

hR=matplotlib.pyplot.hist(img[:,:,0].flatten(),256)[0] hG=matplotlib.pyplot.hist(img[:,:,1].flatten(),256)[0] hB=matplotlib.pyplot.hist(img[:,:,2].flatten(),256)[0] 
       

Podemos representarlos con list_plot.

rx=[0.0,(1.0/255.0)..1] (list_plot(zip(rx,hR),plotjoined=True, color='red',legend_label='Canal rojo')+list_plot(zip(rx,hG),plotjoined=True, color='green',legend_label='Canal verde')+list_plot(zip(rx,hB),plotjoined=True, color='blue',legend_label='Canal azul')).show(ymax=8500) 
       

Dibujaremos también el de la selección que hemos hecho para comparar.

hR2=matplotlib.pyplot.hist(img2[:,:,0].flatten(),256)[0] hG2=matplotlib.pyplot.hist(img2[:,:,1].flatten(),256)[0] hB2=matplotlib.pyplot.hist(img2[:,:,2].flatten(),256)[0] (list_plot(zip(rx,hR2),plotjoined=True, color='red',legend_label='Canal rojo')+list_plot(zip(rx,hG2),plotjoined=True, color='green',legend_label='Canal verde')+list_plot(zip(rx,hB2),plotjoined=True, color='blue',legend_label='Canal azul')).show(ymax=900) 
       

Filtrando a partir de un umbral (thresholding)

Podemos ver que en la selección el canal azul tiene valores relativamente pequeños y que los edificios de la imagen tienen todos un tono azulado. Podemos usar esa información para obtener aplicar un filtro que nos selecciones todos los edificios.

umbralazul=90.0/256. BT=img[:,:,2]>umbralazul matrix_plot(BT) 
       

Podemos usar la matriz BT como una "máscara" de modo que obtengamos una imagen en la que sólo se vea la parte seleccionada.

imgB=array(zeros(img.shape)) imgB[BT]=img[BT] matrix_plot(imgB) 
       

Sin embargo, podemos observar que un filtro del mismo tipo sobre el canal verde para seleccionar las zonas verdes no funciona.

umbralverde=0.5 GT=img[:,:,1]>umbralverde matrix_plot(GT) 
       

El espacio de color HSV

Para obtener un filtro bueno del color verde, usaremos un espacio del color en el que se pueda filtrar por tono. En este caso hemos elegido el espacio HSV. Tenemos que importar el paquete que nos realiza las conversiones de un espacio de color a otro.

import colorsys 
       

Crearemos una matriz que contenga los valores de Hue de la imagen. Para ello creamos una matriz de tipo array y guardamos en cada posición el valor de Hue obtenido de aplicar colorsys.rgb_to_hsv al color rgb original.

H=array(zeros(img.shape[0:2])) for i in range(H.shape[0]): for j in range(H.shape[1]): H[i,j]=colorsys.rgb_to_hsv(img[i,j,0],img[i,j,1],img[i,j,2])[0] 
       

Hacemos lo mismo para la selección, a fin de comparar.

H2=array(zeros(img2.shape[0:2])) for i in range(H2.shape[0]): for j in range(H2.shape[1]): H2[i,j]=colorsys.rgb_to_hsv(img2[i,j,0],img2[i,j,1],img2[i,j,2])[0] 
       

Ahora podemos calcular los histogramas y representarlos. Nótese que se ha dividido entre (1.0*H.shape[0]*H.shape[1]) para obtener los histogramas en frecuencia de aparición (%).

hH=matplotlib.pyplot.hist(H[:,:].flatten(),256)[0]/(1.0*H.shape[0]*H.shape[1]) hH2=matplotlib.pyplot.hist(H2[:,:].flatten(),256)[0]/(1.0*H2.shape[0]*H2.shape[1]) list_plot(zip(rx,hH),plotjoined=True, color='blue',legend_label='Hue imagen')+list_plot(zip(rx,hH2),plotjoined=True, color='green',legend_label='Hue bosque') 
       

Como parece que la mayor parte de la información está entre 0.2 y 0.3, filtramos usando esos valores.

GT=H>0.2 GT[H>0.3]=0 matrix_plot(GT) 
       

Y la utilizamos como máscara.

imgG=array(zeros(img.shape)) imgG[GT]=img[GT] matrix_plot(imgG) 
       

Podemos observar multitud de puntos blancos o grises. Esto se debe a que, cerca de los tonos grises, la crominancia no discrimina bien los tonos. Podemos filtran en términos de saturación o como parece en la imagen que son blancos lo que tenemos, en términos de la intensidad (definida como el valor medio de los tres canales).

Calculamos la intensidad y el histograma de la misma.

I=(img[:,:,0]+img[:,:,1]+img[:,:,2])/3 I2=(img2[:,:,0]+img2[:,:,1]+img2[:,:,2])/3 hI=matplotlib.pyplot.hist(I[:,:].flatten(),256)[0]/(1.0*I.shape[0]*I.shape[1]) hI2=matplotlib.pyplot.hist(I2[:,:].flatten(),256)[0]/(1.0*I2.shape[0]*I2.shape[1]) list_plot(zip(rx,hI),plotjoined=True, color='black',linestyle='--',legend_label='Intensidad imagen')+list_plot(zip(rx,hI2),plotjoined=True, color='gray',legend_label='Intensidad bosque') 
       

Filtramos a partir del valor 0.4.

GT[I>0.4]=0 matrix_plot(GT) 
       

Y lo utilizamos como máscara.

imgG=array(zeros(img.shape)) imgG[GT]=img[GT] matrix_plot(imgG) 
       

Análisis local de la imagen

A continuación utilizaremos las herramientas que no ven la imagen píxel a píxel, sino que consideran la geometría de la imagen.

Conexión

Para calcular las componentes conexas de la imagen utilizamos la función scipy.ndimage.label que recibe una imagen binaria (ceros y unos) y devuelve una imagen donde cada componente conexa tiene asignado un número (etiquetas) y la lista de etiquetas. Si no especificamos nada, considerará los 4-vecinos.

label_im, nb_labels = scipy.ndimage.label(BT) matrix_plot(label_im) 
       

Apertura y cierre

Podemos ver que la imagen tiene multitud de componentes conexas debido a que existen muchos puntos aislados resultado del filtro con umbral. Para eliminar puntos aislados y dar un aspecto más uniforme a la imagen, realizaremos el cierre (para eliminar los puntos blancos) y la apertura (para eliminar los puntos negros) de la imagen, usando las funciones scipy.ndimage.binary_closing y scipy.ndimage.binary_opening

close_img = scipy.ndimage.binary_closing(BT) open_img = scipy.ndimage.binary_opening(close_img) matrix_plot(BT),matrix_plot(open_img) 
       

Podemos de nuevo usar la imagen resultante como máscara para ver qué se ha seleccionado.

imgBB=array(zeros(img.shape)) imgBB[open_img]=img[open_img] matrix_plot(imgBB) 
       

Si ahora calculamos las componentes conexas, vemos que son bastantes menos y tienen una forma más cercana a la buscada.

label_im, nb_labels = scipy.ndimage.label(open_img) matrix_plot(sqrt(2*label_im), cmap='hot') 
       

Seleccionando el tamaño de las componentes conexas

Otro método para obtener menos componentes conexas es eliminar directamente las de menor tamaño. Vamos a calcular las componentes conexas del fitro que hicimos con los tonos verdes.

label_imG, nb_labels = scipy.ndimage.label(GT) matrix_plot(label_imG) 
       

Podemos obtener los tamaños de las componentes conexas sumando el número de posiciones de la imagen que tienen cierta etiqueta.

sizes = scipy.ndimage.sum(GT, label_imG, range(nb_labels + 1)) sizes 
       

Seleccionamos ahora todas las componentes conexas de menos de 1000 pixeles. De las etiquetas, seleccionamos aquellas que no cumplen ese criterio y las ponemos a 0 en la lista.

mask_size = sizes < 1000 remove_pixel = mask_size[label_imG] label_imG[remove_pixel] = 0 
       

Ahora recomponemos la imagen eligiendo los píxeles con las etiquetas que nos han quedado.

labels = unique(label_imG) label_clean = searchsorted(labels, label_imG) matrix_plot(label_clean) 
       

Cambiando los vecinos

Hasta ahora hemos utilizado 4-vecinos, que es lo que venía por defecto. Podemos cambiar los vecinos utilizando el parámetro structure, que recibe un matriz (de tipo array) de vecinos y la utiliza como modelo para escoger los vecinos. Para ello se superpone la matriz de vecinos a la matriz de la imagen, de modo que el píxel de la imagen que estemos estudiando estará en la posición central de la matriz de vecinos y se consideran como vecinos las posiciones marcadas con un 1.

close_imgG = scipy.ndimage.binary_closing(label_clean,structure=ones((5,5))) open_imgG = scipy.ndimage.binary_opening(close_imgG,structure=ones((3,3))) matrix_plot(open_imgG) 
       

Finalmente, utilizamos las regiones obtenidas de máscara para ver lo que se ha seleccionado.

imgGG=array(zeros(img.shape)) imgGG[open_imgG]=img[open_imgG] matrix_plot(imgGG) 
       

Ejercicios

1. Crear un filtro usando los canales H y S para seleccionar las zonas verdes

2. ¿Cuántas zonas verdes de área mayor de 4000 píxeles aparecen en la imagen?

3. A la imagen obtenida en el ejercicio 1, aplicarle la apertura y aplicar cierre al resultado considerando como vecinos los píxeles que disten 3 o menos unidades.

M=zeros([7,7]);M 
       
array([[ 0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.]])
array([[ 0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.]])
for i in range(7): for j in range(7): if (i-3)^2+(j-3)^2<=9: M[i,j]=1; M 
       
array([[ 0.,  0.,  0.,  1.,  0.,  0.,  0.],
       [ 0.,  1.,  1.,  1.,  1.,  1.,  0.],
       [ 0.,  1.,  1.,  1.,  1.,  1.,  0.],
       [ 1.,  1.,  1.,  1.,  1.,  1.,  1.],
       [ 0.,  1.,  1.,  1.,  1.,  1.,  0.],
       [ 0.,  1.,  1.,  1.,  1.,  1.,  0.],
       [ 0.,  0.,  0.,  1.,  0.,  0.,  0.]])
array([[ 0.,  0.,  0.,  1.,  0.,  0.,  0.],
       [ 0.,  1.,  1.,  1.,  1.,  1.,  0.],
       [ 0.,  1.,  1.,  1.,  1.,  1.,  0.],
       [ 1.,  1.,  1.,  1.,  1.,  1.,  1.],
       [ 0.,  1.,  1.,  1.,  1.,  1.,  0.],
       [ 0.,  1.,  1.,  1.,  1.,  1.,  0.],
       [ 0.,  0.,  0.,  1.,  0.,  0.,  0.]])