Segmentazione semantica delle immagini di telerilevamento utilizzando k-Means |  di Aleksei Rozanov |  Marzo 2024

 | Intelligenza-Artificiale

Da zero in Python🐍

Immagine di autore.

OUno dei modelli ML più semplici e geniali, secondo me, è il clustering k-Means. Si riferisce al gruppo di algoritmi di apprendimento non supervisionato ed è in grado di trovare modelli all'interno di un set di dati senza etichetta. La caratteristica più piacevole è che non presenta calcoli complicati e praticamente qualsiasi studente delle scuole superiori può implementare e utilizzare con successo questo metodo. Quindi in questo articolo voglio condividere come è possibile creare da zero l'algoritmo k-Means in Python utilizzando solo insensato E panda librerie e applicarlo a un problema del mondo reale: la segmentazione semantica delle immagini satellitari.

Innanzitutto parliamo dei dati a nostra disposizione.

In uno dei miei articoli precedenti ho parlato del problema del restringimento del lago d'Aral. Di conseguenza, siamo riusciti a ottenere immagini di telerilevamento da MODIS utilizzando Google Earth Engine, che indicano chiaramente che il mare si sta asciugando. Quindi mi sono chiesto: come possiamo stimare il cambiamento della superficie dell'acqua tra il 2000 e il 2023 utilizzando la segmentazione semantica ML? La risposta è k-Mezzi!

Prima di immergerci nella programmazione, diamo un'occhiata ai dati che utilizzeremo in questo tutorial. Si tratta di due immagini RGB della stessa area con un intervallo di 23 anni, tuttavia è chiaro che le proprietà della superficie terrestre e le condizioni atmosferiche (nuvole, aerosol ecc.) sono diverse. Ecco perché ho deciso di addestrare due modelli k-Means separati, uno per ciascuna immagine.

Immagine di autore.

Per seguire il tutorial, è possibile scaricare ed eseguire il notebook Qui.

Innanzitutto importiamo le librerie necessarie e carichiamo i dati sul notebook:

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')

Puoi vedere che l'area coperta dalle immagini è piuttosto ampia, quindi suggerisco di ingrandire un po':

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()

Immagine di autore.

E l'ultimo passaggio prima della fase ML, convertiamo le nostre immagini in panda dataframe (una colonna per ciascun canale immagine). Lo faccio per motivi di visibilità delle mie spiegazioni. Se vuoi ottimizzarlo, è meglio usarlo insensato invece gli array.

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()})

Allora qual è l’idea alla base dell’algoritmo?

Immagina di giudicare il gusto del cibo utilizzando due criteri: dolcezza e prezzo. Tenendo presente questo, ti darò una serie di possibili opzioni per mangiare:

Immagine di autore.

Scommetto che il tuo cervello ha già suddiviso le opzioni in tre gruppi: frutta, bevande e prodotti da forno. Fondamentalmente, hai inconsciamente raggruppato i dati bidimensionali, che sono definiti da una coppia di valori: (dolcezza; prezzo).

Immagine di autore.

In caso di k-Mezzil'obiettivo dell'algoritmo è abbastanza simile: trovare a preimpostato numero di cluster, Knello spazio n-dimensionale (ad esempio, oltre alla dolcezza e al prezzo, vuoi tenere conto della nutrizione, della salute, della presenza del cibo nel tuo frigorifero e, in questo caso, n = 5).

I. Definire il numero di cluster.

Come ho accennato in precedenza, K in k-Means è il numero di cluster che vuoi ottenere alla fine e dovresti impostare questo valore Prima addestrare il modello.

II. Inizializza casualmente i centroidi.

Il centroide è parte integrante di k-Means. Fondamentalmente, il baricentro è un cerchio con un centro, definito da un insieme di coordinate, e ciascun baricentro rappresenta un cluster. Ad esempio, nel nostro esempio precedente ci sono 3 centroidi.

III. Calcolare le distanze e assegnare i cluster.

Ora dobbiamo trovare quanto dista ciascun punto da ciascun baricentro. Sulla base di questi calcoli, assegniamo ciascun punto al baricentro meno distante (cluster).

IV. Calcolare nuovi centroidi.

Ora ciascuno dei nostri cluster contiene almeno un punto, quindi è il momento di ricalcolare i centroidi semplicemente prendendo le coordinate medie di tutti i punti del cluster.

E questo è tutto! Ripetiamo i passaggi 2–4 finché i centroidi non cambiano.

Immagine di autore.

Ora inseriamo questa idea davvero semplice dietro k-Means nel codice Python.

