domingo, 8 de octubre de 2017

Clustering de imágenes en espacio HSL con R

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.

Las fotografías escogidas se explican por sí mismas así que no vamos a descubrir nada nuevo sobre ellas con el 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, lo que constituye una especie de K-means "supervisado".

La primera opción es cargar las imágenes como un array de tres canales RGB. 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: 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 el clustering resulta más efectivo moviéndonos en HSL que en el espacio RGB origen.

~~~

La imagen de La Manga presenta una fuerte discriminación por Matiz, así que desechando las demás variables aplicamos el algoritmo K-means sobre H. Empezaremos con los 3 clústers que se distinguen a simple vista: "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 considerables de la laguna como tierra, y bastantes zonas terrestres como mar.



Logramos discriminar esas áreas incrementando a 5 el número de clústers, que se corresponderían ahora con: "mar", "tierra1", "tierra2", "laguna1" y "laguna2".



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 usada en el cálculo tiene el siguiente aspecto (se muestra con S = 1 y L = 0,5):



Los 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º ("tierra1") 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º ("tierra2"). A ojos de K-means este clúster es el que está a mayor distancia euclídea de H=24º ("tierra1"), cuando por Matiz son clústers muy próximos; de hecho en el histograma queda claro que "tierra2" no es sino la cola izquierda de "tierra1". Al usar tres clústers todas las zonas de "tierra2" quedaban indeseablemente identificadas como "mar". Desplazando cíclicamente el Matiz unos 85º (360 - (250+300)/2) nos habríamos ahorrado el clúster "tierra2".

~~~

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.

Fuente: Infoterra Servicios de Geoinformación S.A.


A priori volvemos a identificar 3 clústers: "bosque", "bosque incendiado" y "poblaciones". Un clustering basado en el Matiz funciona realmente mal, las zonas quemadas se confunden con las de bosque por lo que concluimos que la variable Matiz por sí sola no aporta suficiente información.

Sobre la imagen original nos preguntarnos si las zonas quemadas pudieran compartir un bajo parámetro de Saturación (se observan fuertemente grisáceas). Repetimos el clustering de tres grupos sobre la variable S en solitario, y aunque funciona bastante mejor en la tarea de localizar tierra quemada aún se "equivoca" clasificando como calcinadas amplias áreas no quemadas con poca saturación de color (zonas de monte desplobadas de árboles).

En un tercer intento combinamos la aportación de las dos variables, Saturación y Matiz, sin cambiar el número de clústers. 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).

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

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 siguiente comparación muestra en rojo las áreas clasificadas como "bosque incendiado" en cada una de las ejecuciones.



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. Esto explica porqué el Matiz era una variable insuficiente por sí misma para hacer el clustering.

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.

Puede verse aquí 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.

Fuente: 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.