Segmentación semántica de imágenes de teledetección utilizando k-Means |  de Aleksei Rozanov |  marzo de 2024

Desde cero en Python🐍

Imagen por autor.

ohUno de los modelos de aprendizaje automático más simples y geniales, en mi opinión, es la agrupación de k-medias. Se relaciona con el grupo de algoritmos de aprendizaje no supervisados ​​y es capaz de encontrar patrones dentro de un conjunto de datos sin etiquetar. La característica más agradable es que carece de matemáticas complicadas y, básicamente, cualquier estudiante de secundaria puede implementar y utilizar este método con éxito. Entonces, en este artículo quiero compartir cómo puedes construir el algoritmo k-Means desde cero en Python usando solo engordado y pandas bibliotecas y aplicarlo a un problema del mundo real: la segmentación semántica de imágenes de satélite.

En primer lugar, hablemos de los datos que tenemos.

En uno de mis artículos anteriores hablé del problema de la contracción del Mar de Aral. Como resultado, logramos obtener imágenes de teledetección de MODIS utilizando Google Earth Engine, lo que indica claramente que el mar se está secando. Entonces me pregunté, ¿cómo podemos estimar el cambio de la superficie del agua entre 2000 y 2023 utilizando la segmentación semántica de ML? ¡La respuesta es k-Medios!

Antes de sumergirnos en la codificación, echemos un vistazo a los datos que usaremos en este tutorial. Son dos imágenes RGB de la misma zona con un intervalo de 23 años, sin embargo está claro que las propiedades de la superficie terrestre y las condiciones atmosféricas (nubes, aerosoles, etc.) son diferentes. Por eso decidí entrenar dos modelos k-Means separados, uno para cada imagen.

Imagen por autor.

Para seguir el tutorial, puede descargar y ejecutar el cuaderno. aquí.

En primer lugar, importemos las bibliotecas necesarias y carguemos los datos en el cuaderno:

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.image as mpimg

img = mpimg.imread('MOD_01.jpg')
img2 = mpimg.imread('MOD_24.jpg')

Puedes ver que el área cubierta por las imágenes es bastante grande, por lo que te sugiero hacer un poco de zoom:

img = img[140:600,110:500,:]
img2 = img2[140:600,110:500,:]

fig, ax = plt.subplots(ncols=2, figsize=(16,9))
ax[0].imshow(img)
ax[1].imshow(img2)
for i in range(2):
ax[i].set_facecolor('black')
ax[i].set_xticks([])
ax[i].set_yticks([])
ax[0].set_title('2000-08-01', fontsize=26)
ax[1].set_title('2023-08-01', fontsize=26)
plt.show()

Imagen por autor.

Y el último paso antes de la fase de ML, conviertamos nuestras imágenes a pandas marcos de datos (una columna para cada canal de imagen). Lo hago por el bien de la visibilidad de mis explicaciones. Si desea optimizarlo, es mejor utilizar engordado matrices en su lugar.

df = pd.DataFrame({'R': img[:,:, 0].flatten(), 'G': img[:,:, 1].flatten(), 'B':img[:,:, 2].flatten()})
df2 = pd.DataFrame({'R': img2[:,:, 0].flatten(), 'G': img2[:,:, 1].flatten(), 'B':img2[:,:, 2].flatten()})

Entonces, ¿cuál es la idea detrás del algoritmo?

Imagine que juzga el sabor de la comida utilizando dos criterios: dulzura y precio. Teniendo esto en cuenta, te daré un conjunto de posibles opciones para comer:

Imagen por autor.

Apuesto a que tu cerebro ya ha dividido las opciones en tres grupos: frutas, bebidas y panadería. Básicamente, inconscientemente agrupaste los datos bidimensionales, que están definidos por un par de valores (dulzura; precio).

Imagen por autor.

En el caso de k-mediasel objetivo del algoritmo es bastante similar: encontrar una Preestablecido número de grupos, ken un espacio n-dimensional (por ejemplo, además del dulzor y el precio, desea tener en cuenta la nutrición, la salud, la presencia de los alimentos en su refrigerador y, en este caso, n = 5).

I. Definir el número de conglomerados.

Como mencioné anteriormente, k en k-Means es la cantidad de grupos que desea obtener al final, y se supone que debe establecer este valor antes entrenar el modelo.

II. Inicializa aleatoriamente los centroides.

El centroide es una parte integral de k-Means. Básicamente, el centroide es un círculo con un centro, que se define como un conjunto de coordenadas, y cada centroide representa un grupo. Por ejemplo, en nuestro ejemplo anterior hay 3 centroides.

III. Calcula distancias y asigna grupos.

Ahora necesitamos encontrar qué tan lejos está cada punto de cada centroide. Con base en estos cálculos, asignamos cada punto al centroide (grupo) menos distante.

IV. Calcular nuevos centroides.

Ahora cada uno de nuestros grupos contiene al menos un punto, por lo que es hora de volver a calcular los centroides simplemente tomando las coordenadas medias de todos los puntos del grupo.

¡Y eso es! Repetimos los pasos 2 a 4 hasta que los centroides no cambien.

