domingo, 8 de octubre de 2017

Clustering de imágenes en espacio HSL con R

El procesado de imágenes capturadas por satélite resulta imprescindible para extraer la máxima cantidad de información de las mismas. Vamos a practicar el clustering en R sobre dos fotografías del satélite Pléiades 1, una obtenida sobre La Manga del Mar Menor en diciembre de 2016 y otra de los efectos del incendio que tuvo lugar en La Junquera en julio de 2012.



Dado que las fotografías escogidas se explican por sí mismas, no vamos a descubrir nada nuevo sobre ellas con este análisis. Sin embargo podremos comprobar el rendimiento del algoritmo K-means de R, contrastando el resultado del mismo con lo que resulta obvio inspeccionando las imágenes.

La primera opción es cargar las imágenes como un array de tres canales RGB y analizarlos para tratar de encontrar información discriminante. No obstante he preferido hacer una transformación previa de las coordenadas RGB a un modelo HSL (Hue, Saturation, Lightness), que proporciona tres nuevos ejes más intuitivos sobre los que trabajar: el Tono o Matiz (H), la Saturación (S) y la Luminosidad (L).

Fuente: Wikipedia


La Saturación y la Luminosidad toman valores normalizados entre 0 y 1. El Matiz en cambio sigue una escala cíclica expresada normalmente en grados entre 0 y 360 (360º corresponde al mismo tono que 0º, un rojo puro). Este carácter periódico del Matiz nos causará alguna dificultad como veremos.

Fuente: Adobe


HSL es una transformación clásica en colorimetría, pero ampliamente superada por modelos posteriores más sofisticados desde el punto de vista perceptual (p.ej. CIELAB). Aún así en los dos ejemplos estudiados se demuestra que el clustering resulta más efectivo moviéndonos en estas coordenadas que en el espacio RGB origen de las imágenes.

En concreto la imagen de La Manga presenta una fuerte discriminación por Matiz, por lo que desechando las demás variables empezamos aplicando sobre H el algoritmo K-means con 3 clústers, los que visualmente distinguimos en la imagen: "mar", "tierra" y "laguna".

kmeans(subset(M, select = c("H")), 3, nstart=10, iter.max=10)

El resultado obtenido arroja un parámetro de bondad bastante alto: between_SS / total_SS = 89.1 %, pero la imagen deja patente que usando solo la variable Matiz tres clústers resultan insuficientes ya que se interpretan partes de la laguna como tierra, y algunas zonas terrestres como mar.



Logramos discriminar esas áreas incrementando a 5 el número de clústers, que se corresponden ahora con: "mar", "tierra", "laguna1", "laguna2" y "resto tierra".



Sin ser perfecto, vemos que el clustering ahora sí diferencia muy bien los accidentes geográficos. Con cinco agrupaciones el parámetro de bondad lógicamente mejora: between_SS / total_SS = 97.3 %.

La información de Matiz píxel a píxel usada en el cálculo tiene el siguiente aspecto (se muestra con S = 1 y L = 0,5):



Los cinco centroides caen en los siguientes valores del histograma de frecuencias de Matiz:



Podemos deducir porqué tres o incluso cuatro grupos resultan insuficientes:
  • Por un lado el clúster con centroide H=82º (franja sur de la laguna o "laguna2"), queda más cerca en términos de Matiz del clúster H=24º ("tierra") que del clúster H=144º ("laguna1"). Es el motivo de que gran parte de esa zona de la laguna quedara asimilada al clúster "tierra" cuando se hacían solo tres conjuntos.
  • Pero además nos juega una mala pasada el carácter cíclico de H, que ha forzado introducir el clúster ficticio H=344º ("resto tierra"). A ojos de K-means este clúster es el que está a mayor distancia euclídea de H=24º ("tierra"), cuando por Matiz son clústers muy próximos; de hecho en el histograma queda claro que "resto tierra" no es sino la cola izquierda de "tierra". Al usar tres clústers todas las zonas de "resto tierra" quedaban indeseablemente identificadas como "mar". Desplazando cíclicamente el Matiz unos 85º (360 - (250+300)/2) nos habríamos ahorrado el clúster "resto tierra".

Como segundo ejemplo vamos a tratar de hacer un clustering que ayude a clasificar las zonas de bosque quemadas en julio de 2016 en La Junquera.

© Infoterra Servicios de Geoinformación S.A.


En este caso un clustering basado en el Matiz funciona realmente mal. Aquí se muestra el resultado con los 3 clústers que se perciben al inspeccionar la imagen: "bosque", "bosque incendiado" y "poblaciones", y resaltando en rojo el correspondiente a las zonas arrasadas:



Las zonas quemadas se confunden con las de bosque por lo que en este caso la variable Matiz por sí sola no aporta suficiente información. Mirando la imagen original podemos preguntarnos si esas zonas quemadas pudieran compartir un muy bajo parámetro de Saturación de color (se observan fuertemente grisáceas). Repitiendo el clustering de tres grupos, esta vez sobre la variable S en solitario:



Funciona bastante mejor en la tarea de localizar tierra quemada, pero aún se "equivoca" clasificando como calcinadas amplias áreas no quemadas con poca saturación de color (zonas de monte desplobadas de árboles). Vamos ahora a combinar la aportación de las dos variables, Saturación y Matiz, sin cambiar el número de clústers:

kmeans(subset(M, select = c("S", "H")), 3, nstart=10, iter.max=10)



Esta vez sí tenemos un buen modelado de las áreas arrasadas por el incendio y solo de ellas (salvo las carreteras, seguramente por su color tan neutro, y alguna otra zona puntual).

Como dato curioso este último clustering es el que peor dato de ajuste arroja de los tres: solo Matiz between_SS / total_SS = 86.4 %, solo Saturación between_SS / total_SS = 90.4 %, Matiz + Saturación between_SS / total_SS = 80.4 %. Asumo que el motivo es la inclusión de un mayor número de variables, haciendo más exigente la pertenencia a un clúster.

La ausencia de contenido cercano a H=0º/360º hace que la periodicidad de H no plantee problemas en esta imagen, como podemos ver en el scatterplot de Matiz y Saturación de la fotografía original. Identificamos además los tres clústers y sus centroides.



Resulta interesante ver el efecto del fuego en el color: las zonas verdes (C3) no cambian excesivamente su Matiz tras ser arrasadas por el incendio (C1), lo que sufren principalmente es una fuerte pérdida de Saturación. Los tres clústers con los colores tipo de cada centroide quedan así (se muestran con L = 0,5):



Nos damos cuenta de que el clúster C2 ("poblaciones"), siendo el más pequeño en superficie abarcada, es también el más heterogéneo por la variedad de Matices y Saturaciones que contiene al incluir tanto edificaciones como campos de cultivo. Sería por tanto el más susceptible de ser dividido en subclústers.

En clusteringincendio.gif se puede ver en forma de animación la diferencia entre los tres cluster realizados, y aquí puede verse el resultado del clustering final en alta resolución y falso color.

En la web de Infoterra hay también una imagen de la zona en espectro infrarrojo, en la que se ven con claridad las zonas afectadas por el incendio.

© Infoterra Servicios de Geoinformación S.A.


El código R para realizar el clustering, incluyendo las fórmulas vectorizadas de conversión RGB a HSL, puede encontrarse en clusteringHSL.R.

No hay comentarios:

Publicar un comentario

Por claridad del blog, por favor trata de utilizar una sintaxis lo más correcta posible y no abusar del uso de emoticonos, mayúsculas y similares.