Promemoria: in questo compito abbiamo 3D problema, cioè il nostro X, Y E Z Sono Rosso, Verde E Blu canali di immagini!

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

Nella prima fase creiamo un elenco l per raccogliere tutte le distanze tra i cluster per visualizzarli successivamente e campionare casualmente K punti dal set di dati per renderli i nostri centroidi (o in alternativa, puoi assegnare valori casuali ai centroidi).

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

Ora dobbiamo calcolare le distanze tra centroidi e punti dati. Esistono molti parametri di distanza diversi nella scienza dei dati, ma concentriamoci su quelli seguenti: Euclideo, Manhattan, Chebyshev.

Per la distanza euclidea:

Immagine di autore.

Per Manhattan:

Immagine di autore.

Per Chebyshev:

Immagine di autore.

Per utilizzare queste formule, scriviamo una funzione versatile per Qualunque numero di dimensioni:

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

Quindi ora possiamo semplicemente calcolare le distanze e assegnare un cluster a ciascun punto dati. Pertanto, i nostri nuovi centroidi sono diventati vecchi, quindi li memorizziamo in un'altra variabile e ricalcoliamo quelli nuovi. Per fare ciò iteriamo su ciascun cluster e prendiamo una media su tutte le coordinate (nel nostro caso, sui canali RGB). Pertanto, la variabile new_centroids ha una forma di (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))

Infine, ripetiamo tutti questi passaggi finché le coordinate dei centroidi non cambiano più. Ho espresso questa condizione in questo modo: la differenza tra le coordinate medie del cluster dovrebbe essere inferiore a 0,001. Ma puoi giocare con altri numeri qui.

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))

E questo è tutto. L'algoritmo è pronto per essere addestrato! Quindi impostiamo k = 3 e memorizziamo i risultati nei dizionari.

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')

Ho deciso di confrontare tutti i parametri di distanza per questo particolare compito, come puoi vedere, ed è evidente che qui la distanza di Manhattan è stata la più veloce.

Immagine di autore.

Prima di visualizzare i cluster, convertiamo i nomi dei cluster nel 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))

Il tempo crea le trame finali!

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()

Immagine di autore.
Immagine di autore.
Immagine di autore.

Non è difficile vedere che la distanza euclidea e quella di Manhattan si sono rivelate le più adatte per questo particolare compito. Ma per assicurarci che sia vero, valutiamo i risultati del clustering k-Means utilizzando il coefficiente Silhouette. Questa metrica è perfetta per la valutazione dei risultati dell'addestramento quando non sono presenti valori reali etichettati per i punti raggruppati.

Per calcolarlo useremo imparato funzione (1).

Immagine di imparato.
  • UN – la distanza media tra un campione e tutti gli altri punti della stessa classe.
  • B – la distanza media tra un campione e tutti gli altri punti nel cluster successivo più vicino.

L'intervallo di valori del coefficiente di silhouette è (-1,1). E sì, è computazionalmente costoso, poiché è necessario calcolare più volte le distanze tra migliaia di punti, quindi preparati ad aspettare.

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)}')
Immagine di autore.

Ora puoi vedere che l'abbiamo dimostrato: le distanze euclidee e Manhattan hanno prestazioni altrettanto buone, quindi stimiamo la perdita di superficie dell'acqua, utilizzandole entrambe.

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)}%')

Immagine di autore.

— — — —

Distanza: euclidea
Superficie d'acqua prima: 17125 km²
Area d'acqua dopo: 1960 km²
Variazione: -89%

— — — — —

Distanza: Manhattan
Superficie d'acqua prima: 16244 km²
Area d'acqua dopo: 2003 km²
Variazione: -88%

Come puoi vedere, secondo i nostri risultati di clustering, il cambiamento nella superficie dell'acqua è quasi Perdita d'acqua del 90% (!!!). negli ultimi 23 anni, il che è la prova concreta del fatto che il restringimento del Lago d’Aral è una tragedia planetaria…

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

Riferimento:

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

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

Tutte le mie pubblicazioni su Medium sono gratuite e ad accesso aperto, ecco perché apprezzerei davvero se mi seguissi qui!

Ps Sono estremamente appassionato di (Geo)Data Science, ML/AI e cambiamento climatico. Quindi, se vuoi lavorare insieme su qualche progetto, contattami LinkedIn.

🛰️Segui per saperne di più🛰️

Fonte: towardsdatascience.com

Lascia un commento

Il tuo indirizzo email non sarà pubblicato. I campi obbligatori sono contrassegnati *