Imagen por autor.

Ahora incluyamos esta idea realmente simple detrás de k-Means en código Python.

Recordatorio: en esta tarea tenemos 3D problema, es decir, nuestro X, Y y z son Rojo, Verde y Azul canales de imagen!

def kmeans(data, K, kind):
L = list()
new_centroids = data.sample(K).values

data = distance(data.copy(), new_centroids, kind)
old_centroids = new_centroids.copy()
new_centroids = np.array([data[data.Class == Class][['R', 'G', 'B']].mean().values for Class in data.loc[:,'C1':f'C{K}'].columns])
i = 1
print(f'Iteration: {i}\tDistance: {abs(new_centroids.mean()-old_centroids.mean())}')
while abs(new_centroids.mean()-old_centroids.mean())>0.001:
L.append(abs(new_centroids.mean()-old_centroids.mean()))
data = distance(data, new_centroids, kind)
old_centroids = new_centroids.copy()
new_centroids = np.array([data[data.Class == Class][['R', 'G', 'B']].mean().values for Class in data.loc[:,'C1':f'C{K}'].columns])
i+=1
print(f'Iteration: {i}\tDistance: {abs(new_centroids.mean()-old_centroids.mean())}')
print(f"k-Means has ended with {i} iteratinons")
return data, L

En la primera etapa creamos una lista. l para recopilar todas las distancias entre grupos para visualizarlos después y muestrear aleatoriamente K puntos del conjunto de datos para convertirlos en nuestros centroides (o alternativamente, puede asignar valores aleatorios a los centroides).

L = list()
new_centroids = data.sample(K).values

Ahora necesitamos calcular las distancias entre los centroides y los puntos de datos. Hay muchas métricas de distancia diferentes en ciencia de datos, pero centrémonos en las siguientes: Euclidean, Manhattan, Chebyshev.

Para distancia euclidiana:

Imagen por autor.

Para Manhattan:

Imagen por autor.

Para Chebyshev:

Imagen por autor.

Para usar estas fórmulas, escribamos una función versátil para cualquier número de dimensiones:

def distance(data, centroids, kind):
#kind = euclidean, manhattan, chebyshev
#Here we add to the dataframe as many clusters C-ith as needed
cols=list()
for i in range(1,k+1):
if kind=='euclidean':
data[f'C{i}'] = ((centroids[i-1][0]-data.R)**2+(centroids[i-1][1]-data.G)**2+(centroids[i-1][2]-data.B)**2)**0.5
elif kind=='manhattan':
data[f'C{i}'] = abs(centroids[i-1][0]-data.R)+abs(centroids[i-1][1]-data.G)+abs(centroids[i-1][2]-data.B)
elif kind=='chebyshev':
merged=pd.concat([centroids[i-1][0]-data.R, centroids[i-1][1]-data.G, centroids[i-1][2]-data.B], axis=1)
data[f'C{i}'] = merged.max(axis=1)
cols.append(f'C{i}')
data['Class'] = data[cols].abs().idxmin(axis=1) #assigning clusters to points
return data #returning the dataframe with k cluster columns and one Class column with the final cluster

Ahora podemos simplemente calcular distancias y asignar un grupo a cada punto de datos. Por lo tanto, nuestros nuevos centroides se volvieron viejos, por lo que los almacenamos en otra variable y recalculamos los nuevos. Para hacer eso, iteramos sobre cada grupo y tomamos una media en todas las coordenadas (en nuestro caso, en los canales RGB). Por lo tanto, la variable new_centroids tiene la forma de (k,3).

data = distance(data.copy(), new_centroids, kind)
old_centroids = new_centroids.copy()
new_centroids = np.array([data[data.Class == Class][['R', 'G', 'B']].mean().values for Class in data.loc[:,'C1':f'C{K}'].columns])

Finalmente, repetimos todos estos pasos hasta que las coordenadas de los centroides ya no cambien. Expresé esta condición de la siguiente manera: la diferencia entre las coordenadas promedio del grupo debe ser inferior a 0,001. Pero puedes jugar con otros números aquí.

while abs(new_centroids.mean()-old_centroids.mean())>0.001:
L.append(abs(new_centroids.mean()-old_centroids.mean()))
data = distance(data, new_centroids, kind)
old_centroids = new_centroids.copy()
new_centroids = np.array([data[data.Class == Class][['R', 'G', 'B']].mean().values for Class in data.loc[:,'C1':f'C{K}'].columns])

Y eso es. ¡El algoritmo está listo para ser entrenado! Así que establezcamos k = 3 y almacenemos los resultados en diccionarios.

k = 3
segmented_1, segmented_2, distances_1, distances_2 = {}, {}, {}, {}
segmented_1['euclidean'], distances_1['euclidean'] = kmeans(df, k, 'euclidean')
segmented_2['euclidean'], distances_2['euclidean'] = kmeans(df2, k, 'euclidean')
segmented_1['manhattan'], distances_1['manhattan'] = kmeans(df, k, 'manhattan')
segmented_2['manhattan'], distances_2['manhattan'] = kmeans(df2, k, 'manhattan')
segmented_1['chebyshev'], distances_1['chebyshev'] = kmeans(df, k, 'chebyshev')
segmented_2['chebyshev'], distances_2['chebyshev'] = kmeans(df2, k, 'chebyshev')

Decidí comparar todas las métricas de distancia para esta tarea en particular, como puedes ver, y es evidente que aquí la distancia de Manhattan fue la más rápida.

Imagen por autor.

Antes de visualizar los clústeres, conviertamos los nombres de los clústeres al tipo int:

d = {'C1':0, 'C2': 1, 'C3':2}
for key in segmented_1.keys():
segmented_1[key].Class = segmented_1[key].Class.apply(lambda x: d[x])
segmented_2[key].Class = segmented_2[key].Class.apply(lambda x: d[x])

¡Es hora de hacer las tramas finales!

for key in segmented_1.keys():
fig, ax = plt.subplots(ncols=2, nrows=2, figsize=(10,10))
ax[0, 0].imshow(img)
ax[0, 1].imshow(segmented_1[key].Class.values.reshape(460,390))
ax[0, 0].set_title('MOD09GA RGB', fontsize=18)
ax[0, 1].set_title(f'kMeans\n{key[0].upper()+key[1:]} Distance', fontsize=18)

ax[1, 0].imshow(img2)
ax[1, 1].imshow(segmented_2[key].Class.values.reshape(460,390))
ax[1, 0].set_title('MOD09GA RGB', fontsize=18)
ax[1, 1].set_title(f'kMeans\n{key[0].upper()+key[1:]} Distance', fontsize=18)

for i in range(2):
for j in range(2):
ax[i, j].set_facecolor('black')
ax[i, j].set_xticks([])
ax[i, j].set_yticks([])

plt.savefig(f'{key}.png')
plt.tight_layout()
plt.show()

Imagen por autor.
Imagen por autor.
Imagen por autor.

No es difícil ver que la distancia euclidiana y de Manhattan resultó ser la más adecuada para esta tarea en particular. Pero para asegurarnos de que sea cierto, evaluemos los resultados del agrupamiento de k-Means utilizando el coeficiente de silueta. Esta métrica es perfecta para evaluar los resultados del entrenamiento cuando no hay valores verdaderos etiquetados para los puntos agrupados.

Para calcularlo usaremos aprender función [1].

Imagen por aprender.
  • a – la distancia media entre una muestra y todos los demás puntos de la misma clase.
  • b – la distancia media entre una muestra y todos los demás puntos en la siguiente grupo más cercano.

El rango de valores del coeficiente de silueta es [-1,1]. Y sí, es computacionalmente costoso, ya que es necesario calcular distancias entre miles de puntos varias veces, así que prepárate para esperar.

scores_1, scores_2 = {}, {}
for key in segmented_1.keys(): #key is a metric for the distance estimation
scores_1[key]=round(silhouette_score(segmented_1[key].loc[:, :'C3'], segmented_1[key].Class, metric=key),2)
scores_2[key]=round(silhouette_score(segmented_2[key].loc[:, :'C3'], segmented_2[key].Class, metric=key),2)
print(f'Distance: {key}\t Img 1: {scores_1[key]}\t Img 2: {scores_2[key]}')
Imagen por autor.

Ahora puedes ver que lo hemos demostrado: las distancias euclidiana y de Manhattan tienen un rendimiento igualmente bueno, así que estimemos la pérdida de superficie del agua usando ambas.

for metric, Class in zip(['euclidean', 'manhattan'], [2,1]):
img1_water = np.count_nonzero(segmented_1[metric].Class.values == Class)*500*500*1e-6 #pixel size is 500, so the area is 500*500 and to convert to km2 * 1e-6
img2_water = np.count_nonzero(segmented_2[metric].Class.values == Class)*500*500*1e-6

print(f'Distance: {metric}\tWater Area Before: {round(img1_water)}km\u00b2\tWater Area After: {round(img2_water)}km\u00b2\tChange: -{100-round(img2_water/img1_water*100)}%')

Imagen por autor.

— — — —

Distancia: euclidiana
Área de agua antes: 17125 km²
Área de agua después: 1960 km²
Cambio: -89%

— — — — —

Distancia: Manhattan
Área de agua antes: 16244 km²
Área de agua después: 2003 km²
Cambio: -88%

Como puede ver, según nuestros resultados de agrupación, el cambio en el área de la superficie del agua es casi 90% (!!!) pérdida de agua durante los últimos 23 años, lo que es una prueba real de que la contracción del Mar de Aral es una tragedia planetaria…

=============================================

Referencia:

[1] https://scikit-learn.org/stable/modules/generated/sklearn.metrics.silhouette_score.html#sklearn.metrics.silhouette_score

=============================================

Todas mis publicaciones en Medium son gratuitas y de acceso abierto, ¡por eso te agradecería mucho que me siguieras aquí!

PD: Soy un apasionado de la ciencia (geo)datos, el aprendizaje automático/inteligencia artificial y el cambio climático. Entonces, si quieres trabajar juntos en algún proyecto, por favor contáctame en LinkedIn.

🛰️Sigue para más🛰️